pitakakariki / simr

Power Analysis of Generalised Linear Mixed Models by Simulation
68 stars 19 forks source link

Simr suddenly not replicating previous results #253

Open ajj1988 opened 1 year ago

ajj1988 commented 1 year ago

Hi,

I ran simr ~ 6months ago and got a simulated power calculation no problem (power was ~95%). However, upon returning to my code and rerunning the models I am getting power of 0. I haven't changed the code or data set.

library(lme4) library(lmerTest) library(simr)

options(scipen = 999)

#################power analysis (use data file 'power_calculation_subset.xlsx') power_analysis_data<-readxl::read_excel('power_calculation_subset.xlsx')

table(power_analysis_data$Attentio_CHECK) ###quick look at the numbers of careless responder power_analysis_data$Attentio_CHECK<-ifelse(power_analysis_data$Attentio_CHECK==1,0,1) ####recode as 0 (not careless) and 1 (careless) for analyses table(power_analysis_data$Attentio_CHECK) ###check it worked

final_book<-power_analysis_data

next steps remove any rows with missing data for age,sex and audit (complete cases)

final_book<-final_book[complete.cases(final_book$AGE),] final_book<-final_book[complete.cases(final_book$SEX),] final_book<-final_book[complete.cases(final_book$AUDIT_total),] final_book$x<-final_book$AUDIT_total #####recode AUDIT as 'x' so it's easier to deal with in the simR package

###############################https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12504 ###############################USE SIMR to get power calculations, based on AUDIT SCORE

single_test_model<-glm(final_book$Attentio_CHECK ~ x, family = binomial, data = final_book) summary(single_test_model) #####single level model - no random intercept for study

test_model<-glmer(Attentio_CHECK ~ x + (1|STUDYID), family = binomial, data = final_book) summary(test_model) ####test model with a ranom intercept for study confint(test_model)

fixef(test_model)["x"]<- .07 ####this is the log coefficient for AUDIT ~ careless (see test model) power<-powerSim(test_model) powerSim(test_model2, nsim = 100)

this tells us how much power we had to detect this effect

power #####note, that as this is based on simulations it's possible you won't get the same value(s) as what are reported in the manuscript we have written. However, note that you can run this multiple times and is highly unlikely power < 80%

lowest i've seen so far is 91% - highest 96%

################################# 5 more studies would provide 100% power

test_model_2<-extend(test_model, along = "STUDYID", n = 1000) ###this adds five more studies power2<-powerSim(test_model_2, nsim = 100) summary(power2)

power_calculation_subset.xlsx

pitakakariki commented 1 year ago

Not sure why this used to work, but the problem now is with the missing value handling.

There are a number of cases where Attentio_CHECK is missing that are still in final_book, and if you remove those cases your simulation should work again.

ajj1988 commented 1 year ago

Perfect. Thank you. And thank you for the fantastic R package :)

Andrew Jones Lecturer University of Liverpool Tel: 01517941120 Office Hours: Mondays and Tuesdays 10 – 12.

On 24 Jan 2023, at 04:11, Peter Green @.***> wrote:



Not sure why this used to work, but the problem now is with the missing value handling.

There are a number of cases where Attentio_CHECK is missing that are still in final_book, and if you remove those cases your simulation should work again.

— Reply to this email directly, view it on GitHubhttps://github.com/pitakakariki/simr/issues/253#issuecomment-1401371196, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AR4O2LDD3SZYZ6I67LOEK7DWT5I5ZANCNFSM6AAAAAAUBP6WJE. You are receiving this because you authored the thread.Message ID: @.***>

wiedenhoeft commented 1 year ago

I'm having the same issue. Used to work with R 4.1.2, but all power is 0% with R 4.2.3, without any changes to my script. I do not have missing values though, and I am using fixed seeds for reproducibility.

pitakakariki commented 1 year ago

Hi John, could you please run your analysis in both R 4.1.2 and R 4.2.3+ and show me the results so I can get a sense of what's broken?