혼합 효과 모델에서 공분산 행렬을 지정하기위한 pdBlocked 구문 nlme

Question

혼합 효과 모델이 있고 내 상관 관계를 내 임의 효과 공분산 행렬에 떨어 뜨려 내 dof를 줄입니다. 이렇게하려면 pdBlocked을 사용해야한다고 생각하지만 올바른 구문을 가져 와서 내가 원하는 것을 얻을 수는 없습니다.

예제 코드 : 다음과 같은 공분산 행렬을 제공

library(nlme) 
m3 <- lme(distance ~ age +I(age^2) + I(age^3), data = Orthodont, 
      random = list(Subject = pdBlocked(list(~ age,~0 + I(age^2),~0+I(age^3))))) 

:

getVarCov(m3) 
Random effects variance covariance matrix 
      (Intercept)  age   I(age^2)      I(age^3) 
(Intercept)  5.2217 -0.30418 0.00000000000000 0.00000000000000000000000000 
age    -0.3042 0.04974 0.00000000000000 0.00000000000000000000000000 
I(age^2)   0.0000 0.00000 0.00000000003593 0.00000000000000000000000000 
I(age^3)   0.0000 0.00000 0.00000000000000 0.00000000000000000000002277 
    Standard Deviations: 2.285 0.223 0.000005994 0.000000000004772 

이것은 내가 원하는 가까이는 아니지만 꽤. I(age^3)과 intercept 사이의 상관 관계를 유지하고 싶습니다. age은 0이지만, I(age^2)과의 상관 관계를 허용합니다. 이런 식으로 뭔가 :이 scenrio

getVarCov(m3) 
    Random effects variance covariance matrix 
       (Intercept)  age   I(age^2)      I(age^3) 
    (Intercept)  5.2217 -0.30418 c_value   b_value   
    age    -0.3042 0.04974 d_value   0.00000000000000000000000000 
    I(age^2)   c_value d_value 0.00000000003593 a_value 
    I(age^3)   b_value 0.00000 a_value   0.00000000000000000000002277 
     Standard Deviations: 2.285 0.223 0.000005994 0.000000000004772 

에도

getVarCov(m3) 
Random effects variance covariance matrix 
      (Intercept)  age   I(age^2)      I(age^3) 
(Intercept)  5.2217 -0.30418 0.00000000000000 0.00000000000000000000000000 
age    -0.3042 0.04974 0.00000000000000 0.00000000000000000000000000 
I(age^2)   0.0000 0.00000 0.00000000003593 a_value 
I(age^3)   0.0000 0.00000 a_value   0.00000000000000000000002277 
    Standard Deviations: 2.285 0.223 0.000005994 0.000000000004772 

난 그냥 사람이 제로가되는 선택을 할 수있는 유연한 공분산 행렬을 만드는 방법을 모르겠어요. 이러한 링크는 매우 도움이되었다하지만 여전히 http://rpsychologist.com/r-guide-longitudinal-lme-lmer

https://stats.stackexchange.com/questions/87050/dropping-term-for-correlation-between-random-effects-in-lme-and-interpretting-su?rq=1

감사 어떤 도움을 정확히 을 알아낼 기운 다. 감사합니다

Answer 1

age^2과 age^3 용어를 한 마디로 모으는 것처럼 보입니다.

m4 <- lme(distance ~ age +I(age^2) + I(age^3), data = Orthodont, 
      random = list(Subject = pdBlocked(list(~ age, 
               ~0 + I(age^2)+I(age^3)))), 
      control=lmeControl(opt="optim")) 
getVarCov(m4) 
## Random effects variance covariance matrix 
##    (Intercept)  age I(age^2) I(age^3) 
## (Intercept)  5.00960 -0.225450 0.0000e+00 0.0000e+00 
## age   -0.22545 0.019481 0.0000e+00 0.0000e+00 
## I(age^2)  0.00000 0.000000 4.1676e-04 -1.5164e-05 
## I(age^3)  0.00000 0.000000 -1.5164e-05 5.5376e-07 
## Standard Deviations: 2.2382 0.13957 0.020415 0.00074415 

내가 두 번째 예 ( age 및 age^3 상관, 다른 모든 상관 관계 제로가 아닌) pdBlocked와는 구성 할 수있는 방법이 있다고 생각하지 않습니다

- 용어의 순서를 재 배열 할 수있는 방법이 없습니다 (행/행렬)의 행렬을 사용하여이 행렬이 블록 대각선이되도록합니다. 당신은 원칙적으로에서 자신의 pdMatrix 클래스를 쓸 수 있지만이 모듈 식 디자인을 가지고 lme4에서이 작업을 수행하는 방법을 알아 내기 위해 시작

그 ... 슈퍼 쉽지 않을 것이라고이 작업을 수행 할 것 약간이 더 쉽지만 모델에 추가 문제가 있음을 발견했습니다. 이 데이터 세트에 대해 과다하게 판별됩니다 (실제 데이터 세트인지 여부는 알 수 없음). Orthodont 데이터 세트에는 피험자 당 4 개의 관측치 만 있기 때문에 개인당 4 개의 임의 효과 값 (절편과 3 개의 다항식 값)을 맞추면 임의 효과 효과 편차가 잔차 편차와 혼동되는 모델을 얻을 수 있습니다 이 모델들로부터). 시도 할 경우 lme4 오류가 발생합니다.

그러나 여전히 이것을 원한다면 오류를 무시할 수 있습니다 (위험! 로빈슨!). 먼저 선형 대수학을 수행하고 하위 삼각형의 콜레 스키 팩터 [lme4의 매개 변수화 이 theta[2]*theta[4]+theta[5]*theta[7]과 같음을 확인하기 위해 분산 공분산 행렬을 사용합니다. 여기에서 theta은 콜레 스키 요인 요소의 벡터입니다 (하위 삼각형, 열 우선 순서). 따라서 우리는이 -theta[2]*theta[4]/theta[5]과 같게 설정된 전체 10 개 매개 변수 모델 대신 9 개 매개 변수 모델을 피팅하여이 작업을 수행 할 수 있습니다.

lf <- lFormula(distance ~ age +I(age^2) + I(age^3) + 
       (age+ I(age^2) + I(age^3)|Subject), data = Orthodont, 
       control=lmerControl(check.nobs.vs.nRE="ignore")) 
devfun <- do.call(mkLmerDevfun,lf) 
trans_theta <- function(theta) 
    c(theta[1:6],-theta[2]*theta[4]/theta[5],theta[7:9]) 
devfun2 <- function(theta) { 
    return(devfun(trans_theta(theta))) 
} 
diagval <- (lf$reTrms$lower==0) 
opt <- minqa::bobyqa(fn=devfun2,par=ifelse(diagval,1,0)[-7], 
      lower=lf$reTrms$lower[-7]) 
opt$par <- trans_theta(opt$par) 
opt$conv <- 0 
m1 <- mkMerMod(environment(devfun), opt, lf$reTrms, fr = lf$fr) 
VarCorr(m1) 

그러나이 작업을 수행하는 것이 타당한지를 신중히 생각해 보는 것이 좋습니다. 나는 당신이 실제로 이런 식으로 조건을 떨어 뜨리는 것으로부터 정밀도/힘의 측면에서 많이 얻지는 않을 것이라고 생각한다. (일반적으로, 이후의 모델에서의 가설의 이득은 모델 단순화이다. Harrell 회귀 모델링 전략)이 공분산 구조를 예상 할 수있는 기계 론적 또는 주제 기반의 이유가 없다면 나는 신경 쓰지 않을 것이라고 생각합니다.