Сохранение постоянной дисперсии случайных эффектов в PROC MIXED по сравнению с lmer

Мне было интересно, можно ли поддерживать постоянную дисперсию случайных эффектов в функциях R's lme или lmer (или другой подпрограмме случайных эффектов в R) или, по крайней мере, указать начальные значения.

В SAS это возможно с помощью оператора parms в PROC MIXED. В статье Selya et al. (2012) авторы используют это, чтобы установить параметры дисперсии для модели с более простой структурой с фиксированными эффектами в соответствии с параметрами полной модели.

Конкретный вызов в PROC MIXED, который они используют, таков: parms/parmsdata = fullmodel.AB hold = ... Их цель состоит в том, чтобы поддерживать постоянные оценки дисперсии в моделях с различными фиксированными структурами эффектов (хотя мне интересно, действительно ли это возможно в SAS или R).


sas lme4 mixed-models
person N Brouwer    schedule 16.02.2015    source источник
comment
Поскольку (во флаге) вы характеризуете это как вопрос кодирования и просите перенести его в переполнение стека, я Сделай так. Вполне возможно, что вы не сможете получить желаемый результат напрямую из любого существующего пакета R, и в этом случае потребуется некоторый теоретический анализ, чтобы показать, как выполнять вычисления в первую очередь, что вернет вас прямо сюда, к сайт статистики...   -  person whuber    schedule 18.02.2015

Ответы (1)


arrow_upward
1
arrow_downward

Попробуйте gnlrim на github. Это может сделать максимальное правдоподобие с ограничениями на параметры. Просто установите начальное значение, нижнюю границу и верхнюю границу на одно и то же значение для случайной дисперсии перехвата pmix, и она будет поддерживаться постоянной.

Ниже приведен пример, который показывает, что gnlrim оценивает ту же модель, что и lmer с REML=FALSE. Первые фрагменты предназначены для простого копирования и вставки; последующие фрагменты показывают выполнение соответствующих строк.

Настройте пакеты, данные и подходящие модели (скопируйте и вставьте фрагмент):

library(devtools)
## if you have macOS, grab this version of libstableR:
devtools::install_github("hrbrmstr / libstableR")
devtools::install_github(  "swihart/gnlrim")

## data: 4 individuals with 5 observations
dose <- c(9,12,4,9,11,10,2,11,12,9,9,9,4,9,11,9,14,7,9,8)
y <- c(8.674419, 11.506066, 11.386742, 27.414532, 12.135699,  4.359469,
       1.900681, 17.425948,  4.503345,  2.691792,  5.731100, 10.534971,
       11.220260,  6.968932,  4.094357, 16.393806, 14.656584,  8.786133,
       20.972267, 17.178012)
id <- rep(1:4, each=5)

## fit with lmer
lmer_fit <- lme4::lmer(y~dose + (1|id), REML=FALSE)

## fit with gnlrim
gnlrim_fit <-
gnlrim(y,
       mu=~a+b*dose+rand,
       random="rand",
       nest=id,
       pmu=c(a=8.7,b=0.25),
       pshape = c(shape=1),
       pmix=c(var=3.0938^2),
       p_uppb = c(10,  1, 5, 3.0938^3),
       p_lowb = c( 5, -1, 0, 0)
       )

Запустите эти строки, чтобы увидеть, насколько похожи модели (скопируйте и вставьте фрагмент):

## show fits are the same:

## intercept (a) slope (b)
summary(lmer_fit)$coeff[,1]
gnlrim_fit$coeff[1:2]

## Residuals standard deviation
## sigma_epsilon = 5.58
summary(lmer_fit)$varcor
sqrt(exp(gnlrim_fit$coeff[3]))

## random effects standard deviation
## sigma_id = 3.0938
summary(lmer_fit)$varcor
sqrt(gnlrim_fit$coeff[4])

## likelihood
summary(lmer_fit)$logLik
-gnlrim_fit$maxlike

Возьмите ту же модель, но запустите ее, установив и удерживая константу дисперсии случайных эффектов (скопируйте и вставьте фрагмент):

## Take same model but hold constant
## random effects standard deviation
## sigma_id   :=  9
## sigma^2_id := 81
gnlrim_fit2 <-
  gnlrim(y,
         mu=~a+b*dose+rand,
         random="rand",
         nest=id,
         pmu=c(a=8.7,b=0.25),
         pshape = c(shape=1),
         pmix=c(var=9^2),
         p_uppb = c(10,  1, 5, 9^2),
         p_lowb = c( 5, -1, 0, 9^2)
  )

gnlrim_fit2$coeff
gnlrim_fit2$se

Исполнение строк показа моделей аналогичное:

> ## show fits are the same:
> ## intercept (a) slope (b)
> summary(lmer_fit)$coeff[,1]
(Intercept)        dose 
  8.7117914   0.2488724 
> gnlrim_fit$coeff[1:2]
[1] 8.7118426 0.2488648
> 
> ## Residuals standard deviation
> ## sigma_epsilon = 5.58
> summary(lmer_fit)$varcor
 Groups   Name        Std.Dev.
 id       (Intercept) 3.0938  
 Residual             5.5880  
> sqrt(exp(gnlrim_fit$coeff[3]))
[1] 5.587926
> 
> ## random effects standard deviation
> ## sigma_id = 3.0938
> summary(lmer_fit)$varcor
 Groups   Name        Std.Dev.
 id       (Intercept) 3.0938  
 Residual             5.5880  
> sqrt(gnlrim_fit$coeff[4])
[1] 3.094191
> 
> ## likelihood
> summary(lmer_fit)$logLik
'log Lik.' -64.64964 (df=4)
> -gnlrim_fit$maxlike
[1] -64.64958

Выполнение линий для удержания постоянной дисперсии:

> ## Take same model but hold constant
> ## random effects standard deviation
> ## sigma_id   :=  9
> ## sigma^2_id := 81
> gnlrim_fit2 <-
+   gnlrim(y,
+          mu=~a+b*dose+rand,
+          random="rand",
+          nest=id,
+          pmu=c(a=8.7,b=0.25),
+          pshape = c(shape=1),
+          pmix=c(var=9^2),
+          p_uppb = c(10,  1, 5, 9^2),
+          p_lowb = c( 5, -1, 0, 9^2)
+   )
> 
> gnlrim_fit2$coeff
[1]  9.1349920  0.2012785  3.4258404 81.0000000
> gnlrim_fit2$se
[1] 6.1006729 0.4420228 0.3485940 0.0000000
person swihart    schedule 04.11.2017