Home > Statistics > Calculating power using Monte Carlo simulations, part 4: Multilevel/longitudinal models

## Calculating power using Monte Carlo simulations, part 4: Multilevel/longitudinal models

In my last three posts, I showed you how to calculate power for a t test using Monte Carlo simulations, how to integrate your simulations into Stata’s power command, and how to do this for linear and logistic regression models. In today’s post, I’m going to show you how to estimate power for multilevel/longitudinal models using simulations. You may want to review my earlier post titled “How to simulate multilevel/longitudinal data” before you read this post.

Our goal is to write a program that will calculate power for different values of the model parameters. For example, figure 1 displays power for different numbers of observations at level 1 and 2 in a longitudinal model.

Figure 1: Estimated power for a multilevel/longitudinal model

In my last post, I suggested the following steps for writing such a program, so let’s follow these steps again for a multilevel/longitudinal model.

In this example, let’s imagine that you are planning a longitudinal study of childrens’ weight and you are interested in the interaction of age and sex.

Step 1: Write down the model

The first step toward simulating power is to write down the model.

${\bf weight}_{it} = \beta_0 + \beta_1({\bf age}_{it}) + \beta_2({\bf female}_i) + \beta_3({\bf age}_{it}\times {\bf female}_i) + u_{0i} + u_{1i}({\bf age}) + e_{it}$

Here i indexes children, t indexes time (age), and we assume that $$u_{0i} \sim N(0,\tau_0)$$, $$u_{1i} \sim N(0,\tau_1)$$, $$e_{it} \sim N(0,\sigma)$$, and $${\rm cov}(\tau_0,\tau_1)=0$$.

We will need to create variables for weight, age, female, and the interaction term age$$\times$$female. We will also need to specify reasonable parameter values for $$\beta_0$$, $$\beta_1$$, $$\beta_2$$, $$\beta_3$$, $$\tau_0$$, $$\tau_1$$, and $$\sigma$$.

Step 2: Specify the details of the covariates

Next, we need to think about the covariates in our model. This is a longitudinal study, so we need to specify the starting age, the length of time between measurements, and the total number of measurements. We also need to consider the proportion of males and females in our study. Are we likely to sample 50% females and 50% males?

Let’s assume that we will measure the childrens’ weight every 6 months for 2 years beginning at age 12. And let’s also assume that the sample will be 50% female. The interaction term age$$\times$$female is easy to calculate once we create variables for age and female.

Step 3: Locate or think about reasonable values for the parameters in your model

Next, we need to think about reasonable values for the parameters in our model. We could choose parameter values based on a review of the literature, results from a pilot study, or from publicly available data.

I have chosen to use data from a study of the weight of Asian children. You can download and describe this dataset by typing

. use http://www.stata-press.com/data/r16/childweight.dta, clear
(Weight data on Asian children)

. describe

Contains data from http://www.stata-press.com/data/r16/childweight.dta
obs:           198                          Weight data on Asian children
vars:             5                          23 May 2016 15:12
size:         3,168                          (_dta has notes)
--------------------------------------------------------------------------------
storage   display    value
variable name   type    format     label      variable label
--------------------------------------------------------------------------------
id              int     %8.0g                 child identifier
age             float   %8.0g                 age in years
weight          float   %8.0g                 weight in Kg
brthwt          int     %8.0g                 Birth weight in g
girl            float   %9.0g      bg         gender
--------------------------------------------------------------------------------
Sorted by: id  age


This dataset includes variables for age measured in years (age), weight measured in kilograms (weight), and gender (girl). We can verify that these are longitudinal data by listing some of the observations.

. list in 1/5

+---------------------------------------+
| id       age   weight   brthwt   girl |
|---------------------------------------|
1. | 45   .136893    5.171     4140    boy |
2. | 45   .657084    10.86     4140    boy |
3. | 45   1.21834    13.15     4140    boy |
4. | 45   1.42916     13.2     4140    boy |
5. | 45   2.27242    15.88     4140    boy |
+---------------------------------------+


Child 45 is a boy whose weight was recorded five times over two years. Notice that age is not stored as absolute age. It is stored as decimals representing the number of years since an unknown point in time. This isn’t a problem because we are interested in the change of weight over time. We can fit our model to this dataset to estimate the parameters.

. mixed weight c.age##i.girl || id: age, stddev

Performing EM optimization:

