Mix-effect model
bigcat
April 6, 2017
1. Introduction
+++++++++++++++++++++++++++++++++++++++++
+ Mixed ANOVA for Stress Hormones
+ Response: Stress hormones
+ Factor1: Time (Before, After)
+ Factor2: Sex (Male, Female)
+ Factor 3: Mibyeong2
+++++++++++++++++++++++++++++++++++++++++2. Data visualization
str(dat_cv)## 'data.frame': 316 obs. of 5 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Sex : Factor w/ 2 levels "Men","Women": 1 1 1 2 2 2 2 1 2 1 ...
## $ Mibyeong2: Factor w/ 2 levels "Health2","MB2": 2 2 2 2 2 2 2 2 2 2 ...
## $ Time : Factor w/ 2 levels "Epi.P._B","Epi.P._A": 1 1 1 1 1 1 1 1 1 1 ...
## $ value : num 67.6 133.6 128.8 16 99 ...Hormone_before_and_after
3. Function for mix-effect model by ezANOVA
(using F test)
library(ez)
dat_cv$ID=factor(dat_cv$ID) #ezANOVA requires ID as factor
Mix_ezANOVA=ezANOVA(data = dat_cv, dv = .(value), wid = .(ID),
between = .(Sex,Mibyeong2), within = .(Time), type = 3, detailed = F)## Warning: Data is unbalanced (unequal N per group). Make sure you specified
## a well-considered value for the type argument to ezANOVA().Mix_ezANOVA## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 Sex 1 154 25.38746012 1.299137e-06 * 0.0963139837
## 3 Mibyeong2 1 154 1.44311552 2.314801e-01 0.0060218573
## 5 Time 1 154 118.30006378 8.421180e-21 * 0.2135564122
## 4 Sex:Mibyeong2 1 154 0.58727291 4.446487e-01 0.0024593656
## 6 Sex:Time 1 154 10.53030882 1.440616e-03 * 0.0236008987
## 7 Mibyeong2:Time 1 154 13.00132944 4.199237e-04 * 0.0289785498
## 8 Sex:Mibyeong2:Time 1 154 0.08026927 7.773133e-01 0.0001842169Remark:
- Similar results as performed by GLM
4. Function for mix-effect model by GLM
(using Chisq test for Likelyhood Ratio)
library(nlme)
# Option 1
baseline=lme(value~1,random=~1|ID/Sex/Mibyeong2/Time,data=dat_cv,method="ML") # no factor
TimeM=update(baseline,.~.+Time) # add main effect of Time
MibyeongM=update(TimeM,.~.+Mibyeong2) # add main effect of MB
SexM=update(MibyeongM,.~.+Sex) # add main effect of Sex
Time_Sex=update(SexM,.~.+Time:Sex) # add interaction effect of Time and Sex
MB_Sex=update(Time_Sex,.~.+Mibyeong2:Sex) # add interaction effect of MB and Sex
Time_MB=update(MB_Sex,.~.+Time:Mibyeong2) # add interaction effect of MB and Time
full=update(Time_MB,.~.+Time:Mibyeong2:Sex) # model with all main and interaction effect
anova(baseline,TimeM,MibyeongM,SexM,Time_Sex,MB_Sex,Time_MB,full)## Model df AIC BIC logLik Test L.Ratio p-value
## baseline 1 6 3242.627 3265.162 -1615.314
## TimeM 2 7 3176.720 3203.010 -1581.360 1 vs 2 67.90758 <.0001
## MibyeongM 3 8 3176.959 3207.005 -1580.480 2 vs 3 1.76015 0.1846
## SexM 4 9 3149.888 3183.690 -1565.944 3 vs 4 29.07105 <.0001
## Time_Sex 5 10 3141.407 3178.964 -1560.704 4 vs 5 10.48141 0.0012
## MB_Sex 6 11 3142.806 3184.119 -1560.403 5 vs 6 0.60138 0.4381
## Time_MB 7 12 3131.951 3177.020 -1553.976 6 vs 7 12.85420 0.0003
## full 8 13 3133.869 3182.694 -1553.935 7 vs 8 0.08233 0.7742Remark:
There were significant main effects of the Time (Chisq(2) = 67.9, p < .0001) and Sex (Chisq(2) = 29.07, p < .0001), whereas no main effect of MB significant.
There were significant interaction effects between Time and Sex and between Time and MB.
However, no significant in full model that means interaction between Time and Sex was not different in MB group, and interaction between Time and MB was not different betwen men and women.
#Option 2
full=lme(value~Sex*Mibyeong2*Time,random=~1|ID/Sex/Mibyeong2/Time,data=dat_cv,method="ML")
anova(full)## numDF denDF F-value p-value
## (Intercept) 1 154 691.8156 <.0001
## Sex 1 154 31.4559 <.0001
## Mibyeong2 1 154 1.8537 0.1753
## Time 1 154 95.8997 <.0001
## Sex:Mibyeong2 1 154 0.5873 0.4446
## Sex:Time 1 154 11.4638 0.0009
## Mibyeong2:Time 1 154 13.0593 0.0004
## Sex:Mibyeong2:Time 1 154 0.0803 0.7773Remark (Similar results as in Option 1)
#Summary model full
full=lme(value~Sex*Mibyeong2*Time,random=~1|ID/Sex/Mibyeong2/Time,data=dat_cv,method="ML")
summary(full)## Linear mixed-effects model fit by maximum likelihood
## Data: dat_cv
## AIC BIC logLik
## 3133.869 3182.694 -1553.935
##
## Random effects:
## Formula: ~1 | ID
## (Intercept)
## StdDev: 10.56906
##
## Formula: ~1 | Sex %in% ID
## (Intercept)
## StdDev: 10.56926
##
## Formula: ~1 | Mibyeong2 %in% Sex %in% ID
## (Intercept)
## StdDev: 10.56929
##
## Formula: ~1 | Time %in% Mibyeong2 %in% Sex %in% ID
## (Intercept) Residual
## StdDev: 27.49355 7.259093
##
## Fixed effects: value ~ Sex * Mibyeong2 * Time
## Value Std.Error DF t-value
## (Intercept) 42.45455 7.303239 154 5.813112
## SexWomen -7.84809 9.549323 154 -0.821848
## Mibyeong2MB2 11.51688 9.015314 154 1.277480
## TimeEpi.P._A 62.84091 8.684322 154 7.236133
## SexWomen:Mibyeong2MB2 -9.20905 11.736860 154 -0.784626
## SexWomen:TimeEpi.P._A -24.62155 11.355153 154 -2.168316
## Mibyeong2MB2:TimeEpi.P._A -27.13853 10.720159 154 -2.531541
## SexWomen:Mibyeong2MB2:TimeEpi.P._A 3.95409 13.956365 154 0.283318
## p-value
## (Intercept) 0.0000
## SexWomen 0.4124
## Mibyeong2MB2 0.2034
## TimeEpi.P._A 0.0000
## SexWomen:Mibyeong2MB2 0.4339
## SexWomen:TimeEpi.P._A 0.0317
## Mibyeong2MB2:TimeEpi.P._A 0.0124
## SexWomen:Mibyeong2MB2:TimeEpi.P._A 0.7773
## Correlation:
## (Intr) SexWmn Mb2MB2 TE.P._ SxW:M2MB2
## SexWomen -0.765
## Mibyeong2MB2 -0.810 0.620
## TimeEpi.P._A -0.595 0.455 0.482
## SexWomen:Mibyeong2MB2 0.622 -0.814 -0.768 -0.370
## SexWomen:TimeEpi.P._A 0.455 -0.595 -0.368 -0.765 0.484
## Mibyeong2MB2:TimeEpi.P._A 0.482 -0.368 -0.595 -0.810 0.457
## SexWomen:Mibyeong2MB2:TimeEpi.P._A -0.370 0.484 0.457 0.622 -0.595
## SW:TE. M2MB2:
## SexWomen
## Mibyeong2MB2
## TimeEpi.P._A
## SexWomen:Mibyeong2MB2
## SexWomen:TimeEpi.P._A
## Mibyeong2MB2:TimeEpi.P._A 0.620
## SexWomen:Mibyeong2MB2:TimeEpi.P._A -0.814 -0.768
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.63931548 -0.10864864 -0.02450976 0.05184494 1.97669971
##
## Number of Observations: 316
## Number of Groups:
## ID Sex %in% ID
## 158 158
## Mibyeong2 %in% Sex %in% ID Time %in% Mibyeong2 %in% Sex %in% ID
## 158 316Remark: Interpret summary(full) carefully
Ref: Discovering Statistic using R