set.seed(123) n <- 25 (x <- sort(round(runif(n, 1,10), 2))) N <- 9 Y <- matrix(NA, n, N) LM <- as.list(1:N) for(i in 1:N) { Y[,i] <- 2 + x + rnorm(n, s = 2) LM[[i]] <- lm(Y[,i] ~ x) } setwd("/u/maechler/Vorl/comput-statist") ps.latex("ex-lm-rand-0.eps", width=7.6, height=8) ## FIXME: rather use a lattice xyplot() ! mult.fig(N, main="Regression Line is Random", mar = c(1,1,1,1)) for(i in 1:N) { plot(x, Y[,i], cex = 0.5, ylim = range(Y), axes = FALSE, frame = TRUE, xlab="", ylab="") # "Y" abline(2,1, col = "blue") abline(LM[[i]], lwd = 1.5, col = "red") } ps.end() ## Much nicer: using lattice xyplot() : library(lattice) str( d.ex <- data.frame(y = c(Y), ## one long vector run= as.factor(col(Y)), # coding the 9 different "runs" x = rep(x,N))) ps.latex("ex-lm-rand.eps", width=7.6, height=8) ##^^ uses "black & white" (black FG, white BG) xyplot(y ~ x | run, data = d.ex, type = c("p","r"), cex = 0.5, col.line = "red", lwd = 1.5, panel = function(x,y, ...) { panel.abline(2,1, col = "dark blue") panel.xyplot(x,y, ...) }, main = "Regression Line is Random") ps.end() ## Random ness shown in beta - space tit <- expression("True"~~ beta ~~ " and " ~~ hat(beta)~"s") plot(t(sapply(LM, coef)), main = tit, xlim = c(-1, 5), ylim = c(0.5, 1.5)) points(c(2, 1), col = "red", pch=8, cex=3, lwd=2) text (c(2, 1), expression(beta), col = "red", cex=2, adj = c(-1.5, 1)) ## More rand.regr <- function(N, n = 25, beta = c(2,1), sigma = 2) { x <- sort(round(runif(n, 1,10), 2)) Y <- matrix(NA, n, N) LM <- as.list(1:N) for(i in 1:N) { Y[,i] <- beta[1] + beta[2]* x + rnorm(n, sd = sigma) LM[[i]] <- lm(Y[,i] ~ x) } list(lm = LM, Y = Y) } dim(beta.hat <- sapply(rand.regr(100)\$lm, coef)) ## add these new beta hats to the plot points(t(beta.hat)) ## and if you want, repeat the above a few times ...