Iteration 0:   log likelihood = -339.79909
Iteration 1:   log likelihood = -339.41232
Iteration 2:   log likelihood = -339.41033
Iteration 3:   log likelihood = -339.41032

Computing standard errors:

Mixed-effects ML regression                     Number of obs     =        198
Group variable: id                              Number of groups  =         68

Obs per group:
min =          1
avg =        2.9
max =          5

Wald chi2(3)      =     680.51
Log likelihood = -339.41032                     Prob > chi2       =     0.0000

------------------------------------------------------------------------------
weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
age |   3.588488   .1892333    18.96   0.000     3.217598    3.959379
|
girl |
girl  |   -.468204   .2954284    -1.58   0.113    -1.047233     .110825
|
girl#c.age |
girl  |  -.2397793    .267607    -0.90   0.370    -.7642793    .2847208
|
_cons |   5.351429   .2076606    25.77   0.000     4.944421    5.758436
------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Independent              |
sd(age) |   .5662947   .1219764      .3712777    .8637461
sd(_cons) |   .2384136   .4034962      .0086445    6.575403
-----------------------------+------------------------------------------------
sd(Residual) |   1.165676   .0762599      1.025395    1.325148
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 20.38                 Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.


We can use the parameter estimates from this model to help us choose the parameters for our simulation. For example, we could begin our simulation by setting the simulation parameters equal to the parameter estimates above: $$\beta_0=5.35$$, $$\beta_1=3.59$$, $$\beta_2=-0.47$$, $$\beta_3=-0.24$$, $$\tau_0=0.24$$, $$\tau_1=0.57$$, and $$\sigma=1.17$$.

Step 4: Simulate a single dataset assuming the alternative hypothesis, and fit the model

Next, we create a simulated dataset based on our assumptions about the model under the alternative hypothesis. We will simulate 5 observations at 6-month increments for 200 children. You may wish to read my previous blog post titled “How to simulate multilevel/longitudinal data” if you have not done this before.

The code block below is similar to the code we used to create the data for our linear regression model, but there are two features that are unique to simulating multilevel data. First, we generate the child-level variables child and female as well as the random deviations u_0i and u_1i. Second, we use expand to make five copies of the observations for each child. Third, we use bysort child: to generate age within each child. Fourth, we generate the observation-level variable interaction and the observation-level error (e_ij). Finally, we generate the variable weight using a linear combination of parameter values and the other variables that we generated.

set seed 16
clear
set obs 200
generate child = _n
generate female = rbinomial(1,0.5)
generate u_0i = rnormal(0,0.25)
generate u_1i = rnormal(0,0.60)
expand 5
bysort child: generate age = (_n-1)*0.5
generate interaction = age*female
generate e_ij = rnormal(0,1.2)
generate weight = 5.35 + 3.6*age + (-0.5)*female + (-0.25)*interaction ///
+ u_0i + age*u_1i + e_ij


Let’s list the data for child 1 to check our work. The simulated data includes five observations for weight, age, and female. Our dataset also includes the random deviations that we would not observe in a real dataset.

. list child weight age female u_0i u_1i e_ij if child==1

+------------------------------------------------------------------+
| child     weight   age   female       u_0i       u_1i       e_ij |
|------------------------------------------------------------------|
1. |     1   6.957936     0        1   .0530039   .6878794   2.054933 |
2. |     1   8.517486    .5        1   .0530039   .6878794   1.595542 |
3. |     1   10.28054     1        1   .0530039   .6878794   1.339654 |
4. |     1   11.48091   1.5        1   .0530039   .6878794   .5210837 |
5. |     1   14.25277     2        1   .0530039   .6878794   1.274008 |
+------------------------------------------------------------------+


We can then use mixed to fit a model to our simulated data.

. mixed weight age i.female c.age#i.female || child: age , stddev nolog noheader

------------------------------------------------------------------------------
weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
age |   3.654038   .1029402    35.50   0.000     3.452279    3.855797
1.female |  -.3765285   .1389315    -2.71   0.007    -.6488291   -.1042279
|
female#c.age |
1  |  -.4101041    .142071    -2.89   0.004    -.6885581   -.1316501
|
_cons |   5.309637   .1006654    52.75   0.000     5.112336    5.506937
------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
child: Independent           |
sd(age) |   .6628522   .0544682      .5642497    .7786855
sd(_cons) |   .3342406   .1071213       .178343    .6264151
-----------------------------+------------------------------------------------
sd(Residual) |   1.190916   .0316823      1.130411    1.254659
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 157.47                Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.


