# Simulate data, estimate model, plot model set.seed(1234) x <- sample(x = 100, size = 30, replace = TRUE) beta0 <- 10 beta1 <- 2 sigmae <- 20 E.y <- beta0 + beta1 * x y <- rnorm(n = 30, mean = E.y, sd = sigmae) plot(x = x, y = y, main = "Population and sample regression model", xlim = range(x), ylim = c(-20, 210), pch = 1, cex = 1.0, panel.first = grid(col = "gray", lty = "dotted")) mod.fit <- lm(formula = y ~ x) mod.fit curve(expr = mod.fit$coefficients[1] + mod.fit$coefficients[2]*x, xlim = c(min(x), max(x)), col = "blue", add = TRUE, n = 1000, lwd = 1) curve(expr = beta0 + beta1 * x, xlim = c(min(x), max(x)), col = "red", add = TRUE, n = 1000, lwd = 1) ################################################################################ # Repeat 10 times sample.model <- function(y) { mod.fit <- lm(formula = y ~ x) beta0 <- mod.fit$coefficients[1] beta1 <- mod.fit$coefficients[2] curve(expr = beta0 + beta1 * x, xlim = c(min(x), max(x)), col = "blue", add = TRUE, n = 1000, lwd = 1) c(beta0, beta1) } set.seed(2574) y <- matrix(rnorm(n = 30*10, mean = E.y, sd = sigmae), nrow = 10, byrow = TRUE) plot(x = x, y = y[1,], type = "n", ylab="y", main = "Population and sample regression model", col = "red", xlim = range(x), ylim = c(0, 210), panel.first = grid(col = "gray", lty = "dotted")) est <- apply(X = y, MARGIN = 1, FUN = sample.model) curve(expr = beta0 + beta1 * x, xlim = c(min(x), max(x)), col = "red", add = TRUE, n = 1000, lwd = 1) options(width = 60) est apply(X = est, MARGIN = 1, FUN = mean) ################################################################################ # Repeat 100 times set.seed(2574) y <- matrix(rnorm(n=30*100, mean=E.y, sd=sigmae), nrow=100, byrow=T) plot(x = x, y = y[1,], type = "n", ylab="y", main = "Population and sample regression model", col = "red", xlim = range(x), ylim = c(0, 210), panel.first = grid(col = "gray", lty = "dotted")) est <- apply(X = y, MARGIN = 1, FUN = sample.model) curve(expr = beta0 + beta1 * x, xlim = c(min(x), max(x)), col = "red", add = TRUE, n = 1000, lwd = 1) apply(X = est, MARGIN = 1, FUN = mean) apply(X = est, MARGIN = 1, FUN = var) sigmae^2/sum((x-mean(x))^2) sigmae^2*(1/30+mean(x)^2/sum((x-mean(x))^2)) mod.fit <- lm(formula = y[1,] ~ x) summary(mod.fit) mod.fit <- lm(formula = y[10,] ~ x) summary(mod.fit)