######################################################################
# Regression for GPA data                                            #
#                                                                    #
######################################################################

#Read in the data - This assumes that the file is in what R recognizes
#  as the current directory (folder)
gpa <- read.csv(file = "gpa.csv")
head(gpa)
# On my computer, I stored the file at c:\data. To specify this folder,
#   I could use the following code to read in the file. Notice the use of
#  "\\" rather than "\".
# gpa <- read.csv(file = "c:\\data\\gpa.csv")
# head(gpa)

# Other ways to read in data
# 1. Use file.choose() in the file argument. This will open a window allowing
#    the user to navigate to the folder/directory location.
# gpa <- read.csv(file = file.choose())
#
# 2. Read the file directly from the book's website
# gpa <- read.csv(file = "http://www.chrisbilder.com/stat705/section1/gpa.csv")
#
# 3. Within RStudio, select IMPORT DATASET from the ENVIRONMENT tab (usually in
#    the upper right window). Select FROM TEXT (BASE) and navigate to the location
#    of the file. Select import. RStudio will put the read.csv() code into the CONSOLE tab
#    (usually in the lower left window) when import is complete. This code can be
#    copied/pasted into a program for future use.

# Examples of accessing parts of a data frame
  # See names of variables
  names(gpa)

  # One observation for one variable
  gpa[1,1] # First observation of first variable [row #1, column #1]

  gpa[1,] # First observation, all variables

  # One variable
  gpa$HS.GPA
  gpa[,1]
  gpa[,"HS.GPA"]
  gpa["HS.GPA"]

  # First observation of first variable again
  gpa$HS.GPA[1]

  # Use of :
  1:10
  # First 10 observations of first variable
  gpa[1:10,1]
  # First 10 observations
  gpa[1:10,]


#Scatter plot
plot(x = gpa$HS.GPA, y = gpa$College.GPA, xlab = "HS GPA", ylab = "College GPA",
  main = "College GPA vs. HS GPA", xlim = c(0,4.5), ylim = c(0,4.5), col = "red",
  pch = 1, cex = 1.0, panel.first = grid(col = "gray", lty = "dotted"))

 
###########################################################################
# Regression model
    
  #Fit the simple linear regression model and save the results in mod.fit
  mod.fit <- lm(formula = College.GPA ~ HS.GPA, data = gpa)

  #A very brief look of what is inside of mod.fit - see the summary function for a better way
  mod.fit

  #See the names of all of the object components
  names(mod.fit)
  mod.fit$coefficients 
  options(width = 60)
  mod.fit$residuals
  mod.fit$fitted.values

  #Put some of the components into a data.frame object
  save.fit <- data.frame(gpa, College.GPA.hat = round(mod.fit$fitted.values,2), residuals = round(mod.fit$residuals,2))

  #Print contents of save.fit
  head(save.fit)

  #Summarize the information stored in mod.fit
  summary(object = mod.fit)

  # Calculate betahat0 and betahat1 using the equations directly
  x <- gpa$HS.GPA
  y <- gpa$College.GPA
  n <- nrow(gpa)
  betahat1 <- (sum(x*y) - 1/n * sum(y) * sum(x)) / (sum(x^2) - 1/n * sum(x)^2)
  betahat0 <- mean(y) - betahat1*mean(x)
  data.frame(betahat0, betahat1)

  # Show how to find an estimated college GPA
  mod.fit$coefficients # vector of length 2
  mod.fit$coefficients[1] # beta_hat0
  mod.fit$coefficients[2] # beta_hat1
  x <- 3.5 # A potential HS GPA
  mod.fit$coefficients[1] + mod.fit$coefficients[2]*x
  as.numeric(mod.fit$coefficients[1] + mod.fit$coefficients[2]*x) # Remove label

  pred.set <- data.frame(HS.GPA = 3.5)
  predict(object = mod.fit, newdata = pred.set)

  
  # Examine object oriented language aspects of R
  class(mod.fit)
  class(gpa)
  class(gpa$HS.GPA)
  methods(class = "lm")
  options(width = 70)
  methods(generic.function = "summary")


  plot(x = gpa$HS.GPA, y = gpa$College.GPA, xlab = "HS GPA", ylab = "College GPA", main = "College GPA vs. HS GPA", 
     xlim = c(0,4.5), ylim = c(0,4.5), col = "red", pch = 1, cex = 1.0, panel.first = grid(col = "gray", lty = "dotted"))
  curve(expr = mod.fit$coefficients[1] + mod.fit$coefficients[2]*x, xlim = c(min(gpa$HS.GPA),
    max(gpa$HS.GPA)), col = "blue", add = TRUE, n = 1000, lwd = 2)
  # Alternative way to use curve()
  # curve(expr = predict(object = mod.fit, newdata = data.frame(HS.GPA = x)), xlim = c(min(gpa$HS.GPA),
  #   max(gpa$HS.GPA)), col = "blue", add = TRUE, n = 1000, lwd = 2)