We can then test the null hypothesis that the interaction term equals zero using a likelihood-ratio test. Recall from my last post that we can calculate a likelihood-ratio test by fitting the model with and without the intereaction term female#c.age and use lrtest to calculate the test statistic.

mixed weight age i.female c.age#i.female || child: age , stddev nolog noheader
estimates store full
mixed weight age i.female || child: age , stddev nolog noheader
estimates store reduced
lrtest full reduced


The p-value for our test is 0.0041, so we would reject the null hypothesis that the interaction term equals zero.

. lrtest full reduced

Likelihood-ratio test                                 LR chi2(1)  =      8.23
(Assumption: reduced nested in full)                  Prob > chi2 =    0.0041


In this example, both models converged. But sometimes models fail to converge for some of the datasets when we fit models to many simulated datasets. This is often because the process of simulating data can result in data that are incompatible with the model. But there may be other reasons, and you should consider them carefully if you observe many models that fail to converge.

Lack of model convergence can also cause our simulation to run far longer than necessary. By default, mixed will run 16,000 maximum-likelihood iterations before it stops. In the code block below, I have added options and lines to handle models that fail to converge. The option iter(200) in the mixed commands will stop them from running if a model fails to converge in the first 200 iterations. mixed will store a value of 1 to e(converged) if the model converges and 0 otherwise. The second line of the code block stores the value of e(converged) from the first model to a local macro named conv1. The fifth line of the code block stores the value of e(converged) from the second model to a local macro named conv2. The last line of the code block stores a value to the local macro reject. The value that is stored in reject is determined by the conditional function cond(). If both mixed models converged and conv1’+conv2’==2, then r(p)<alpha’ will be evaluated. A value of 1 will be stored to reject if the likelihood-ratio p-value, r(p), is less than the alpha level specified in alpha’ and 0 otherwise. If either of the mixed models fails to converge and conv1’+conv2′ does not equal 2, then a missing value will be stored to reject.

