Post-Hoc Power Analysis for Multilevel Models in R
Examples from B & L Chapter 10
Overview
This tutorial illustrates how to do power analysis for multilevel
models using the simr package in R. The examples follow
Chapter 10 (“Statistical Power for Intensive Longitudinal Designs”) from
Bolger & Laurenceau’s book Intensive longitudinal methods: An
introduction to diary and experience sampling research (referred to
as “B & L” from here on). Please see the chapter for additional
context about the analyses.
This tutorial covers:
A. Power Analysis for Time Course
This example fits a simple multilevel model, i.e., growth model, that examines how variables change over time.
B. Power Analysis for Within-Person Process
This example fits a multilevel process model that investigates causal processes between within-person covariates (rather than time itself) and their outcomes.
C. Power Analysis for Categorical Outcomes
This example fits a multilevel model for categorical data, i.e., generalized linear mixed effect model.
In all three cases, we load example data from B & L and conduct
post-hoc power analysis for the multilevel models. After loading the
data, we fit a multilevel model using the lme4 package.
Based on the fitted lme4 object, we then perform a
simulation-based power analysis for the specified multilevel model using
functions from the simr package. Finally, we illustrate how
statistical power changes as the number of persons or the number of
occasions increases, thus enabling us to determine how many
persons/occasions are necessary for ensuring adequate power.
Preliminaries
Loading libraries used in this script.
A. Power Analysis for Time Course
This example makes use of the Intimacy data set from Chapter 4 in B & L. The data was contributed by 50 wives from 50 married heterosexual couples, with each person randomly assigned to a 16-week marital therapy treatment (N = 25) or a control (N = 25) group.
We first load the Intimacy data (N = 50, T = 16).
#set filepath for data file
filepath <- "https://raw.githubusercontent.com/The-Change-Lab/collaborations/main/Bolger_Laurenceau_2013/BL2013_data_intimacy.csv"
#read in the .csv file using the url() function
data_intimacy <- read.csv(file=url(filepath), header=TRUE)We next describe the data. The outcome variable is
intimacy scores, the time variable is time01
(scaled time variable; range = 0 to 1), and the treatment indicator is
treatment.
## vars n mean sd median trimmed mad min max range skew
## id 1 800 25.50 14.44 25.50 25.50 18.53 1 50.00 49.00 0.00
## time 2 800 7.50 4.61 7.50 7.50 5.93 0 15.00 15.00 0.00
## time01 3 800 0.50 0.31 0.50 0.50 0.40 0 1.00 1.00 0.00
## intimacy 4 800 3.47 1.63 3.33 3.43 1.61 0 9.41 9.41 0.28
## treatment 5 800 0.50 0.50 0.50 0.50 0.74 0 1.00 1.00 0.00
## kurtosis se
## id -1.21 0.51
## time -1.21 0.16
## time01 -1.21 0.01
## intimacy -0.06 0.06
## treatment -2.00 0.02## [1] 50## [1] 16## time
## id 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 7 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 8 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 9 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 10 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 11 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 12 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 13 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 14 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 15 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 16 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 17 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 18 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 19 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 20 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 21 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 22 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 23 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 24 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 25 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 26 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 27 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 28 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 29 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 30 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 31 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 32 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 33 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 34 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 35 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 36 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 37 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 38 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 39 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 40 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 41 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 42 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 43 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 44 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 45 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 46 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 47 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 48 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 49 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 50 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1We see that there are 50 persons with 16 occasions each (coded 0 to 15). The data are fully complete.
1) Running the Multilevel with the Available Data
We run the multilevel using the lmer() function from the
lme4 library because the simr power analysis
package was built to work with lmer().
#Run linear growth model without the AR(1) errors
fit_lgmodel <- lme4::lmer(intimacy ~ 1 + time01 + treatment + time01:treatment +
(1 + time01 | id),
data=data_intimacy,
na.action=na.exclude)
#examine model summary
summary(fit_lgmodel)## Linear mixed model fit by REML ['lmerMod']
## Formula: intimacy ~ 1 + time01 + treatment + time01:treatment + (1 + time01 |
## id)
## Data: data_intimacy
##
## REML criterion at convergence: 2834.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6242 -0.6798 -0.0190 0.6433 2.4463
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 0.6858 0.8281
## time01 1.8936 1.3761 -0.45
## Residual 1.6925 1.3010
## Number of obs: 800, groups: id, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.89897 0.20703 14.003
## time01 0.73520 0.34720 2.118
## treatment -0.05644 0.29279 -0.193
## time01:treatment 0.92143 0.49101 1.877
##
## Correlation of Fixed Effects:
## (Intr) time01 trtmnt
## time01 -0.599
## treatment -0.707 0.423
## tm01:trtmnt 0.423 -0.707 -0.599The results are very similar to those in the book even through the
model does not include the auto-correlation of residuals.
lmer() does not provide for alternative error
structures.
The object obtained from lmer(), formally referred to as
a “merMod” object, is used in the power analysis. Note that this
object must be obtained from the lme4::lmer() function and
NOT the lmerTest::lmer() function. If you use the
lmerTest version, the subsequent steps will
fail.
2) Examining Statistical Power
We can simulate power directly from the fitted model using the
powerSim() function. We supply the function with the fitted
model object (fit_lgmodel), the specific parameter we would
like to examine ("time01"), the type of test to use
("lr"), the number of simulations to run
(nsim=100), and the alpha level to use in the evaluations
of significance (alpha=.05). Typically, power analysis
simulations use 1000 or more simulations. To save time here, we only use
100.
Testing Fixed Effect Parameters
For fixed effect parameters, we evaluate all three parameters in the model: (1) time01, (2) treatment, (3) time01:treatment interaction effect.
Simulating power directly from fitted model for the time01 effect.
#Power for the effect of interest
#Note "time01" must exactly match what appears in the lmer output.
SimPower_Fixed <- simr::powerSim(fit_lgmodel,
test=simr::fixed("time01", method="lr"),
seed=1234, #set for replication
nsim=100, #set low for time or high for precision
alpha=.05, progress=TRUE) #options## Simulating: | |Simulating: |= |Simulating: |== |Simulating: |=== |Simulating: |==== |Simulating: |===== |Simulating: |====== |Simulating: |======= |Simulating: |======== |Simulating: |========= |Simulating: |========== |Simulating: |=========== |Simulating: |============ |Simulating: |============= |Simulating: |============== |Simulating: |=============== |Simulating: |================ |Simulating: |================= |Simulating: |================== |Simulating: |=================== |Simulating: |==================== |Simulating: |===================== |Simulating: |====================== |Simulating: |======================= |Simulating: |======================== |Simulating: |========================= |Simulating: |========================== |Simulating: |=========================== |Simulating: |============================ |Simulating: |============================= |Simulating: |============================== |Simulating: |=============================== |Simulating: |================================ |Simulating: |================================= |Simulating: |================================== |Simulating: |=================================== |Simulating: |==================================== |Simulating: |===================================== |Simulating: |====================================== |Simulating: |======================================= |Simulating: |======================================== |Simulating: |========================================= |Simulating: |========================================== |Simulating: |=========================================== |Simulating: |============================================ |Simulating: |============================================= |Simulating: |============================================== |Simulating: |=============================================== |Simulating: |================================================ |Simulating: |================================================= |Simulating: |================================================== |Simulating: |=================================================== |Simulating: |==================================================== |Simulating: |===================================================== |Simulating: |====================================================== |Simulating: |======================================================= |Simulating: |======================================================== |Simulating: |========================================================= |Simulating: |========================================================== |Simulating: |=========================================================== |Simulating: |============================================================ |Simulating: |============================================================= |Simulating: |============================================================== |Simulating: |=============================================================== |Simulating: |================================================================ |Simulating: |================================================================= |Simulating: |==================================================================|## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## Power for predictor 'time01', (95% confidence interval):
## 62.00% (51.75, 71.52)
##
## Test: Likelihood ratio
## Effect size for time01 is 0.74
##
## Based on 100 simulations, (100 warnings, 0 errors)
## alpha = 0.05, nrow = 800
##
## Time elapsed: 0 h 0 m 5 s
##
## nb: result might be an observed power calculationAs we specified in the function call, powerSim()
calculated the power for predictor time01 based on 100
simulations. In each simulation, the fixed effect is tested using the
likelihood ratio test (method="lr") at significance level
\(\alpha = .05\)
(alpha=.05). For other types of testing, such as z-test
(method="z") or t-test (method="t"), please
refer to ?simr::tests.
Power is usually set as 80%. We see that with an effect size
(or, estimate) of 0.74, power for the time01 predictor is
62%, which is less than ideal. (Note that effect size here is in raw
units of the data – NOT standardized units.) Results also include a 95%
confidence interval for the estimated power: between 51.75% and
71.52%.
Alternatively, instead of using the fixed() function,
one can specify a (smaller) model formula without the parameter of
interest using the compare() function. If one is testing a
single fixed effect, fixed() function will be enough, but
for more complicated tests, compare() function can be
useful to specify the exact model formula.
#Power for the effect of interest
#Note "time01" must exactly match what appears in the lmer output.
SimPower_Fixed <- simr::powerSim(fit_lgmodel,
test=simr::compare(.~1 + treatment + time01:treatment +
(1 + time01 | id)),
seed=1234, #set for replication
nsim=100, #set low for time or high for real
alpha=.05, progress=TRUE) #options## Simulating: | |Simulating: |= |Simulating: |== |Simulating: |=== |Simulating: |==== |Simulating: |===== |Simulating: |====== |Simulating: |======= |Simulating: |======== |Simulating: |========= |Simulating: |========== |Simulating: |=========== |Simulating: |============ |Simulating: |============= |Simulating: |============== |Simulating: |=============== |Simulating: |================ |Simulating: |================= |Simulating: |================== |Simulating: |=================== |Simulating: |==================== |Simulating: |===================== |Simulating: |====================== |Simulating: |======================= |Simulating: |======================== |Simulating: |========================= |Simulating: |========================== |Simulating: |=========================== |Simulating: |============================ |Simulating: |============================= |Simulating: |============================== |Simulating: |=============================== |Simulating: |================================ |Simulating: |================================= |Simulating: |================================== |Simulating: |=================================== |Simulating: |==================================== |Simulating: |===================================== |Simulating: |====================================== |Simulating: |======================================= |Simulating: |======================================== |Simulating: |========================================= |Simulating: |========================================== |Simulating: |=========================================== |Simulating: |============================================ |Simulating: |============================================= |Simulating: |============================================== |Simulating: |=============================================== |Simulating: |================================================ |Simulating: |================================================= |Simulating: |================================================== |Simulating: |=================================================== |Simulating: |==================================================== |Simulating: |===================================================== |Simulating: |====================================================== |Simulating: |======================================================= |Simulating: |======================================================== |Simulating: |========================================================= |Simulating: |========================================================== |Simulating: |=========================================================== |Simulating: |============================================================ |Simulating: |============================================================= |Simulating: |============================================================== |Simulating: |=============================================================== |Simulating: |================================================================ |Simulating: |================================================================= |Simulating: |==================================================================|## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## Power for model comparison, (95% confidence interval):
## 62.00% (51.75, 71.52)
##
## Test: Likelihood ratio
## Comparison to . ~ 1 + treatment + time01:treatment + (1 + time01 | id)
##
## Based on 100 simulations, (1 warning, 0 errors)
## alpha = 0.05, nrow = 800
##
## Time elapsed: 0 h 0 m 5 s
##
## nb: result might be an observed power calculationWe see that the results are identical. For simplicity going forward,
we will use the fixed() function.
Simulating power directly from fitted model for the treatment effect.
#Power for the effect of interest
SimPower_Fixed <- simr::powerSim(fit_lgmodel,
test=simr::fixed("treatment", method="lr"),
seed=1234, #set for replication
nsim=100, #set low for time or high for real
alpha=.05, progress=FALSE) #options## boundary (singular) fit: see help('isSingular')## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## Power for predictor 'treatment', (95% confidence interval):
## 5.00% ( 1.64, 11.28)
##
## Test: Likelihood ratio
## Effect size for treatment is -0.056
##
## Based on 100 simulations, (100 warnings, 0 errors)
## alpha = 0.05, nrow = 800
##
## Time elapsed: 0 h 0 m 5 s
##
## nb: result might be an observed power calculationWe see that with an effect size of -0.06, power for the
treatment predictor is 5%, which is very low.
Simulating power directly from fitted model for the time01:treatment interaction effect.
#Power for the effect of interest
SimPower_Fixed <- simr::powerSim(fit_lgmodel,
test=simr::fixed("time01:treatment", method="lr"),
seed=1234, #set for replication
nsim=100, #set low for time or high for real
alpha=.05, progress=FALSE) #options## boundary (singular) fit: see help('isSingular')## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## Power for predictor 'time01:treatment', (95% confidence interval):
## 45.00% (35.03, 55.27)
##
## Test: Likelihood ratio
## Effect size for time01:treatment is 0.92
##
## Based on 100 simulations, (1 warning, 0 errors)
## alpha = 0.05, nrow = 800
##
## Time elapsed: 0 h 0 m 5 s
##
## nb: result might be an observed power calculationWe see that with an effect size of 0.92, power for the
time01:treatment interaction is 45%, which is lower than
ideal.
All three of these evaluations are very similar to what was obtained using Mplus and presented in the book: the powers are lower than the ideal threshold of 80%. As for the next steps, we will extend the size of the data to examine how statistical power changes as the number of persons or the number of occasions increases.
Testing Random Effect Parameters
Additionally, we can conduct power analysis for random effects. To do so, we evaluate the covariance parameter between the random intercept and random slope (the random intercept/slope cannot be separated from the covariance parameter).
Specifically, we compare the full model (fit_lgmodel)
with a reduced model without the covariance parameter between random
intercept/slope using the compare() function. Note that
|| specifies that the random intercept and random slope are
independent of each other.
#Power for the effect of interest
SimPower_Random <- simr::powerSim(fit_lgmodel,
test=simr::compare(.~1 + time01 + treatment + time01:treatment +
(1 + time01 || id)),
seed=1234, #set for replication
nsim=100, #set low for time or high for real
alpha=.05, progress=TRUE) #options## Simulating: | |Simulating: |= |Simulating: |== |Simulating: |=== |Simulating: |==== |Simulating: |===== |Simulating: |====== |Simulating: |======= |Simulating: |======== |Simulating: |========= |Simulating: |========== |Simulating: |=========== |Simulating: |============ |Simulating: |============= |Simulating: |============== |Simulating: |=============== |Simulating: |================ |Simulating: |================= |Simulating: |================== |Simulating: |=================== |Simulating: |==================== |Simulating: |===================== |Simulating: |====================== |Simulating: |======================= |Simulating: |======================== |Simulating: |========================= |Simulating: |========================== |Simulating: |=========================== |Simulating: |============================ |Simulating: |============================= |Simulating: |============================== |Simulating: |=============================== |Simulating: |================================ |Simulating: |================================= |Simulating: |================================== |Simulating: |=================================== |Simulating: |==================================== |Simulating: |===================================== |Simulating: |====================================== |Simulating: |======================================= |Simulating: |======================================== |Simulating: |========================================= |Simulating: |========================================== |Simulating: |=========================================== |Simulating: |============================================ |Simulating: |============================================= |Simulating: |============================================== |Simulating: |=============================================== |Simulating: |================================================ |Simulating: |================================================= |Simulating: |================================================== |Simulating: |=================================================== |Simulating: |==================================================== |Simulating: |===================================================== |Simulating: |====================================================== |Simulating: |======================================================= |Simulating: |======================================================== |Simulating: |========================================================= |Simulating: |========================================================== |Simulating: |=========================================================== |Simulating: |============================================================ |Simulating: |============================================================= |Simulating: |============================================================== |Simulating: |=============================================================== |Simulating: |================================================================ |Simulating: |================================================================= |Simulating: |==================================================================|## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## Power for model comparison, (95% confidence interval):
## 56.00% (45.72, 65.92)
##
## Test: Likelihood ratio
## Comparison to . ~ 1 + time01 + treatment + time01:treatment + (1 + time01 || id)
##
## Based on 100 simulations, (1 warning, 0 errors)
## alpha = 0.05, nrow = 800
##
## Time elapsed: 0 h 0 m 5 s
##
## nb: result might be an observed power calculationWe see that with an effect size of -0.45, power for the covariance parameter is 56%, lower than ideal. This is very similar to what was presented in the book.
3) Extending the Size of the Data
Having found that our analysis is underpowered, we would like to evaluate what power is necessary for detecting the interaction at different sample sizes – whether by adding more persons or adding more occasions within person.
To evaluate sample sizes that are larger than the original data, we
first “extend” the data by creating a new lmer() object.
Specifically, we use the extend() function in the simr
package.
We do that by extending the data that is inside the
lmer() model object to include more participants or more
occasions. Let’s look at size of the data within the lmer()
object.
## [1] 50## [1] 16In a model entitled fit_lgmodelpersons, we extend the
number of persons from 50 to 200.
#extending number of persons
fit_lgmodelpersons <- simr::extend(fit_lgmodel, along="id", n=200)
#checking length of data, number of persons
length(unique(getData(fit_lgmodelpersons)$id))## [1] 200## [1] 16In a separate model entitled fit_lgmodeltimes, we extend
the number of occasions from 16 to 60.
#extending number of occasions
fit_lgmodeltimes <- simr::extend(fit_lgmodel, along="time", n=60)
#checking length of data, number of persons
length(unique(getData(fit_lgmodeltimes)$id))## [1] 50## [1] 60Using the new model objects with extended data, we can examine what power looks like with N = 200 persons and/or T = 60 occasions.
We first run a power curve analysis for a range of sample sizes (or,
persons) by specifying along="id" – each with
T = 16 occasions. Specific sample sizes ("breaks")
were selected so that the total number of observations, N *
T, is near 500, 1000, 1500, 2000, 2500, and 3000.
#running power curve
powerCurve_persons <- simr::powerCurve(fit_lgmodelpersons,
simr::fixed("time01:treatment", method="lr"),
along="id",
breaks=c(32,63,94,125,157,188),
nsim=100, alpha=.05, progress=FALSE)## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## nrow nlevels successes trials mean lower upper
## 1 512 32 28 100 0.28 0.1947936 0.3786670
## 2 1008 63 54 100 0.54 0.4374116 0.6401566
## 3 1504 94 76 100 0.76 0.6642645 0.8397754
## 4 2000 125 87 100 0.87 0.7879593 0.9289270
## 5 2512 157 92 100 0.92 0.8484424 0.9648284
## 6 3008 188 98 100 0.98 0.9296161 0.9975687The power estimates at each sample size can be found in the
mean column. From these, we can see how power increases as
number of persons increases: from 0.28 with 512 observations to 0.98
with 3008 observations.
We first run a power curve analysis for a range of sample sizes (or,
occasions) by specifying within="id" – each with
N = 50 persons. Specific sample sizes ("breaks")
were selected so that the total number of observations, N *
T, is near 500, 1000, 1500, 2000, 2500, and 3000.
As in the book, we examine the time01:treatment interaction parameter.
#running power curve
powerCurve_times <- simr::powerCurve(fit_lgmodeltimes,
simr::fixed("time01:treatment", method="lr"),
within="id",
breaks=c(10,20,30,40,50,60),
nsim=100, alpha=.05, progress=FALSE)## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## nrow nlevels successes trials mean lower upper
## 1 500 10 29 100 0.29 0.2035742 0.3892660
## 2 1000 20 49 100 0.49 0.3886442 0.5919637
## 3 1500 30 51 100 0.51 0.4080363 0.6113558
## 4 2000 40 54 100 0.54 0.4374116 0.6401566
## 5 2500 50 59 100 0.59 0.4871442 0.6873800
## 6 3000 60 60 100 0.60 0.4972092 0.6967052The power estimates at each sample size can be found in the
mean column. In this case, we can see how power does not
increase very much as the number of occasions increases: from 0.29 with
500 observations to 0.60 with 3000 observations.
4) Plotting Power Curves
Of course, power curves are often useful in visual form. It takes a
few steps to get them into a form useful for ggplot().
We first prepare the summary outputs for plotting. We will display
nlevels on the x-axis (persons/times) and power estimates
(mean with 95% confidence interval [lower,
upper]) on the y-axis.
#preparing persons power curve data
powerdata_persons <- as.data.frame(summary(powerCurve_persons))
powerdata_persons$persons <- powerdata_persons$nlevels
powerdata_persons$times <- 16
powerdata_persons$design <- "persons"
#preparing times power curve data
powerdata_times <- as.data.frame(summary(powerCurve_times))
powerdata_times$persons <- 50
powerdata_times$times <- powerdata_times$nlevels
powerdata_times$design <- "times"
#merging together
powerdata_all <- rbind(powerdata_persons, powerdata_times)Now we can plot the power curves. Specify design as the
group variable to show the difference between the effects of adding more
participants/occasions.
#Plotting curves
powerdata_all %>%
ggplot(aes(y=mean, x=nrow, group=factor(design))) +
#power goal
geom_hline(yintercept=0.8, color="gray40",lty = 2) +
#persons curve
geom_line(color="black") +
geom_errorbar(aes(ymin=lower, ymax=upper), width=50, color="black") +
geom_point(aes(shape=design), size=3, color="black", show.legend=FALSE) +
scale_shape_manual(values=c(15,17)) +
scale_x_continuous(limits=c(500,3100),
breaks=c(500,1000,1500,2000,2500,3000),
name="Number of Total Observations") +
scale_y_continuous(limits=c(0,1.0),
breaks=c(0,.2,.4,.6,.8,1.0),
name="Power") +
ggtitle("Power Curves for time01:treatment Interaction") +
theme_classic()
Nice! We can see that power increases with both larger sample size (squares) and additional occasions (triangles), but at a much slower rate for occasions. Our plot mirrors Figure 10.1 in B & L. Please refer to the book for more in-depth discussion of the implications, but the general takeaway here is that more level-2 units (here, persons) is better.
B. Power Analysis for the Within-Person Process
This example makes use of the Process data set from Chapter 5 of B & L. The dataset contains daily relationship conflicts and daily relational intimacy from 66 women for 28 days.
We first load the Process data (N = 66, T = 28).
#set filepath for data file
filepath <- "https://raw.githubusercontent.com/The-Change-Lab/collaborations/main/Bolger_Laurenceau_2013/BL2013_data_process.csv"
#read in the .csv file using the url() function
data_process <- read.csv(file=url(filepath),header=TRUE)We next describe the data. The outcome variable is
intimacy scores and the parameter of interest is
relqual:confcw, or the interaction between
relqual (relationship quality; range = 0 to 1) and
confcw (raw scores of conflict; within-person
mean-centered). Please see Chapter 5 of B & L for further
details.
## vars n mean sd median trimmed mad min max range skew
## id 1 1848 33.50 19.06 33.50 33.50 24.46 1.00 66.00 65.00 0.00
## time 2 1848 13.50 8.08 13.50 13.50 10.38 0.00 27.00 27.00 0.00
## time7c 3 1848 0.00 1.15 0.00 0.00 1.48 -1.93 1.93 3.86 0.00
## intimacy 4 1848 4.82 2.28 4.84 4.83 2.34 0.00 10.00 10.00 -0.02
## conflict 5 1848 0.22 0.42 0.00 0.15 0.00 0.00 1.00 1.00 1.34
## confc 6 1848 0.00 0.42 -0.22 -0.07 0.00 -0.22 0.78 1.00 1.34
## confcb 7 1848 0.00 0.16 -0.04 -0.02 0.16 -0.22 0.42 0.64 0.75
## confcw 8 1848 0.00 0.38 -0.11 -0.04 0.16 -0.64 0.96 1.61 1.04
## relqual 9 1848 0.61 0.49 1.00 0.63 0.00 0.00 1.00 1.00 -0.43
## kurtosis se
## id -1.20 0.44
## time -1.21 0.19
## time7c -1.21 0.03
## intimacy -0.50 0.05
## conflict -0.20 0.01
## confc -0.20 0.01
## confcb -0.35 0.00
## confcw 0.04 0.01
## relqual -1.81 0.01## [1] 66## [1] 28## time
## id 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 7 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 8 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 9 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 10 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 11 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 12 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 13 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 14 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 15 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 16 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 17 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 18 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 19 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 20 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 21 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 22 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 23 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 24 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 25 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 26 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 27 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 28 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 29 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 30 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 31 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 32 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 33 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 34 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 35 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 36 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 37 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 38 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 39 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 40 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 41 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 42 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 43 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 44 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 45 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 46 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 47 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 48 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 49 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 50 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 51 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 52 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 53 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 54 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 55 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 56 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 57 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 58 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 59 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 60 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 61 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 62 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 63 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 64 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 65 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 66 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1We see that there are 66 persons with 28 occasions each (coded 0 to 27). The data are fully complete.
1) Running the Multilevel Model with the Available Data
We run a multilevel using lmer(). Remember, the
simr package is built to work with lmer().
Along with intimacy, relqual, and
confcw, this model formula also includes
confcb (raw scores of conflict; between-person
mean-centered) and time7c (time variable, scaled that
1-unit difference corresponds to the passage of 1-week).
#Run multilevel model
fit_cpmodel <- lme4::lmer(intimacy ~ 1 + relqual + confcb + relqual:confcb + confcw +
confcw:relqual + time7c +
(1 + confcw | id),
data=data_process,
na.action=na.exclude)
#examine model summary
summary(fit_cpmodel)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## intimacy ~ 1 + relqual + confcb + relqual:confcb + confcw + confcw:relqual +
## time7c + (1 + confcw | id)
## Data: data_process
##
## REML criterion at convergence: 7817.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.00949 -0.67665 -0.00451 0.67142 2.91605
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 0.8007 0.8948
## confcw 2.6990 1.6429 0.27
## Residual 3.5908 1.8949
## Number of obs: 1848, groups: id, 66
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.53208 0.22131 20.478
## relqual 0.64735 0.28212 2.295
## confcb -0.84166 1.10431 -0.762
## confcw -2.02443 0.37061 -5.462
## time7c -0.02820 0.03865 -0.730
## relqual:confcb 2.52633 1.68192 1.502
## relqual:confcw 1.03381 0.49232 2.100
##
## Correlation of Fixed Effects:
## (Intr) relqul confcb confcw time7c rlql:cnfcb
## relqual -0.784
## confcb -0.520 0.408
## confcw 0.166 -0.131 0.036
## time7c 0.001 -0.001 -0.001 -0.008
## relql:cnfcb 0.341 -0.038 -0.657 -0.024 -0.002
## relql:cnfcw -0.125 0.178 -0.027 -0.753 0.005 0.040The results are similar to those in the book.
2) Examining Statistical Power
As before, we can simulate power directly from the fitted model using
the powerSim() function. We supply the function with the
fitted model object (fit_cpmodel), the specific parameter
we would like to examine ("relqual:confcw"), the type of
test to use ("lr"), the number of simulations to run
(nsim=100), and the alpha level to use in the evaluations
of significance (alpha=.05). Typically, power analysis
simulations use 1000 or more simulations. To save time here, we only use
100.
Testing Fixed Effect Parameters
This time, we evaluate the main parameter of interest (the relqual:confcw interaction effect) by simulating power directly from the fitted model.
#Power for the effect of interest
SimPower_Fixed <- simr::powerSim(fit_cpmodel,
test=simr::fixed("relqual:confcw", method="lr"),
seed=1234, #set for replication
nsim=100, #set low for time or high for accuracy
alpha=.05, progress=FALSE) #options## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## Power for predictor 'relqual:confcw', (95% confidence interval):
## 58.00% (47.71, 67.80)
##
## Test: Likelihood ratio
## Effect size for relqual:confcw is 1.0
##
## Based on 100 simulations, (1 warning, 0 errors)
## alpha = 0.05, nrow = 1848
##
## Time elapsed: 0 h 0 m 10 s
##
## nb: result might be an observed power calculationWe see that with an effect size of 1.03, power for the relqual:confcw interaction is 58%, just as was found in the book.
Testing Random Effect Parameters
Additionally, we can conduct power analysis for random effects by evaluating the covariance parameter between the random intercept and random slope (the random intercept/slope cannot be separated from the covariance parameter).
Specifically, we compare the full model (fit_cpmodel)
with a reduced model without the covariance parameter between random
intercept/slope using the compare() function. Note that
|| specifies that the random intercept and random slope are
independent of each other.
#Power for the effect of interest
#Note "time01" must exactly match what appears in the lmer output.
SimPower_Random <- simr::powerSim(fit_cpmodel,
test=simr::compare(.~1 + relqual + confcb + relqual:confcb + confcw +
confcw:relqual + time7c +
(1 + confcw || id)),
seed=1234, #set for replication
nsim=100, #set low for time or high for real
alpha=.05, progress=TRUE) #options## Simulating: | |Simulating: |= |Simulating: |== |Simulating: |=== |Simulating: |==== |Simulating: |===== |Simulating: |====== |Simulating: |======= |Simulating: |======== |Simulating: |========= |Simulating: |========== |Simulating: |=========== |Simulating: |============ |Simulating: |============= |Simulating: |============== |Simulating: |=============== |Simulating: |================ |Simulating: |================= |Simulating: |================== |Simulating: |=================== |Simulating: |==================== |Simulating: |===================== |Simulating: |====================== |Simulating: |======================= |Simulating: |======================== |Simulating: |========================= |Simulating: |========================== |Simulating: |=========================== |Simulating: |============================ |Simulating: |============================= |Simulating: |============================== |Simulating: |=============================== |Simulating: |================================ |Simulating: |================================= |Simulating: |================================== |Simulating: |=================================== |Simulating: |==================================== |Simulating: |===================================== |Simulating: |====================================== |Simulating: |======================================= |Simulating: |======================================== |Simulating: |========================================= |Simulating: |========================================== |Simulating: |=========================================== |Simulating: |============================================ |Simulating: |============================================= |Simulating: |============================================== |Simulating: |=============================================== |Simulating: |================================================ |Simulating: |================================================= |Simulating: |================================================== |Simulating: |=================================================== |Simulating: |==================================================== |Simulating: |===================================================== |Simulating: |====================================================== |Simulating: |======================================================= |Simulating: |======================================================== |Simulating: |========================================================= |Simulating: |========================================================== |Simulating: |=========================================================== |Simulating: |============================================================ |Simulating: |============================================================= |Simulating: |============================================================== |Simulating: |=============================================================== |Simulating: |================================================================ |Simulating: |================================================================= |Simulating: |==================================================================|## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## Power for model comparison, (95% confidence interval):
## 43.00% (33.14, 53.29)
##
## Test: Likelihood ratio
## Comparison to . ~ 1 + relqual + confcb + relqual:confcb + confcw + confcw:relqual + time7c + (1 + confcw || id)
##
## Based on 100 simulations, (0 warnings, 0 errors)
## alpha = 0.05, nrow = 1848
##
## Time elapsed: 0 h 0 m 9 s
##
## nb: result might be an observed power calculationWe see that with an effect size of 0.27, power for the covariance parameter is 43%, which is lower than ideal.
3) Extending the Size of the Data
Given this low power, we would like to evaluate what power is necessary to detect the interaction at different sample sizes –whether by either adding more persons or adding more occasions within person.
To evaluate sample sizes that are larger than the original data, we
first “extend” the data by creating a new lmer() object
that includes either more participants or more occasions. Specifically,
we use the extend() function from the simr
package.
Extending the number of persons from 66 to 200.
#extending number of persons
fit_cpmodelpersons <- simr::extend(fit_cpmodel, along="id", n=200)
#checking length of data, number of persons
length(unique(getData(fit_cpmodelpersons)$id))## [1] 200## [1] 28Extending the number of occasions from 28 to 70.
#extending number of occasions
fit_cpmodeltimes <- simr::extend(fit_cpmodel, along="time", n=70)
#checking length of data, number of persons
length(unique(getData(fit_cpmodeltimes)$id))## [1] 66## [1] 70Using the new model objects with extended data, we can examine what power looks like all the way up to N = 200 persons and/or T = 28 occasions.
To start, we run a power curve analysis for a range of sample sizes
(or, persons) by specifying along="id" – each with
T = 28 occasions.
#running power curve
powerCurve_persons <- simr::powerCurve(fit_cpmodelpersons,
simr::fixed("relqual:confcw", method="lr"),
along="id",
breaks=c(35,66,83,95,118,142,165),
nsim=100, alpha=.05, progress=FALSE)## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## nrow nlevels successes trials mean lower upper
## 1 980 35 31 100 0.31 0.2212888 0.4103146
## 2 1848 66 40 100 0.40 0.3032948 0.5027908
## 3 2324 83 52 100 0.52 0.4177898 0.6209945
## 4 2660 95 57 100 0.57 0.4671337 0.6686090
## 5 3304 118 69 100 0.69 0.5896854 0.7787112
## 6 3976 142 73 100 0.73 0.6319837 0.8139336
## 7 4620 165 84 100 0.84 0.7532124 0.9056897Based on the power estimates for each sample size
(mean), we see how power increases as the number of persons
increases: from 0.31 with 980 observations to 0.84 with 4620
observations.
We next run a power curve analysis for a range of sample sizes (or,
occasions) by specifying within="id" – each with
N = 66 persons. Specific sample sizes ("breaks")
were selected so that the total number of observations, N *
T, is near 500, 1000, 1500, 2000, 2500, and 3000.
#running power curve
powerCurve_times <- simr::powerCurve(fit_cpmodeltimes,
simr::fixed("relqual:confcw", method="lr"),
within="id",
breaks=c(15,28,35,40,50,60,70),
nsim=100, alpha=.05, progress=FALSE)## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## nrow nlevels successes trials mean lower upper
## 1 990 15 53 100 0.53 0.4275815 0.6305948
## 2 1848 28 57 100 0.57 0.4671337 0.6686090
## 3 2310 35 59 100 0.59 0.4871442 0.6873800
## 4 2640 40 58 100 0.58 0.4771192 0.6780145
## 5 3300 50 62 100 0.62 0.5174607 0.7152325
## 6 3960 60 65 100 0.65 0.5481506 0.7427062
## 7 4620 70 65 100 0.65 0.5481506 0.7427062Based on the power estimates for each sample size
(mean), we see how power does not increase very much as the
number of occasions increases: from 0.53 with 990 observations to 0.65
with 4620 observations.
4) Plotting Power Curves
Of course, power curves are often useful in visual form. It takes a
few steps to get them into a form useful for ggplot().
Prepare the summary outputs for plotting. We will display
nlevels on the x-axis (persons/times) and power estimates
(mean with 95% confidence interval [lower,
upper]) on the y-axis.
#preparing persons power curve data
powerdata_persons <- as.data.frame(summary(powerCurve_persons))
powerdata_persons$persons <- powerdata_persons$nlevels
powerdata_persons$times <- 28
powerdata_persons$design <- "persons"
#preparing times power curve data
powerdata_times <- as.data.frame(summary(powerCurve_times))
powerdata_times$persons <- 66
powerdata_times$times <- powerdata_times$nlevels
powerdata_times$design <- "times"
#merging together
powerdata_all <- rbind(powerdata_persons, powerdata_times)We next plot the power curves. Specify design as the
group variable so that we can see the differences in effects of adding
more participants (squares) or more occasions (triangles).
#Plotting curves
powerdata_all %>%
ggplot(aes(y=mean, x=nrow, group=factor(design))) +
#power goal
geom_hline(yintercept=0.8, color="gray40",lty = 2) +
#persons curve
geom_line(color="black") +
geom_errorbar(aes(ymin=lower, ymax=upper), width=50, color="black") +
geom_point(aes(shape=design), size=3, color="black", show.legend=FALSE) +
scale_shape_manual(values=c(15,17)) +
scale_x_continuous(limits=c(950,5000),
breaks=c(990,1848,2320,2640,3300,3960,4620),
name="Number of Total Observations") +
scale_y_continuous(limits=c(0,1.0),
breaks=c(0,.2,.4,.6,.8,1.0),
name="Power") +
theme_classic() +
ggtitle("Power Curves for relqual:confcw Interaction")
Nice! We can see that power increases with both adding more persons (squares) and adding more occasions (triangles) – but less so with additional occasions. This is the same as the Figure 10.2 in B & L. Please refer to B & L for more in-depth discussion of the implications, but the general takeaway is that more level-2 units (here, persons) is better.
C. Power Analysis for the Categorical Outcomes
This example makes use of the Categorical data set from B & L Chapter 6. The dataset features 61 heterosexual couples for a 4-week period (27 days).
We first load the Categorical data (N = 61, T = 27).
#set filepath for data file
filepath <- "https://raw.githubusercontent.com/The-Change-Lab/collaborations/main/Bolger_Laurenceau_2013/BL2013_data_categorical.csv"
#read in the .csv file using the url() function
data_categorical <- read.csv(file=url(filepath), header=TRUE)We next describe the data. The outcome variable is the male partner’s
evening report of the occurrence of a conflict that day
(pconf; 1 if reported else 0) and the predictor variable is
the female partner’s morning report of anger/irritability
(amang).
## vars n mean sd median trimmed mad min max range skew
## id 1 1345 45.61 30.19 38.00 44.44 37.06 1.00 102.00 101.00 0.28
## time 2 1345 14.66 7.76 14.00 14.60 10.38 2.00 28.00 26.00 0.06
## time7c 3 1345 0.00 1.15 -0.10 -0.01 1.54 -1.88 1.98 3.85 0.06
## pconf 4 1345 0.14 0.35 0.00 0.06 0.00 0.00 1.00 1.00 2.02
## lpconf 5 1345 0.16 0.36 0.00 0.07 0.00 0.00 1.00 1.00 1.88
## lpconfc 6 1345 0.00 0.36 -0.16 -0.09 0.00 -0.16 0.84 1.00 1.88
## amang 7 1345 0.49 1.11 0.00 0.22 0.00 0.00 10.00 10.00 3.96
## amangc 8 1345 0.00 1.11 -0.49 -0.27 0.00 -0.49 9.51 10.00 3.96
## amangcb 9 1345 0.00 0.49 -0.18 -0.10 0.27 -0.47 1.61 2.08 1.72
## amangcw 10 1345 0.00 1.00 -0.17 -0.15 0.27 -2.10 9.40 11.50 3.78
## kurtosis se
## id -1.29 0.82
## time -1.18 0.21
## time7c -1.18 0.03
## pconf 2.09 0.01
## lpconf 1.55 0.01
## lpconfc 1.55 0.01
## amang 21.05 0.03
## amangc 21.05 0.03
## amangcb 2.32 0.01
## amangcw 22.51 0.03## [1] 61## [1] 27#crossing of persons and occasions
head(xtabs(~id+time, data=data_categorical), 10) #just top of output## time
## id 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0
## 6 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 9 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 10 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 11 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1
## 12 1 1 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 13 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1
## 14 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0We see that there are 61 persons with 27 occasions (coded 2 to 28). Notably, this time the data are not fully complete.
We can check the number of observations contributed by each person.
#calculating intraindividual stats (counts)
istats <- data_categorical %>%
group_by(id) %>%
summarize(icount_pconf = sum(!is.na(pconf)), #count of observations
icount_amang = sum(!is.na(amang)))
describe(istats)## vars n mean sd median trimmed mad min max range skew
## id 1 61 46.30 29.86 42 45.24 40.03 1 102 101 0.26
## icount_pconf 2 61 22.05 6.52 25 23.35 2.97 5 27 22 -1.43
## icount_amang 3 61 22.05 6.52 25 23.35 2.97 5 27 22 -1.43
## kurtosis se
## id -1.32 3.82
## icount_pconf 0.88 0.84
## icount_amang 0.88 0.84## icount_pconf
## 5 6 11 12 13 14 15 17 19 20 21 22 23 24 25 26 27
## 2 3 1 1 2 1 1 1 2 1 2 2 7 2 10 2 21We see that on average individual completed 22 occasions, and that many completed 25 or 27.
1) Running the Multilevel Model with the Available Data
We run multilevel using lmer() from the
lme4 package. The model formula includes the male partner’s
evening report of a conflict that day (pconf) and yesterday
(lpconfc), the female partner’s morning report of anger
after within- (amangcw) and between-subject
(amangcb) mean-centering, and time7c (time
variable, scaled that 1-unit difference corresponds to the passage of
1-week). The predictor variable of interest is amangcw.
#Run multilevel model
fit_catmodel <- lme4::glmer(pconf ~ 1 + amangcw + amangcb + lpconfc + time7c +
(1 | id),
family=binomial,
data=data_categorical,
na.action=na.exclude)
#examine model summary
summary(fit_catmodel)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: pconf ~ 1 + amangcw + amangcb + lpconfc + time7c + (1 | id)
## Data: data_categorical
##
## AIC BIC logLik deviance df.resid
## 1084.8 1116.0 -536.4 1072.8 1339
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.7981 -0.4224 -0.3533 -0.2893 4.5619
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.2475 0.4975
## Number of obs: 1345, groups: id, 61
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.90161 0.10894 -17.456 < 2e-16 ***
## amangcw 0.21565 0.06738 3.200 0.00137 **
## amangcb -0.22437 0.22398 -1.002 0.31649
## lpconfc 0.31826 0.20892 1.523 0.12768
## time7c -0.19222 0.07073 -2.718 0.00657 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) amngcw amngcb lpcnfc
## amangcw -0.108
## amangcb 0.050 -0.103
## lpconfc 0.006 -0.083 0.052
## time7c 0.132 -0.004 -0.027 0.102The results are very similar to the results reported in the book using MPlus.
2) Examining Statistical Power
As before, we can simulate power directly from the fitted model using
the powerSim() function. We supply the function with the
fitted model object (fit_catmodel), the specific parameter
we would like to examine ("amangcw"), the type of test to
use ("lr"), the number of simulations to run
(nsim=100), and the alpha to use in the evaluations of
significance (alpha=.05). Typically, power analysis
simulations use 1000 or more simulations. To save time here, we only use
100.
Testing Fixed Effect Parameters
This time, we just evaluate the main parameter of interest.
We simulate power directly from fitted model for the amangcw effect.
#Power for the effect of interest
SimPower_Direct <- simr::powerSim(fit_catmodel,
test=simr::fixed("amangcw","lr"),
seed=1234, #set for replication
nsim=100, #set low for time or high for real
alpha=.05, progress=FALSE) #options## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## Power for predictor 'amangcw', (95% confidence interval):
## 79.00% (69.71, 86.51)
##
## Test: Likelihood ratio
## Effect size for amangcw is 0.22
##
## Based on 100 simulations, (0 warnings, 0 errors)
## alpha = 0.05, nrow = 1345
##
## Time elapsed: 0 h 0 m 29 s
##
## nb: result might be an observed power calculationWe see that with an effect size of 0.22, power for the amangcw effect is 79%. Since the data object we are using has the missing data in it, we are doing the power analysis with, on average, 22 occasions (i.e., we do not need to reduce the data to get that analysis, as B & L do in Chapter 10.5).
Testing Random Effect Parameters
Following the book, we can also test power for the variance term for
the intercept. We use the powerSim() function and specify
the model to compare against, which in this case is a model without any
random effects (simr::random()).
#Power for the single random effect
SimPower_variance <- simr::powerSim(fit_catmodel,
test=simr::random(), #empty because only one
seed=1234, #set for replication
nsim=100, #set low for time or high for real
alpha=.05, progress=FALSE) #options## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## Power for a single random effect, (95% confidence interval):
## 0.00% ( 0.00, 3.62)
##
## Test: Exact restricted LRT (package RLRsim)
##
## Based on 100 simulations, (0 warnings, 100 errors)
## alpha = 0.05, nrow = 1345
##
## Time elapsed: 0 h 0 m 17 s
##
## nb: result might be an observed power calculation## stage index message
## 1 Testing 1 Invalid <m> specified. \n
## 2 Testing 2 Invalid <m> specified. \n
## 3 Testing 3 Invalid <m> specified. \n
## 4 Testing 4 Invalid <m> specified. \n
## 5 Testing 5 Invalid <m> specified. \n
## 6 Testing 6 Invalid <m> specified. \nUnfortunately, this power analysis did not work because the
simr package only supports linear mixed models when testing
single random effects (c.f., this code would work with
lmer() fits). For more information, see the vignette here:
https://cran.r-project.org/web/packages/simr/vignettes/examples.html
CAUTION! At this point, random effects can be tested
for linear mixed models (i.e., lmer() fits) with continuous
dependent variables, not for generalized mixed models (i.e.,
glmer() fits) with categorical dependent variables. Future
work is necessary for the simr package, but in this
tutorial, we will proceed by treating the categorical model fit as a
lmer() fit.
#Run multilevel model with lmer()
fit_catmodel <- lme4::lmer(pconf ~ 1 + amangcw + amangcb + lpconfc + time7c +
(1 | id),
data=data_categorical,
na.action=na.exclude)
#Power for the single random effect
SimPower_variance <- simr::powerSim(fit_catmodel,
test=simr::random(), #empty because only one
seed=1234, #set for replication
nsim=100, #set low for time or high for real
alpha=.05, progress=FALSE) #options## Warning in observedPowerWarning(sim): This appears to be an "observed power"
## calculation## Power for a single random effect, (95% confidence interval):
## 95.00% (88.72, 98.36)
##
## Test: Exact restricted LRT (package RLRsim)
##
## Based on 100 simulations, (0 warnings, 0 errors)
## alpha = 0.05, nrow = 1345
##
## Time elapsed: 0 h 0 m 10 s
##
## nb: result might be an observed power calculationNote that the lmer() function was used only for
illustrative purposes and NOT for analytical purposes.
Conclusion
This tutorial provided an introduction for conducting power analysis
for multilevel model in R using the simr package. We fit
models using “pilot” data and then extended those data in order to
examine power across a range of sample sizes (both persons and
occasions).
There is also the option to focus on effect sizes, rather than sample sizes, to achieve appropriate power called the Smallest Effect Size of Interest (SESOI). Please check out our next tutorial if this is of interest.
Thanks for powering through with us!