0

It's a bit of a long question so thanks for bearing with me.

Here's my data

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

I constructed a mixed effect model

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)

All the predictors are standardised and I am interested in knowing how does y changes with changes in x5while keeping other variables at their mean values (equal to 0 since all the variables are standardised).

This is how I do it.

# 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)

Now I want to plot the fitted as a quadratic function of x5. I do this:

i<-order(data$x5)  

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

I expected this to produce a plot of y as a quadratic function of x5. But As you can see, I get the following plot which does not have any quadratic curve. Can anyone tell me what I am doing wrong here?

This is what I get

user53020
  • 889
  • 2
  • 10
  • 33
  • Yes. I have it. – user53020 Apr 07 '17 at 19:14
  • I have just modified my question to reflect what I used to predict using the merTools – user53020 Apr 07 '17 at 19:17
  • It's a bit confusing for me. If I extract the coefficient of x5 from mod and draw the quadratic myself, will it account for the fact that other predictors are held constant. Sorry my maths/statistics is a bit weak. hence the confusion. – user53020 Apr 07 '17 at 19:22
  • Will it be possible for you to post this as a solution, so that I can mark it and see if I have got it right. Thanks – user53020 Apr 07 '17 at 19:26
  • 1
    `plot(data$x5, 49.59507 + data$x5 * -1.09963 + data$x5^2 * -1.28696)` – erc Apr 07 '17 at 19:38

1 Answers1

2

I'm not sure where predictInterval comes from, but you can do this with predict. The trick is just to make sure you set your random effects to 0. Here's how you can do that

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

The re.form=NA part drops out the random effect

enter image description here

MrFlick
  • 195,160
  • 17
  • 277
  • 295