Covariance pattern models_1019
yuwenchen
10/21/2020
Data were mean electromyographic (EMG) amplitudes in the left brow region, measured over 90 seconds in response to different kind of music presented to 22 subjects. Measurements were taken while the subjects listened to each of four kinds of music: a relaxing piece, followed by pieces designed to elicit positive affect, agitation, and sadness.
Source: Vasey, M.W., & Thayer, J.F. (1987). The continuing problem of false positives in repeated measures ANOVA in psychophysiology: A multivariate solution. Psychophysiology, 24, 479-486.
pacman::p_load(tidyverse, nlme, car)
dta<- read.table("left_brow_emg.txt", h=T)
head(dta)## Patient Relax Positive Agitate Sad
## 1 S01 143 368 345 772
## 2 S02 142 155 161 178
## 3 S03 109 167 356 956
## 4 S04 123 135 137 187
## 5 S05 276 216 232 307
## 6 S06 235 386 398 425# means & variances
colMeans(dta[, -1])## Relax Positive Agitate Sad
## 217.6364 260.7273 394.3636 559.1364var(dta[, -1]) #variance matrice## Relax Positive Agitate Sad
## Relax 9929.385 8100.229 9589.853 8600.814
## Positive 8100.229 14768.113 13425.056 11809.420
## Agitate 9589.853 13425.056 45621.481 39267.710
## Sad 8600.814 11809.420 39267.710 170963.838cor(dta[, -1])## Relax Positive Agitate Sad
## Relax 1.0000000 0.6689192 0.4505739 0.2087499
## Positive 0.6689192 1.0000000 0.5172124 0.2350249
## Agitate 0.4505739 0.5172124 1.0000000 0.4446300
## Sad 0.2087499 0.2350249 0.4446300 1.0000000# wide format to long format
dtaL <- gather(dta, key = "Emotion", value = "EMG", 2:5) %>%
mutate(Emotion = as.factor(Emotion))
str(dtaL)## 'data.frame': 88 obs. of 3 variables:
## $ Patient: chr "S01" "S02" "S03" "S04" ...
## $ Emotion: Factor w/ 4 levels "Agitate","Positive",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ EMG : int 143 142 109 123 276 235 208 267 183 245 ...# set plot theme to black and white
# this removes the default gray background
ot <- theme_set(theme_bw())
ggplot(data = dtaL, aes(reorder(Emotion, EMG, mean), EMG)) +
geom_point(pch = 20, col="gray") +
stat_summary(fun.data = "mean_cl_boot", size = rel(1.1)) +
geom_hline(yintercept = mean(dtaL$EMG), lty = 2) +
coord_flip()+
labs(x = "Emotion", y = "Mean EMG Amplitude") 
univariate analysis
- aov should only be used for balanced data
- change the contrasts option to from default treatment type to sum type
- this is to ensure type III SSQ is used (contrast=list(Emotion=“contr.sum”))
summary(m0 <- aov(EMG ~ Emotion + Error(Patient/Emotion),
contrasts = list(Emotion = "contr.sum"), data = dtaL))##
## Error: Patient
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 21 2220062 105717
##
## Error: Patient:Emotion
## Df Sum Sq Mean Sq F value Pr(>F)
## Emotion 3 1560726 520242 11.51 4.1e-06 ***
## Residuals 63 2846877 45189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# same if you do not test for the patient effects here
summary(m0a <- lm(EMG ~ Patient + Emotion, data = dtaL))##
## Call:
## lm(formula = EMG ~ Patient + Emotion, data = dtaL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -366.51 -127.98 -33.65 103.93 749.08
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 443.40 113.30 3.913 0.000226 ***
## PatientS02 -248.00 150.31 -1.650 0.103945
## PatientS03 -10.00 150.31 -0.067 0.947169
## PatientS04 -261.50 150.31 -1.740 0.086797 .
## PatientS05 -149.25 150.31 -0.993 0.324546
## PatientS06 -46.00 150.31 -0.306 0.760594
## PatientS07 -186.25 150.31 -1.239 0.219915
## PatientS08 116.50 150.31 0.775 0.441213
## PatientS09 -54.50 150.31 -0.363 0.718136
## PatientS10 210.00 150.31 1.397 0.167292
## PatientS11 65.75 150.31 0.437 0.663304
## PatientS12 9.00 150.31 0.060 0.952445
## PatientS13 -258.25 150.31 -1.718 0.090695 .
## PatientS14 14.25 150.31 0.095 0.924773
## PatientS15 -317.50 150.31 -2.112 0.038636 *
## PatientS16 -154.25 150.31 -1.026 0.308728
## PatientS17 -55.25 150.31 -0.368 0.714430
## PatientS18 -194.50 150.31 -1.294 0.200403
## PatientS19 125.50 150.31 0.835 0.406918
## PatientS20 -86.00 150.31 -0.572 0.569266
## PatientS21 294.75 150.31 1.961 0.054315 .
## PatientS22 106.75 150.31 0.710 0.480213
## EmotionPositive -133.64 64.09 -2.085 0.041127 *
## EmotionRelax -176.73 64.09 -2.757 0.007616 **
## EmotionSad 164.77 64.09 2.571 0.012521 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 212.6 on 63 degrees of freedom
## Multiple R-squared: 0.5705, Adjusted R-squared: 0.4068
## F-statistic: 3.486 on 24 and 63 DF, p-value: 3.773e-05multivariate approach
# make ellipses in pairwise plot
scatterplotMatrix(~ Relax + Positive + Agitate + Sad, data = dta, smooth = F,
reg.line = F, ellipse = T, diag = "none", col = "darkgray") 
- repeated measures within-factor
- make sure you do not include the ID as part of the variables to be analyzed
emo <- as.factor(colnames(dta[, -1]))- run the usual linear model first and save the output object
- this is a one-sample design, Intercept is the only predictor
m1 <- lm(as.matrix(dta[, -1]) ~ 1)-idata, idesign inform Anova the repeated measures design
m1m <- Anova(m1, idata = data.frame(emo), idesign = ~ emo)# report only univariate output
summary(m1m, multivariate = F)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 11276284 1 2220062 21 106.665 1.098e-09 ***
## emo 1560726 3 2846877 63 11.513 4.104e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## emo 0.10299 1.6884e-08
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## emo 0.48021 0.0006406 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## emo 0.5060518 0.0004967266# report only multivariate output
summary(m1m, univariate = F)##
## Type III Repeated Measures MANOVA Tests:
##
## ------------------------------------------
##
## Term: (Intercept)
##
## Response transformation matrix:
## (Intercept)
## Relax 1
## Positive 1
## Agitate 1
## Sad 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 45105136
##
## Multivariate Tests: (Intercept)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.835506 106.6646 1 21 1.098e-09 ***
## Wilks 1 0.164494 106.6646 1 21 1.098e-09 ***
## Hotelling-Lawley 1 5.079265 106.6646 1 21 1.098e-09 ***
## Roy 1 5.079265 106.6646 1 21 1.098e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: emo
##
## Response transformation matrix:
## emo1 emo2 emo3
## Relax 0 0 1
## Positive 0 1 0
## Agitate 1 0 0
## Sad -1 -1 -1
##
## Sum of squares and products for the hypothesis:
## emo1 emo2 emo3
## emo1 597301.1 1081733 1237938
## emo2 1081733.0 1959056 2241948
## emo3 1237937.5 2241948 2565690
##
## Multivariate Tests: emo
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.5474483 7.661383 3 19 0.0014864 **
## Wilks 1 0.4525517 7.661383 3 19 0.0014864 **
## Hotelling-Lawley 1 1.2096921 7.661383 3 19 0.0014864 **
## Roy 1 1.2096921 7.661383 3 19 0.0014864 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mixed-effects analysis
covariance patterns - same as using lm except for estimation methods
summary(m0b <- gls(EMG ~ Patient + Emotion, data = dtaL))## Generalized least squares fit by REML
## Model: EMG ~ Patient + Emotion
## Data: dtaL
## AIC BIC logLik
## 944.4433 1000.165 -446.2216
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 443.3977 113.30335 3.913368 0.0002
## PatientS02 -248.0000 150.31388 -1.649881 0.1039
## PatientS03 -10.0000 150.31388 -0.066527 0.9472
## PatientS04 -261.5000 150.31388 -1.739693 0.0868
## PatientS05 -149.2500 150.31388 -0.992922 0.3245
## PatientS06 -46.0000 150.31388 -0.306026 0.7606
## PatientS07 -186.2500 150.31388 -1.239074 0.2199
## PatientS08 116.5000 150.31388 0.775045 0.4412
## PatientS09 -54.5000 150.31388 -0.362575 0.7181
## PatientS10 210.0000 150.31388 1.397077 0.1673
## PatientS11 65.7500 150.31388 0.437418 0.6633
## PatientS12 9.0000 150.31388 0.059875 0.9524
## PatientS13 -258.2500 150.31388 -1.718072 0.0907
## PatientS14 14.2500 150.31388 0.094802 0.9248
## PatientS15 -317.5000 150.31388 -2.112247 0.0386
## PatientS16 -154.2500 150.31388 -1.026186 0.3087
## PatientS17 -55.2500 150.31388 -0.367564 0.7144
## PatientS18 -194.5000 150.31388 -1.293959 0.2004
## PatientS19 125.5000 150.31388 0.834920 0.4069
## PatientS20 -86.0000 150.31388 -0.572136 0.5693
## PatientS21 294.7500 150.31388 1.960897 0.0543
## PatientS22 106.7500 150.31388 0.710181 0.4802
## EmotionPositive -133.6364 64.09405 -2.085004 0.0411
## EmotionRelax -176.7273 64.09405 -2.757312 0.0076
## EmotionSad 164.7727 64.09405 2.570796 0.0125
##
## Correlation:
## (Intr) PtnS02 PtnS03 PtnS04 PtnS05 PtnS06 PtnS07 PtnS08 PtnS09
## PatientS02 -0.663
## PatientS03 -0.663 0.500
## PatientS04 -0.663 0.500 0.500
## PatientS05 -0.663 0.500 0.500 0.500
## PatientS06 -0.663 0.500 0.500 0.500 0.500
## PatientS07 -0.663 0.500 0.500 0.500 0.500 0.500
## PatientS08 -0.663 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS09 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS10 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS11 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS12 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS13 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS14 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS15 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS16 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS17 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS18 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS19 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS20 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS21 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS22 -0.663 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## EmotionPositive -0.283 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## EmotionRelax -0.283 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## EmotionSad -0.283 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## PtnS10 PtnS11 PtnS12 PtnS13 PtnS14 PtnS15 PtnS16 PtnS17 PtnS18
## PatientS02
## PatientS03
## PatientS04
## PatientS05
## PatientS06
## PatientS07
## PatientS08
## PatientS09
## PatientS10
## PatientS11 0.500
## PatientS12 0.500 0.500
## PatientS13 0.500 0.500 0.500
## PatientS14 0.500 0.500 0.500 0.500
## PatientS15 0.500 0.500 0.500 0.500 0.500
## PatientS16 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS17 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS18 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS19 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS20 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS21 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## PatientS22 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## EmotionPositive 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## EmotionRelax 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## EmotionSad 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## PtnS19 PtnS20 PtnS21 PtnS22 EmtnPs EmtnRl
## PatientS02
## PatientS03
## PatientS04
## PatientS05
## PatientS06
## PatientS07
## PatientS08
## PatientS09
## PatientS10
## PatientS11
## PatientS12
## PatientS13
## PatientS14
## PatientS15
## PatientS16
## PatientS17
## PatientS18
## PatientS19
## PatientS20 0.500
## PatientS21 0.500 0.500
## PatientS22 0.500 0.500 0.500
## EmotionPositive 0.000 0.000 0.000 0.000
## EmotionRelax 0.000 0.000 0.000 0.000 0.500
## EmotionSad 0.000 0.000 0.000 0.000 0.500 0.500
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.7241433 -0.6020576 -0.1582857 0.4888895 3.5238212
##
## Residual standard error: 212.5759
## Degrees of freedom: 88 total; 63 residualunequal variances and same correlation(gls)
summary(m2 <- gls(EMG ~ Emotion,
weights = varIdent(form = ~ 1 | Emotion),
correlation = corCompSymm(form = ~ 1 | Patient),
data = dtaL))## Generalized least squares fit by REML
## Model: EMG ~ Emotion
## Data: dtaL
## AIC BIC logLik
## 1124.511 1146.388 -553.2555
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | Patient
## Parameter estimate(s):
## Rho
## 0.4245937
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Emotion
## Parameter estimates:
## Relax Positive Agitate Sad
## 1.000000 1.205723 2.124736 4.401306
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 394.3636 44.74068 8.814431 0.0000
## EmotionPositive -133.6364 41.00877 -3.258726 0.0016
## EmotionRelax -176.7273 40.55987 -4.357196 0.0000
## EmotionSad 164.7727 84.08261 1.959653 0.0534
##
## Correlation:
## (Intr) EmtnPs EmtnRl
## EmotionPositive -0.828
## EmotionRelax -0.883 0.809
## EmotionSad -0.064 0.220 0.196
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.4360744 -0.7807650 -0.2423782 0.6585354 3.1019030
##
## Residual standard error: 98.76635
## Degrees of freedom: 88 total; 84 residualunequal variances and unstructured correlations(gls)
summary(m3 <- gls(EMG ~ Emotion,
weights = varIdent(form = ~ 1 | Emotion),
correlation = corSymm(form = ~ 1 | Patient),
data = dtaL))## Generalized least squares fit by REML
## Model: EMG ~ Emotion
## Data: dtaL
## AIC BIC logLik
## 1127.745 1161.776 -549.8723
##
## Correlation Structure: General
## Formula: ~1 | Patient
## Parameter estimate(s):
## Correlation:
## 1 2 3
## 2 0.669
## 3 0.451 0.517
## 4 0.209 0.235 0.445
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Emotion
## Parameter estimates:
## Relax Positive Agitate Sad
## 1.000000 1.219558 2.143505 4.149460
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 394.3636 45.53790 8.660120 0.0000
## EmotionPositive -133.6364 39.04513 -3.422613 0.0010
## EmotionRelax -176.7273 40.66006 -4.346459 0.0000
## EmotionSad 164.7727 79.21513 2.080066 0.0406
##
## Correlation:
## (Intr) EmtnPs EmtnRl
## EmotionPositive -0.823
## EmotionRelax -0.885 0.879
## EmotionSad -0.080 0.070 0.076
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.4109335 -0.7650935 -0.2381349 0.6510508 3.0745201
##
## Residual standard error: 99.646
## Degrees of freedom: 88 total; 84 residualunequal variances and unstructured correlations (lme)
summary(m4 <- lme(EMG ~ Emotion, random = ~ 1 | Patient,
weights = varIdent(form = ~ 1 | Emotion),
correlation = corCompSymm(form = ~ 1 | Patient),
data = dtaL))## Linear mixed-effects model fit by REML
## Data: dtaL
## AIC BIC logLik
## 1121.332 1145.64 -550.6659
##
## Random effects:
## Formula: ~1 | Patient
## (Intercept) Residual
## StdDev: 83.51979 57.91493
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | Patient
## Parameter estimate(s):
## Rho
## 0.3390305
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Emotion
## Parameter estimates:
## Relax Positive Agitate Sad
## 1.000000 1.588484 3.404611 7.209578
## Fixed effects: EMG ~ Emotion
## Value Std.Error DF t-value p-value
## (Intercept) 394.3636 45.65416 63 8.638066 0.0000
## EmotionPositive -133.6364 39.91051 63 -3.348400 0.0014
## EmotionRelax -176.7273 39.59459 63 -4.463420 0.0000
## EmotionSad 164.7727 84.58341 63 1.948050 0.0559
## Correlation:
## (Intr) EmtnPs EmtnRl
## EmotionPositive -0.816
## EmotionRelax -0.880 0.882
## EmotionSad -0.129 0.240 0.208
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.1316206 -0.6892089 -0.2082434 0.4716320 2.5200430
##
## Number of Observations: 88
## Number of Groups: 22comparing models
anova(m2, m3, m4)## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 9 1124.511 1146.388 -553.2555
## m3 2 14 1127.745 1161.776 -549.8723 1 vs 2 6.766257 0.2386
## m4 3 10 1121.332 1145.640 -550.6659 2 vs 3 1.587102 0.8111# show estimated covariance
intervals(m4, which= "var-cov")## Approximate 95% confidence intervals
##
## Random Effects:
## Level: Patient
## lower est. upper
## sd((Intercept)) 50.67523 83.51979 137.6522
##
## Correlation structure:
## lower est. upper
## Rho 0.004190219 0.3390305 0.6710206
## attr(,"label")
## [1] "Correlation structure:"
##
## Variance function:
## lower est. upper
## Positive 0.6064687 1.588484 4.160614
## Agitate 1.4249696 3.404611 8.134474
## Sad 3.0221699 7.209578 17.198904
## attr(,"label")
## [1] "Variance function:"
##
## Within-group standard error:
## lower est. upper
## 24.65004 57.91493 136.07034residual plots
plot(m0b, resid(., type = "pearson") ~ fitted(.) | Emotion,
layout = c(4, 1), aspect = 2, abline = 0, pch = 20, cex = .8,
ylim = c(-2, 3.3),
xlab = "Ftted values", ylab = "Pearson residuals", main = "Linear model")
plot(m3, resid(., type = "pearson") ~ fitted(.) | Emotion,
layout = c(4, 1), aspect = 2, abline = 0, pch = 20, cex = .8,
ylim = c(-2, 3.3),
xlab = "Ftted values", ylab = "Pearson residuals", main = "Covariance patterns")
plot(m4, resid(., type = "pearson") ~ fitted(.) | Emotion,
layout = c(4, 1), aspect = 2, abline = 0, pch = 20, cex = .8,
ylim = c(-2, 3.3),
xlab = "Ftted values", ylab = "Pearson residuals",
main = "Covariance patterns + random effects")
dta_m4 <- data.frame(dtaL, M4 = fitted(m4), M0 = fitted(m0b))
dtaL_m4 <- dta_m4 %>%
unite(PE, Patient, Emotion) %>%
gather(key = Model, value = Estimate, 2:4) %>%
separate(PE, c("Patient", "Emotion"))ggplot(dtaL_m4, aes(reorder(Emotion, Estimate, mean), Estimate, color = Model)) +
stat_summary(fun.data = mean_se, geom = "errorbar",
width = .2, position = position_dodge(.2)) +
coord_flip()+
labs(x = "Emotion", y = "Mean EMG") +
theme(legend.position = c(.9, .2))