Some adjustments to make power analysis work

2 files changed, 15 insertions, 6 deletions

diff --git a/R/power.R b/R/power.R
index ad1558a..2dc5675 100644
--- a/R/power.R
+++ b/R/power.R
@@ -13,6 +13,7 @@ powerAnalysis <- function(object, coefs=NULL, design=NULL, alpha=0.05, nulls=lis
 
   submodels <- names(object@estimates@estimates)
   coefs <- check_coefs(coefs, object)
+  coefs <- generate_random_effects(coefs, object)
   fit_temp <- replace_estimates(object, coefs)
 
   T <- 1
@@ -151,11 +152,13 @@ bootstrap_data <- function(data, nsims, design){
   lapply(1:nsims, function(i) data[M_samps[[i]], J_samps[[i]]])
 }
 
-check_coefs <- function(coefs, fit, formulas, template=FALSE){
+check_coefs <- function(coefs, fit, template=FALSE){
   required_subs <- names(fit@estimates@estimates)
   required_coefs <- lapply(fit@estimates@estimates, function(x) names(x@estimates))
   required_lens <- lapply(required_coefs, length)
 
+  formulas <- sapply(names(fit), function(x) get_formula(fit, x))
+
   # If there are random effects, adjust the expected coefficient names
   # to remove the b vector and add the grouping covariate name
   rand <- lapply(formulas, lme4::findbars)
@@ -274,13 +277,15 @@ setMethod("summary", "unmarkedPower", function(object, ...){
     x
   })
 
+  coefs_no_rand <- unlist(object@coefs)[!grepl("b_", names(unlist(object@coefs)))]
+
   pow <- sapply(1:npar, function(ind){
     submod <- sum_dfs[[1]]$submodel[ind]
     param <- sum_dfs[[1]]$param[ind]
     ni <- nulls[[submod]][param]
 
     pcrit <- sapply(sum_dfs, function(x) wald(x$Estimate[ind], x$SE[ind], ni)) < object@alpha
-    direct <- sapply(sum_dfs, function(x) diff_dir(x$Estimate[ind], unlist(object@coefs)[ind], ni))
+    direct <- sapply(sum_dfs, function(x) diff_dir(x$Estimate[ind], coefs_no_rand[ind], ni))
     mean(pcrit & direct, na.rm=T)
   })
 
@@ -292,11 +297,14 @@ setMethod("summary", "unmarkedPower", function(object, ...){
     ni
   })
 
-  out <- cbind(sum_dfs[[1]][,1:2], effect=unlist(object@coefs), null=all_nulls,  power=pow)
+  effect_no_random <- unlist(object@coefs)[!grepl("b_",names(unlist(object@coefs)))]
+
+  out <- cbind(sum_dfs[[1]][,1:2], effect=effect_no_random, null=all_nulls,  power=pow)
   rownames(out) <- NULL
   names(out) <- c("Submodel", "Parameter", "Effect", "Null", "Power")
   out
 })
+
 setMethod("show", "unmarkedPower", function(object){
   cat("\nModel:\n")
   print(object@call)
diff --git a/R/simulate.R b/R/simulate.R
index 7793820..db10552 100644
--- a/R/simulate.R
+++ b/R/simulate.R
@@ -63,8 +63,8 @@ setMethod("simulate", "character",
   function(object, nsim=1, seed=NULL, formulas, coefs=NULL, design, guide=NULL, ...){
   model <- blank_umFit(object)
   fit <- suppressWarnings(simulate_fit(model, formulas, guide, design, ...))
-  coefs <- check_coefs(coefs, fit, formulas)
-  coefs <- generate_random_effects(coefs, fit, formulas)
+  coefs <- check_coefs(coefs, fit)
+  coefs <- generate_random_effects(coefs, fit)
   fit <- replace_estimates(fit, coefs)
   ysims <- suppressWarnings(simulate(fit, nsim))
   umf <- fit@data
@@ -85,8 +85,9 @@ setMethod("simulate", "character",
 })
 
 
-generate_random_effects <- function(coefs, fit, formulas){
+generate_random_effects <- function(coefs, fit){
   required_subs <- names(fit@estimates@estimates)
+  formulas <- sapply(names(fit), function(x) get_formula(fit, x))
   rand <- lapply(formulas, lme4::findbars)
   if(!all(sapply(rand, is.null))){
     rvar <- lapply(rand, function(x) unlist(lapply(x, all.vars)))