我正在使用nlme软件包中的gls函数。您可以复制并粘贴以下代码以重现我的分析。
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))
真正的分析从这里开始：
# 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）中，以下结果部分对我来说很重要，我很乐意检索它们。
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
我能得到的最接近的方法是使用以下方法。
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），但是它们都不起作用。我无法将五个数字作为一个干净的数值向量。
提前致谢！
mod1$modelStruct$varStruct在R控制台中运行时，R首先检查它的类
> class(mod1$modelStruct$varStruct)
[1] "varIdent" "varFunc"
然后调度相应的print功能。在这种情况下，它是nlme:::print.varFunc。即，实际运行的命令是nlme:::print.varFunc(mod1$modelStruct$varStruct)。
如果运行nlme:::print.varFunc，则可以看到它的功能主体
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不可见地返回。
因此，为了获得cor / var，您需要
coef(mod1$modelStruct$corStruct, uncons = FALSE, allCoef = TRUE)
coef(mod1$modelStruct$varStruct, uncons = FALSE, allCoef = TRUE)
本文收集自互联网，转载请注明来源。
如有侵权，请联系[email protected] 删除。
我来说两句