in class

Makefile

@@ -129,10 +129,16 @@ Sources += $(wildcard docs/code/*.*)
 
 ######################################################################
 
+sims.Rout: code/sims.R
+	$(pipeR)
+
+simlist.Rout: simlist.R sims.Rout
+	$(pipeR)
+
 testsims.Rout: code/testsims.R
 	$(pipeR)
 
-sims.Rout: code/sims.R
+simCIs.Rout: simCIs.R testsims.rds
 	$(pipeR)
 
 australia.Rout: code/australia.R

docs/code/sims.R

@@ -1,42 +1,89 @@
 library(glmmTMB)
 library(dplyr)
+library(purrr)
+library(broom.mixed)
 
-set.seed(101)
+## These ...'s aren't the best practice; should figure out more about map
+sim <- function(groupSize, days, β0, β_treat, β_day
+	, sdint, sdslope, corr, sdres=2, ...
+){
 
-nGroup <- 15
-nBats <- nGroup*2
-nGroup=3; days=seq(3, 60, by=3)
-β0=10; β_treat=10; β_day=2
-sdint=3; sdslope=0.1; corr=0; sdres=2
+	nBats <- groupSize*2 ## Treatment and control
 
-batID <- as.factor(1:nBats) ## Make this look less like a number with paste()
-treatment <- as.factor(rep(c("control", "exercise"), each=nGroup))
+	batID <- as.factor(1:nBats) ## Make this look less like a number with paste()
+	treatment <- as.factor(rep(c("control", "exercise"), each=groupSize))
 
-tdat <- data.frame(batID, treatment)
+	tdat <- data.frame(batID, treatment)
 
-dat <- expand.grid(batID=batID, day=days)
+	dat <- expand.grid(batID=batID, day=days)
 
-dat <- full_join(dat, tdat, by="batID")
+	dat <- full_join(dat, tdat, by="batID")
 
-y <- simulate_new( ~ 1 + treatment + day + (1 + day | batID)
-	, nsim = 1
-	, family = "gaussian"
-	, newdata = dat
-	, newparams = list(
-		beta = c(β0, β_treat, β_day)
-		, theta = c(log(sdint), log(sdslope), corr)
-		, betad = log(sdres)
+	y <- simulate_new( ~ 1 + treatment + day + (1 + day | batID)
+		, nsim = 1
+		, family = "gaussian"
+		, newdata = dat
+		, newparams = list(
+			beta = c(β0, β_treat, β_day)
+			, theta = c(log(sdint), log(sdslope), corr)
+			, betad = log(sdres)
+		)
 	)
-)
 
-dat$flightTime <- y[[1]]
+	dat$flightTime <- y[[1]]
+	return(dat)
+}
+
+fit <- function(dat){
+	return(glmmTMB(flightTime ~ 1 + treatment + day + (1 + day | batID)
+		, data=dat
+		, family = "gaussian"
+	))
+}
+
+simCIs <- function(simfun, fitfun, ...){
+	dat <- simfun(...)
+	fit <- fitfun(dat)
+	c <- as.data.frame(confint(fit))
+	return(data.frame(vname=rownames(c)
+		, est = c$Est
+		, low = c$"2.5"
+		, high = c$"97.5"
+	))
+}
+
+set.seed(101)
+nReps <- 1000
 
-fit <- glmmTMB(flightTime ~ 1 + treatment + day + (1 + day | batID)
-	, data=dat
-	, family = "gaussian"
+print(simCIs(simfun=sim, fitfun=fit
+	, groupSize = 15, days=seq(3, 60, by=3)
+	, β0=10, β_treat=10, β_day=2
+	, sdint=3, sdslope=0.1, corr=0, sdres=2
+))
+
+cis <- map_dfr(1:nReps, simCIs, .id="sim", simfun=sim, fitfun=fit
+	, groupSize = 15, days=seq(3, 60, by=3)
+	, β0=10, β_treat=10, β_day=2
+	, sdint=3, sdslope=0.1, corr=0, sdres=2
 )
 
-summary(fit)
+summary(cis)
+
+## saveRDS(cis, "sims.rds")
+## cis <- readRDS("sims.rds")
+
+## summary(cis |> mutate_if(is.character, as.factor))
 
-confint(fit)
+## Check above
+treatmentTrue <- 10
 
+## We are only interested in “power” to detect the true effect direction
+print(cis
+	|> filter(vname=="treatmentexercise")
+	|> summarize(
+		toohigh=mean(low>treatmentTrue)
+		, toolow=mean(high<treatmentTrue)
+		, ci_width=mean(high-low)
+		, power = mean(low>0)
+	)
+)

docs/code/testsims.R

@@ -1,10 +1,11 @@
 library(glmmTMB)
 library(dplyr)
+library(purrr)
 
-set.seed(101)
+library(shellpipes)
 
 sim <- function(nGroup, days, β0, β_treat, β_day
-	, sdint, sdslope, corr, sdres
+	, sdint, sdslope, corr, sdres, ...
 ){
 	nBats <- nGroup*2
 	batID <- as.factor(1:nBats)
@@ -40,16 +41,40 @@ fit <- function(dat){
 		, data=dat
 		, family = "gaussian"
 	)
-	summary(fit)
-	confint(fit)
+	## summary(fit)
+	return(confint(fit))
+}
+
+simCIs <- function(simfun, fitfun, ...)
+{
+	dat <- simfun(...)
+	fit <- as.data.frame(fitfun(dat))
+	return(data.frame(vname = row.names(fit)
+		, est = fit$Est
+		, low = fit$"2.5"
+	))
 }
 
 ######################################################################
 
+numSims <- 10
+
 set.seed(2134)
 s0 <- sim(nGroup=3, days=seq(3, 60, by=3)
 	, β0=10, β_treat=10, β_day=2
 	, sdint=3, sdslope=0.1, corr=0, sdres=2
 )
 
-fit(s0)
+simCIs(sim, fit, nGroup=3, days=seq(3, 60, by=3)
+	, β0=10, β_treat=10, β_day=2
+	, sdint=3, sdslope=0.1, corr=0, sdres=2
+)
+
+ciList <- map(.x=1:numSims, .f=simCIs
+	, simfun=sim, fitfun=fit
+	, nGroup=3, days=seq(3, 60, by=3)
+	, β0=10, β_treat=10, β_day=2
+	, sdint=3, sdslope=0.1, corr=0, sdres=2
+)
+
+rdsSave(ciList)