R의 glm 모델에서 상관 행렬을 검색하는 방법
문제
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)
요약(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)를 시도했지만 어느 것도 작동하지 않습니다.5개의 숫자를 깨끗한 숫자 벡터로 얻을 수 있는 방법은 없습니다.
미리 감사드립니다!
해결책
당신이 달릴 때 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)
제휴하지 않습니다 StackOverflow