Modérateur : Groupe des modérateurs
Tillard a écrit :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
Tillard a écrit :est -il possible d'obtenir - de maniere simple - l'erreur standard des estimateurs des effets aleatoires (variance, covariance) apres une utilisation de lme ?
Code : Tout sélectionner
> fm <- lme(distance ~ age + Sex, data = Orthodont, random = ~ age)
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
Code : Tout sélectionner
> nu <- log(c(1.3626967, 2.7966864, 5.7396886))
> nu
[1] 0.3094656 1.0284353 1.7474050
Code : Tout sélectionner
> c(nu[2] - nu[1], nu[3] - nu[2])
[1] 0.7189697 0.7189697
Code : Tout sélectionner
> ((nu[2] - nu[1]) / 1.96)^2
[1] 0.1345578
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
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
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
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
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
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]
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 Lesnoff a écrit :Tillard a écrit :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
Tillard a écrit : j'essaierai de faire une comparaison avec MlWin et/ou stata des que possible
Retourner vers « Archives : Fonctions statistiques »
Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité