Почему введение эффекта случайного уклона увеличило SE наклона?

Я пытаюсь проанализировать влияние года на переменную logInd для конкретной группы лиц (у меня есть 3 группы). Самая простая модель:

> fix1 = lm(logInd ~ 0 + Group + Year:Group, data = mydata)
> summary(fix1)

Call:
lm(formula = logInd ~ 0 + Group + Year:Group, data = mydata)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.5835 -0.3543 -0.0024  0.3944  4.7294 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
Group1       4.6395740  0.0466217  99.515  < 2e-16 ***
Group2       4.8094268  0.0534118  90.044  < 2e-16 ***
Group3       4.5607287  0.0561066  81.287  < 2e-16 ***
Group1:Year -0.0084165  0.0027144  -3.101  0.00195 ** 
Group2:Year  0.0032369  0.0031098   1.041  0.29802    
Group3:Year  0.0006081  0.0032666   0.186  0.85235    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.7926 on 2981 degrees of freedom
Multiple R-squared: 0.9717,     Adjusted R-squared: 0.9716 
F-statistic: 1.705e+04 on 6 and 2981 DF,  p-value: < 2.2e-16

Мы видим, что Group1 значительно снижается, а Groups2 и 3 увеличиваются, но незначительно.

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

> mix1a = lmer(logInd ~ 0 + Group + Year:Group + (1|Individual), data = mydata)
> summary(mix1a)
Linear mixed model fit by REML 
Formula: logInd ~ 0 + Group + Year:Group + (1 | Individual) 
   Data: mydata 
  AIC  BIC logLik deviance REMLdev
 4727 4775  -2356     4671    4711
Random effects:
 Groups     Name        Variance Std.Dev.
 Individual (Intercept) 0.39357  0.62735 
 Residual               0.24532  0.49530 
Number of obs: 2987, groups: Individual, 103

Fixed effects:
              Estimate Std. Error t value
Group1       4.6395740  0.1010868   45.90
Group2       4.8094268  0.1158095   41.53
Group3       4.5607287  0.1216522   37.49
Group1:Year -0.0084165  0.0016963   -4.96
Group2:Year  0.0032369  0.0019433    1.67
Group3:Year  0.0006081  0.0020414    0.30

Correlation of Fixed Effects:
            Group1 Group2 Group3 Grp1:Y Grp2:Y
Group2       0.000                            
Group3       0.000  0.000                     
Group1:Year -0.252  0.000  0.000              
Group2:Year  0.000 -0.252  0.000  0.000       
Group3:Year  0.000  0.000 -0.252  0.000  0.000

Это дало ожидаемый эффект - SE уклонов (коэффициенты Group1-3: Year) теперь ниже, а остаточный SE также ниже.

Люди также отличаются по наклону, поэтому я также ввел эффект случайного наклона:

> mix1c = lmer(logInd ~ 0 + Group + Year:Group + (1 + Year|Individual), data = mydata)
> summary(mix1c)
Linear mixed model fit by REML 
Formula: logInd ~ 0 + Group + Year:Group + (1 + Year | Individual) 
   Data: mydata 
  AIC  BIC logLik deviance REMLdev
 2941 3001  -1461     2885    2921
Random effects:
 Groups     Name        Variance  Std.Dev. Corr   
 Individual (Intercept) 0.1054790 0.324775        
            Year        0.0017447 0.041769 -0.246 
 Residual               0.1223920 0.349846        
Number of obs: 2987, groups: Individual, 103

Fixed effects:
              Estimate Std. Error t value
Group1       4.6395740  0.0541746   85.64
Group2       4.8094268  0.0620648   77.49
Group3       4.5607287  0.0651960   69.95
Group1:Year -0.0084165  0.0065557   -1.28
Group2:Year  0.0032369  0.0075105    0.43
Group3:Year  0.0006081  0.0078894    0.08

Correlation of Fixed Effects:
            Group1 Group2 Group3 Grp1:Y Grp2:Y
