This vignette details how you can set up and execute a basic power
analysis for a bivariate random intercept crosslagged panel model
(RICLPM) using the powRICLPM
package. Throughout, an
illustrating example will be used in which we wish to detect a small
crosslagged effect \(\beta_{2}\)
(defined here as the effect of \(a_{1}^{*}\) to \(b_{2}^{*}\), where \(a_{1}^{*}\) and \(b_{2}^{*}\) denote the latent withinunit
components of \(a_{1}\) and \(b_{2}\), respectively) of 0.2
(standardized). For the design of our power analysis we follow the steps
in the strategy as described in Mulder (2022).
Various extensions are available for this basic power analysis, and are
described in the Vignette Extensions.
Step 1: Determine experimental conditions
Before performing the power analysis, you must first determine the experimental conditions of interest. Experimental conditions (or: simulation conditions) are defined by characteristics of the study design that can impact statistical power. This includes, among others, characteristics like the sample size and the number of repeated measures. Decide on the number of repeated measures that will be used in the simulations, as well as the range of sample sizes over which you want to simulate the power.
For this example, we take a sample size range from 100 to 1000 first, increasing with steps of 100. Let the numbers of repeated measures range from 3 to 5. If these experimental conditions do not lead to the desired amount of power for detecting the small crosslagged effect, the ranges can be extended later.
Step 2: Choose population parameter values
Next, determine population parameter values for generating data from the RICLPM. This requires the specification of:

Phi
: Standardized autoregressive and crosslagged effects for the withinunit components of the model. These values are collected in a matrix, with columns representing predictors and rows representing outcomes. 
within_cor
: A correlation for the withinunit components. 
ICC
: The proportion of variance at the betweenunit level (relative to the total variance). 
RI_cor
: The correlation between the random intercepts.
For our example, the parameter values are set to:
Phi < matrix(c(.4, .1, .2, .3), ncol = 2, byrow = T)
# The .2 refers to our standardized crosslagged effect of interest
within_cor < 0.3
ICC < 0.5
RI_cor < 0.3
If you are unsure if you have specified the Phi
matrix
as intended, you can use the check_Phi()
function to give
you a summary of how the effects in your Phi
are
interpreted.
Steps 35: Perform the power analysis
Steps 3 to 5 are automated by the powRICLPM()
function.
As input, you must provide:
 the desired power level using the
target_power
argument,  the range of sample sizes to simulate the power for using the
search_lower
,search_upper
, andsearch_step
arguments (alternatively, you can specify this directly by providing a vector of sample sizes to thesample_size
argument),  the number of time points for the simulated data using the
time_points
argument,  the population values
Phi
,within_cor
,ICC
, andRI_cor
, and  the number of Monte Carlo replications we want to perform per
experimental condition in the
reps
argument.
Optionally, you can specify:

skewness
andkurtosis
: An integer (vector) that determines the skewness and kurtosis for the simulated observed variables, respectively. Suppose we have reason to believe the \(A\) and \(B\) variables are positively skewed, and have heavy tails (i.e., a higher kurtosis) we can include the argumentsskewness = 1
andkurtosis = 0.5
(default: 0). 
alpha
: A numeric value denoting the significance criterion (default: 0.05). 
seed
: An integer to control the starting point of the random number generator. This is important to use if you want to replicate the results. 
reliability
: A numeric value representing the reliability of the indicators (i.e., the proportion of true score variance) (default: 1). 
bootstrap_reps
: A numeric value, denoting the number of bootstrap replications are used to simulate the uncertainty around the power analysis results (default: 1000). 
estimator
: A character string, denoting the estimator that is used for estimating models (default: “ML”; whenskewness
and/orkurtosis
values are set to nonzero values, it defaults to “MLR”).
The constraints
, bounds
, and
estimate_ME
arguments can be set as well to extend the
basic power analysis setup. These options are further described in the
Vignette Extensions.
For our example, we can perform the power analysis by running:
output < powRICLPM(
target_power = 0.8,
search_lower = 100,
search_upper = 1000,
search_step = 50,
time_points = c(3, 4, 5),
ICC = ICC,
RI_cor = RI_cor,
Phi = Phi,
within_cor = 0.3,
reps = 1000
)
Parallel processing using furrr
Performing a Monte Carlo power analysis with a large number of
replications, and across multiple simulation conditions can be
timeconsuming. To speed up the process, it is recommended to perform
the power analysis across simulation conditions in parallel
(i.e., on multiple cores). The powRICLPM()
function has
implemented future
’s parallel processing capabilities using
the furrr
package.
Load the furrr
package, and use its plan()
function to change the power analysis execution from sequential
(i.e., singlecore, the default), to multisession (i.e.,
multicore). Use the workers
argument to specify how many
cores you want to use. Next, run the powRICLPM
analysis,
and the power analysis will run on the specified number of cores. This
can result in a significant reduction of computing time. For more
information on other parallel execution strategies, see
?furrr::plan()
.
Progress bar using progressr
It can be useful to get an approximation of the progress of the
powRICLPM
analysis while running the code, especially when
running the analysis in parallel. powRICLPM()
has
implemented progress notifications using the progressr
package. Simply put, there are two options through which you can get
progress notifications:
 You can subscribe to progress updates from a specific expression by
