problème avec lme, comment utiliser le summary

Questions sur les fonctions statistiques de R

Modérateur : Groupe des modérateurs

Gaëlle Nogues
Messages : 2
Enregistré le : 14 Juil 2006, 11:29

problème avec lme, comment utiliser le summary

Messagepar Gaëlle Nogues » 27 Juil 2006, 15:15

Bonjour,

C'est mon premier petit projet en R... Et j'ai deja des soucis...

J'ai des données sur 31 individus qui sont mesurées 7 fois dans le temps. J'ai ajusté un modèle avec la fonction "lme" à mes données, et quand je trace un graphique avec les valeurs ajustées par rapport a ma réponse, ça a l'air juste. Je cherche à retrouver la fonction correspondant a mon ajustement, et pour ça j'ai demandé un "summary" de mon modèle :

Code : Tout sélectionner

Linear mixed-effects model fit by maximum likelihood
 Data: glucose2
        AIC      BIC   logLik
  -763.5905 -687.423 400.7952

Random effects:
 Formula: ~TimePoint | SubjectID
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev      Corr 
(Intercept) 0.085272828 (Intr)
TimePoint   0.001960009 -0.694

 Formula: ~TimePoint | Period %in% SubjectID
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev      Corr 
(Intercept) 0.116220940 (Intr)
TimePoint   0.001457469 -0.843
Residual    0.066005142       

Fixed effects: log(BloodGlucose) ~ poly(TimePoint, 2) * Grade
                                             Value     Std.Error  DF   t-value p-value
(Intercept)                               5.879948 0.0318297 341 184.73177  0.0000
poly(TimePoint, 2)1                  -0.753267 0.4282114 341  -1.75910  0.0795
poly(TimePoint, 2)2                  -3.949686 0.1361768 341 -29.00410  0.0000
GradeDiabGrade 1                       0.260718 0.0445524  27   5.85194  0.0000
GradeDiabGrade 2                      0.253268 0.0445524  27   5.68473  0.0000
GradeDiabGrade 3                      0.252461 0.0456082  27   5.53544  0.0000
poly(TimePoint, 2)1:GradeDiabGrade 1  3.165915 0.6015212 341   5.26318  0.0000
poly(TimePoint, 2)2:GradeDiabGrade 1  1.368242 0.1895393 341   7.21878  0.0000
poly(TimePoint, 2)1:GradeDiabGrade 2  4.095052 0.6015212 341   6.80783  0.0000
poly(TimePoint, 2)2:GradeDiabGrade 2  1.292958 0.1895393 341   6.82158  0.0000
poly(TimePoint, 2)1:GradeDiabGrade 3  4.121715 0.6179009 341   6.67051  0.0000
poly(TimePoint, 2)2:GradeDiabGrade 3  1.472031 0.1927865 341   7.63555  0.0000


Je voudrais tracer la courbe correspondant a mon modèle, donc :

y(t)= [b0 + b01 grade_1 + b02 grade_2 + b03 grade_3]
+ [b1 + b11 grade_1 + b12 grade_2 + b13 grade_3] t
+ [b2 + b21 grade_1 + b22 grade_2 + b23 grade_3] t²

avec grade_i =0 si l'individu appartient au groupe i, 1 sinon,
b0 <- 5.889098
b1 <- -0.790891
b2 <- -3.949542
b01 <- 0.258151
b02 <- 0.243972
b03 <- 0.243311
b11 <- 3.168266
b21 <- 1.368098
b12 <- 4.130165
b22 <- 1.292814
b13 <- 4.159338
b23 <- 1.471887

je trouve quelque chose qui n'a rien a voir avec mes données. Apparement, les échelles ne sont pas bonnes du tout (il faudrait diviser par environ 2 pour l'ordonnée et environ 60 pour l'abscisse..)...

Est-ce que quelqu'un comprend mon problème et si oui, pourrait m'aider ?

merci
Gaëlle

Pierre Bady
Messages : 405
Enregistré le : 02 Mai 2006, 07:46

