Рисование квадратичной функции в модели смешанного эффекта

Это немного длинный вопрос, так что спасибо, что выдержали меня.

Вот мои данные

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

Я построил модель смешанного эффекта

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)

Все предикторы стандартизированы, и мне интересно узнать, как y изменяется с изменениями в x5, сохраняя при этом средние значения других переменных (равные 0, поскольку все переменные стандартизированы).

Вот как я это делаю.

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

Теперь я хочу построить подгонку как квадратичную функцию x5. Я делаю это:

i<-order(data$x5)  

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

Я ожидал, что это даст график y как квадратичную функцию x5. Но, как видите, я получаю следующий график, в котором нет квадратичной кривой. Может ли кто-нибудь сказать мне, что я делаю неправильно здесь?

Вот что я получаю


r regression lme4 quadratic
person user53020    schedule 07.04.2017    source источник
comment
Да. У меня есть это.   -  person user53020    schedule 07.04.2017
comment
Я только что изменил свой вопрос, чтобы отразить то, что я использовал для прогнозирования с помощью merTools.   -  person user53020    schedule 07.04.2017
comment
Это немного сбивает меня с толку. Если я извлеку коэффициент x5 из mod и самостоятельно нарисую квадрат, будет ли это учитывать тот факт, что другие предикторы остаются постоянными. Извините, моя математика/статистика немного слаба. отсюда путаница.   -  person user53020    schedule 07.04.2017
comment
Не могли бы вы опубликовать это как решение, чтобы я мог отметить его и посмотреть, правильно ли я понял. Спасибо   -  person user53020    schedule 07.04.2017

Ответы (1)


arrow_upward
2
arrow_downward

Я не уверен, откуда взялось predictInterval, но вы можете сделать это с помощью predict. Хитрость заключается в том, чтобы убедиться, что вы установили свои случайные эффекты на 0. Вот как вы можете это сделать.

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

Часть re.form=NA выпадает из случайного эффекта

введите здесь описание изображения

person MrFlick    schedule 07.04.2017
comment
Спасибо. predictInterval исходит из библиотеки (merTools) - person user53020; 07.04.2017