Как изменить случайные эффекты в lmer

Есть ли способ изменить (перезаписать) случайные эффекты в lmer-модели? Для фиксированных эффектов есть слот под названием my_lmer@beta, и я могу изменить фиксированные эффекты, используя:

my_lmer@beta[1] <- 0.5

Есть ли аналогичный способ сделать это для случайных эффектов? Содержит ли lmer-объект уже случайные эффекты или рассчитываются позже, например, ranef().


person winwin    schedule 30.08.2016    source источник
comment
Хорошо, я нашел случайные эффекты в my_lmer@pp$b(1). Но похоже, что изменить этот объект с помощью my_lmer@pp$b(1) <- 0 невозможно.   -  person winwin    schedule 30.08.2016
comment
зачем вам это?   -  person Roland    schedule 30.08.2016


Ответы (1)


Было бы действительно хорошо узнать, что вы пытаетесь сделать. Изменение внутренних компонентов объектов ссылочного класса особенно опасно - вы можете случайно изменить копии объектов или вызывают ошибки сегментации ... из здесь,

library(lme4)
fm1 <- lmer(Reaction~Days+(1|Subject),sleepstudy) ## just for example
fm1@pp$getRefClass()$methods()

покажет вам методы ... однако вы должны пойти немного глубже ... оказывается (l. 85 из src/predModule.cpp), что то, что b на самом деле делает, чтобы взять внутренний u

 VectorXd merPredD::b(const double& f) const {return d_Lambdat.adjoint() * u(f);}

который, в свою очередь, вызывает

VectorXd merPredD::u(const double& f) const {return d_u0 + f * d_delu;}

это означает, что для изменения b вам нужно будет изменить соответствующие значения u0; в настоящее время я не думаю, что это возможно.

Для справки, это код (из здесь) который оценивает отклонение, когда случайные эффекты смещаются на (вектор) z от их оценочных значений ...

rr <- m@resp                        ## extract response module
u0 <- getME(m,"u")                  ## conditional modes
L <- getME(m,"L")
## sd <- 1/getME(pp,"L")@x
## filled elements of L matrix==diag for simple case
## for more general case need the following -- still efficient
sd <- sqrt(diag(chol2inv(L)))
## fixed-effects contribution to linear predictor
fc <- getME(m,"X") %*% getME(m,"beta")
ZL <- t(getME(m,"Lambdat") %*% getME(m,"Zt"))
## evaluate the unscaled conditional density on the deviance scale
dc <- function(z) {
    uu <- u0 + z * sd    ## displace conditional modes
    ##  should still work if z is a vector (by recycling, because u values
    ##  applying to each group are stored adjacent to each other)
    rr$updateMu(fc + ZL %*% uu)     ## update linear predictor
    drc <- unname(as.vector(tapply(rr$devResid(), ff, sum)))
    uuc <- colSums(matrix(uu * uu,nrow=nvar))
    (drc + uuc)[grps]
}
person Ben Bolker    schedule 31.08.2016
comment
Спасибо за ответы! Намерение состояло в том, чтобы сравнить прогнозируемые значения от модели к гипотетической второй модели с заданными оценками. В основном я хотел перейти к оценкам на данный априор и посмотреть, насколько вероятны предсказанные значения. Я подумал, что проще отредактировать lmer-объект и использовать метод predict, а не умножать его на себя. - person winwin; 31.08.2016
comment
на самом деле, я думаю, что было бы проще размножить его самому; вы можете извлечь матрицы X и (транспонированные) Z через getME(fitted_model,c("X","Zt")); тогда X %*% beta + t(Zt) %*% b должен дать вам прогнозы ... - person Ben Bolker; 01.09.2016
comment
Спасибо за совет! - person winwin; 02.09.2016