Como recuperar a matriz de correlação de modelos glm em R

Question

Estou usando a função gls do pacote nlme.Você pode copiar e colar o código a seguir para reproduzir minha análise.

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))

A verdadeira análise começa aqui:

# 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)

No resumo (mod1), as seguintes partes dos resultados são de meu interesse e eu adoraria recuperar.

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 

O mais próximo que posso chegar é usar o seguinte método.

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

No entanto, seja o que for que você armazenou com temperatura, não consigo extrair os cinco números.Tentei as.numeric(temp) e unclass(temp), mas nenhum deles funciona.Não há como obter os cinco números como um vetor numérico limpo.

Desde já, obrigado!

Answer 1

Quando você corre mod1$modelStruct$varStruct no console R, R primeiro inspeciona sua classe

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

e então despachar o correspondente print função.Neste caso, é nlme:::print.varFunc.ou seja, o comando real em execução é nlme:::print.varFunc(mod1$modelStruct$varStruct).

Se você correr nlme:::print.varFunc, você pode ver o corpo da função dele

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>

O que ele faz é avaliar coef e imprimi-lo, e o não avaliado x é retornado de forma invisível.

Portanto, para obter o cor/var, você precisa

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