Messagepar Pierre Bady » 28 Juil 2006, 12:40

Bonjour,

Pour répondre à ta question, je vais utiliser l’exemple qui est utilisé dans l’aide de la fonction lme.


Code : Tout sélectionner

> require(nlme)
[1] TRUE
> ?lme
> Orthodont[1:3,]
Grouped Data: distance ~ age | Subject
  distance age Subject  Sex
1       26   8     M01 Male
2       25  10     M01 Male
3       29  12     M01 Male
...


108 relevés pour pour 27 individus différents
on construit le modèle suivant :

distance = a + b*age + A(subject) + B(subject)*age + erreur

avec :
partie fixe : intercept (a) et age (b)
partie aléatoire : intercept (A(subject)) et age (B(subject))

on construit le modèle dans R avec la fonction lme :

Code : Tout sélectionner

> fm1 <- lme(distance ~ age, data = Orthodont) # random is ~ age

# on utilise le summary de lme
> summary(fm1)
Linear mixed-effects model fit by REML
 Data: Orthodont
       AIC      BIC    logLik
  454.6367 470.6173 -221.3183

Random effects:
 Formula: ~age | Subject
 Structure: General positive-definite
            StdDev    Corr 
(Intercept) 2.3270339 (Intr)
age         0.2264276 -0.609
Residual    1.3100399       

Fixed effects: distance ~ age
                Value Std.Error DF   t-value p-value
(Intercept) 16.761111 0.7752461 80 21.620375       0
age          0.660185 0.0712533 80  9.265334       0
 Correlation:
    (Intr)
age -0.848

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max
-3.223106023 -0.493760866  0.007316632  0.472151091  3.916032740

Number of Observations: 108
Number of Groups: 27
>


les coefficients sont contenus dans l'objet "coefficient" du modèle:

Code : Tout sélectionner

> fm1$coefficient
$fixed
(Intercept)         age
 16.7611111   0.6601852

$random
$random$Subject
    (Intercept)          age
M16  -0.1877576 -0.068853674
M05  -1.1766673  0.025600294
M02  -0.7275013  0.014507807
M11   0.8904889 -0.118825808
M07  -0.6079718  0.034899975
M08   0.8602964 -0.094736144
M03  -0.1739021  0.035852320
M12  -0.9979934  0.114563960
M13  -4.1295406  0.413668190
M14   0.9043445 -0.014119814
M09  -0.4443942  0.135908474
M15  -0.5349719  0.208177467
M06   1.2176440  0.083191189
M04   3.0004525 -0.065884754
M01   1.0515851  0.215684343 .. .


Pour obtenir les valeurs « fittées » (ajustées,prédites, .. .):

Pour calculer les valeurs prédites fixées, c’est comme d’habitude :
modèle : distance = a + b*age

Code : Tout sélectionner

> 0.6601852 *Orthodont$age[1] +16.7611111
[1] 22.04259


si on rajoute l’effet aléatoire :
modèle : distance = a + b*age + A(subject) + B(subject)*age

Code : Tout sélectionner

> 16.7611111 + 0.6601852 *Orthodont$age[1] + 1.0515851  + 0.215684343* Orthodont$age[1]
[1] 24.81965

Il faut faire attention à l’ordre des individus et pas se tromper de sujet ! .. . je l'ai fait pour mon premier essai et j’ai cherché mon erreur pendant un p’tit moment :)

ces valeurs sont aussi contenues dans l’objet « fitted » du modèle:

Code : Tout sélectionner

> fm1$fitted[1:3,]
     fixed  Subject
1 22.04259 24.81965
2 23.36296 26.57139
3 24.68333 28.32313


et on peut également les extraire directement avec la fonction « fitted » (fitted.lme) :

Code : Tout sélectionner

> ?fitted.lme
> # fixe + aléatoire
> fitted(fm1,level=1)[1:3]
     M01      M01      M01