Group2       0.000                            
Group3       0.000  0.000                     
Group1:Year -0.285  0.000  0.000              
Group2:Year  0.000 -0.285  0.000  0.000       
Group3:Year  0.000  0.000 -0.285  0.000  0.000

Но теперь, вопреки ожиданиям, SE склонов (коэффициенты Group1-3: Year) теперь намного выше, даже выше, чем без какого-либо случайного эффекта вообще!

Как это возможно? Я ожидаю, что случайный эффект «съест» необъяснимую изменчивость и увеличит «достоверность» оценки!

Однако остаточный SE ведет себя как ожидалось - он ниже, чем в модели случайного перехвата.

Вот данные, если это необходимо.

редактировать

Теперь я осознал удивительный факт. Если я делаю линейную регрессию для каждого отдельного человека и затем запускаю ANOVA для результирующих уклонов, я получаю точно такой же результат, как и модель случайных уклонов ! Знаете ли вы, почему?

indivSlope = c()
for (indiv in 1:103) {
    mod1 = lm(logInd ~ Year, data = mydata[mydata$Individual == indiv,])
    indivSlope[indiv] = coef(mod1)['Year']
}

indivGroup = unique(mydata[,c("Individual", "Group")])[,"Group"]


anova1 = lm(indivSlope ~ 0 + indivGroup)
summary(anova1)

Call:
lm(formula = indivSlope ~ 0 + indivGroup)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.176288 -0.016502  0.004692  0.020316  0.153086 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
indivGroup1 -0.0084165  0.0065555  -1.284    0.202
indivGroup2  0.0032369  0.0075103   0.431    0.667
indivGroup3  0.0006081  0.0078892   0.077    0.939

Residual standard error: 0.04248 on 100 degrees of freedom
Multiple R-squared: 0.01807,    Adjusted R-squared: -0.01139 
F-statistic: 0.6133 on 3 and 100 DF,  p-value: 0.6079

Вот данные, если это необходимо.

r mixed-model lme4-nlme random-effects-model любознательный
источник

Вам нужен годовой фиксированный эффект, если у вас будет годовой фиксированный эффект взаимодействия. В общем, вы не можете включить термин взаимодействия, не включая также основные эффекты. Вы действительно думаете, что нет фиксированного компонента эффекта года? И, если это так, как может быть фиксированный год: групповое взаимодействие?

Джон

И почему нет фиксированного перехвата? Вы можете иметь как фиксированные, так и случайные.

Джон

@ Джон, эта модель полностью действительна. Это только вопрос желаемого кодирования категориальной переменной. Таким образом, - это перехват группы , а - это наклон внутри группы . Если включить основное влияние Года и пересечения, то оценки будут представлять собой разницу между пересечением Группы и Группы 1, а также аналогично наклонам. Group

i

$i$

i

$i$ Group

i

$i$ :Year

i

$i$

i

$i$

Анико

@ Джон, тем не менее, это не по теме моего вопроса: поверь мне, это нормально, я много экспериментировал с этим. Моя первая модель lm полностью эквивалентна logInd ~ Year*Group, только коэффициенты имеют различную форму, не более того. Зависит от вашего вкуса и того, какая форма коэффициентов вам нравится, не более того. В моей 1-ой модели не исключено «годовой основной эффект», поскольку вы пишете ... logInd ~ Year*Groupделает то же самое, тогда Yearкоэффициент является не основным эффектом, а Group1: Year.

Любопытно

Хорошо, аккуратно, не считал и перехват 0, и группу категоричными.

Джон

Ответы:

Я думаю, что проблема с вашими ожиданиями :) Обратите внимание, что, когда вы добавили случайный перехват для каждого человека, стандартная ошибка перехватов увеличилась. Поскольку каждый человек может иметь свой собственный перехват, средний показатель по группе менее определен. То же самое произошло со случайным уклоном: вы больше не оцениваете один общий (внутри группы) уклон, а усредняете переменные уклоны.

РЕДАКТИРОВАТЬ: Почему лучшая модель не дает более точную оценку?

