diff --git a/R/powerAnalysis.R b/R/powerAnalysis.R index 639f3e0..724205c 100644 --- a/R/powerAnalysis.R +++ b/R/powerAnalysis.R @@ -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