如何从R中的glm模型中检索相关矩阵 [英] How to retrieve correlation matrix from glm models in R

问题描述

我正在使用nlme软件包中的gls函数.您可以复制并粘贴以下代码以重现我的分析.

I am using the gls function from the nlme package. You can copy and paste the following code to reproduce my analysis.

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

真正的分析从这里开始:

The real analysis starts from here:

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

在summary(mod1)中，结果的以下部分对我很感兴趣，我很乐意检索.

In the summary(mod1), the following parts of the results are of interest to me that I would love to retrieve.

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 

我能得到的最接近的方法是使用以下方法.

The closest I can get is to use the following method.

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

但是，无论您存储的是什么温度，我都无法获得这五个数字.我尝试了as.numeric(temp)和unclass(temp)，但是它们都不起作用.我无法将五个数字作为一个干净的数值向量.

However, whatever you stored with temp, I cannot get the five numbers out. I tried as.numeric(temp) and unclass(temp), but none of them works. There is no way I can just get the five numbers as a clean numeric vector.

提前谢谢！

推荐答案

在R控制台中运行mod1$modelStruct$varStruct时，R首先检查它的类

When you run mod1$modelStruct$varStruct in R console, R first inspects the class of it

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

，然后分派相应的print函数.在这种情况下，它是nlme:::print.varFunc.也就是说，实际运行的命令是nlme:::print.varFunc(mod1$modelStruct$varStruct).

and then dispatch the corresponding print function. In this case, it is nlme:::print.varFunc. i.e., the actual command running is nlme:::print.varFunc(mod1$modelStruct$varStruct).

如果运行nlme:::print.varFunc，则可以看到它的功能主体

If you run nlme:::print.varFunc, you can see the function body of it

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>

它所做的是评估coef并将其打印出来，并且未评估的x会以不可见的方式返回.

What it does is evaluating the coef and print it, and the unevaluated x is returned invisibly.

因此，要获取cor/var，您需要

Therefore, in order to get the cor/var, you need

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

这篇关于如何从R中的glm模型中检索相关矩阵的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！