2017-04-07 1 views
0

C'est une question un peu longue, merci de me supporter.Dessiner une fonction quadratique dans un modèle à effets mixtes

Voilà mes données

https://www.dropbox.com/s/jo22d68a8vxwg63/data.csv?dl=0

Je construit un modèle d'effet mixte

library(lme4) 
mod <- lmer(sqrt(y) ~ x1 + I(x1^2) + x2 + I(x2^2) + x3 + I(x3^2) + x4 + I(x4^2) + x5 + I(x5^2) + 
     x6 + I(x6^2) + x7 + I(x7^2) + x8 + I(x8^2) + (1|loc) + (1|year), data = data) 

Tous les prédicteurs sont standardisés et je suis intéressé à savoir comment le fait y changements avec les changements dans x5 tout en conservant autres variables à leurs valeurs moyennes (égal à 0 puisque toutes les variables sont normalisées).

Voici comment je le fais.

# make all predictors except x5 equal to zero 

data$x1<-0 
data$x2<-0 
data$x3<-0 
data$x4<-0 
data$x6<-0 
data$x7<-0 
data$x8<-0 

# Use the predict function 
library(merTools) 
fitted <- predictInterval(merMod = mod, newdata = data, level = 0.95, n.sims = 1000,stat = "median",include.resid.var = TRUE) 

Maintenant, je veux tracer la équipée en fonction du second degré de x5. Je fais ceci:

i<-order(data$x5) 

plot(data$x5[i],fitted$fit[i],type="l") 

Je me attendais à ce produire une parcelle de y en fonction quadratique de x5. Mais comme vous pouvez le voir, j'obtiens le graphique suivant qui n'a pas de courbe quadratique. Quelqu'un peut-il me dire ce que je fais mal ici?

This is what I get

+0

Oui. Je l'ai. – user53020

+0

Je viens de modifier ma question pour refléter ce que j'avais l'habitude de prédire en utilisant les merTools – user53020

+0

C'est un peu déroutant pour moi. Si j'extrais le coefficient de x5 de mod et dessine moi-même le quadratique, cela rendra-t-il compte du fait que les autres prédicteurs sont maintenus constants. Désolé mes maths/statistiques est un peu faible. d'où la confusion. – user53020

Répondre

2

Je ne sais pas où predictInterval vient, mais vous pouvez le faire avec predict. L'astuce est juste pour vous assurer que vous définissez vos effets aléatoires à 0. Voici comment vous pouvez le faire

newdata <- data 
newdata[,paste0("x", setdiff(1:8,5))] <- 0 
y <- predict(mod, newdata=newdata, re.form=NA) 
plot(data$x5, y) 

La partie re.form=NA abandonne l'effet aléatoire

enter image description here

+0

Merci. 'predictInterval' vient de la bibliothèque (merTools) – user53020