mgaughan/24_deb_pkg_gov
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

milestones now sim at the right freq

Browse Source 
corrects the coef(summary(model)) notation
This commit is contained in:
Kaylea Champion 2023-11-14 14:29:36 -08:00
parent 0b1467c823
commit 03fd8dd891
1 changed files with 17 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
R/powerAnalysis.R
View File

@ -17,24 +17,27 @@ l2p <- function(b) {
#Matt:
makeDataNew <- function(n) {
tDF <- data.frame(
Group0=rnorm(n=n, mean=-0.1296376, sd=1.479847), # up.fac.mean
Group1=rlnorm(n=n, mean=1.685715, sd = 0.2532059), # mmt
Group2=rbinom(n=n, size=1, prob=c(0.247, 0.753)), #milestones
Group3=rnorm(n=n, mean=4351.578, sd=1408.811) # age
up.fac.mean=rnorm(n=n, mean=-0.1296376, sd=1.479847), # up.fac.mean
mmt=rlnorm(n=n, mean=1.685715, sd = 0.2532059), # mmt
## this generates a 50-50 split of milestones --v
#milestones=rbinom(n=n, size=1, prob=c(0.247, 0.753)), #milestones
milestones=rbinom(n=n, size=1, prob=.247), #milestones
age=rnorm(n=n, mean=4351.578, sd=1408.811) # age
)
#sDF <- melt(tDF, id.vars = 0) #AKA the index is the unique id, as far as that goes
colnames(tDF) <- c('up.fac.mean', 'mmt', 'milestones', 'age')
## can name these in the data.frame constructor method directly
#colnames(tDF) <- c('up.fac.mean', 'mmt', 'milestones', 'age')
return(tDF)
}
makeDataNew2 <- function(n) {
tDF <- data.frame(
Group0=rnorm(n=n, mean=-0.1296376, sd=1.479847), # up.fac.mean
Group1=rlnorm(n=n, mean=6.220282, sd = 2.544058) # formal.score
up.fac.mean=rnorm(n=n, mean=-0.1296376, sd=1.479847), # up.fac.mean
formal.score=rlnorm(n=n, mean=6.220282, sd = 2.544058) # formal.score
)
tDF[is.na(tDF) | tDF=="Inf"] = NA
#sDF <- melt(tDF, id.vars = 0) #AKA the index is the unique id, as far as that goes
colnames(tDF) <- c('up.fac.mean', 'formal.score')
##colnames(tDF) <- c('up.fac.mean', 'formal.score')
return(tDF)
}
@ -51,10 +54,10 @@ powerCheck <- function(n, nSims) { #run a power calculation on the dataset given
#have updated for kkex through here, now need to look at the underproduction work
#m1.sim <- lm(up.fac.mean ~ ((mmt)/ (milestones/age)), data=simData)
m1.sim <- lm(up.fac.mean ~ mmt + milestones + age, data=simData)
p0 <- coef(summary(m1.sim))[1,4]
p1 <- coef(summary(m1.sim))[1,4]
p2 <- coef(summary(m1.sim))[1,4]
p3 <- coef(summary(m1.sim))[1,4]
p0 <- coef(summary(m1.sim))[1,4] #intercept
p1 <- coef(summary(m1.sim))[2,4] #mmt
p2 <- coef(summary(m1.sim))[3,4] #milestones
p3 <- coef(summary(m1.sim))[4,4] #age
signif0[s] <- p0 <=.05
signif1[s] <- p1 <=.05
signif2[s] <- p2 <=.05
@ -69,16 +72,14 @@ powerCheck2 <- function(n, nSims) { #run a power calculation on the dataset give
#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 <- makeDataNew2(n) # make some data
#have updated for kkex through here, now need to look at the underproduction work
#m1.sim <- lm(up.fac.mean ~ ((mmt)/ (milestones/age)), data=simData)
m1.sim <- lm(up.fac.mean ~ formal.score, data=simData)
p0 <- coef(summary(m1.sim))[1,2]
p1 <- coef(summary(m1.sim))[1,2]
p0 <- coef(summary(m1.sim))[1,4]
p1 <- coef(summary(m1.sim))[2,4]
signif0[s] <- p0 <=.05
signif1[s] <- p1 <=.05
signifM[s] <- p0 <=.05 & p1 <=.05