#Get packages library(WebPower) library(Superpower) library(pwr) library(simr) library(dplyr) library(longpower) library(lme4) library(nlme) library(gee) library(powerSurvEpi) library(survival) library(semPower) #reset graphs dev.off() #------------------------------------------------------------------------------------------ #Repeated Measures ANOVA: WebPower c1 <- sqrt(4/(1-0.5)); c1 wp.rmanova(ng = 3, nm = 4, f = 0.25*c1, nscor = 1, alpha = 0.05, power = 0.8, type = 1) c0 <- sqrt(4/(1+(4-1)*0.5)); c0 wp.rmanova(ng = 3, nm = 4, f = 0.25*c0, nscor = 1, alpha = 0.05, power = 0.8, type = 0) #------------------------------------------------------------------------------------------ #REPEATED MEASURES ANOVA #Repeated Measures ANOVA: Superpower #Stock Example design <- "3b*3b" n <- 20 mu <- c(20, 20, 20, 20, 20, 20, 20, 20, 25) # Enter means in the order that matches the labels below. sd <- 5 labelnames <- c("Factor_A", "a1", "a2", "a3", "Factor_B", "b1", "b2", "b3") # # the label names should be in the order of the means specified above. design_result <- ANOVA_design(design = design, n = n, mu = mu, sd = sd, labelnames = labelnames) plot(design_result) res = ANOVA_exact2(design_result, alpha_level = 0.05) plot_power(design_result, desired_power = 80, min_n = 15, max_n = 100) #My example design <- "4b*3b" n <- 15 mu <- c(20, 20, 20, 20, 19, 21, 20, 18, 22, 20, 17, 23) sd <- 5 labelnames <- c("Repeats", "T0", "T1", "T2", "T3", "Groups", "g1", "g2", "g3") design_result <- ANOVA_design(design = design, n = n, mu = mu, sd = sd, labelnames = labelnames) plot(design_result) res = ANOVA_exact2(design_result, alpha_level = 0.05) plot_power(design_result, desired_power = 80, min_n = 15, max_n = 100) #------------------------------------------------------------------------------------------ #GENERALIZED AND MIXED MODELS #LM and GLM #(https://slcladal.github.io/pwr.html#Advanced_Power_Analysis) pwrglm1 <- pwr.f2.test(u = 1, f2 = .02, sig.level = 0.05, power=0.80) pwrglm1 #power calculation ssgglm1 <-round((pwrglm1$v+2)/2,0) ssgglm1 #sample size per group #LMM #get data regdat <- base::readRDS(url("https://slcladal.github.io/data/regdat.rda", "rb")) head(regdat) #generate model fixed <-c(0.52, 0.52, 0.52, 0.52, 0.52); rand <-list(0.5, 0.1); res<-2 m1 <- makeLmer(y ~ (1|Sentence) + (1|ID) + Group * SentenceType + WordOrder, fixef=fixed, VarCorr=rand, sigma=res, data=regdat) #power analysis sim_wo <- powerSim(m1, nsim=20, test = fcompare(y ~ Group * SentenceType)) sim_wo sim_gst <- powerSim(m1, nsim=20, test = fcompare(y ~ WordOrder + Group + SentenceType)) sim_gst #power curves m1_as <- extend(m1,along="Sentence",n=120) pcurve_m1_as_gst <- powerCurve(m1_as, test=fcompare(y ~ WordOrder + Group + SentenceType), along="Sentence", nsim = 20, breaks=seq(20, 120, 20)) plot(pcurve_m1_as_gst) m1_ap <- extend(m1, along="ID", n=120) pcurve_m1_ap_gst <- powerCurve(m1_ap, test=fcompare(y ~ Group + SentenceType + WordOrder), along = "ID", nsim=20) plot(pcurve_m1_ap_gst) #GLMM #get data simdat <- data.frame(sub <- rep(paste0("Sub", 1:10), each = 20), cond <- rep(c(rep("Control", 10),rep("Test", 10)), 10), itm <- as.character(rep(1:10, 20))) %>% rename(Subject = 1,Condition = 2,Item = 3) %>% mutate_if(is.character, factor) head(simdat,20) #generate model fixed <-c(0.52, 0.52); rand<-list(0.5, 0.1) m2 <- makeGlmer(y ~ (1|Subject) + (1|Item) + Condition, family = "binomial", fixef=fixed, VarCorr=rand, data=simdat) #power analysis rsim_m2_c <- powerSim(m2, fixed("ConditionTest", "z"), nsim=20) rsim_m2_c #power curves m2_ai <- extend(m2, along="Item", n=40) pcurve_m2_ai_c <- powerCurve(m2_ai, fixed("ConditionTest", "z"), along = "Item", nsim = 20, breaks=seq(10, 40, 5)) plot(pcurve_m2_ai_c) m2_as <- extend(m2, along="Subject", n=40) pcurve_m2_as_c <- powerCurve(m2_as, fixed("ConditionTest", "z"), along = "Subject", nsim = 20, breaks=seq(10, 40, 5)) plot(pcurve_m2_as_c) m2_as <- extend(m2, along="Subject", n=30) m2_asi <- extend(m2_as, along="Item", n=30) pcurve_m2_asi_c <- powerCurve(m2_asi, fixed("ConditionTest", "z"), along = "Subject", breaks = seq(5, 30, 5), nsim = 20) plot(pcurve_m2_asi_c) #------------------------------------------------------------------------------------------ #LONGITUDINAL AND SURVIVAL ANALYSIS #Longitudinal # General lmmpower(delta=1.5, t = seq(0,1.5,0.25), sig2.i = 55, sig2.s = 24, sig2.e = 10, cov.s.i=0.8*sqrt(55)*sqrt(24), power = 0.80) lmmpower(n=208, t = seq(0,1.5,0.25), sig2.i = 55, sig2.s = 24, sig2.e = 10, cov.s.i=0.8*sqrt(55)*sqrt(24), power = 0.80) lmmpower(beta = 5, pct.change = 0.30, t = seq(0,1.5,0.25), sig2.i = 55, sig2.s = 24, sig2.e = 10, cov.s.i=0.8*sqrt(55)*sqrt(24), power = 0.80) # General Continued fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) lmmpower(fm1, pct.change = 0.30, t = seq(0,9,1), power = 0.80) fm2 <- lme(Reaction ~ Days, random=~Days|Subject, sleepstudy) lmmpower(fm2, pct.change = 0.30, t = seq(0,9,1), power = 0.80) # random intercept only fm3 <- lme(Reaction ~ Days, random=~1|Subject, sleepstudy) lmmpower(fm3, pct.change = 0.30, t = seq(0,9,1), power = 0.80) fm4 <- gee(Reaction ~ Days, id = Subject, data = sleepstudy, corstr = "exchangeable") lmmpower(fm4, pct.change = 0.30, t = seq(0,9,1), power = 0.80) # An Alzheimer's Disease example using ADAS-cog pilot estimates sig2.i <- 55 # var of random intercept sig2.s <- 24 # var of random slope sig2.e <- 10 # residual var cov.s.i <- 0.8*sqrt(sig2.i)*sqrt(sig2.s) # covariance of slope and intercept cov.t <- function(t1, t2, sig2.i, sig2.s, cov.s.i){ sig2.i + t1*t2*sig2.s + (t1+t2)*cov.s.i} t <- seq(0,1.5,0.25) n <- length(t) R <- outer(t, t, function(x,y){cov.t(x,y, sig2.i, sig2.s, cov.s.i)}) R <- R + diag(sig2.e, n, n) diggle.linear.power(d=1.5, t=t, R=R, sig.level=0.05, power=0.80) edland.linear.power(delta=1.5, t=t, sig2.s = 24, sig2.e = 10, sig.level=0.05, power = 0.80) #Survival #Deaths required for Cox Propo. Hazards regression w/ 2 covariates # generate a toy pilot data set X1 <- c(rep(1, 39), rep(0, 61)) set.seed(123456) X2 <- sample(c(0, 1), 100, replace = TRUE) res <- numDEpi(X1 = X1, X2 = X2, power = 0.8, theta = 2, alpha = 0.05) print(res) psi <- 0.505 # proportion of subjects died of the disease of interest. ceiling(res$D / psi) # total number of subjects required to achieve the desired power #Power Calc. in Analysis of Survival Data for Clinical Trials data(Oph) resA <- powerCT(formula=Surv(times, status)~group, dat=Oph, nE=200, nC=200, RR=0.70, alpha=0.05) resA resB <- powerCT(formula=Surv(times, status)~group, dat=Oph, nE=250, nC=250, RR=0.70, alpha=0.05) resB resC <- powerCT(formula=Surv(times, status)~group, dat=Oph, nE=300, nC=300, RR=0.70, alpha=0.05) resC #Sample size calc. in Analysis of Survival Data for Clinical Trials resD <- ssizeCT(formula = Surv(times, status) ~ group, dat = Oph, power = 0.8, k = 1, RR = 0.7, alpha = 0.05) resD #------------------------------------------------------------------------------------------ #STRUCTURAL EQUATION MODELS AND LATENT MIXED MODELS #SEM ap <- semPower.aPriori(effect = .10, effect.measure = 'RMSEA', alpha = .05, power = .80, df = 100) summary(ap) ap2 <- semPower.aPriori(effect = .10, effect.measure = "F0", alpha = .05, power = .80, df = 100) summary(ap2) ap3 <- semPower.aPriori(effect = 1.0, effect.measure = "F0", alpha = .05, power = .80, df = 100) summary(ap3) semPower.powerPlot.byN(effect = .10, effect.measure = 'RMSEA', alpha = .05, df = 100, power.min = .05, power.max = .99) semPower.powerPlot.byN(effect = 1.0, effect.measure = 'F0', alpha = .05, df = 100, power.min = .05, power.max = .99) #Latent #N/A