# Load libraries
if (!require(simr)) {install.packages('simr')}
## Warning: package 'simr' was built under R version 3.4.4
## Warning: package 'lme4' was built under R version 3.4.4
if (!require(lmerTest)) {devtools::install_github('cran/lmerTest')}
if (!require(r2glmm)) {install.packages('r2glmm')}
if (!require(effects)) {install.packages('effects')}
## Warning: package 'carData' was built under R version 3.4.4
# Set seed and define contrasts
set.seed(123)
options(contrasts=c("contr.sum", "contr.poly"))
# Read in the data
data <- read.csv("simdat.csv")
# Scale CS
data$stCs_tot <- as.numeric(scale(data$Cs_tot))
# For the toy dataset, do listwise deletion, subset the Partner pic condition,
# & compute a dummy for sad partner vs the rest
data <- data[!is.na(data$FacExp) & !is.na(data$stCs_tot) & !is.na(data$duration),]
Pdata <- data[data$Person == "P",]
data$Psad.vs.not <- ifelse(data$Person == "P" & data$FacExp == "S", 1, 0) #not used so far
# Null (without the primary fixed effect of Face expression) and full model
# Including CS and duration as covariates, with participants treated as a random factor
model0 <- lmer(Temperature ~ stCs_tot + duration + (1|PPNR), data = Pdata)
model <- lmer(Temperature ~ FacExp + stCs_tot + duration + (1|PPNR), data = Pdata)
#model
summary(model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Temperature ~ FacExp + stCs_tot + duration + (1 | PPNR)
## Data: Pdata
##
## REML criterion at convergence: 13612.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5907 -0.5047 0.0600 0.5257 3.5892
##
## Random effects:
## Groups Name Variance Std.Dev.
## PPNR (Intercept) 13.9602 3.7363
## Residual 0.4334 0.6583
## Number of obs: 6510, groups: PPNR, 70
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.878e+01 4.470e-01 6.820e+01 64.395 <2e-16 ***
## FacExp1 -2.092e-02 9.078e-03 6.438e+03 -2.304 0.0213 *
## stCs_tot -3.843e-01 4.467e-01 6.800e+01 -0.860 0.3926
## duration 2.943e-03 5.847e-05 6.438e+03 50.340 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) FcExp1 stCs_t
## FacExp1 -0.011
## stCs_tot 0.000 0.000
## duration -0.038 0.371 0.000
# Show and plot treatment contrast plot
allEffects(model)
## model: Temperature ~ FacExp + stCs_tot + duration
##
## FacExp effect
## FacExp
## N S
## 29.64400 29.68584
##
## stCs_tot effect
## stCs_tot
## -2 -1 -0.03 1 2
## 30.43747 30.05322 29.68049 29.28470 28.90045
##
## duration effect
## duration
## 90 200 300 400 500
## 29.05188 29.37562 29.66993 29.96424 30.25855
plot(allEffects(model)[1])
# Set the smallest effect size of interest (SESOI)
# Please note that the model term represents deviation contrast.
# To get simple treatment contrast (difference btw neutral and sad, multiply by 2)
SESOI <- -.025 #SESOI set at a 1/20 degree of C difference between sad and neutral
fixef(model)["FacExp1"] <- SESOI
# Compute R^2 for the model terms (after changing to SESOI)
r2beta(model)
## Effect Rsq upper.CL lower.CL
## 1 Model 0.427 0.443 0.411
## 4 duration 0.265 0.283 0.248
## 3 stCs_tot 0.238 0.255 0.221
## 2 FacExp1 0.001 0.003 0.000
# As a benchmark, approximate Cohen's d for the effect of FacExp
r <- sqrt(r2beta(model)$Rsq[4])
paste("d = ", round((2*r)/sqrt(1 - r^2), 3))
## [1] "d = 0.067"
# POWER - Main effect
# Power to detect SESOI with 93 temperature measurements (~15 per picture) and 70 participants (original experiment)
pwr1 <- powerSim(model, test = compare(model0, method = "lr"), nsim = 300, progress = F)
pwr1
## Power for model comparison, (95% confidence interval):
## 73.67% (68.30, 78.56)
##
## Test: Likelihood ratio
## Comparison to Temperature ~ stCs_tot + duration + (1 | PPNR)
##
## Based on 300 simulations, (0 warnings, 0 errors)
## alpha = 0.05, nrow = 6510
##
## Time elapsed: 0 h 3 m 7 s
# Power curve, keeping SESOI, for varying N
pwr2 <- extend(model, along = "PPNR", n = 150)
power.curve <- powerCurve(pwr2, test = compare(model0, method = "lr"), nsim = 300, along = "PPNR", progress = F)
power.curve
## Power for model comparison, (95% confidence interval),
## by largest value of PPNR:
## 3: 8.67% ( 5.74, 12.44) - 279 rows
## 19: 34.33% (28.97, 40.01) - 1767 rows
## 36: 52.33% (46.52, 58.10) - 3348 rows
## 52: 66.00% (60.33, 71.35) - 4836 rows
## 68: 77.67% (72.53, 82.25) - 6324 rows
## 85: 85.00% (80.45, 88.84) - 7905 rows
## 101: 92.33% (88.72, 95.08) - 9393 rows
## 117: 95.00% (91.89, 97.17) - 10881 rows
## 134: 97.33% (94.81, 98.84) - 12462 rows
## 150: 98.33% (96.15, 99.46) - 13950 rows
##
## Time elapsed: 0 h 31 m 8 s
plot(power.curve)
#POWER - Interaction of Face Expression and Communal Strength
model0.int <- lmer(Temperature ~ FacExp + stCs_tot + duration + (1|PPNR), data = Pdata)
model.int <- lmer(Temperature ~ FacExp * stCs_tot + duration + (1|PPNR), data = Pdata)
#model.int
summary(model.int)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Temperature ~ FacExp * stCs_tot + duration + (1 | PPNR)
## Data: Pdata
##
## REML criterion at convergence: 13553.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5330 -0.5180 0.0409 0.5150 3.5623
##
## Random effects:
## Groups Name Variance Std.Dev.
## PPNR (Intercept) 13.962 3.737
## Residual 0.429 0.655
## Number of obs: 6510, groups: PPNR, 70
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.878e+01 4.470e-01 6.819e+01 64.392 < 2e-16 ***
## FacExp1 -2.082e-02 9.032e-03 6.437e+03 -2.305 0.0212 *
## stCs_tot -3.711e-01 4.467e-01 6.800e+01 -0.831 0.4091
## duration 2.944e-03 5.817e-05 6.437e+03 50.614 < 2e-16 ***
## FacExp1:stCs_tot 6.850e-02 8.387e-03 6.437e+03 8.168 3.75e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) FcExp1 stCs_t duratn
## FacExp1 -0.011
## stCs_tot 0.000 0.000
## duration -0.038 0.371 0.000
## FcExp1:stC_ 0.000 0.001 0.004 0.002
# Sensitivity power analysis for the fixed N.
# Estimating the power to detect ES from .01 to .05 for the interaction
N <- 120 # Enter the value determined by the previous power estimation for the focal hypothesis
es.vect <- seq(.01, .04, .005)
pwr3 <- NA
for(i in 1:length(es.vect)){
fixef(model.int)["FacExp1:stCs_tot"] <- es.vect[i]
pwr3.model <- extend(model.int, along = "PPNR", n = N)
pwr3[i] <- round(as.numeric(summary(powerSim(pwr3.model, test = compare(model0.int, method = "lr"), nsim = 300, progress = F)))[3],3)
}
cbind("Unstd ES" = es.vect, Power = pwr3)
## Unstd ES Power
## [1,] 0.010 0.363
## [2,] 0.015 0.660
## [3,] 0.020 0.873
## [4,] 0.025 0.957
## [5,] 0.030 1.000
## [6,] 0.035 1.000
## [7,] 0.040 1.000
# Check how does this unstandardized ES translates to eta squared
int.ES <- .025 # this is to be changed to the resulting magnitude of the interaction term
fixef(model.int)["FacExp1:stCs_tot"] <- int.ES
r2beta(model.int)
## Effect Rsq upper.CL lower.CL
## 1 Model 0.425 0.441 0.409
## 4 duration 0.268 0.285 0.251
## 3 stCs_tot 0.221 0.238 0.204
## 5 FacExp1:stCs_tot 0.001 0.004 0.000
## 2 FacExp1 0.001 0.003 0.000
# As a benchmark, approximate Cohen's d for the effect of FacExp
paste("n^2 = ", round(r2beta(model.int)$Rsq[4], 3))
## [1] "n^2 = 0.001"