Давайте подумаем об этом наоборот: почему исходная модель недооценивает стандартную ошибку? Это предполагает независимость наблюдений, которые не являются независимыми. Вторая модель ослабляет это предположение (таким образом, что влияет на перехваты), а третья ослабляет его еще больше.

РЕДАКТИРОВАТЬ 2: отношения со многими конкретными пациентами моделей

Ваше наблюдение является известным свойством (и если бы у вас было всего два года, то модель случайных эффектов была бы эквивалентна парному t-критерию). Я не думаю, что смогу справиться с реальными доказательствами, но, возможно, выписав две модели, отношения станут яснее. Давайте проигнорируем группирующую переменную, так как это только усложнит запись. Я буду использовать греческие буквы для случайных эффектов и латинские буквы для фиксированных эффектов.

Модель случайных эффектов ( - субъект, - репликация внутри субъекта): где и . $i$ $j$

Y_{i j} = a + α_{i} + (b + β_{i}) x_{i j} + ϵ_{i j},

$Y_{ij} = a + \alpha_i + (b+\beta_i)x_{ij} + \epsilon_{ij},$

(α_{i}, β_{i})^{'} \sim N (0, Σ)

$(\alpha_i,\beta_i)'\sim N(0,\Sigma)$

ϵ_{i j} \sim N (0, σ^{2})

$\epsilon_{ij}\sim N(0,\sigma^2)$

Когда вы подгоняете отдельные модели для каждого предмета, тогда где .

Y_{i j} = a_{i} + b_{i} x_{i j} + ϵ_{i j},

$Y_{ij} = a_i + b_i x_{ij}+ \epsilon_{ij},$

ϵ_{i j} \sim N (0, σ_{i}^{2})

$\epsilon_{ij}\sim N(0,\sigma_i^2)$

[Примечание: следующее на самом деле просто махнул рукой:]

Вы можете увидеть много общего между этими двумя моделями: соответствует а - . Среднее значение соответствует , потому что случайные эффекты усредняются до 0. Неограниченная корреляция случайного пересечения и наклона приводит к тому, что модели просто можно подобрать отдельно. Я не уверен, как единственное предположение согласуется с предметно- , но я бы предположил, что улавливает разницу. $a_i$ $a+\alpha_i$ $b_i$ $b+\beta_i$ $b_i$ $b$ $\sigma$ $\sigma_i$ $\alpha_i$

Анико
источник

Спасибо Анико. Вы правы, мои вычисления подтверждают это, но я хотел бы понять, почему ... Это кажется нелогичным. Я улучшил модели - введя случайные эффекты, я лучше описал структуру ошибок. Остаточная ошибка это подтверждает - все ниже и ниже. Так что с этими лучшими, более точными моделями я бы ожидал более точного наклона ... Я знаю, что где-то не так, пожалуйста, помогите мне увидеть это.

Любопытно

Спасибо Анико, это интересная точка зрения! Меня интересуют только уклоны (группа *: год), а не перехват здесь .. поэтому мой первый шаг по введению случайного эффекта itcept ослабил это предположение о независимости и привел к более низкой SE .. (уклона ..), а затем следующий шаг было, вероятно, слишком много (??) и сделал наоборот (еще хуже SE ..) .. может быть, мне нужно подумать об этом, спасибо.

любопытно

Теперь я также удивлен очень интересным фактом - пожалуйста, смотрите мое редактирование. Знаете ли вы, почему это?

любопытно

Я не думаю, что предположение о независимости было слишком ослаблено! Это было неправильно с самого начала.

Анико

Томас, «точная» модель не означает, что оценки будут более точными. В качестве крайнего примера возьмем любую модель без данных, которая вам нравится, например, модель, в которой все ответы равны нулю. Эта модель абсолютно уверена в своей оценке нуля. Поэтому он настолько точен, насколько это возможно, но, возможно, он также и неправдоподобен. Поэтому предоставление модели большего объема для соответствия параметрам обычно означает, что эти параметры соответствуют меньшей точности, а не большей. Лучшая модель, потому что она может количественно определить неопределенность, не отраженную в худшей модели, часто имеет большие стандартные ошибки.

whuber