I am working on designing a study based on some initial pilot data and would like to conduct a sample size estimation. I have previously used simr with linear mixed effects models. In this case, I think a logistic (glm) or logistic mixed model (glmer) are more appropriate, but am unsure if I am approaching things correctly/defensibly.
I have data from 29 subjects with two continuous predictors (scores on two assessment), a binary outcome of event occurrence in the last 12 months, and the count of events in the last 12 months.
Example data:
Data <- structure(list(Subject = 1:29, Score = c(25, 11.7, 11.7, 43.3,
25, 38.3, 31.7, 35, 30, 36.7, 16.7, 36.7, 11.7, 36.7, 13.3, 96.7,
10, 20, 25, 36.7, 45, 51.7, 68.3, 38.3, 18.3, 5, 18.3, 20, 33.3),
Score1 = c(100L, 100L, 100L, 80L, 100L, 100L, 80L, 90L, 100L,
60L, 100L, 90L, 90L, 100L, 100L, 70L, 100L, 90L, 100L, 40L, 100L,
40L, 10L, 80L, 100L, 100L, 100L, 90L, 100L), Count = c(0, 0,
2, 0, 0, 1, 0, 2, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 4, 3, 2, 12,
1, 0, 1, 0, 1, 0), Bin = c(0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1,
0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0)),
class = "data.frame", row.names = c(NA,
-29L), .Names = c("Subject", "Score", "Score1", "Count", "Bin"))
I have tried I have tried poisson and binomial logistic and logistic mixed models:
# Logistic models
Mod1 <- glm(Count ~ Score + Score1, family = "poisson", data = Data)
Mod2 <- glm(Bin ~ Score + Score1, family = "binomial"(link = logit), data = Data)
# Logistic mixed models
Mod3 <- glmer(Count ~ Score + Score1 + (1|Subject), family = "poisson", data = Data)
Mod4 <- glmer(Bin ~ Score + Score1 + (1|Subject), family = "binomial"(link = logit), data = Data)
I have tried some estimating power using powerSim in the simr package:
powerSim(Mod1, nsim = 100)
Power for predictor 'Score', (95% confidence interval):
5.10% ( 3.82, 6.65)
Test: z-test
Based on 1000 simulations, (0 warnings, 0 errors)
alpha = 0.05, nrow = 29
powerSim(Mod2)
Power for predictor 'Score', (95% confidence interval)
3.20% ( 2.20, 4.49)
Test: z-test
Based on 1000 simulations, (0 warnings, 0 errors)
alpha = 0.05, nrow = 29
powerSim(Mod3)
Power for predictor 'Score', (95% confidence interval):
4.60% ( 3.39, 6.09)
Test: z-test
Effect size for Score is 9.1e-05
Based on 1000 simulations, (11 warnings, 0 errors)
alpha = 0.05, nrow = 29
powerSim(Mod4)
Power for predictor 'Score', (95% confidence interval):
1.70% ( 0.99, 2.71)
Test: z-test
Effect size for Score is -0.0038
Based on 1000 simulations, (715 warnings, 1 error)
alpha = 0.05, nrow = 29
Power is low, so I'd like to know how much larger the sample should be to get closer to 80%.
My first question is whether a standard (glm) is more appropriate than a mixed (glmer) approach?
I've seen some conflicting views on the appropriateness and risk of over-fitting with single observation groupings.
My second questions is how to go about extending the models with the extend function in simr.
How do I extend the model to estimate the effects of more subjects?
I have tried the following for the standard (glm) models:
powerCurve(extend(Mod1, along = "Score", n = 100))
Power for predictor 'Score', (95% confidence interval),
by largest value of Score:
3: 0.00% ( 0.00, 0.37) - 7 rows
14: 2.70% ( 1.79, 3.90) - 23 rows
25: 5.40% ( 4.08, 6.99) - 40 rows
35: 4.70% ( 3.47, 6.20) - 54 rows
46: 4.50% ( 3.30, 5.98) - 74 rows
57: 3.50% ( 2.45, 4.83) - 87 rows
68: 4.20% ( 3.04, 5.64) - 106 rows
78: 3.60% ( 2.53, 4.95) - 122 rows
89: 3.70% ( 2.62, 5.06) - 138 rows
100: 4.50% ( 3.30, 5.98) - 155 rows
powerCurve(extend(Mod2, along = "Score", n = 100))
Power for predictor 'Score', (95% confidence interval),
by largest value of Score:
3: 0.00% ( 0.00, 0.37) - 7 rows
14: 3.80% ( 2.70, 5.18) - 23 rows
25: 5.80% ( 4.43, 7.43) - 40 rows
35: 5.60% ( 4.26, 7.21) - 54 rows
46: 4.90% ( 3.65, 6.43) - 74 rows
57: 5.80% ( 4.43, 7.43) - 87 rows
68: 7.90% ( 6.30, 9.75) - 106 rows
78: 8.40% ( 6.76, 10.29) - 122 rows
89: 10.50% ( 8.67, 12.57) - 138 rows
100: 10.40% ( 8.58, 12.46) - 155 rows
These results don't make sense to me, since each row represents an observation and the model has one observation per subject. The output goes to 155 rows rather than 100. How do I simulate the effects of adding more subjects? Am I going about this all wrong?
For the mixed (glmer) models, I got rank-deficiency errors during simulation, so I dropped "Score1" and used:
Mod5 <- glmer(Count ~ Score + (1|Subject), family = "poisson", data = Data)
powerCurve(extend(Mod5, along = "Subject", n = 100), along = "Subject)
Power for predictor 'Score', (95% confidence interval),
by largest value of Subject:
3: 0.60% ( 0.22, 1.30) - 3 rows
14: 9.00% ( 7.30, 10.95) - 14 rows
25: 60.70% (57.59, 63.74) - 25 rows
35: 64.90% (61.85, 67.86) - 35 rows
46: 81.10% (78.53, 83.48) - 46 rows
57: 86.50% (84.22, 88.56) - 57 rows
68: 87.50% (85.29, 89.49) - 68 rows
78: 94.30% (92.68, 95.65) - 78 rows
89: 96.60% (95.28, 97.63) - 89 rows
100: 96.60% (95.28, 97.63) - 100 rows
Mod6 <- glmer(Bin ~ Score + (1|Subject), family = "binomial"(link = logit, data = Data)
powerCurve(extend(Mod6, along = "Subject", n = 100), along = "Subject")
Power for predictor 'Score', (95% confidence interval),
by largest value of Subject:
3: 0.00% ( 0.00, 0.37) - 3 rows
14: 1.10% ( 0.55, 1.96) - 14 rows
25: 2.10% ( 1.30, 3.19) - 25 rows
35: 2.60% ( 1.71, 3.79) - 35 rows
46: 6.20% ( 4.79, 7.88) - 46 rows
57: 12.70% (10.70, 14.92) - 57 rows
68: 16.10% (13.87, 18.53) - 68 rows
78: 20.50% (18.04, 23.14) - 78 rows
89: 27.00% (24.27, 29.87) - 89 rows
100: 29.90% (27.08, 32.84) - 100 rows
These look more like what I'd expect. But, I had to drop one of the predictors to get it run and more importantly I'm unsure about whether or not a mixed model is appropriate and/or defensible.
