Comment récupérer la matrice de corrélation de glm modèles R

Question

Je suis l'aide de la gls fonction du package nlme.Vous pouvez copier et coller le code suivant pour reproduire mon analyse.

library(nlme)  # Needed for gls function

# Read in wide format
tlc = read.table("http://www.hsph.harvard.edu/fitzmaur/ala2e/tlc.dat",header=FALSE)
names(tlc) = c("id","trt","y0","y1","y4","y6")
tlc$trt = factor(tlc$trt, levels=c("P","A"), labels=c("Placebo","Succimer"))

# Convert to long format
tlc.long = reshape(tlc, idvar="id", varying=c("y0","y1","y4","y6"), v.names="y", timevar="time", direction="long")

# Create week numerical variable
tlc.long$week = tlc.long$time-1
tlc.long$week[tlc.long$week==2] = 4
tlc.long$week[tlc.long$week==3] = 6

tlc.long$week.f = factor(tlc.long$week, levels=c(0,1,4,6))

L'analyse réelle commence ici:

# Including group main effect assuming unstructured covariance:
mod1 = gls(y ~ trt*week.f, corr=corSymm(, form= ~ time | id), 
       weights = varIdent(form = ~1 | time), method = "REML", data=tlc.long)
summary(mod1)

Dans le résumé(mod1), les parties suivantes des résultats sont d'intérêt pour moi que j'aimerais récupérer.

Correlation Structure: General
 Formula: ~time | id 
 Parameter estimate(s):
 Correlation: 
  1     2     3    
2 0.571            
3 0.570 0.775      
4 0.577 0.582 0.581
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | time 
 Parameter estimates:
       1        2        3        4 
1.000000 1.325880 1.370442 1.524813 

Le plus proche que je peux obtenir est d'utiliser la méthode suivante.

temp = mod1$modelStruct$varStruct
Variance function structure of class varIdent representing
       1        2        3        4 
1.000000 1.325880 1.370442 1.524813 

Cependant, ce que vous avez stocké avec le temp, je ne peux pas obtenir les cinq numéros.J'ai essayé aussi.numérique(temp) et sans classification(temp), mais aucun d'entre eux travaille.Il n'y a aucun moyen que je peux obtenir les cinq nombres propre numériques vectorielles.

Merci à l'avance!

Answer 1

Lorsque vous exécutez mod1$modelStruct$varStruct dans la R de la console, R inspecte d'abord la classe de ce

> class(mod1$modelStruct$varStruct)
[1] "varIdent" "varFunc" 

et puis l'envoi de l'correspondant print fonction.Dans ce cas, il est nlme:::print.varFunc.c'est à dire, la commande en cours d'exécution est nlme:::print.varFunc(mod1$modelStruct$varStruct).

Si vous exécutez nlme:::print.varFunc, vous pouvez voir le corps de la fonction de

function (x, ...) 
{
    if (length(aux <- coef(x, uncons = FALSE, allCoef = TRUE)) > 
        0) {
        cat("Variance function structure of class", class(x)[1], 
            "representing\n")
        print(aux, ...)
    }
    else {
        cat("Variance function structure of class", class(x)[1], 
            "with no parameters, or uninitialized\n")
    }
    invisible(x)
}
<bytecode: 0x7ff4bf688df0>
<environment: namespace:nlme>

Ce qu'il fait est de l'évaluation de la coef et l'imprimer, et le non évaluée x est retourné de manière invisible.

Par conséquent, afin d'obtenir le cor/var, vous avez besoin

coef(mod1$modelStruct$corStruct, uncons = FALSE, allCoef = TRUE)
coef(mod1$modelStruct$varStruct, uncons = FALSE, allCoef = TRUE)