From 59fd70c34d6f7b88f7c546fa34890cea1b4b6727 Mon Sep 17 00:00:00 2001
From: Kaylea Champion <khascall@gmail.com>
Date: Tue, 14 Nov 2023 08:39:43 -0800
Subject: [PATCH] adds this to view

---
 R/powerAnalysis.orig.R | 54 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 54 insertions(+)
 create mode 100644 R/powerAnalysis.orig.R

diff --git a/R/powerAnalysis.orig.R b/R/powerAnalysis.orig.R
new file mode 100644
index 0000000..222032b
--- /dev/null
+++ b/R/powerAnalysis.orig.R
@@ -0,0 +1,54 @@
+# This is semi-generic code for doing a power analysis of a logistic regression with 4
+# levels in a factor
+# when there's some pilot values already available and defined
+#modelled heavily the simulation example explained in:
+#http://meeting.spsp.org/2016/sites/default/files/Lane%2C%20Hennes%2C%20West%20SPSP%20Power%20Workshop%202016.pdf
+
+library('batman')
+library('reshape')
+
+l2p <- function(b) {
+ odds <- exp(b) 
+ prob <- odds/(1+odds)
+ return(prob)
+}
+
+
+makeData <- function(n) { #make a random dataset of size n
+ #4 group IDs 
+ tDF <- data.frame(
+     Group0=rbinom(n=n, size=1, prob=l2p(pilot.b0)),        #ASK: what about se in pilot data?
+     Group1=rbinom(n=n, size=1, prob=l2p(pilot.b0 + pilot.b1)), # shouldn't my probs 
+     Group2=rbinom(n=n, size=1, prob=l2p(pilot.b0 + pilot.b2)), # include se?
+     Group3=rbinom(n=n, size=1, prob=l2p(pilot.b0 + pilot.b3)))
+ sDF <- melt(tDF, id.vars = 0) #AKA the index is the unique id, as far as that goes
+ colnames(sDF) <- c('source', 'nd')
+ 
+ return(sDF) 
+}
+
+powerCheck <- function(n, nSims) { #run a power calculation on the dataset given
+  #set up some empty arrays b/c R
+  signif0 <- rep(NA, nSims)
+  signif1 <- rep(NA, nSims)
+  signif2 <- rep(NA, nSims)
+  signif3 <- rep(NA, nSims)
+  signifM <- rep(NA, nSims)
+  for (s in 1:nSims) {           # repeatedly we will....
+    simData <- makeData(n)       # make some data
+    m1.sim <- glm(nd ~ source,   # give the anticipated regression a try
+                  family=binomial(link="logit"), data=simData)
+    p0 <- coef(summary(m1.sim))[1,4]
+    p1 <- coef(summary(m1.sim))[2,4]
+    p2 <- coef(summary(m1.sim))[3,4]
+    p3 <- coef(summary(m1.sim))[4,4]
+    signif0[s] <- p0 <=.05
+    signif1[s] <- p1 <=.05
+    signif2[s] <- p2 <=.05
+    signif3[s] <- p3 <=.05
+    signifM[s] <- p0 <=.05 & p1 <=.05 & p2 <=.05 & p3 <=.05
+  }  
+  power <- c(mean(signif0), mean(signif1), mean(signif2), mean(signif3), mean(signifM))
+  return(power)
+}
+