R/myboot.R

In reza-niazi/MLRpack: Introductory Pacakge

#' Title
#'
#' @param iter The amount of times to iterate
#' @param alpha The level of significance for our confidence interval
#'
#' @return Histograms of distribution of regression coefficients and confidnce intervals
#' @export
#'
#' @examples
myboot = function(iter,alpha)
{
  b0 = numeric(iter)
  b1=numeric(iter)
  b2=numeric(iter)
  b3 =numeric(iter)
  n=18  #Number of rows of data
  Y=cast$CASTINGS
  counter = 0
  for(i in 1:iter){
    indices = sample(x = 1:n, size = n, replace = TRUE)
    X = cbind(rep(1,18),cast$INCENTIVE[indices],cast$PDUMMY[indices],cast$INC_PDUM[indices])
    beta_hat = try(solve(t(X)%*%X)%*%t(X)%*%Y,silent=TRUE)

    if(inherits(beta_hat,"try-error")){
      counter = counter+1
    }
    else{
      b0[i] = beta_hat[1]
      b1[i] = beta_hat[2]
      b2[i] = beta_hat[3]
      b3[i] = beta_hat[4]
    }
  }
  hist(x = b0, probability = TRUE,
       main = "Bootstrapped beta0",
       xlab = "beta0 estimates",col="blue")
  hist(x = b1, probability = TRUE,
       main = "Bootstrapped beta1",
       xlab = "beta1 estimates",col="blue")
  hist(x = b2, probability = TRUE,
       main = "Bootstrapped beta2",
       xlab = "beta2 estimates",col="blue")
  hist(x = b3, probability = TRUE,
       main = "Bootstrapped beta3",
       xlab = "beta3 estimates",col="blue")

  beta0ave = sum(b0)/iter
  beta1ave = sum(b1)/iter
  beta2ave = sum(b2)/iter
  beta3ave = sum(b3)/iter

  ci_b0 = 1366*c(1,1)+c(-1,1)*qnorm(1-(1-alpha)/2)*sqrt(n/(n-1))*sqrt(sum((b0-beta0ave)**2)/(iter-1))
  ci_b1 = 6.217*c(1,1)+c(-1,1)*qnorm(1-(1-alpha)/2)*sqrt(n/(n-1))*sqrt(sum((b1-beta1ave)**2)/(iter-1))
  ci_b2 = 47.78*c(1,1)+c(-1,1)*qnorm(1-(1-alpha)/2)*sqrt(n/(n-1))*sqrt(sum((b2-beta2ave)**2)/(iter-1))
  ci_b3 = .03333*c(1,1)+c(-1,1)*qnorm(1-(1-alpha)/2)*sqrt(n/(n-1))*sqrt(sum((b3-beta3ave)**2)/(iter-1))

  model = lm(cast$CASTINGS~cast$INCENTIVE+cast$PDUMMY+cast$INC_PDUM)

  list = list(Num_singularities = counter,ci0 = ci_b0,ci1=ci_b1,
              ci2=ci_b2,ci3=ci_b3,summ=summary(model))
  names(list) = c("Number of Singularities","CI for beta0","CI for beta1","CI for beta2","CI for beta3","Linear Model Summary")
  print(list)
}