############################################################################
# Plot with y^ +- 2*sigma_epsilon

  # Examine first observation
  sum.fit <- summary(mod.fit)
  names(sum.fit)
  sum.fit$sigma^2  # sigma^2_epsilon
  low <- mod.fit$fitted.values[1] - 2*sum.fit$sigma
  up <- mod.fit$fitted.values[1] + 2*sum.fit$sigma
  data.frame(low,up)
  gpa[1,]

  plot(x = gpa$HS.GPA, y = gpa$College.GPA, xlab = "HS GPA", ylab = "College GPA", main = "College GPA vs. HS GPA", 
     xlim = c(0,4.5), ylim = c(0,4.5), col = "red", pch = 1, cex = 1.0, panel.first = grid(col = "gray", lty = "dotted"))
  curve(expr = mod.fit$coefficients[1] + mod.fit$coefficients[2]*x, from = min(gpa$HS.GPA), to = max(gpa$HS.GPA), col = "blue",
    add = TRUE, n = 1000, lwd = 2)

  low <- mod.fit$fitted.values - 2*sum.fit$sigma
  up <- mod.fit$fitted.values + 2*sum.fit$sigma

  segments(x0 = gpa$HS.GPA, y0 = low, x1 = gpa$HS.GPA, y1 = up, col = "darkgreen") 

  #Also can get sigma^2_epsilon with sum(mod.fit$residuals^2)/mod.fit$df.residual


############################################################################
# Find the estimated variances of the regression parameters

  names(sum.fit)
  sum.fit$coefficients
  class(sum.fit$coefficients)
  sum.fit$coefficients[,2]  # Standard deviation
  sum.fit$coefficients[,2]^2  # Variance

  # Variance-covariance matrix
  methods(generic.function = "vcov")
  vcov(mod.fit)


############################################################################
# Confidence intervals for beta's

  confint(object = mod.fit, level = 0.95)
   

############################################################################
# CI for E(y) and PI for y
     
  more.gpa <- data.frame(HS.GPA = c(2, 3, 4))
  predict(object = mod.fit, newdata = more.gpa, se.fit = TRUE, interval = "confidence", level = 0.95)
  predict(object = mod.fit, newdata = more.gpa, se.fit = TRUE, interval = "prediction", level = 0.95) 
  
  #Notice no se.fit or interval arguments
  save.pred1 <- predict(object = mod.fit, newdata = data.frame(HS.GPA = 2))
  save.pred1
  #Notice no se.fit argument
  save.pred2 <- predict(object = mod.fit, newdata = data.frame(HS.GPA = 2), interval = "confidence", level = 0.95)
  save.pred2
  save.pred2[,2]
             
  #CI and PI bands plot
  plot(x = gpa$HS.GPA, y = gpa$College.GPA, xlab = "HS GPA", ylab = "College GPA", main = "College GPA vs. HS GPA", 
       xlim = c(0,4.5), ylim = c(0,4.5), col = "black", pch = 1, cex = 1.0, panel.first = grid(col = "gray", lty = "dotted"))
  curve(expr = predict(object = mod.fit, newdata = data.frame(HS.GPA = x)), 
        col = "red", lty = "solid", lwd = 1, add = TRUE, xlim = c(min(gpa$HS.GPA), max(gpa$HS.GPA)))
  curve(expr =  predict(object = mod.fit, newdata = data.frame(HS.GPA = x), interval = "confidence", level = 0.95)[,2], 
        col = "darkgreen", lty = "dashed", lwd = 1, add = TRUE, xlim = c(min(gpa$HS.GPA), max(gpa$HS.GPA)))
  curve(expr =  predict(object = mod.fit, newdata = data.frame(HS.GPA = x), interval = "confidence", level = 0.95)[,3], 
        col = "darkgreen", lty = "dashed", lwd = 1, add = TRUE, xlim = c(min(gpa$HS.GPA), max(gpa$HS.GPA)))
  curve(expr =  predict(object = mod.fit, newdata = data.frame(HS.GPA = x), interval = "prediction", level = 0.95)[,2], 
        col = "blue", lty = "dashed", lwd = 1, add = TRUE, xlim = c(min(gpa$HS.GPA), max(gpa$HS.GPA)))
  curve(expr =  predict(object = mod.fit, newdata = data.frame(HS.GPA = x), interval = "prediction", level = 0.95)[,3], 
        col = "blue", lty = "dashed", lwd = 1, add = TRUE, xlim = c(min(gpa$HS.GPA), max(gpa$HS.GPA)))
  legend(locator(1), legend = c("Sample model", "95% CI", "95% PI"), col = c("red", "darkgreen", "blue"),
         lty = c("solid", "dashed", "dashed"), bty = "n", cex = 0.75)

 
#