diff --git a/R/calculatePower.R b/R/calculatePower.R index e48afc6..c69389a 100644 --- a/R/calculatePower.R +++ b/R/calculatePower.R @@ -31,9 +31,11 @@ data1$milestones <- as.numeric(data1$milestones > 0) + 1 # (2) - Run the model on the pilot data data1$formal.score <- data1$mmt / (data1$milestones/data1$age) table(data1$milestones) -hist(data1$old_mmt) #inequality of participation +hist(data1$old_mmt, prob=TRUE) #inequality of participation hist(data1$formal.score) hist(data1$age/365) +data1$new_mmt <- data1$mmt - 1 +hist(data1$new_mmt, prob=TRUE) data1$new.age <- as.numeric(cut(data1$age/365, breaks=c(0,9,12,15,17), labels=c(1,2,3,4))) data1$formal.score <- data1$mmt / (data1$milestones/data1$new.age) @@ -42,7 +44,7 @@ table(data1$formal.score) hist(data1$formal.score) kmodel1 <- lm(up.fac.mean ~ mmt, data=data1) summary(kmodel1) -kmodel1 <- lm(up.fac.mean ~ old_mmt, data=data1) +kmodel1 <- lm(up.fac.mean ~ new_mmt, data=data1) summary(kmodel1) kmodel1 <- lm(up.fac.mean ~ new.age, data=data1) summary(kmodel1) @@ -54,14 +56,14 @@ cor.test(data1$mmt, data1$up.fac.mean) cor.test(data1$milestones, data1$up.fac.mean) cor.test(data1$age, data1$up.fac.mean) -g <- ggplot(data1, aes(x=mmt, y=up.fac.mean)) + +g <- ggplot(data1, aes(x=new_mmt, y=up.fac.mean)) + geom_point() + geom_smooth() g data2 <- subset(data1, (data1$age / 365) < 14 ) hist(floor(data2$age)) -g <- ggplot(data2, aes(x=formal.score, y=up.fac.mean)) + +g <- ggplot(data2, aes(x=mmt, y=up.fac.mean)) + geom_point() + geom_smooth() g @@ -109,9 +111,9 @@ powerCheck(n, nSims) #powerCheck2(n, nSims) like doesn't really work #Sample values -powerCheck(50, 1000) +powerCheck(100, 1000) powerCheck(200, 1000) -powerCheck(500, 1000) +powerCheck(300, 1000) powerCheck2(50, 1000) powerCheck2(200, 1000) diff --git a/R/powerAnalysis.R b/R/powerAnalysis.R index e4a0ccd..2c6f0ce 100644 --- a/R/powerAnalysis.R +++ b/R/powerAnalysis.R @@ -19,7 +19,8 @@ makeDataNew <- function(n) { tDF <- data.frame( ## don't sim the outcome #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 + #mmt=rlnorm(n=n, mean=1.685715, sd = 0.2532059), # mmt + new_mmt=rbeta(n=n, 5, 1), #mmt=rlogis(n=n, location = 1.685715), ## this generates a 50-50 split of milestones --v #milestones=rbinom(n=n, size=1, prob=c(0.247, 0.753)), #milestones @@ -57,21 +58,21 @@ powerCheck <- function(n, nSims) { #run a power calculation on the dataset given ## outcome goes here --v # e.g. simData$up.fac.mean <- (usefuleffsizeA * mmt) + (usefuleffsizeB * milestones) + rnorm(n=1, mean=0, sd=1) ##plus some noise #simData$up.fac.mean <- (-2.075 * simData$mmt) + (0.4284 * simData$milestones) + rnorm(n=1, mean=0, sd=1) - simData$up.fac.mean <- (0.4284 * simData$milestones) + rnorm(n=1, mean=0, sd=1) + simData$up.fac.mean <- (-2.0745 * simData$new_mmt) + (0.4284 * simData$milestones) + rnorm(n=n, mean=0, sd=1) #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) ## could leave age out for now? #m1.sim <- lm(up.fac.mean ~ mmt + milestones + age, data=simData) - m1.sim <- lm(up.fac.mean ~ mmt + milestones, data=simData) + m1.sim <- lm(up.fac.mean ~ new_mmt + milestones, data=simData) p0 <- coef(summary(m1.sim))[1,4] #intercept p1 <- coef(summary(m1.sim))[2,4] #mmt - #p2 <- coef(summary(m1.sim))[3,4] #milestones + 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 + signif2[s] <- p2 <=.05 #signif3[s] <- p3 <=.05 - signifM[s] <- p0 <=.05 & p1 <=.05 #& p2 <=.05 #& 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)