wrapping this expression with
with_progress({...})
.  You can subscribe to progress updates from everywhere by running
handlers(global = T)
.
The second option is not fully developed yet for the
furrr
package, so instead I focus on the first.
Implementing the with_progress({...})
option, as well as
parallel execution of the powRICLPM
analysis, results in
the below code for the example:
# Load the furrr and progressr packages
library(furrr)
library(progressr)
# Check how many cores are available
future::availableCores()
# Plan powRICLPM analysis to run on 1 core less than number of available cores
plan(multisession, workers = 7) # For the case of 8 available cores
# Run the powRICLPM analysis
with_progress({ # Subscribe to progress updates
output < powRICLPM(
target_power = 0.8,
search_lower = 100,
search_upper = 1000,
search_step = 50,
time_points = c(3, 4, 5),
ICC = ICC,
RI_cor = RI_cor,
Phi = Phi,
within_cor = 0.3,
reps = 1000
)
})
# Revert back to sequential execution of code
plan(sequential)
For more information about progress notification options using
progressr
for endusers, including auditory and email
updates, see https://progressr.futureverse.org.
Step 6: Summarize results
The powRICLPM()
function creates a
powRICLPM
object: A list with results, upon which we can
call print()
, summary()
, give()
,
and plot()
functions to print, summarize, extract results,
and visualize the results, respectively.
print()
outputs a textual summary of the power analysis
design contained within the object it was called upon. It does not
output any performance metrics computed by the power analysis.
summary()
can be used in one of four ways. First,
summary can be used simply like print()
to get information
about the design of the power analysis (the different experimental
conditions), as well as the number of problems the occurred per
condition (e.g., nonconvergence, fatal estimation errors, or
inadmissible results). Second, by specifying the
parameter = "..."
argument in summary()
, the
function will print the results specifically for that parameter across
all experimental conditions. Third, if you specify a specific
experimental condition using summary()
’s
sample_size
, time_points
, and ICC
arguments, performance measures are outputted for all parameters in that
experimental condition.
# Summary of study design
summary(output)
# Summary of results for a specific parameter, across simulation conditions
summary(output, parameter = "wB2~wA1")
# Summary of all parameter for a specific simulation condition
summary(output, sample_size = 400, time_points = 4, ICC = 0.5)
give()
extracts various bits of information from an
powRICLPM
object. The exact information to be extracted is
given by the what = "..."
argument:

what = "conditions"
gives the different experimental conditions per row, where each condition is defined by a unique combination of sample size, number of time points and ICC. 
what = "estimation_problems"
gives the proportion of fatal errors, inadmissible values, or nonconverged estimations (columns) per experimental conditions (row). 
what = "results"
gives the average estimateaverage
, minimum estimateminimum
, standard deviation of parameter estimatesSD
, the average standard errorSEavg
, the mean square errorMSE
, the average width of the confidence intervalaccuracy
, the coverage ratecoverage
, and the proportion of times the pvalue was lower than the significance criterionpower
. It requires setting theparameter = "..."
argument. 
what = "names"
gives the parameter names contained within thepowRICLPM
object.
# Extract experimental conditions
give(output, what = "conditions")
# Extract estimation problems
give(output, what = "estimation_problems")
# Extract results for crosslagged effect "wB2~wA1"
give(output, what = "results", parameter = "wB2~wA1")
# Extract parameter names
give(output, what = "names")
Finally, plot()
creates a ggplot2
plot for
a specific parameter (specified using the parameter = "..."
argument) with sample size on the xaxis, the simulated power on the
yaxis, lines grouped by number of timepoints, and plots wrapped by
proportion of betweenunit variance. plot()
returns a
ggplot2
object that can be fully customized using
ggplot2
functionality. For example, you can change the
scales, add titles, change geoms, etc. More information about options in
the ggplot2
framework can be found at https://ggplot2book.org/index.html.
In the below example, I add a title and change the labels on the
xaxis:
# Create basic plot of powRICLPM object
p < plot(output, parameter = "wB2~wA1")
p
# Adjust plot aesthetics
p2 < p +
ggplot2::labs(
title = "Power analysis for RICLPM",
caption = "Based on 1000 replications."
) +
ggplot2::scale_x_continuous(
name = "Sample size",
breaks = seq(100, 1000, 100),
guide = ggplot2::guide_axis(n.dodge = 2)
)
p2