24.81965 26.57139 28.32313
> # fixe
> fitted(fm1,level=0)[1:3]
     M01      M01      M01
22.04259 23.36296 24.68333



ensuite pour le graphique ... on ne ce casse pas la tête et on utilise l’objet fitted du modèle ou la fonction "fitted" ... :)

Code : Tout sélectionner

> plot(Orthodont$age,Orthodont$distance,ylim=c(16,35),xlim=c(6,15),pch=19,cex=1.5)
> points(Orthodont$age,fm1$fitted[,2],pch=18,cex=1.5,col="green")
> points(Orthodont$age,fm1$fitted[,1],pch=17,cex=1.5,col="orange")
>
> legend(7,35,legend=c("observed","aléatoire+fixe","fixe"),col=c("black","green","orange"),
+   lty=0,pch=c(19,18,17),pt.cex=1.5)
>



Si tu as plusieurs niveaux, ... eh bien ... c'est toujours le même principe.


en espérant avoir aidé un peu :)

@+++

pierre
=@===--------¬-------¬------¬-----¬
liens utiles :
http://www.gnurou.org/Writing/SmartQuestionsFr
http://neogrifter.free.fr/welcomeOnInternet.jpg
]<((((*< -------------------------------

Gaëlle Nogues
Messages : 2
Enregistré le : 14 Juil 2006, 11:29

Messagepar Gaëlle Nogues » 31 Juil 2006, 07:53

Merci pour cette aide !

Une grande partie fonctionne bien en effet, je n'ai aucun problème avec le graphique, qui me donne exactement ce que j'attend. Mais si je cherche à calculer les valeurs prédites en me servant de mon modèle, j'obtiens quelque chose de totalement différent... :

Code : Tout sélectionner

> glucose2[1:3,]
Grouped Data: BloodGlucose ~ TimePoint | SubjectID/Period
  Period SubjectID       Grade TimePoint log(BloodGlucose)
1      1       101   DiabGrade 1         5         5.55
2      1       101   DiabGrade 1        10         5.84
3      1       101   DiabGrade 1        15         5.90


2 périodes, 31 individus différents évalués 7fois sur chaque période.


modèle :

effets fixes :

log(BloodGlucose) = a + b*TimePoint + c*TimePoint²+ A(Grade) + B(Grade)*TimePoint + C(Grade)*TimePoint² + erreur

effets fixes + aléatoires :

log(BloodGlucose) = a + b*TimePoint + c*TimePoint²+ A(Grade) + B(Grade)*TimePoint + C(Grade)*TimePoint² + D(Subject/Period)+E(Subject/Period)*TimePoint + erreur


Code : Tout sélectionner

> gluc.lme3 <- lme(fixed=log(BloodGlucose) ~ poly(TimePoint,2) *Grade, random= ~ TimePoint | SubjectID/Period, data = glucose2, method="ML", correlation=corCompSymm(form=~TimePoint|SubjectID/Period))
> summary(gluc.lme3)
Linear mixed-effects model fit by maximum likelihood
 Data: glucose2
        AIC       BIC   logLik
  -761.5905 -681.4142 400.7952

Random effects:
 Formula: ~TimePoint | SubjectID
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev      Corr 
(Intercept) 0.085272850 (Intr)
TimePoint   0.001960010 -0.694

 Formula: ~TimePoint | Period %in% SubjectID
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev      Corr 
(Intercept) 0.116509124 (Intr)
TimePoint   0.001457468 -0.841
Residual    0.065495100       

Correlation Structure: Compound symmetry
 Formula: ~TimePoint | SubjectID/Period
 Parameter estimate(s):
        Rho
-0.01563569
Fixed effects: log(BloodGlucose) ~ poly(TimePoint, 2) * Grade
                                         Value Std.Error  DF   t-value p-value
(Intercept)                           5.879948 0.0318297 341 184.73176  0.0000
poly(TimePoint, 2)1                  -0.753267 0.4282114 341  -1.75910  0.0795
poly(TimePoint, 2)2                  -3.949686 0.1361768 341 -29.00410  0.0000
GradeDiabGrade 1                      0.260718 0.0445524  27   5.85194  0.0000
GradeDiabGrade 2                      0.253268 0.0445524  27   5.68473  0.0000
GradeDiabGrade 3                      0.252461 0.0456082  27   5.53544  0.0000
poly(TimePoint, 2)1:GradeDiabGrade 1  3.165915 0.6015212 341   5.26318  0.0000
poly(TimePoint, 2)2:GradeDiabGrade 1  1.368242 0.1895393 341   7.21878  0.0000
poly(TimePoint, 2)1:GradeDiabGrade 2  4.095052 0.6015212 341   6.80783  0.0000
poly(TimePoint, 2)2:GradeDiabGrade 2  1.292958 0.1895393 341   6.82158  0.0000
poly(TimePoint, 2)1:GradeDiabGrade 3  4.121715 0.6179009 341   6.67051  0.0000
poly(TimePoint, 2)2:GradeDiabGrade 3  1.472031 0.1927865 341   7.63555  0.0000     
...


Code : Tout sélectionner

> gluc.lme3$coefficient$fixed
                         (Intercept)                  poly(TimePoint, 2)1
                           5.8799478                           -0.7532672
                 poly(TimePoint, 2)2                     GradeDiabGrade 1
                          -3.9496861                            0.2607178
                    GradeDiabGrade 2                     GradeDiabGrade 3
                           0.2532682                            0.2524612
poly(TimePoint, 2)1:GradeDiabGrade 1 poly(TimePoint, 2)2:GradeDiabGrade 1
                           3.1659145                            1.3682416
poly(TimePoint, 2)1:GradeDiabGrade 2 poly(TimePoint, 2)2:GradeDiabGrade 2
                           4.0950520                            1.2929576
poly(TimePoint, 2)1:GradeDiabGrade 3 poly(TimePoint, 2)2:GradeDiabGrade 3
                           4.1217146                            1.4720312

> gluc.lme3$coefficient$random$Period
        (Intercept)     TimePoint
101/1 -1.636298e-01  9.457504e-04
101/2  5.684887e-02  1.909495e-05
102/1  1.601810e-01 -1.217918e-03
103/1 -1.732456e-01  1.623104e-03
103/2  1.753465e-01 -1.738190e-03



mais pourtant, si on calcule les valeurs ajustées pour l'individu 101, lors de la période 1, au temps 5 :

Code : Tout sélectionner

> 5.8799478 + 0.2607178+ (-0.7532672 + 3.1659145)*5+ (-3.9496861 + 1.3682416)*5^2 -1.636298e-01+9.457504e-04*5
[1] -46.49111


et :

Code : Tout sélectionner

> gluc.lme3$fitted[1:3,]
     fixed SubjectID   Period
1 5.853959  5.808029 5.649128
2 5.966183  5.921487 5.767315
3 6.064377  6.020915 5.871471

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 04 Aoû 2006, 10:24

Bonjour,

Sans avoir eu le temps de lire l'exemple de façon détaillée, il faut savoir que les objects stockés dans les résultats de la fct lme sont souvent sous une forme particulière (log, ou divisé par l'écart-type résiduel, etc). Il est donc prudent, voire indispensable, d'utiliser les fonctions pré-programmées pour extraire les coefficients et calculer les valeurs prédites (ou autres) plutôt que de manipuler les objets "bruts". Voir en particulier les fonctions:
* fixef (coefficients des effets fixes),
* ranef (coefficients des effets aléatoires),
* fitted (valeurs ajustées),
* predict (valeurs prédites), etc.

Pour cette dernière (predict), regarder en détail les arguments newdata (permet de prédire à des valeurs particulières des variables explicatives) et level (permet de prédire en tenant compte ou pas des effets aléatoires).

Bonne chance,

Renaud


Retourner vers « Archives : Fonctions statistiques »

Qui est en ligne

Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité

cron