mixed weight age i.female c.age#i.female || child: age, iter(200)
local conv1 = e(converged)
estimates store full
mixed weight age i.female || child: age, iter(200)
local conv2 = e(converged)
estimates store reduced
lrtest full reduced
local reject = cond(conv1'+conv2'==2, r(p)<alpha', .)


Step 5: Write a program to create the datasets, fit the models, and use simulate to test the program

Next, let’s write a program that creates datasets under the alternative hypothesis, fits mixed models, tests the null hypothesis of interest, and uses simulate to run many iterations of the program.

The code block below contains the syntax for a program named simmixed. The default parameter values in the syntax command are similar to the parameters that we estimated using the Asian children dataset. And we use lrtest to test the null hypothesis that the parameter for the age$$\times$$female interaction equals zero.

capture program drop simmixed
program simmixed, rclass
version 16
// PARSE INPUT
syntax, n1(integer)             ///
n(integer)              ///
[ alpha(real 0.05)        ///
intercept(real 5.35)    ///
age(real 3.6)           ///
female(real -0.5)       ///
interact(real -0.25)    ///
u0i(real 0.25)          ///
u1i(real 0.60)          ///
eij(real 1.2) ]

// COMPUTE POWER
quietly {
drop _all
set obs n'
generate child = _n
generate female = rbinomial(1,0.5)
generate u_0i = rnormal(0,u0i')
generate u_1i = rnormal(0,u1i')
expand n1'
bysort child: generate age = (_n-1)*0.5
generate interaction = age*female
generate e_ij = rnormal(0,eij')
generate weight = intercept' + age'*age + female'*female + ///
interact'*interaction  + u_0i + age*u_1i + e_ij

mixed weight age i.female c.age#i.female || child: age, iter(200)
local conv1 = e(converged)
estimates store full
mixed weight age i.female || child: age, iter(200)
local conv2 = e(converged)
estimates store reduced
lrtest full reduced
local reject = cond(conv1' + conv2'==2, (r(p)<alpha'), .)
}
// RETURN RESULTS
return scalar reject = reject'
return scalar conv = conv1'+conv2'
end


We then use simulate to run simmixed 10 times using the default parameter values for 5 observations on each of 200 children.

. simulate reject=r(reject) converged=r(conv), reps(10) seed(12345):
>              simmixed, n1(5) n(200)

command:  simmixed, n1(5) n(200)
reject:  r(reject)
converged:  r(conv)

Simulations (10)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
..........


simulate saved the results of the hypothesis tests to a variable named reject. The mean of reject is our estimate of the power to test the null hypothesis that the age$$\times$$sex interaction term equals zero, assuming that the weight of 200 children is measured 5 times.

. summarize reject

Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
reject |         10          .3    .4830459          0          1


Step 6: Write a program called power_cmd_mymethod, which allows you to run your simulations with power

We could stop with our quick simulation if we were interested only in a specific set of assumptions. But it’s easy to write an additional program named power_cmd_simmixed that will allow us to use Stata’s power command to create tables and graphs for a range of sample sizes.

capture program drop power_cmd_simmixed
program power_cmd_simmixed, rclass
version 16
// PARSE INPUT
syntax, n1(integer)             ///
n(integer)              ///
[ alpha(real 0.05)        ///
intercept(real 5.35)    ///
age(real 3.6)           ///
female(real -0.5)       ///
interact(real -0.25)    ///
u0i(real 0.25)          ///
u1i(real 0.60)          ///
eij(real 1.2)           ///
reps(integer 1000) ]

// COMPUTE POWER
quietly {
simulate reject=r(reject), reps(reps'):                            ///
simmixed, n1(n1') n(n') alpha(alpha') intercept(intercept')     ///
age(age') female(female') interact(interact')          ///
u0i(u0i') u1i(u1i') eij(eij')
summarize reject
}

// RETURN RESULTS
return scalar power = r(mean)
return scalar n1 = n1'
return scalar N = n'
return scalar alpha = alpha'
return scalar intercept = intercept'
return scalar age = age'
return scalar female = female'
return scalar interact = interact'
return scalar u0i = u0i'
return scalar u1i = u1i'
return scalar eij = eij'
end


Step 7: Write a program called power_cmd_mymethod_init so that you can use numlists for all parameters

It’s also easy to write a program named power_cmd_simmixed_init that will allow us to simulate power for a range of values for the parameters in our model.

capture program drop power_cmd_simmixed_init
program power_cmd_simmixed_init, sclass
version 16
sreturn clear
// ADD COLUMNS TO THE OUTPUT TABLE
sreturn local pss_colnames "n1 intercept age female interact u0i u1i eij"
// ALLOW NUMLISTS FOR ALL PARAMETERS
sreturn local pss_numopts  "n1 intercept age female interact u0i u1i eij"
end


Using power simmixed

Now, we can use power simmixed to simulate power for a variety of assumptions. The example below simulates power for a range of sample sizes at both levels 1 and 2. Level 2 sample sizes range from 100 to 500 children in increments of 100. At level 1, we consider 5 and 6 observations per child.

. power simmixed, n1(5 6) n(100(100)500) reps(1000)
>               table(n1 N power)
>               graph(ydimension(power) xdimension(N) plotdimension(n1)
>               xtitle(Level 2 Sample Size) legend(title(Level 1 Sample Size)))
xxxxxxxxxxxxxxxxxxxxxxxxxxx
Estimated power
Two-sided test

+-------------------------+
|      n1       N   power |
|-------------------------|
|       5     100   .2629 |
|       6     100    .313 |
|       5     200    .397 |
|       6     200    .569 |
|       5     300    .621 |
|       6     300    .735 |
|       5     400    .734 |
|       6     400    .855 |
|       5     500    .828 |
|       6     500    .917 |
+-------------------------+
`

Figure 2: Estimated power for a multilevel/longitudinal model

The table and graph above indicate that 80% power is achieved with three combinations of sample sizes. Given our assumptions, we estimate that we will have at least 80% power to detect an interaction parameter of -0.25 with 400 children measured 6 times each and 500 children measured 5 or 6 times each.

Running our simulations with power gives us the flexibility to use power simmixed to estimate power for different values of any model parameter we wish. For example, we could consider different values of the interaction term along with different numbers of observations per child.

In this blog post, I showed you how to simulate statistical power for the interaction term in a multilevel regression model. You may be interested in simulating power for variations of these models, and you can modify this example for your own purposes. In my next post, I will show you how to simulate power for path coefficients in structural equation models.

Categories: Statistics Tags: