2013-05-20 2 views
2

nlme 패키지의 lme 모델에서 임의의 용어의 분산을 구할 수있는 방법이 있습니까?nlme의 lme에서 무작위 효과 분산 추정에 액세스

Random effects: 
Formula: ~t | UID 
Structure: General positive-definite, Log-Cholesky parametrization 
      StdDev  Corr 
(Intercept) 520.310397 (Intr) 
t    3.468834 0.273 
Residual  31.071987 

다른 말로하면, 나는 3.468834에 도착하고 싶습니다.

+0

예, 그렇지만 어느 요소에 있어야하며 어떤 숫자 여야합니까? 나는 lme 객체의 구조에서 길을 잃는다. – drw

답변

4

이것은 하지 어려운 그; VarCorr 접근 자 메서드는이 정보를 정확하게 복구하도록 설계되었습니다. VarCorr 메서드가 숫자가 아닌 문자 매트릭스로 분산 공분산을 반환하기 때문에 구조가 손실되지 않고 숫자로 변환하려면 storage.mode을 사용하고 구조에 대해서는 경고를 무시하기 위해 suppressWarnings을 사용하기 때문에 조금 더 어려웠습니다.

library(nlme) 
fit <- lme(distance ~ Sex, data = Orthodont, random = ~ age|Subject) 
vc <- VarCorr(fit) 
suppressWarnings(storage.mode(vc) <- "numeric") 
vc[1:2,"StdDev"] 
## (Intercept)   age 
## 7.3913363 0.6942889 

귀하의 경우에는 vc["t","StdDev"]을 사용하십시오.

+0

+1 내 시도보다 훨씬 튼튼하다. 'VarCorr'을 몰랐습니다. – Roland

1

이것은 인쇄 방법 중 하나에서 계산됩니다 (나는 print.summary.pdMat으로 추정 됨). 가장 쉬운 방법은 출력을 캡처하는 것입니다.

library(nlme) 

fit <- lme(distance ~ Sex, data = Orthodont, random = ~ age|Subject) 
summary(fit) 

# Linear mixed-effects model fit by REML 
# Data: Orthodont 
# AIC  BIC logLik 
# 483.1635 499.1442 -235.5818 
# 
# Random effects: 
# Formula: ~age | Subject 
# Structure: General positive-definite, Log-Cholesky parametrization 
#    StdDev Corr 
# (Intercept) 7.3913363 (Intr) 
# age   0.6942889 -0.97 
# Residual 1.3100396 
# <snip/> 

ttt <- capture.output(print(summary(fit$modelStruct), sigma = fit$sigma)) 
as.numeric(unlist(strsplit(ttt[[6]],"\\s+"))[[2]]) 
#[1] 0.6942889 

그리고 여기에 그것을 계산하는 방법입니다 :

fit$sigma * attr(corMatrix(fit$modelStruct[[1]])[[1]],"stdDev") 
#(Intercept)   age 
# 7.3913363 0.6942889 
1
> fit <- lme(distance ~ Sex, data = Orthodont, random = ~ age|Subject) 
> getVarCov(fit) 
Random effects variance covariance matrix 
      (Intercept)  age 
(Intercept)  54.6320 -4.97540 
age    -4.9754 0.48204 
    Standard Deviations: 7.3913 0.69429 
> # In contrast to VarCorr(), this returns a numeric matrix: 
> str(getVarCov(fit)) 
random.effects [1:2, 1:2] 54.632 -4.975 -4.975 0.482 
- attr(*, "dimnames")=List of 2 
    ..$ : chr [1:2] "(Intercept)" "age" 
    ..$ : chr [1:2] "(Intercept)" "age" 
- attr(*, "class")= chr [1:2] "random.effects" "VarCov" 
- attr(*, "group.levels")= chr "Subject" 
> unclass(getVarCov(fit)) 
      (Intercept)  age 
(Intercept) 54.631852 -4.975417 
age   -4.975417 0.482037 
attr(,"group.levels") 
[1] "Subject" 
관련 문제