2

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 を取得したいと考えています。

4

3 に答える 3

4

This is not that difficult; the VarCorr accessor method is designed precisely to recover this information. It's a little bit harder than it should be since the VarCorr method returns the variance-covariance as a character matrix rather than as numeric (I use storage.mode to convert to numeric without losing the structure, and suppressWarnings to ignore the warnings about NAs)

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

In your case, you would use vc["t","StdDev"].

于 2013-05-20T14:59:50.797 に答える
1

これは、印刷方法の 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
于 2013-05-20T14:00:50.733 に答える