앞에서 질문했던 데이터에 대한 추가 질문입니다.

[박상일]

1)앞서 질문 드렸던것처럼

본 실험이 2개의 그룹(FUT/HPR 그룹당 4마리)각 개체마다 4가지 test에서 6가지 검사 항목을 측정하는데 1시간 간격으로 6회 측정합니다.

데이터는 longitudinal data 형식으로 변환하였으며

앞서 질문드렸던 것 처럼 각 test당 각 검사 항목의 변이(분산)을 그룹간 차이가 있는지를 test 하려고 합니다.

var.test 함수를이용해서(안재형 선생님의 책을 참고하였습니다)

test 하려고 우선은 첫 1시간째 검사 결과를 그룹간 차이가 있는지 test 하려고 하였습니다.

종속 변수는 6가지 검사 항목이고

여기서 각 시간마다 4가지 test가 있기 때문에(intem,extem,heptem, fibtem)

우선은 각 시간대 별로 데이터 프레임을 만들었습니다.

data1=ecmo[ecmo$times==1,]

data1

이렇게 6개를 만들었고

각 데이터 프레임에서 각각의 test에 대한 var.test를 하려고 하였는데

var.test(alpha~group,data1[data1$test==intem,])를 실행하면

다음에 오류`[.data.frame`(data1, data1$test == intem, ) : 

  개체 'intem'이 없습니다

이렇게 출력이 됩니다.

무슨 문제가 있을까요?

2) mixed model을 만들어 보았습니다. random intercept만 있는 모델(mod4hg)보다는 random slope이 있는 모델(mod3hg)가 유의한 차이가 있는것으로 나와 mod3hg를 선택하였습니다.

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

> # fit Model 3 

> mod3hg = lmer(hg~group+as.numeric(times)+(as.numeric(times)|ID), ecmo) 

> # fit Model 4 

> mod4hg = lmer(hg~group+as.numeric(times)+(1|ID), ecmo) 

> # model comparison 

> anova(mod3hg, mod4hg)

Data: ecmo

Models:

mod4hg: hg ~ group + as.numeric(times) + (1 | ID)

mod3hg: hg ~ group + as.numeric(times) + (as.numeric(times) | ID)

       Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)    

mod4hg  5 265.10 281.17 -127.55                             

mod3hg  7 222.47 244.98 -104.24 46.626      2  7.503e-11 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

mod3hg를 요약하면

> summary(mod3hg)

Linear mixed model fit by REML 

Formula: hg ~ group + as.numeric(times) + (as.numeric(times) | ID) 

   Data: ecmo 

   AIC   BIC logLik deviance REMLdev

 227.6 250.1 -106.8    208.5   213.6

Random effects:

 Groups   Name              Variance Std.Dev. Corr  

 ID       (Intercept)       0.288567 0.53718        

          as.numeric(times) 0.031485 0.17744  0.027 

 Residual                   0.138543 0.37221        

Number of obs: 184, groups: ID, 8

Fixed effects:

                  Estimate Std. Error t value

(Intercept)       10.62589    0.27687   38.38

groupHPR          -0.70019    0.39186   -1.79

as.numeric(times) -0.53685    0.06502   -8.26

Correlation of Fixed Effects:

            (Intr) grpHPR

groupHPR    -0.706       

as.nmrc(tm) -0.017 -0.004

해석은 FUT 그룹의 시간에 따른 측정치(hemoglobin 농도)는 0.53685씩 감소하며

HPR 그룹의 그것은 FUT 그룹보다 0.70019 더 감소한다. 라고 결론을 내려고 합니다.

맞는 결론 일까요?(여기서 염려가 되는 점은 t value가 -1.79로 너무 작다는 것입니다. 이 값이 의미가 있으려면 단측검정이라 하더라도 자유도가 12 이상이어야 하는데...ㅠㅠ)

긴 글 읽어 주셔서 감사합니다.

참고로 헤모글로빈 농도에 대한 그래프를 보여 드리면

[박상일]

아직도 mixed model 을 충분히 이해하지 못하고 있어 질문을 드리면,랜덤효과를 모델에 추가한다는 것은 랜덤효과에 의한 분산팽창효과를 제거하고 고정효과에 의한 측정변수의 차이를 비교한다는 것인가요?일요일 하루를 연구실에서 그생각만 하다가 이제야 집에 왔네요^^

[안재형]

1) interm대신 "interm"이라고 ""를 해줘야할겁니다.
2) 이건 여기서 설명하기 힘든데요. time과 group의 interaction이 유의한가 봐야할겁니다. 그림에서 두 group의 slop이 차이가 있는지 검정하는건데 londitudinal study에서는 시간이 지남에 따라 어느 그룹의 효과가 더 있는지 보기때문에 slop의 차이를 봅니다. groupHPR은 단순히 두 그룹의 시작점의 차이입니다. 즉 baseline의 차이가 있는지 보는거여서 별 의미가 없습니다.

[박상일]

매번 감사합니다. 좀더 공부를.....

[박상일]

맞는 개념인지는 몰라도 혼합모형이라는게 결국은 회귀분석에서의 계수(여기서는 고정효과죠?)의 유의성을 찾는 것인데 모델이라는게 오차항이 있기 마련이어서 이 오차항을 분해하여 랜덤 절편/기울기로 해석을해 결국에는 마직막에 남은 오차항을 최소화하는데,다변량을 다루기 때문에 분산/공분산 행렬의 구조를 따지고, 모델(?)의 적합성 비교를 위해REMLE/MLE를 사용한다. 모델간 비교를 위해서 AIC/BIC(아직 이건모름), log lik ratio를 구하고 카이제곱분포를 이용해 통계량을 구한다. 개념은 어느정도 잡은건가요?

[안재형]

회귀계수, 분산등 모든 모수의 추정은 ML/REML로 합니다. 모든 모수가 추정된 최종 ML/REML 값은 모형이 커지면 증가합니다. AIC, BIC는 ML/REML 값에 모수의 수를 penalty로 준 값으로, 모수의 수가 증가할때 AIC, BIC가 무조건 증가하지않아 모형간 비교가 가능합니다.