## ----include = FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- knitr::opts_chunk$set( collapse = FALSE, comment = '#>' ) old = options(width = 200) ## ----warning=F, message=F----------------------------------------------------------------------------------------------------------------------------------------------------------------------------- library(dplyr) library(tidyr) library(glmmTMB) library(emmeans) library(powerutilities) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- power.t.test(n = 8, delta = 3, sd = 2, sig.level = 0.05, power = NULL, type = 'two.sample', alternative = 'two.sided') ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- nrep = 8 # Number of replicates nfac = 2 # Number of treatments # Expand grid creates a data set # with all combinations of treatments ex1 = expand.grid(Rep = factor(1:nrep), Trt = factor(LETTERS[1:nfac])) # Set the response values to correspond # to the expected treatment differences ex1$Y = ifelse(ex1$Trt == 'A', 0, 3) head(ex1, n = 10) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex1_mod = set_glmm(Y ~ Trt, data = ex1, disp = 2) # Confirm it worked: sigma(ex1_mod) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- power_ftest(ex1_mod) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- power_ftest(ex1_mod, ddf = 'asymptotic') ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- (ex1_emm = emmeans(ex1_mod, ~Trt)) # Estimated treatment means ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- contr_list = list('Trt A - Trt B' = c(1, -1)) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- power_contrast(ex1_emm, contr_list) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- (ex1_emm = emmeans(ex1_mod, ~Trt)) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- power_contrast(ex1_emm, contr_list) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex1_emm = emmeans(ex1_mod, ~Trt) power_contrast(ex1_emm, contr_list, ddf = 15, infer = c(T, T)) |> print(digits = 3) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex1_mod2 = set_glmm(Y ~ Trt, data = ex1, disp = c(1.75, 2.25), dispformula = ~ Trt - 1) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- (ex1_emm2 = emmeans(ex1_mod2, ~Trt)) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- power_contrast(ex1_emm2, contr_list) |> print(digits = 3) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex2_trt_dict = expand.grid(gen = c('GoldenRain', 'Marvellous', 'Victory'), nitro = factor(c(0, 20, 40, 60))) |> arrange(gen) ex2_trt_dict ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex2_trt_dict$yield = c(1400 + 150*0, 1400 + 150*1, 1400 + 150*2, 1400 + 150*3, 1500 + 150*0, 1500 + 150*1, 1500 + 150*2, 1500 + 150*3, 1400 + 200*0, 1400 + 200*1, 1400 + 200*2, 1400 + 200*3) ex2_trt_dict ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- nblk = 6 ex2_des_df = data.frame(block = factor(1:nblk)) ex2_des_df ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex2 = merge(ex2_trt_dict, ex2_des_df) head(ex2, 10) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- theta_finder(yield ~ gen*nitro + (1|block/gen), data = ex2) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex2_mod = set_glmm(yield ~ gen*nitro + (1|block/gen), data = ex2, re_terms = c(200, 250), disp = 225) VarCorr(ex2_mod, sigma = extract_disp(ex2_mod)) # It worked ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- power_ftest(ex2_mod) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- (ex2_emm1 = emmeans(ex2_mod, ~ gen|nitro, ddf = 'kenward-roger')) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex2_contr1 = list('Golden - Victory' = c(1, 0, -1), 'Marvellous - Victory' = c(0, 1, -1)) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- power_contrast(ex2_emm1, ex2_contr1) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex2_mod2 = set_glmm(yield ~ gen*nitro + (1|block), data = ex2, re_terms = 250, disp = 225) power_ftest(ex2_mod2) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex3_trt_dict = expand.grid(gen = c('GoldenRain', 'Marvellous', 'Victory'), nitro = c(0, 20, 40, 60)) |> arrange(gen) ex3_trt_dict$yield = c(1400 + 150*0, 1400 + 150*1, 1400 + 150*2, 1400 + 150*3, 1500 + 150*0, 1500 + 150*1, 1500 + 150*2, 1500 + 150*3, 1400 + 200*0, 1400 + 200*1, 1400 + 200*2, 1400 + 200*3) nblk = 6 ex3_des_df = data.frame(block = factor(1:nblk)) ex3 = merge(ex3_trt_dict, ex3_des_df) head(ex3, 10) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex3_mod = set_glmm(yield ~ nitro*gen + (1|block/gen), data = ex3, re_terms = c(250, 200), disp = 225) VarCorr(ex3_mod, sigma = extract_disp(ex3_mod)) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- power_ftest(ex3_mod, cov.reduce = F) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- (ex3_emt = emtrends(ex3_mod, ~ gen, var = 'nitro', ddf = 'kenward-roger')) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex3_contr1 = list('Slope Victory - Golden' = c(-1, 0, 1)) power_contrast(ex3_emt, ex3_contr1, ddf = 'kenward-roger') ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex3_mod2 = set_glmm(yield ~ factor(nitro)*gen + (1|block/gen), data = ex3, re_terms = c(250, 200), disp = 225) (ex3_emm2 = emmeans(ex3_mod2, ~ nitro:gen)) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- poly(c(0, 0.2, 0.4, 0.6), degree = 1)[, 1] ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ex3_contr2 = list('Slope Vic ≠ 0' = c(-0.671, -0.224, 0.224, 0.671, 0, 0, 0, 0, 0, 0, 0, 0), 'Slope Vic - Slope Golden' = c(+0.671, +0.224, -0.224, -0.671, 0, 0, 0, 0, -0.671, -0.224, +0.224, +0.671)) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- power_contrast(ex3_emm2, ex3_contr2, ddf = 'kenward-roger') ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- (ex3_emm3 = emmeans(ex3_mod2, ~ nitro|gen)) ## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # Naturally, since we created data with zero curvature, there is no power # for detecting quadratic or cubic trends power_contrast(ex3_emm3, poly.emmc) ## ----include = FALSE---------------------------------------------------------- ## Notes for future development of the vignette # Example 4: A Repeated Measures Study with AR(1) Correlation Structure # Forthcoming # Example 5: Data with a Gaussian Spatial Correlation Structure # Forthcoming # Example 6: A Random Slopes Model with Correlated Slope and Intercept # Forthcoming options(old)