Class 13

Simulation Study 1

library(magrittr)
library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.4.0      ✔ purrr   1.0.0 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.5.0 
✔ readr   2.1.3      ✔ forcats 0.5.2 ── Conflicts ───────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract()   masks magrittr::extract()
✖ dplyr::filter()    masks stats::filter()
✖ dplyr::lag()       masks stats::lag()
✖ purrr::set_names() masks magrittr::set_names()

sim1 <- function(i){
  x1 <- rnorm(1000, 2)
  x2 <- rnorm(1000, -4)
  y <- 3 + 2 * x1 + 4 * x2 + rnorm(1000, sd = sqrt(2))
  df <- tibble(x1, x2, y)
  res <- df %$% lm(y ~ x1 + x2)
  return(c(coef(res),sigma(res)))
}
sim2 <- function(i){
  x1 <- rnorm(1000, 2)
  x2 <- rnorm(1000, -4)
  y <- 3 + 2 * x1 + 4 * x2 + rnorm(1000, sd = sqrt(2))
  df <- tibble(x1, x2, y)
  res <- df %$% lm(y ~ x1)
  return(c(coef(res),sigma(res)))
}

# sim1()
print("Correct Model")

[1] "Correct Model"

print("Average estimate")

[1] "Average estimate"

results1 <- sapply(1:1000, sim1)
rowMeans(results1)

(Intercept)          x1          x2             
   3.000889    1.996928    3.998753    1.413139 

print("Theoretical Standard Error Estimates")

[1] "Theoretical Standard Error Estimates"

apply(results1, 1, sd)

(Intercept)          x1          x2             
 0.20244465  0.04521681  0.04406489  0.03066468 

print("Incorrect Model")

[1] "Incorrect Model"

print("Average estimate")

[1] "Average estimate"

results2 <- sapply(1:1000, sim2)
rowMeans(results2)

(Intercept)          x1             
 -12.994301    1.998836    4.247681 

print("Theoretical Standard Error Estimates")

[1] "Theoretical Standard Error Estimates"

apply(results2, 1, sd)

(Intercept)          x1             
 0.28824676  0.13099029  0.09396826 

Simulation Study 2

sim3 <- function(i){
  x1 <- rnorm(1000, 8)
  # x2 <- rnorm(1000, -4)
  y <- 3 + 2 * log(x1)  + rnorm(1000, sd = sqrt(2))
  df <- tibble(x1, y)
  res <- df %$% lm(y ~ log(x1))
  return(c(coef(res),sigma(res)))
}
sim4 <- function(i){
  x1 <- rnorm(1000, 8)
  # x2 <- rnorm(1000, -4)
  y <- 3 + 2 * log(x1) + rnorm(1000, sd = sqrt(2))
  df <- tibble(x1, y)
  res <- df %$% lm(y ~ x1)
  return(c(coef(res),sigma(res)))
}


print("Correct Model")

[1] "Correct Model"

print("Average estimate")

[1] "Average estimate"

results1 <- sapply(1:1000, sim3)
rowMeans(results1)

(Intercept)     log(x1)             
   3.020812    1.989513    1.412930 

print("Theoretical Standard Error Estimates")

[1] "Theoretical Standard Error Estimates"

apply(results1, 1, sd)

(Intercept)     log(x1)             
 0.73301536  0.35310754  0.03167802 

print("Incorrect Model")

[1] "Incorrect Model"

print("Average estimate")

[1] "Average estimate"

results2 <- sapply(1:1000, sim4)
rowMeans(results2)

(Intercept)          x1             
  5.1104472   0.2539425   1.4129081 

print("Theoretical Standard Error Estimates")

[1] "Theoretical Standard Error Estimates"

apply(results2, 1, sd)

(Intercept)          x1             
 0.36313921  0.04505431  0.03127642 

Matrix Formulation

library(palmerpenguins)
penguins %<>% drop_na
ymat <- penguins$body_mass_g
print("flipper length and bill length")

[1] "flipper length and bill length"

xmat <- cbind(1, penguins$flipper_length_mm ,penguins$bill_length_mm)
solve(t(xmat)%*%xmat)%*%t(xmat)%*%ymat

             [,1]
[1,] -5836.298732
[2,]    48.889692
[3,]     4.958601

print("lm results")

[1] "lm results"

penguins %$% lm(body_mass_g ~ flipper_length_mm + bill_length_mm)

Call:
lm(formula = body_mass_g ~ flipper_length_mm + bill_length_mm)

Coefficients:
      (Intercept)  flipper_length_mm     bill_length_mm  
        -5836.299             48.890              4.959  

print("flipper length and island")

[1] "flipper length and island"

penguins %<>% mutate(island_torg = relevel(island, ref = "Torgersen"))
xmat <- penguins %$% model.matrix(~flipper_length_mm+island_torg)
solve(t(xmat)%*%xmat)%*%t(xmat)%*%ymat

                         [,1]
(Intercept)       -4889.32199
flipper_length_mm    44.88982
island_torgBiscoe   201.46081
island_torgDream    -63.90430

penguins %$% lm(body_mass_g ~ flipper_length_mm + island_torg)

Call:
lm(formula = body_mass_g ~ flipper_length_mm + island_torg)

Coefficients:
      (Intercept)  flipper_length_mm  island_torgBiscoe   island_torgDream  
         -4889.32              44.89             201.46             -63.90