Talk about power analtsis for the Sevilla R users meeting
#SevillaR talk | |
#The problem: | |
time <- c(2000:2015) | |
abundance <- rnorm(16, 150, 50) #poison?? | |
plot(abundance ~ time, t = "l") | |
#can I detect a trend? | |
m <- lm(abundance ~ time) | |
summary(m) | |
abline(m, col = "red") | |
#but, can I detect a trend? Power analysis. | |
library(pwr) | |
pwr.r.test(n = 16, r = cor(time, abundance), sig.level = 0.05, power = NULL) | |
#how many more years would I need? | |
(sample_needed <- pwr.r.test(n = NULL, r = cor(time, abundance), sig.level = 0.05, power = 0.80)) | |
sample_needed[1] < length(time) | |
#same for anovas (pwr.anova.test) | |
# simple, right? but data is complex | |
abundance2 <- jitter(abundance, 1000) | |
abundance3 <- jitter(abundance, 1000) | |
plot(abundance ~ time, t = "p") | |
points(abundance2 ~ time, col = "red") | |
points(abundance3 ~ time, col = "blue") | |
#build a data frame | |
data <- data.frame(time = rep(time, 3), abundance = c(abundance, abundance2, abundance3), | |
site = c(rep("uno", 16), rep("dos", 16), rep("tres", 16))) | |
head(data) | |
library(lme4) | |
m <- lmer(abundance ~ time + (1|site), data) | |
summary(m) | |
abline(fixef(m)) | |
#so, what is my power? how many years? #do simulations! | |
library(devtools) | |
install_github("pitakakariki/simr") | |
library(simr) | |
#Specify the effect size | |
fixef(m) | |
fixef(m)["time"] | |
fixef(m)["time"] <- -2 #a 10% change in 10 years, i.e. 20 individuos sobre 200 en 10 años | |
fixef(m) | |
abline(fixef(m), col = "red") | |
#Calculate power | |
powerSim(m) #slow pre-calculated | |
ps <- lastResult() | |
#Use extend to increase the sample size. | |
model2 <- extend(m, along = "time", n=30) | |
powerSim(model2) | |
#Power curve over range of sample sizes. | |
pc2 <- powerCurve(model2) | |
print(pc2) | |
plot(pc2) | |
#Increase the number or size of groups. | |
#Power curve with more groups. | |
model3 <- extend(m, along="site", n=5) #this is equivalent at saying more sites, equal number of years | |
pc3 <- powerCurve(model3, along="site") | |
print(pc3) | |
plot(pc3) | |
#Power curve with larger groups. | |
model4 <- extend(m, within="time+site", n=5) #this is like taking more than one value per site (same years and sites) | |
pc4 <- powerCurve(model4, within="time+site", nsim = 99) | |
print(pc4) | |
plot(pc4) | |
#If pilot data is not available, simr can be used to create lme4 objects from scratch as a starting point. | |
#https://github.com/pitakakariki/simr/blob/master/vignettes/fromscratch.Rmd | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment