Skip to content

Instantly share code, notes, and snippets.

@AkselA
Last active August 21, 2018 09:46
Show Gist options
  • Save AkselA/f4516b9df572ad2cc66bf0c8c8d5eb29 to your computer and use it in GitHub Desktop.
Save AkselA/f4516b9df572ad2cc66bf0c8c8d5eb29 to your computer and use it in GitHub Desktop.
shootings <- structure(list(Year = 1982:2017, Fatalities = c(8,
0, 28, 0, 15, 6, 7, 15, 10, 35, 9, 23, 5, 6, 6, 9, 14, 42, 7,
5, 0, 7, 5, 17, 21, 53, 17, 39, 9, 19, 71, 35, 15, 37, 65, 99
)), row.names = c(NA, 36L), class = "data.frame")
library(caret)
set.seed(1)
k <- 5
yeargroups <- cut(shootings[,1], nrow(shootings)/k)
shootings$Fold <- createFolds(yeargroups, k, list=FALSE)
head(shootings)
# Year Fatalities Fold
# 1 1982 8 3
# 2 1983 0 2
# 3 1984 28 4
# 4 1985 0 1
# 5 1986 15 2
# 6 1987 6 5
par(mar=c(3, 3, 1, 1), mgp=c(1.8, 0.6, 0))
plot(shootings[1:2], col=rainbow(k)[shootings$Fold], pch=16, cex=1.5)
aggregate(Year ~ Fold, shootings, mean)
# Fold Year
# 1 1 2000.286
# 2 2 1998.125
# 3 3 1999.571
# 4 4 1999.286
# 5 5 2000.429
ll <- list()
ngroups <- 1:15
for (i in ngroups) {
lli <- replicate(200, {
if (i > 1) {
yeargroups <- cut(shootings[,1], i)
} else {
yeargroups <- shootings[,1]
}
shootings$Fold <- createFolds(yeargroups, k, list=FALSE)
a <- aggregate(Year ~ Fold, shootings, mean)[,2]
round(diff(range(a)), 4)
})
ll[[i]] <- lli
}
ll.mat <- simplify2array(ll)
par(xaxs="i")
matplot(ngroups, t(ll.mat), type="l", lty=1, lwd=0.2, col="#00000044",
xaxt="n", xlab="Number of strata", ylab="Max absolute difference")
lines(ngroups, apply(ll.mat, 2, mean), col="#2078FA", lwd=4)
axis(1, at=ngroups)
legend("topleft", legend="Mean", bty="n", lwd=4, col="#2078FA")
quartz.save("nstrat.png", dpi=110)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment