##### Demonstration Script #5a
##### MATH322
#####
##### OLS Regression in R
##### just the basics




##### Preamble

# Load two important packages
library(lawstat)
library(lmtest)


# Load additional functionality
source("http://rfs.kvasaheim.com/rfs.R")



# Load data
dt = read.csv("http://rfs.kvasaheim.com/data/crime.csv")
attach(dt)



### Get to know the data

# Look at data
dt

# Look at variable names
names(dt)

# Look at data summary
summary(dt)



##### The First Model: 
#         Violent Crime in 1990 vs in 2000


mod1 = lm(vcrime00~vcrime90)   ## depVar ~ indepVar
summary(mod1)                  ## important output!

### Check the assumptions
e = residuals(mod1)      ## Calculate the residuals

normoverlay(e)           ## Do they look Normal-enough?
shapiro.test(e)          ## They are definitely Not Normal :(

### Close enough????
# Are e-bars Normally distributed???


ts = numeric()
for(i in 1:1e3) {
  x = sample(e1, size=51, replace=TRUE)
  ts[i] = sum(x)
}

normoverlay(ts)
shapiro.test(ts)   ### Yay!

# Now, we have evidence that the sample size is large enough
# that the sample sums are very close to Normal. Because of
# our formulas for b0, b1, etc., this is all we need.



### Residuals plot for the other three requirements

plot(vcrime90,e)
abline(h=0, lty="dotted")   ## horizontal line at y=0


### Statistical tests for the other three requirements

runs.test(e, orderBy=vcrime90)  ## constant expected value
bptest(mod1)                    ## constant variance
runs.test(e, orderBy=vcrime90)  ## independence of residuals


# Since all four requirements were met, this model is 
# appropriate. =)



### Results

summary(mod1)








##### The Second Model:
#         Property Crime in 1990 vs in 2000


### Fit the property crime model
mod2 = lm(pcrime00~pcrime90)

# calculate the residuals
e2 = residuals(mod2)

### Check requirements

normoverlay(e2)   ## Looks close enough
shapiro.test(e2)  ## Cannot reject the Normality hypothesis

plot(pcrime90,e2) ## Looks like constant expected value and variance
abline(h=0, lty=3)

runs.test(e2, orderBy=pcrime90)
bptest(mod2)



# Since all four requirements were met, this model is 
# appropriate. =)



### Results

summary(mod2)




### The End
#
#
#
#####