Erreurs standards des estimateurs des effets aléatoires

[Tillard]

est -il possible d'obtenir - de maniere simple - l'erreur standard des estimateurs des effets aleatoires (variance, covariance) apres une utilisation de lme ou glmmPQL, sans passer par la fonction intervals?

la fonction intervals.lme permet d'obtenir un intervalle de confiance des effets aleatoires a partir d'un objet crée avec lme ou glmmPQL; mais ne donne pas directement l'erreur standard (comme le font d'autres packages comme stata ou mlwin - il y a d'ailleurs peut etre une raison a cela, en rapport avec la difficulte de tester directement des effets aleatoires via le test de Wald).

J'en profite pour vous adresser mes meilleurs voeux pour 2005 ain si qu'une longue vie au forum.

[Matthieu Lesnoff]

est -il possible d'obtenir - de maniere simple - l'erreur standard des estimateurs des effets aleatoires (variance, covariance) apres une utilisation de lme ou glmmPQL, sans passer par la fonction intervals?

la fonction intervals.lme permet d'obtenir un intervalle de confiance des effets aleatoires a partir d'un objet crée avec lme ou glmmPQL; mais ne donne pas directement l'erreur standard

Tu cherches donc les variances-covariances des composantes de la variance. En utilisant "search" du site R, j'ai retrouve un email de Bates (2003) :

"Other software for fitting mixed-effects models, such as SAS PROC
MIXED and HLM, return standard errors along with the estimates of the
variances and covariances of the random effects. We don't return
standard errors of estimated variances because we don't think they are
useful. A standard error for a parameter estimate is most useful when
the distribution of the estimator is approximately symmetric, and
these are not. "

Donc pas evident qu'il existe qq chose sous lme. En remarque, les intervalles de confiance des composantes de la variance donnés par intervals.lme sont aussi informatifs (si ce n'est plus) que les variances, car les lois des composantes sont connues pour être dissymétriques en général.

Matthieu

[Matthieu Lesnoff]

est -il possible d'obtenir - de maniere simple - l'erreur standard des estimateurs des effets aleatoires (variance, covariance) apres une utilisation de lme ?

Je ne pense pas qu'il existe de fonction permettant de récupérer **directement** les variances-covariances des composantes de la variance d'un modèle mixte ajusté avec lme. Par contre, voici une approche qui permet de les calculer assez facilement (si je ne me suis pas planté !).

En fait, d'après le bouquin de Pinheiro et BAtes (2002, chap. 2), lme n'estime pas directement les composantes de la variance, mais passe par des transformées de ces composantes. Les estimations de ces transformées (et non les composantes elle-mêmes), utilisant la méthode ML ou REML, sont supposées asymptotiquement normales avec comme matrice de variances-covariances l'inverse de la matrice d'information de Fischer issue du processus d'optimisation.

Prenons l'exemple du modèle à 2 effets aléatoires suivant (exemple de l'aide lme) :

Code : Tout sélectionner

> fm <- lme(distance ~ age + Sex, data = Orthodont, random = ~ age)

La fonction intervals.lme donne les IC des effets fixes, des ecart-types des effets aléatoires, des coefficients de corrélation entre les effets aléatoires, et enfin de sigma (sigma = Within-group standard error) :

Code : Tout sélectionner

> intervals(fm)
Approximate 95% confidence intervals

 Fixed effects:
                 lower       est.      upper
(Intercept) 15.8715571 17.6351776 19.3987981
age             0.5183895  0.6601852  0.8019808
SexFemale  -3.7053901 -2.1454359 -0.5854816
attr(,"label")
[1] "Fixed effects:"

 Random Effects:
  Level: Subject
                                    lower       est.      upper
sd((Intercept))          1.3626967  2.7966864  5.7396886
sd(age)                    0.1023954  0.2264063  0.5006066
cor((Intercept),age) -0.9542751 -0.7658093 -0.1417022

 Within-group standard error:
   lower     est.    upper
1.084774 1.310072 1.582162

D'après les résultats ci-dessus, les estimations des deux composantes du modèle sont :

sd_1 = 2.7966864,
sd_2 = 0.2264063
(attention, intervals.lme fournit des écart-types et non des variances)

et leur corrélation est :
cor[sd_1, cor_sd2] = -0.7658093.

