# Use SRS samples and make SRS folds
srssrsds <- data.frame(df = c(), MSE = c())
# Use SRS samples and evaluate on full pop
srspopds <- data.frame(df = c(), MSE = c())
# Use Cluster samples and make Cluster folds
clusclusds <- data.frame(df = c(), MSE = c())
# Use Cluster samples and make SRS folds
clussrsds <- data.frame(df = c(), MSE = c())
# Use Cluster samples and evaluate on full pop
cluspopds <- data.frame(df = c(), MSE = c())
# Use Strat samples and make Strat folds
stratstratds <- data.frame(df = c(), MSE = c())
# Use Strat samples and make SRS folds
stratsrsds <- data.frame(df = c(), MSE = c())
# Use Strat samples and evaluate on full pop
stratpopds <- data.frame(df = c(), MSE = c())
modelsToFit <- c("Response~splines::ns(Predictor, df=1)",
"Response~splines::ns(Predictor, df=2)",
"Response~splines::ns(Predictor, df=3)",
"Response~splines::ns(Predictor, df=4)",
"Response~splines::ns(Predictor, df=5)",
"Response~splines::ns(Predictor, df=6)")
# Making as many simple random samples as we specify for 'loops'
for (i in 1:loops) {
# Take a SRS
sim.srs <- dplyr::sample_n(splinepop.df, n)
# Using our SRS function on SRS samples
srs.CV.out <- cv.svy(sim.srs, modelsToFit,
nfolds = 5, fpcID = "fpcSRS") %>% as.vector()
srssrsds2 <- data.frame(df = 1:6, MSE = srs.CV.out)
srssrsds <- rbind(srssrsds, srssrsds2)
# Eval SRS samples on full pop
sub.srs <- dplyr::sample_n(sim.srs, n*.8)
srs.des <- svydesign(ids = ~1, fpc = ~fpcSRS, data = sub.srs)
srs.pop.out <- sapply(1:6,
function(ii) {
yhat <- predict(svyglm(modelsToFit[ii], srs.des),
newdata = splinepop.df)
return(mean((yhat - splinepop.df$Response)^2))
})
srspopds2 <- data.frame(df = 1:6, MSE = srs.pop.out)
srspopds <- rbind(srspopds, srspopds2)
}
# Making as many cluster samples as we specify for 'loops'
for (i in 1:loops) {
# Take a Cluster sample
c <- unique(splinepop.df[["Cluster"]])
sim.clus <- splinepop.df[splinepop.df[["Cluster"]] %in% sample(c, n/clusSizes),]
# Using our Cluster function on Cluster samples
clus.CV.out <- cv.svy(sim.clus, modelsToFit,
clusterID = "Cluster", nfolds = 5, fpcID = "fpcCluster") %>% as.vector()
clusclusds2 <- data.frame(df = 1:6, MSE = clus.CV.out)
clusclusds <- rbind(clusclusds, clusclusds2)
# Using our SRS function on Cluster samples
srs.CV.out <- cv.svy(sim.clus, modelsToFit,
nfolds = 5, fpcID = "fpcSRS") %>% as.vector()
clussrsds2 <- data.frame(df = 1:6, MSE = srs.CV.out)
clussrsds <- rbind(clussrsds, clussrsds2)
# Eval 80% Cluster samples on full pop
sub.c <- unique(sim.clus[["Cluster"]])
sub.clus <- sim.clus[sim.clus[["Cluster"]] %in% sample(sub.c, (n/clusSizes)*.8),]
clus.des <- svydesign(ids = ~Cluster, fpc = ~fpcCluster, data = sub.clus)
clus.pop.out <- sapply(1:6,
function(ii) {
yhat <- predict(svyglm(modelsToFit[ii], clus.des),
newdata = splinepop.df)
return(mean((yhat - splinepop.df$Response)^2))
})
cluspopds2 <- data.frame(df = 1:6, MSE = clus.pop.out)
cluspopds <- rbind(cluspopds, cluspopds2)
}
# Making as many stratified samples as we specify for 'loops'
for (i in 1:loops) {
# Take a Strat sample
stratcounts <- rep(n/nrStrata, nrStrata)
names(stratcounts) <- 1:nrStrata
s <- stratsample(splinepop.df$Stratum, stratcounts)
sim.strat <- splinepop.df[s,]
# Using our Strat function on Strat samples
strat.CV.out <- cv.svy(sim.strat, modelsToFit,
strataID = "Stratum", nfolds = 5, fpcID = "fpcStratum") %>% as.vector()
stratstratds2 <- data.frame(df = 1:6, MSE = strat.CV.out)
stratstratds <- rbind(stratstratds, stratstratds2)
# Using our SRS function on Strat samples
srs.CV.out <- cv.svy(sim.strat, modelsToFit,
nfolds = 5, fpcID = "fpcSRS") %>% as.vector()
stratsrsds2 <- data.frame(df = 1:6, MSE = srs.CV.out)
stratsrsds <- rbind(stratsrsds, stratsrsds2)
# Eval 80% Strat samples on full pop
sub.stratcounts <- rep((n/nrStrata)*.8, nrStrata)
names(sub.stratcounts) <- 1:nrStrata
sub.s <- survey::stratsample(sim.strat$Stratum, sub.stratcounts)
sub.strat <- sim.strat[sub.s,]
strat.des <- svydesign(ids = ~1, strata = ~Stratum, fpc = ~fpcStratum, data = sub.strat)
strat.pop.out <- sapply(1:6,
function(ii) {
yhat <- predict(svyglm(modelsToFit[ii], strat.des),
newdata = splinepop.df)
return(mean((yhat - splinepop.df$Response)^2))
})
stratpopds2 <- data.frame(df = 1:6, MSE = strat.pop.out)
stratpopds <- rbind(stratpopds, stratpopds2)
}
# Making the degrees of freedom variable a factor variable
srssrsds$df <- as.factor(srssrsds$df)
srspopds$df <- as.factor(srspopds$df)
clusclusds$df <- as.factor(clusclusds$df)
clussrsds$df <- as.factor(clussrsds$df)
cluspopds$df <- as.factor(cluspopds$df)
stratstratds$df <- as.factor(stratstratds$df)
stratsrsds$df <- as.factor(stratsrsds$df)
stratpopds$df <- as.factor(stratpopds$df)
usethis::use_data(srssrsds, srspopds,
clusclusds, clussrsds, cluspopds,
stratstratds, stratsrsds, stratpopds,
internal = FALSE, overwrite = TRUE)
# TODO: for now setting internal=FALSE
# so that we don't overwrite the sysdata.Rda for intro.Rmd;
# but eventually when we're ready to replace that old vignette,
# we can return to internal=TRUE here
# if that's indeed better for R packages
# ...
# OK, now we've written those dataframes,
# AND the 4 below, all into R/sysdata.Rda together
if(!rerunSims) {
srssrsds <- surveyCV:::srssrsds
clusclusds <- surveyCV:::clusclusds
clussrsds <- surveyCV:::clussrsds
stratstratds <- surveyCV:::stratstratds
stratsrsds <- surveyCV:::stratsrsds
srspopds <- surveyCV:::srspopds
cluspopds <- surveyCV:::cluspopds
stratpopds <- surveyCV:::stratpopds
}
# Find the y-ranges of MSEs
yminMSE <- min(min(srssrsds$MSE),
min(clusclusds$MSE), min(clussrsds$MSE),
min(stratstratds$MSE), min(stratsrsds$MSE))
ymaxMSE <- max(max(srssrsds$MSE),
max(clusclusds$MSE), max(clussrsds$MSE),
max(stratstratds$MSE), max(stratsrsds$MSE))
srspopds <- dplyr::group_by(srspopds, df) %>%
dplyr::summarise(MSE = mean(MSE)) %>%
dplyr::mutate(df = as.numeric(df))
cluspopds <- dplyr::group_by(cluspopds, df) %>%
dplyr::summarise(MSE = mean(MSE)) %>%
dplyr::mutate(df = as.numeric(df))
stratpopds <- dplyr::group_by(stratpopds, df) %>%
dplyr::summarise(MSE = mean(MSE)) %>%
dplyr::mutate(df = as.numeric(df))
# Making ggplot objects for the MSEs and AICs
plot.srssrs <- ggplot(data = srssrsds, mapping = aes(x = df, y = MSE/1e4)) +
geom_boxplot() +
geom_line(data = srspopds, mapping = aes(x = df, y = MSE/1e4)) +
ggtitle("CV: SRS folds,\nSRS sample") +
scale_y_log10(limits = c(yminMSE, ymaxMSE)/1e4)
plot.blank <- rectGrob(width = 0, height = 0) # empty spot here -- no corresponding sim
plot.clusclus <- ggplot(data = clusclusds, mapping = aes(x = df, y = MSE/1e4)) +
geom_boxplot() +
geom_line(data = cluspopds) +
ggtitle("CV: Cluster folds,\nCluster sample") +
scale_y_log10(limits = c(yminMSE, ymaxMSE)/1e4)
plot.clussrs <- ggplot(data = clussrsds, mapping = aes(x = df, y = MSE/1e4)) +
geom_boxplot() +
geom_line(data = cluspopds) +
ggtitle("CV: SRS folds,\nCluster sample") +
scale_y_log10(limits = c(yminMSE, ymaxMSE)/1e4)
plot.stratstrat <- ggplot(data = stratstratds, mapping = aes(x = df, y = MSE/1e4)) +
geom_boxplot() +
geom_line(data = stratpopds) +
ggtitle("CV: Strat folds,\nStrat sample") +
scale_y_log10(limits = c(yminMSE, ymaxMSE)/1e4)
plot.stratsrs <- ggplot(data = stratsrsds, mapping = aes(x = df, y = MSE/1e4)) +
geom_boxplot() +
geom_line(data = stratpopds) +
ggtitle("CV: SRS folds,\nStrat sample") +
scale_y_log10(limits = c(yminMSE, ymaxMSE)/1e4)
gridExtra::grid.arrange(plot.srssrs, plot.clussrs, plot.stratsrs,
plot.blank, plot.clusclus, plot.stratstrat,
ncol = 3,
top = paste0("Simulated Data (Sample Sizes = ", n,
", Clusters = ", n/clusSizes,
" or Strata = ", nrStrata,
", Loops = ", loops, ", Folds = 5)"))
