Pourquoi l'introduction d'un effet de pente aléatoire a élargi le SE de la pente?

J'essaie d'analyser l'effet de l'année sur logInd variable pour un groupe particulier d'individus (j'ai 3 groupes). Le modèle le plus simple:

> 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

Nous pouvons voir que le Groupe 1 est en baisse significative, les Groupes 2 et 3 en augmentation mais pas de manière significative.

De toute évidence, l'individu doit être un effet aléatoire, donc j'introduis un effet d'interception aléatoire pour chaque individu:

> 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

Il a eu un effet attendu - le SE des pentes (coefficients Groupe1-3: Année) est désormais plus faible et le SE résiduel est également plus faible.

Les individus sont également différents en pente, j'ai donc également introduit l'effet de pente aléatoire:

> 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

Mais maintenant, contrairement à ce que l'on attendait, les SE des pentes (coefficients Groupe1-3: Année) sont désormais bien plus élevés, voire plus que sans effet aléatoire du tout!

Comment est-ce possible? Je m'attendrais à ce que l'effet aléatoire «mange» la variabilité inexpliquée et augmente la «certitude» de l'estimation!

Cependant, le SE résiduel se comporte comme prévu - il est inférieur à celui du modèle d'interception aléatoire.

Voici les données si besoin.

Éditer

Maintenant, j'ai réalisé un fait étonnant. Si je fais la régression linéaire pour chaque individu séparément et que j'exécute ensuite l'ANOVA sur les pentes résultantes, j'obtiens exactement le même résultat que le modèle de pente aléatoire! Sais-tu pourquoi?

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

Voici les données si besoin.

r mixed-model lme4-nlme random-effects-model Curieuse
la source

Vous avez besoin d'un effet fixe d'une année si vous voulez avoir un effet fixe d'interaction année: groupe. En général, vous ne pouvez pas inclure un terme d'interaction sans inclure également les effets principaux. Pensez-vous vraiment qu'il n'y a pas de composante fixe à l'effet année? Et si oui, comment pourrait-il y avoir une année fixe: interaction de groupe?

John

Et pourquoi aucune interception fixe? Vous pouvez avoir les deux, fixes et aléatoires.

John

Group

i

$i$

i

$i$ Group

i

$i$ :Year

i

$i$

i

$i$

@John, c'est hors sujet à ma question, néanmoins: croyez-moi, c'est OK, j'ai fait beaucoup d'expériences avec ça. Mon premier modèle lm est entièrement équivalent à logInd ~ Year*Group, seuls les coefficients sont de forme différente, rien de plus. Cela dépend de votre goût et de la forme des coefficients que vous aimez, rien de plus. Il n'y a aucune exclusion de "Effet principal de l'année" dans mon 1er modèle pendant que vous écrivez ... logInd ~ Year*Groupfait exactement la même chose, le Yearcoefficient n'est donc pas l'effet principal, mais le Groupe1: Année.

Curieux

OK, soigné, n'avait pas considéré à la fois l'interception 0 et le groupe comme catégoriques.

John

Je pense que le problème est avec vos attentes :) Notez que lorsque vous avez ajouté une interception aléatoire pour chaque individu, l'erreur standard des interceptions a augmenté. Étant donné que chaque individu peut avoir sa propre interception, la moyenne du groupe est moins certaine. La même chose s'est produite avec la pente aléatoire: vous n'évaluez plus une pente commune (au sein d'un groupe), mais la moyenne de différentes pentes.

EDIT: Pourquoi un meilleur modèle ne donne-t-il pas une estimation plus précise?

Réfléchissons à l'inverse: pourquoi le modèle initial sous-estime-t-il l'erreur standard? Il suppose l'indépendance des observations qui ne sont pas indépendantes. Le deuxième modèle assouplit cette hypothèse (d'une manière qui affecte les interceptions), et le troisième l'assouplit davantage.

EDIT 2: relation avec de nombreux modèles spécifiques au patient

Votre observation est une propriété connue (et si vous n'aviez que deux ans, le modèle à effets aléatoires serait équivalent à un test t apparié). Je ne pense pas pouvoir gérer une vraie preuve, mais peut-être que l'écriture des deux modèles rendra la relation plus claire. Ignorons la variable de regroupement, car cela compliquerait simplement la notation. J'utiliserai des lettres grecques pour les effets aléatoires et des lettres latines pour les effets fixes.

$i$ $j$

{Oui}_{je j} = une + α_{je} + (b + β_{je}) X_{je j} + ϵ_{je 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)$

{Oui}_{je j} = {une}_{je} + b_{je} X_{je j} + ϵ_{je 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)$

[Remarque: ce qui suit est vraiment juste de la main:]

$a_i$ $a+\alpha_i$ $b_i$ $b+\beta_i$ $b_i$ $b$ $\sigma$ $\sigma_i$ $\alpha_i$

Aniko
la source

Merci Aniko. Vous avez raison, mes calculs le confirment, mais j'aimerais voir pourquoi ... Cela semble contre-intuitif. J'ai amélioré les modèles - en introduisant des effets aléatoires, j'ai mieux décrit la structure d'erreur. Une erreur résiduelle le confirme - elle est de plus en plus basse. Donc, avec ces modèles meilleurs et plus précis, je m'attendrais à une pente plus précise ... Je sais que je me trompe quelque part, aidez-moi à le voir.

Curieux

Merci Aniko, c'est un point de vue intéressant! Je ne suis intéressé que par les pentes (groupe *: année), pas par interception ici .. donc ma première étape d'introduction d'un effet itcept aléatoire a assoupli cette hypothèse d'indépendance et conduit à la SE inférieure .. (de pente ..) puis à l'étape suivante était probablement trop (??) et a fait le contraire (encore pire SE ..) .. peut-être que je dois y penser, merci.

Curieux

Maintenant, je suis également étonné par un fait très intéressant - veuillez consulter mon montage. Pourriez-vous savoir pourquoi?

Curieux

Je ne pense pas que l'hypothèse d'indépendance ait été trop assouplie! C'était mal au départ.

Aniko

Tomas, un modèle «précis» ne signifie pas que les estimations seront plus précises. À titre d'exemple extrême, prenez n'importe quel modèle sans données que vous aimez, comme celui qui prédit que toutes les réponses sont nulles. Ce modèle est absolument certain dans son estimation de zéro. Il est donc aussi précis que possible - mais il est probablement aussi erroné que possible. Donner à un modèle une plus grande portée pour ajuster les paramètres signifie donc généralement que ces paramètres sont ajustés avec moins de précision, pas plus. Un meilleur modèle, car il peut quantifier l'incertitude non captée par un modèle pire, a souvent des erreurs standard plus importantes.

whuber

Pourquoi l'introduction d'un effet de pente aléatoire a élargi le SE de la pente?

Mais maintenant, contrairement à ce que l'on attendait, les SE des pentes (coefficients Groupe1-3: Année) sont désormais bien plus élevés, voire plus que sans effet aléatoire du tout!

Éditer

Réponses: