# Demo code for Bootstrap 5


# --- Multiple regression

  Duncan <- read.table(":Data:duncan.dat",  header=T) 

  names(Duncan)
  labels(Duncan)
  attach(Duncan)

  pairs(Duncan)
  
  plot(income, prestige)
  identify(income, prestige, labels=labels(Duncan)[[1]])

  regr <- lm(prestige ~ income + education, data=Duncan)
  regr

  summary(regr)

  plot(regr)

  # now specify model differently to reveal difference of effects
  regr <- lm(prestige ~ I(income+education) + income, data=Duncan)
  regr
  summary(regr)  # difference is not significant
  plot(regr)     # role for outliers?

  # set aside those pesky outliers and the difference *is* signif
  regr.noout <- update(regr, subset=c(-6,-16,-27))
  summary(regr.noout)

  
  # make sure MASS still on search path
  search()
  library(MASS)
 
  # rob regr with nominal SE thinks not significant
  rob.regr <- rlm(prestige ~ I(income+education) + income, data=Duncan)
  summary(rob.regr)
 
  # check which observations are down-weighted
  n <- length(prestige)
  plot(1:n, rob.regr$w)
  identify(1:n, rob.regr$w, labels=labels(Duncan)[[1]])
  
  # each robust regression is a sequence of weighted regressions
  my.rregr <- function (data,indices) {
    data <- data[indices,]
    fit <- rlm(prestige ~ I(income+education) + income, data=data, maxit=50)
    coefficients(fit)
  }
  boot.rregr <- boot(Duncan, my.rregr, 500)
  boot.rregr

  # look close to normal
  plot(boot.rregr, index=2) 
  macintosh()  # opens another window on mac
  plot(boot.rregr, index=3) 

  # clearly *not* significant
  boot.ci(boot.rregr, index=2, type=c("norm", "perc", "bca"))
 
  boot.ci(boot.rregr, index=3, type=c("norm", "perc", "bca"))