L'estimation des IC de sd_1, sd_2 et cor[sd_1, cor_sd2] (= tableau "Random effect" ci-dessus) passe par l'estimation d'une transformation log pour les sd_1 et sd_2 et d'une transformation logit généralisée pour cor[sd_1, cor_sd2] (tout est expliqué p. 93 dans le bouquin de P&B). Les paramètres transformés sont appelés "nu", et ce sont eux qui sont estimés, ainsi que leur variance. Les "nu" sont supposés asymptotiquement normaux, leurs IC sont donc symétriques, du type nu +- 1.96 s.e.(nu) pour un risque 5%. Les IC de sd_1, sd_2 et cor[sd_1, cor_sd2] sont calculés a posteriori à partir de ceux des nu
en prenant la transformée réciproque (= exp(x) pour sd_1 et sd_2, = exp(x - 1) / (exp x + 1) pour cor[sd_1, cor_sd2]).

Par exemple, l'estimation et l'IC de nu_sd1 est :

Code : Tout sélectionner

> nu <- log(c(1.3626967, 2.7966864, 5.7396886))
> nu
[1] 0.3094656 1.0284353 1.7474050

On peut vérifier que cet intervalle est bien symétrique :

Code : Tout sélectionner

> c(nu[2] - nu[1], nu[3] - nu[2])
[1] 0.7189697 0.7189697


La variance Var[nu_sd1] est :

Code : Tout sélectionner

> ((nu[2] - nu[1]) / 1.96)^2
[1] 0.1345578


De la même manière, l'estimation, l'IC et la variance de nu_sd2 sont :

Code : Tout sélectionner

> nu <- log(c(0.1023954, 0.2264063, 0.5006066)) # estimation et IC de nu_sd2
> nu
[1] -2.2789135 -1.4854241 -0.6919347
> ((nu[2] - nu[1]) / 1.96)^2 # variance de nu_sd2
[1] 0.1638967

et pour nu_cor[sd_1, cor_sd2] :

Code : Tout sélectionner

> rho <- c(-0.9542751, -0.7658093, -0.1417022)
> nu <- log((1 + rho) / (1 - rho)) # estimation et IC de nu_cor[sd_1, cor_sd2]
> nu
[1] -3.7551316 -2.0202287 -0.2853245
> ((nu[2] - nu[1]) / 1.96)^2 # variance de nu_cor[sd_1, cor_sd2]
[1] 0.7834986


En résumé :
Var[nu_sd_1] = 0.1345578,
Var[nu_sd_2] = 0.1638967,
Var[nu_cor[sd_1, cor_sd2]] = 0.7834986.

L'objet "apVAr" donne directment ces valeurs (et en + les covariances !) :

Code : Tout sélectionner

> fm$apVar
                  reStruct.Subject1 reStruct.Subject2 reStruct.Subject3       lSigma
reStruct.Subject1        0.13456279         0.1197552       -0.28755048 -0.012850260
reStruct.Subject2        0.11975517         0.1639028       -0.29949463 -0.015558300
reStruct.Subject3       -0.28755048        -0.2994946        0.78352838  0.028812250
lSigma                      -0.01285026        -0.0155583        0.02881225  0.009270341

En fait, fm$apVar est la matrice de var-cov des estimations des transformées de tous ce qui n'est pas "effet fixe". Ici il s'agit des transformées log pour sd_1, sd_2 et sigma, et de la transformée logistique generalisée pour cor[sd_1, cor_sd2]. Cela varie bien sûr en focntion des modèles spécifiés.

Si on ne s'interesse qu'aux composantes sd_1 et sd_2, la matrice de var-cov du vecteur [nu_sd_1 nu_sd_2]' est :

Code : Tout sélectionner

> V <- fm$apVar[1:2, 1:2]
> V
                  reStruct.Subject1 reStruct.Subject2
reStruct.Subject1         0.1343769         0.1195592
reStruct.Subject2         0.1195592         0.1636969

A partir de cette matrice V, on peut calculer la matrice de var-cov du vecteur [sd_1 sd_2]' (ce que l'on recherche) en utilisant la "Delta method".
Pour "x" un paramètre unidimensionnel, cette methode est basée sur la formule d'approximation suivante :

Var[f(x)] = Var[x] * (df(x) / dx)^2.

Ici, les formules d'approximation sont donc :

Var[sd_i] = Var[nu_i] * exp(2 * nu_i)
Sd[sd_i] = Sd[nu_i] * exp(nu_i)
Var[var_i] = 4 * Var[nu_i] * exp(4 * nu_i)
Sd[var_i] = 2 * Sd[nu_i] * exp(2 * nu_i)

Cela donne :

Code : Tout sélectionner

