Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
# Time-stamp: <2010-12-10 14:28:19 chl>
#
# A toy example to illustrate SS idea in the Linear Model.
#
do.it <- function(x, g) {
stopifnot(is.factor(g))
stripchart(x ~ g, method="jitter", pch=19, col="darkgrey",
vertical=TRUE, ylab="")
points(1:nlevels(g), tapply(x, g, mean), pch="x", cex=1.2, col=4)
lines(1:nlevels(g), tapply(x, g, mean), col=4)
arrows(1:nlevels(g), tapply(x, g, mean)-tapply(x, g, sd)/sqrt(n),
1:nlevels(g), tapply(x, g, mean)+tapply(x, g, sd)/sqrt(n),
angle=90, length=.1, code=3, col=4)
}
pdf("1.pdf", width=8, height=5)
op <- par(mfrow=c(1,3), las=1)
set.seed(58)
# Case 1: Two-sample t-test
n <- c(20, 20)
x <- c(rnorm(n[1], mean=7.5, sd=2), rnorm(n[2], mean=10, sd=2))
g <- factor(rep(letters[1:2], each=n))
do.it(x, g)
# total variation
arrows(1.5, min(x), 1.5, max(x), code=3, angle=30, lwd=2,
length=.2, col="darkgrey")
# group-specific variation
for (i in seq_along(levels(g)))
arrows(abs(i+1.5)/2, min(x[as.numeric(g)==i]),
abs(i+1.5)/2, max(x[as.numeric(g)==i]),
code=3, angle=30, lwd=1.5, length=.1)
# Case 2: One-way ANOVA
n <- 10
k <- 3
es <- c(70,-80,0) # deviation /mean
g <- gl(3, n, , n*k, labels=letters[1:3])
x <- NULL
for (i in seq_along(levels(g)))
x <- append(x, rnorm(n, mean=200+es[i], sd=75))
do.it(x, g)
# total variation
arrows(2.5, min(x), 2.5, max(x), code=3, angle=30, lwd=2,
length=.2, col="darkgrey")
# variation between means (BSS)
arrows(1.5, min(tapply(x, g, mean)), 1.5, max(tapply(x, g, mean)),
code=3, angle=30, lwd=2, length=.2, col="darkgrey")
# group-specific total variation (WSS)
for (i in seq_along(levels(g)))
arrows(i+.15, min(x[as.numeric(g)==i]),
i+.15, max(x[as.numeric(g)==i]),
code=3, angle=30, lwd=1.5, length=.1)
# Case 3: Simple linear regression
n <- 30
# x <- runif(n, 5, 20)
x <- sample(seq(5, 20), n, replace=TRUE)
y <- 2*x/10 + rnorm(n)
mod.lm <- lm(y ~ x)
new <- data.frame(x = seq(5, 20, by=.5))
cb <- predict(lm(y ~ x), new, interval="confidence")
plot(y ~ x, xlab="", ylab="", cex=.7, pch=19, col="darkgrey")
abline(mod.lm, col=4)
lines(new$x, cb[,"lwr"], lty=3, col=4)
lines(new$x, cb[,"upr"], lty=3, col=4)
# total variation
arrows((5+20)/2, min(y), (5+20)/2, max(y), code=3, angle=30, lwd=2,
length=.2, col="darkgrey")
# x-specific variation
x0 <- as.numeric(names(table(x)[table(x)>1][1]))
arrows(x0, min(y[x==x0]), x0, max(y[x==x0]),
code=3, angle=30, lwd=1.5, length=.1)
par(op)
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment