PowerSim on interaction term and second order modelling

[fowcai]

Hi all, I'm new to GitHub so apologies if my entry is not formatted correcting for you to help me appropriately. Happy to edit/add info if needed!

My study design is 4 timepoints with 30, 30, 16, and 11 subjects at timepoints 1 through 4, with some split of sex between them. There is a heavy bias towards females at the last 2 timepoints. The data I'm looking to model is change in metabolite levels (we test 35 of them so each model has to be re-run for each metabolite, yay!) with age and/or sex, while controlling for water.linewidth and random effect per subject, and weighting based on the reliability of the metabolite measurement.

These are the ideal models that I would be running if I am powered to do so:

    model.int = lmer(metabolite ~ ns(age,1)*sex + water.lw + (1|subject), data=LME_MRS_abs_WT, weights=1/((stdev.absolute), center=FALSE). 

model2.int = lmer(metabolite ~ ns(age,2)*sex + water.lw + (1|subject), data=LME_MRS_abs_WT, weights = 1/((stdev.absolute), center=FALSE)) 

ns indicates a natural spline modelling the age, either examining change in metabolite level via first order (i.e. linear) or second order effects of age.

I would like to test the power for the interaction terms as well as the main effects for both models. I am 99.9 % certain that I am underpowered for the interaction terms, and possible even the sex effect on it's own, but would like to justify this.

I tested the power of the AGE term via the following comparison of model0 (no age term) against model1. I am focusing on only main effects and no interaction right now to ensure that in its simplest form, I can use the model to detect age, and possibly sex effects if I'm lucky.

model0 = lmer(metabolite ~ sex + water.lw +  (1|subject),etc
model1 = lmer(metabolite ~ ns(age,1) + sex + water.lw + (1|subject), etc

Set smallest effect size of interest:

SESOI = 0.03
fixef(model1)["ns(age, 1)"] = SESOI
power_age=powerSim(model1, test=compare(model0, method="lr"), nsim=300)

This gives me a power of 94%, which is great. Just to confirm, this can be described as "the power afforded to this model by a linear age variable to detect 3% changes in metabolite level, keeping all other variables in the equation constant" , correct?

I then do the exact same process for the main effect of sex, using the following model0 and the same model1.

model0 = lmer(metabolite ~ ns(age,1) + water.lw + (1|subject), etc
model1 = lmer(metabolite ~ ns(age,1) + sex + water.lw + (1|subject), etc 

SESOI = 0.03
fixef(model1)["sexF"] = SESOI

power_sex = powerSim(model1, test=compare(model0, method="lr"), nsim=300))

This gives me 75% power, which is also good.

Now I want to know if I have similar power when I go to work these terms into an itneraction - I assume I won't be want to see.

model0 = lmer(metabolite ~ ns(age,1) + sex + water.lw + (1|subject), etc
model.int = lmer(metabolite ~ ns(age,1)*sex + water.lw + (1|subject), etc
SESOI = 0.03
fixef(lme.int)["ns(age, 1)"] = SESOI
power.int = powerSim(lme.int, test=compare(lme0.int, method="lr"), nsim=300)

which gives me a power of 53% - I believe the way to interpret this is that the inclusion of the interaction term, as opposed to the two variables separately, only provides a 53% power to detect a 3% change in metabolite level with age, and I am therefore underpowered to use the interaction term in the model. Is that right?

AND afterwards, I do the same comparison for detecting a 3% change in metabolite level between sexes, again, assuming I'm interpreting this correctly. The following gives me 49% power.

SESOI = 0.03
fixef(lme.int)["sexF"] = SESOI
power.int=powerSim(lme.int, test=compare(lme0.int, method="lr"),nsim=300)

So, what I am understanding from this is that I am underpowered to detect a 3% changes in metabolite level with age and sex when the terms are interacting, but as separate main effects I am adequately powered, correct? I attempted to run powerCurve to see if more subjects would improve the power of the interaction term but I don't think it works with "longitudinal data" or at least unbalanced subject numbers. Am I correct in assuming that?

Now, moving into the second order effects, this is the model that I want to test. Knowing that my interaction term is underpowered I am keeping age and sex as separate effects.

model0 = lmer(metabolite ~ ns(age,1) + sex + water.lw + (1|subject)
model1 = lmer(metabolite ~ ns(age,2) + sex + water.lw + (1|subject)
SESOI = 0.03
fixef(model1)["ns(age, 2)1"] = SESOI
power_age_1=powerSim(model1, test=compare(model0, method="lr"),nsim=300)

which gives a power of 35%.

same idea for the second order effect, which gives a power of 37%.

fixef(model1)["ns(age, 2)2"] = SESOI
power_age_2 = powerSim(model1, test=compare(model0, method="lr"),nsim=300)

So, by including a second order modelling of age (which simultaneously evaluates first and second order effects of age) I have 43% power to detect linear change in metabolite level with age, and 42% power to detect a second order change in metabolite level with age, yes?

Finally, to test how this clearly underpowered second order equation will evaluate metabolite differences between sex, I did the following, which surprisingly, gives a power of 69%, which is a good 20% higher than that given for the effect of age with this model.

model0 = lmer(metabolite ~ ns(age,2) + water.lw + (1|subject)
model1 = lmer(metabolite ~ ns(age,2) + sex + water.lw + (1|subject)
SESOI = 0.03
fixef(model1)["sexF"] = SESOI
power_age_1 = powerSim(model1, test=compare(model0, method="lr"),nsim=300)

So, after this painfully long novel, I have some questions:

1. Do the above methods make sense? Am I testing the "main effects" of age and sex within the interaction and within the second order equation, correctly? No one I know works with rsim and there's only so far the documentation can take me, given my specific data format and model.
2. As a general question, does it always require more power to detect second order effects as opposed to first order? Even if the data is better fit using a second order model (via AIC comparison testing)?
3. Should the powerCurve tool work on this data or given the longitudinal nature with unbalanced subjects per timepoint remove that functionality? Can I use it to see what my power would be if I had 30 subjects per timepoint rather than 30, 30, 16, and 11?
4. I will be running these analysis in 35 different metabolites, each of which will give me different power for the main effects. If I am adequately powered to detect a sex effect in one metabolite but not another, I thought it would make sense to present all of the p-values and FDR-corrected values next to one another, with the power of each test next to the effect so the reader knows immediately in some cases sex acted as a simple covariate, and in other cases we were powered to actually make definitely conclusions about differences between sex. Would that be an appropriate way of dealing with this or should I simply not show/discuss main effects of sex unless I am powered in all 35 metabolites?

Thank you so much for your help and your time. I know this is super long but I wanted to make sure all of the information was here for you and others to use and understand!

[pitakakariki]

Question 1

I haven't worked with ns before, but most of what you're doing with powerSim look sensible. I think your interpretation of the effect size is mistaken though. 0.03 is a 3% change in a Poisson GLMM, but in a Gaussian LMM (i.e. lmer) effects are additive rather than multiplicative.

So e.g. when you set the fixed effect for sexF to 0.03, you're saying that the expected value for metabolite in observations coded female is 0.03 higher than in males, in whatever units you measure metabolite in. So if a typical value for metabolite was 10 units, that's only a 0.3% increase. If a typical value is 0.1 units, 0.03 is a 30% increase.

Unfortunately I can't help you with interpretation of the ns effect sizes. But if the effects aren't easily interpretable you could use predict to get a sense of what sort of difference is implied.

[pitakakariki]

Question 2

That's a good question. I'd expect that second order effects would usually be harder to detect since they would tend to be smaller and their estimates would tend to have higher standard errors. But you could have exceptions where the first order effects are small and the second order effects are large.

[pitakakariki]

Question 3

At this point I would not recommend using powerCurve on this kind of unbalanced data. It will "work", sort of, but some of the decisions it will make in the background will a) opaque, and b) outright bad.

What you can do is replace your data with a balanced design like this:

getData(model1) <- balanced_data_frame

Then you can e.g. increase your subjects per time point with extend. If you use extend I would also recommend using getData and running some crosstabs to make sure it did what you asked.

[pitakakariki]

Question 4

I think this mostly depends on norms within your field, but I think your plan for this is good. I would also include confidence intervals for your estimates.

[FelipePeg]

Hi, I am using a previous study data to predict the sample size needed to have a power of 0.8 in a new study but I am struggling to find out the error: model1 <- glmer (correct ~ repet * lex + (1|subj) + (1|item), data= ff4,subset=(ff4$repet!='repeat'), family = "binomial")

test<-powerSim(model1, test=fixed("repettranspose:lexword", method="z"), nsim=20)

with the following warning message: Warning message:=====================================================================| In observedPowerWarning(sim) : This appears to be an "observed power" calculation

And then:

head(test$errors) stage index message 1 Simulating 1 length mismatch in beta (4!=6) 2 Simulating 2 length mismatch in beta (4!=6) 3 Simulating 3 length mismatch in beta (4!=6) 4 Simulating 4 length mismatch in beta (4!=6) 5 Simulating 5 length mismatch in beta (4!=6) 6 Simulating 6 length mismatch in beta (4!=6)

print(test) Power for predictor 'repettranspose:lexword', (95% confidence interval): 0.00% ( 0.00, 16.84)

Test: z-test Effect size for repettranspose:lexword is -1.0

Based on 20 simulations, (0 warnings, 20 errors) alpha = 0.05, nrow = NA

Time elapsed: 0 h 0 m 2 s

nb: result might be an observed power calculation

Any idea would be very helpful! thx Felipe

[pitakakariki]

You'll need to subset the data directly rather than passing it as a parameter to glmer.

[FelipePeg]

Thanks a lot for your prompt answer pitakakariki !! It is a bit weird for me that the syntax 'subset' runs smoothly on the lme4 but when we use the obtained model in simR power analysis it won't work. Anyway, I followed your advice and created a new dataframe with the selected data, run the LME model, then the simR and ... It works!!

However, the results obtained with powerCurve(model1, along='subj', test=fixed("repettranspose:lexword"), nsim=200) have a bizarre "all or none" pattern (below). Any thoughts ? (ps. I had 28 subjects in my original study so I suspect that it has to do with it):

(...) by number of levels in subj: 3: 0.00% ( 0.00, 16.84) - 463 rows 6: 0.00% ( 0.00, 16.84) - 937 rows 9: 0.00% ( 0.00, 16.84) - 1414 rows 11: 0.00% ( 0.00, 16.84) - 1720 rows 14: 0.00% ( 0.00, 16.84) - 2195 rows 17: 0.00% ( 0.00, 16.84) - 2662 rows 20: 0.00% ( 0.00, 16.84) - 3137 rows 22: 0.00% ( 0.00, 16.84) - 3449 rows 25: 0.00% ( 0.00, 16.84) - 3921 rows 28: 100.0% (83.16, 100.0) - 4358 rows

[pitakakariki]

Looks like it was running into some errors maybe?

Try head(lastResult()$errors) to see what happened.

[FelipePeg]

It gives that:

head(lastResult()$errors) stage index message 1 Fitting 1 variable lengths differ (found for 'subj') 2 Fitting 2 variable lengths differ (found for 'subj') 3 Fitting 3 variable lengths differ (found for 'subj') 4 Fitting 4 variable lengths differ (found for 'subj') 5 Fitting 5 variable lengths differ (found for 'subj') 6 Fitting 6 variable lengths differ (found for 'subj')

[FelipePeg]

One additional information: when I try to extend the number of subjects, I have the following error:

model2 <- extend(model1, along='subj', n=40)

Error in oldValues[] <- as.character(unique(object[[along]])) : replacement has length zero

    |    

[pitakakariki]

Is subj in your data frame or outside of it?

[FelipePeg]

You found it !! It was outside of it. So I used the dataframe$subject_nr as the factor and it worked!! Thanks a lot for your help!