> diag(varnu) * exp(2 * nu) # Var[sd_i]
reStruct.Subject1 reStruct.Subject2
      1.051022944       0.008391072

> diag(varnu)^0.5 * exp(nu)  # Sd[sd_i]
reStruct.Subject1 reStruct.Subject2
        1.0251941         0.0916028

> 4 * diag(varnu) * exp(4 * nu) # Var[var_i]
reStruct.Subject1 reStruct.Subject2
     32.882112978       0.001720498

> 2 * diag(varnu)^0.5 * exp(2 * nu) # Sd[var_i]
reStruct.Subject1 reStruct.Subject2
       5.73429272        0.04147889


J'ai vérifié cette approche sur un autre exemple, détaillé dans la doc SAS pour les modèles mixtes (Littel et al. 1996, p.236 Output 6.1). SAS fournit Sd[var_i]. Avec la formule "Sd[var_i] = 2 * Sd[nu_i] * exp(2 * nu_i)", j'obtiens bien les mêmes resultats que ceux donnés par SAS. J'essaierai (...) de faire une fiche à partir de cet exemple SAS en utilisant lme.

De manière +compacte, l'approche se résume à quelques lignes :

Code : Tout sélectionner

> fm <- lme(distance ~ age + Sex, data = Orthodont, random = ~ age)
> nu <- log(intervals(fm)$reStruct$Subject[1:2, 2])
> varnu <- fm$apVar[1:2, 1:2]
> diag(varnu) * exp(2 * nu) # Var[sd_i]
> diag(varnu)^0.5 * exp(nu)  # Sd[sd_i]
> 4 * diag(varnu) * exp(4 * nu) # Var[var_i]
> 2 * diag(varnu)^0.5 * exp(2 * nu) # Sd[var_i]


Il existe sûrement d'autres approches ...

En remarque, la dernière ligne du print(apVar) corresponds aux nu :

Code : Tout sélectionner

> fm$apVar
                  reStruct.Subject1 reStruct.Subject2 reStruct.Subject3       lSigma
reStruct.Subject1        0.13437692        0.11955915       -0.28709610 -0.012819791
reStruct.Subject2        0.11955915        0.16369695       -0.29901893 -0.015526035
reStruct.Subject3       -0.28709610       -0.29901893        0.78240885  0.028745251
lSigma                  -0.01281979       -0.01552603        0.02874525  0.009262611
attr(,"Pars")
reStruct.Subject1 reStruct.Subject2 reStruct.Subject3            lSigma
        1.0284353        -1.4854243        -2.0202285         0.2700821
attr(,"natural")
[1] TRUE

> z <- log(c(2.7966864, 0.2264063)) # vecteur[nu_sd_1 nu_sd_2]
> z
[1]  1.028435 -1.485424


Matthieu

[Tillard]

est -il possible d'obtenir - de maniere simple - l'erreur standard des estimateurs des effets aleatoires (variance, covariance) apres une utilisation de lme ?

Je ne pense pas qu'il existe de fonction permettant de récupérer **directement** les variances-covariances des composantes de la variance d'un modèle mixte ajusté avec lme. Par contre, voici une approche qui permet de les calculer assez facilement (si je ne me suis pas planté !).

...
De manière +compacte, l'approche se résume à quelques lignes :

> fm <- lme(distance ~ age + Sex, data = Orthodont, random = ~ age)
> nu <- log(intervals(fm)$reStruct$Subject[1:2, 2])
> library(MASS)
> res <- exp(mvrnorm(n = 50000, mu = nu, Sigma = fm$apVar[1:2, 1:2]))
> cov(res)
reStruct.Subject1 reStruct.Subject2
reStruct.Subject1 1.29699987 0.09367221
reStruct.Subject2 0.09367221 0.01074878

Il existe sûrement d'autres approches ... Ca vaut peut-être le coup de poser la question sur le forum officiel R (en précisant que ton pb est de comparer les resultats avec MlWin qui fournit les variances des composantes).

Matthieu

Merci matthieu pour cette reponse, précise et claire: j'essaierai de faire une comparaison avec MlWin et/ou stata des que possible

[Matthieu Lesnoff]

j'essaierai de faire une comparaison avec MlWin et/ou stata des que possible

J'ai modifié la fin de ma reponse précédente (cf+haut). J'ai utilisé la "delta method" à la place du Monte Carlo. Avec cette méthode, on obtient les mêmes résultats de variance des composantes que ceux donnés par SAS.

Matthieu