Доверительный интервал случайных эффектов с lmer

Я использую lmer из пакета lme4 для расчета доверительного интервала для компонента дисперсии.

Когда я подбираю модель, появляются предупреждающие сообщения:

fit <- lmer(Y~X+Z+X:Z+(X|group),data=sim_data)
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

Я много искал, чтобы понять, почему возникает ошибка, и, наконец, пришел к выводу, что there is difference between error and warning in the R world .

Я хочу вычислить доверительный интервал для параметров модели и запустить код, который показывает ошибку:

 confint.merMod(fit,oldNames=FALSE)
 Computing profile confidence intervals ...
 Error in if (all(x[iu <- upper.tri(x)] == 0)) t(x[!iu]) else t(x)[!iu] : 
 missing value where TRUE/FALSE needed

Есть ли другой способ получить CI случайных эффектов с lmer?

РЕДАКТИРОВАТЬ :

simfun <- function(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1){
    N <- sum(rep(n_j,J))  
    x <- rnorm(N)         
    z <- rnorm(J)        

    mu <- c(0,0)
    sig <- matrix(c(sig2_0,sig01,sig01,sig2_1),ncol=2)
    u   <- rmvnorm(J,mean=mu,sigma=sig)

    b_0j <- g00 + g01*z + u[,1]
    b_1j <- g10 + g11*z + u[,2]

    y <- rep(b_0j,each=n_j)+rep(b_1j,each=n_j)*x + rnorm(N,0,0.5)
    data <- data.frame(Y=y,X=x,Z=rep(z,each=n_j),group=rep(1:J,each=n_j))
  } 

noncoverage <- function(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1){
    dat <- simfun(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1)
    fit <- lmer(Y~X+Z+X:Z+(X|group),data=dat)
 }

comb1 = replicate(1000,noncoverage(10,5,1,.3,.3,.3,(1/18),0,(1/18)))
comb26 = replicate(1000,noncoverage(100,50,1,.3,.3,.3,(1/8),0,(1/8)))

r statistics lme4 mixed-models lmer
person user81411    schedule 27.07.2015    source источник
comment
Если вы получаете предупреждения о сходимости, вы можете ожидать проблем в некоторых последующих анализах. confint — это одна из функций, которые обычно не работают в таком случае. Вы должны быть осторожны и не должны доверять этой модели. Исследуйте, почему он не сошелся.   -  person Roland    schedule 27.07.2015
comment
нам, вероятно, нужно выяснить основную причину предупреждений (хотя это полезно для вас, чтобы понять разницу между ошибками и предупреждениями). Вы читали ?convergence страницу руководства, особенно ниже. Если вы видите предупреждения о конвергенции?) Можете ли вы опубликовать дополнительную информацию о своей проблеме (например, воспроизводимый пример)?   -  person Ben Bolker    schedule 27.07.2015

Ответы (4)


arrow_upward
3
arrow_downward

Это зависит от того, что именно вы ищете из доверительных интервалов, но функция sim в пакете arm предоставляет отличный способ получить повторные выборки из апостериорного значения объекта lmer или glmer, чтобы получить представление о изменчивости коэффициентов как фиксированные, так и случайные условия.

В пакете merTools мы написали обертку, которая упрощает процесс извлечения этих значений и взаимодействия с ними:

library(merTools)
randomSims <- REsim(fit, n.sims = 500)
# and to plot it
plotREsim(REsim(fit, n.sims = 500))

В merTools есть ряд других инструментов для их изучения. Однако, если вам нужны фактические результирующие симуляции, вам лучше использовать arm::sim.

person jknowles    schedule 13.08.2015
comment
только что наткнулся на ваш пакет merTools - опция plotREsim (REsim (fit)) просто фантастическая - у вас есть возможность преобразовать шкалу из z-распределения обратно в вероятности, сохраняя более темный цвет, если группы отличаются знаком. из среднего? - person cath's_lab; 30.10.2019
comment
Вам следует открыть новый вопрос на SO или опубликовать вопрос на GitHub (похоже, нам не помешало бы улучшить документацию), и кто-то поможет вам получить то, что вам нужно. - person jknowles; 30.10.2019

arrow_upward
1
arrow_downward

Более простая модель показывает результаты в виде матрицы. Вы можете получить доступ к оценкам и стандартным ошибкам модели, чтобы вычислить доверительные интервалы.

Поскольку первая строка — это оценка точки пересечения модели, вторая строка — это та, которая показывает оценку управляемой вами переменной. Первый столбец — это оценка эффекта, а второй столбец — стандартная ошибка.

Таким образом, заменив MyModel на имя вашей модели и округлив число до двух, round(coef(summary(MyModel))[2,1],2)-round(coef(summary(MyModel))[2,2],2)*2 даст вам нижний предел доверительного интервала, в то время как простое изменение вычитания для сложения в предыдущей формуле даст вам верхний предел оценки: round(coef(summary(MyModel))[2,1],2)+round(coef(summary(MyModel))[2,2],2)*2

person HernanLG    schedule 07.11.2017

arrow_upward
1
arrow_downward

Вы можете использовать parameters-package, который предлагает вам методы вычисления CI для фиксированных эффектов, включая Kenward -Аппроксимация Роджера для небольших размеров выборки, или стандартные ошибки для фиксированных и случайных эффектов, или полные сводки модели. См. примеры ниже.

library(parameters)
library(lme4)
#> Loading required package: Matrix
model <- lmer(mpg ~ wt + (1 | gear), data = mtcars)

ci(model)
#>     Parameter CI   CI_low   CI_high
#> 1 (Intercept) 95 31.89749 40.482616
#> 2          wt 95 -6.29877 -3.791242

ci_kenward(model)
#>     Parameter CI    CI_low   CI_high
#> 1 (Intercept) 95 31.208589 41.171513
#> 2          wt 95 -6.573142 -3.516869

standard_error(model)
#>     Parameter        SE
#> 1 (Intercept) 2.1901246
#> 2          wt 0.6396873

standard_error(model, effects = "random")
#> $gear
#>   (Intercept)
#> 3   0.6457169
#> 4   0.6994964
#> 5   0.9067223

model_parameters(model, df_method = "satterthwaite")
#> Parameter   | Coefficient |   SE |         95% CI |     t |    df |      p
#> --------------------------------------------------------------------------
#> (Intercept) |       36.19 | 2.19 | [31.90, 40.48] | 16.52 | 13.85 | < .001
#> wt          |       -5.05 | 0.64 | [-6.30, -3.79] | -7.89 | 21.92 | < .001

Создано 07 февраля 2020 г. с помощью пакета reprex (v0.3.0)

person Daniel    schedule 07.02.2020

arrow_upward
0
arrow_downward

Проблема здесь в том, что ваш случайный эффект равен (X|группа), и я предполагаю, что он должен быть (1|группа), модель случайного перехвата или включать и то, и другое как (1+X|группа). Наличие только (X | group) вызовет проблемы, поскольку вы допускаете различные наклоны, но не позволяете варьировать точку пересечения.

person Ken Beath    schedule 07.02.2020