# Simple Regression in R

data = read.csv("http://stat.wharton.upenn.edu/~khyuns/stat431/PennSexSurvey.txt")
attach(data) # Automatically extracts the individual measurements from data

x = HEIGHT
y = WEIGHT

# Fit regression model
model = lm(y ~ x)

coefficients(model) # coefficients of the model
plot(x,y,ylab="Weight",xlab="Height") # Gives a scatterplot
abline(model)       # And the fitted line

summary(model)      # summary of the fitted model
plot(model)         # Various plots about the model 

resid(model) # Residuals, for the current x
predict(model) # predicted values, based on the current x
predict(model,list(x=95)) # predicted values, for x = 95
confint(model) # confidence interval for the slopes of the model
predict(model,interval='confidence',level=0.95) # 95% confidence interval for the fitted value, based on current x
predict(model,interval='prediction',level=0.95) # 95% prediction interval for the fitted value, based on current x

# Volcano Demonstration
persp(volcano)
title("Maunga Whau\nOne of 50 Volcanoes in the Auckland Region.",font.main=4)

# Take a slice of the volcano
slice = 4
plot(volcano[slice,],type="l",main="Cross-Section of Maunga Whau Volcano",xlab="Distance",ylab="Altitude",col="green")

# Manually set up stations
nStations = 25
stations = sample(length(volcano[slice,]),nStations)
points(stations,volcano[slice,stations],pch=17,col="black")

# Take measurements each day and fit regression
maxSeismic = 5
baseSeismic = 2
precision = 1

par(mfrow = c(1,3))
nTotalDays = 365
slope.fit = rep(0,nTotalDays)
intercept.fit = rep(0,nTotalDays)
for (i in 1:nTotalDays) {
  # Volcano
  plot(volcano[slice,],type="l",main="Cross-Section of Maunga Whau Volcano",xlab="Distance",ylab="Altitude",col="green",lwd=2)
  points(stations,volcano[slice,stations],pch=17,col="black")
  
  # Measurements
  measurements = baseSeismic + maxSeismic*volcano[slice,stations]/max(volcano[slice,]) + rnorm(length(stations),0,sd=precision)
  plot(c(1,length(volcano[slice,])),c(4.5,8),type="n",main=paste("Seismic Activity on Day",i),xlab="Station Location",ylab="Seismic Activity")
  points(stations,rep(4.5,length(stations)), pch=17)
  points(stations,measurements)
  model = lm(measurements ~ stations)
  abline(model,lwd=2)
  
  slope.fit[i] = coefficients(model)[2]
  intercept.fit[i] = coefficients(model)[1]
  
  # Overall fits
   plot(c(1,length(volcano[slice,])),c(4.5,7.5),type="n",main=paste("Seismic Activity from Day 1 to Day",i),xlab="Station Location",ylab="Seismic Activity")
   for (j in 1:i) abline(a=intercept.fit[j],b=slope.fit[j])
}




# Cool plots
z = - 2 * volcano        # Exaggerate the relief
x = 10 * (1:nrow(z))   # 10 meter spacing (S to N)
y = 10 * (1:ncol(z))   # 10 meter spacing (E to W)
persp(x, y, z, theta = 120, phi = 15, scale = FALSE, axes = FALSE)

z0 <- min(z) - 20
z <- rbind(z0, cbind(z0, z, z0), z0)
x <- c(min(x) - 1e-10, x, max(x) + 1e-10)
y <- c(min(y) - 1e-10, y, max(y) + 1e-10)
fill <- matrix("green3", nrow = nrow(z)-1, ncol = ncol(z)-1)
fill[ , i2 <- c(1,ncol(fill))] <- "gray"
fill[i1 <- c(1,nrow(fill)) , ] <- "gray"
par(bg = "lightblue")
persp(x, y, z, theta = 120, phi = 15, col = fill, scale = FALSE, axes = FALSE)
title(main = "Maunga Whau\nOne of 50 Volcanoes in the Auckland Region.",
       font.main = 4)

fcol <- fill

zi <- volcano[ -1,-1] + volcano[ -1,-61] + volcano[-87,-1] + volcano[-87,-61]  ## / 4
fcol[-i1,-i2] <-terrain.colors(20)[cut(zi, stats::quantile(zi, seq(0,1, length.out = 21)), include.lowest = TRUE)]
persp(x, y, 2*z, theta = 110, phi = 40, col = fcol, scale = FALSE,ltheta = -120, shade = 0.4, border = NA, box = FALSE)

par(oldpar)