Фиксирование значений дисперсии в lme4

Я использую lme4 Пакет R для создания линейной смешанной модели с использованием lmer() функция. В этой модели у меня есть четыре случайных эффекта и один фиксированный эффект (перехват). Мой вопрос о предполагаемых дисперсиях случайных эффектов. Можно ли задавать начальные значения параметров ковариации аналогично тому, как это можно сделать в SAS с PARMS аргумент.

В следующем примере предполагаемые отклонения:

c(0.00000, 0.03716, 0.00000, 0.02306)

Я хотел бы исправить это (например)

c(0.09902947, 0.02460464, 0.05848691, 0.06093686)

так что не оценили.

    > summary(mod1)
    Linear mixed model fit by maximum likelihood  ['lmerMod']
    Formula: log_cumcover_mod ~ (1 | kildestationsnavn) + (1 | year) + (1 |  
        kildestationsnavn:year) + (1 | proevetager)
       Data: res

         AIC      BIC   logLik deviance df.resid 
       109.9    122.9    -48.9     97.9       59 

    Scaled residuals: 
        Min      1Q  Median      3Q     Max 
    -2.1056 -0.6831  0.2094  0.8204  1.7574 

    Random effects:
     Groups                 Name        Variance Std.Dev.
     kildestationsnavn:year (Intercept) 0.00000  0.0000  
     kildestationsnavn      (Intercept) 0.03716  0.1928  
     proevetager            (Intercept) 0.00000  0.0000  
     year                   (Intercept) 0.02306  0.1518  
     Residual                           0.23975  0.4896  
    Number of obs: 65, groups:  
    kildestationsnavn:year, 6; kildestationsnavn, 3; proevetager, 2; year, 2

    Fixed effects:
                Estimate Std. Error t value
    (Intercept)   4.9379     0.1672   29.54
Источник

1 ответ

Решение

Это возможно, если немного взломать. Вот воспроизводимый пример:

Подходит оригинальная модель:

library(lme4)
set.seed(101)
ss <- sleepstudy[sample(nrow(sleepstudy),size=round(0.9*nrow(sleepstudy))),]
m1 <- lmer(Reaction~Days+(1|Subject)+(0+Days|Subject),ss)
fixef(m1)
## (Intercept)        Days 
##   251.55172    10.37874

Восстановите функцию отклонения (в данном случае REML-критерий):

dd <- as.function(m1)

Я собираюсь установить стандартные отклонения на ноль, чтобы мне было с чем сравнивать, то есть коэффициенты регулярной линейной модели. (Вектор параметров для dd представляет собой вектор, содержащий столбчатые, нижнетреугольные, сцепленные факторы Холецкого для масштабированных членов случайных эффектов в модели. К счастью, если у вас есть только скалярные / перехватывающие случайные эффекты (например, (1|x)), то они соответствуют стандартным отклонениям случайных эффектов, масштабированным по стандартному отклонению модели).

(ff <- dd(c(0,0)))  ## new REML: 1704.708
environment(dd)$pp$beta(1)  ## new parameters
##    [1] 251.11920  10.56979

Матчи:

coef(lm(Reaction~Days,ss))
## (Intercept)        Days 
##   251.11920    10.56979

Если вы хотите построить новый merMod объект вы можете сделать это следующим образом...

opt <- list(par=c(0,0),fval=ff,conv=0)
lmod <- lFormula(Reaction~Days+(1|Subject)+(0+Days|Subject),ss)
m1X <- mkMerMod(environment(dd), opt, lmod$reTrms, fr = lmod$fr,
         mc = quote(hacked_lmer()))

Теперь предположим, что мы хотим установить дисперсии для определенного ненулевого значения (скажем, (700,30)). Это будет немного сложнее из-за масштабирования по остаточному стандартному отклонению...

newvar <- c(700,30)
ff2 <- dd(sqrt(newvar)/sigma(m1))
opt2 <- list(par=c(0,0),fval=ff,conv=0)
m2X <- mkMerMod(environment(dd), opt, lmod$reTrms, fr = lmod$fr,
         mc = quote(hacked_lmer()))
VarCorr(m2X)
unlist(VarCorr(m2X))
##   Subject Subject.1 
## 710.89304  30.46684

Так что это не дает нам того, чего мы хотели (потому что остаточная дисперсия меняется...)

buildMM <- function(theta) {
   dd <- as.function(m1)
   ff <- dd(theta)
   opt <- list(par=c(0,0),fval=ff,conv=0)
   mm <- mkMerMod(environment(dd), opt, lmod$reTrms, fr = lmod$fr,
         mc = quote(hacked_lmer()))
   return(mm)
}

objfun <- function(x,target=c(700,30)) {
   mm <- buildMM(sqrt(x))
   return(sum((unlist(VarCorr(mm))-target)^2))
}
s0 <- c(700,30)/sigma(m1)^2
opt <- optim(fn=objfun,par=s0)
mm_final <- buildMM(sqrt(opt$par))
summary(mm_final)
##  Random effects:
##  Groups    Name        Variance Std.Dev.
##  Subject   (Intercept) 700      26.458  
##  Subject.1 Days         30       5.477  
##  Residual              700      26.458  
## Number of obs: 162, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  251.580      7.330   34.32
## Days          10.378      1.479    7.02

Между прочим, обычно не рекомендуется использовать случайные эффекты, когда группирующие переменные имеют очень небольшое количество (например, <5 или 6) уровней: см. Здесь...

Источник

Для больших наборов данных, возможно, стоит отметить, что шаг оптимизации в ответе Бена Болкера можно упростить до одномерной задачи оптимизации. Финалtheta всегда будет скейлером c(700, 30).

Особенно, если есть больше скалярных / случайных эффектов только для перехвата, стоит изменить objfun следующим образом:

objfun <- function(x,target=c(700,30)) {
  scaled_theta <- s0*x
  mm <- buildMM(scaled_theta)
  return(sum((unlist(VarCorr(mm))-target)^2))
}
s0 <- sqrt(c(700,30)/sigma(m1)^2)
opt <- optim(fn=objfun,par=1, method = "L-BFGS-B", lower = 0)
Источник
Другие вопросы по тегам