########### R code for Chapter 2

### Load and plot new versions of Amherst and Mount Airy data

# this will upload directly from the website

par(cex=1.3)
Y1=read.csv('http://rls.sites.oasis.unc.edu/faculty/rs/source/Data/amh1.csv')
plot(Y1$Year,Y1$Temp,pch=20,xlab='Year',ylab='Temperature',
main='Annual Mean Temperature in Amherst, MA')
lines(Y1$Year,lm(Temp~Year,data=Y1)$fitted)

Y2=read.csv('http://rls.sites.oasis.unc.edu/faculty/rs/source/Data/mta1.csv')
plot(Y2$MtAiry,Y2$Charl,pch=20,xlab='Mount Airy',ylab='Charleston',
main='Summer Mean Temp. in Mount Airy and Charleston')
lines(Y2$MtAiry[order(Y2$MtAiry)],lm(Charl~MtAiry,data=Y2)$fitted[order(Y2$MtAiry)])

# The rest of the code uses the old versions of these datasets

amh=read.table('http://rls.sites.oasis.unc.edu/faculty/rs/source/Data/amherst.s',header=T)
cha=read.table('http://rls.sites.oasis.unc.edu/faculty/rs/source/Data/charles.dat',header=F)
names(cha)=c('Year','MtAiry','Charleston')

# initial code to produce figures in Ch2
par(mfrow=c(1,2),cex=1.3)
plot(amh$year,amh$temp,xlab='Year',ylab='Temperature',pch=20)
lines(amh$year,lm(temp~year,amh)$fitted,lw=3)
plot(cha$MtAiry,cha$Charleston,xlab='Mount Airy',ylab='Charleston',pch=20)
lines(cha$MtAiry,lm(Charleston~MtAiry,cha)$fitted,lw=3)

x=amh$year
y=amh$temp

# calculate initial statistics
n=length(x)
sumx=sum(x)
sumy=sum(y)
sumxx=sum(x*x)
sumxy=sum(x*y)
sumyy=sum(y*y)
xbar=sumx/n
ybar=sumy/n
sxx=sumxx-n*xbar*xbar
sxy=sumxy-n*xbar*ybar
syy=sumyy-n*ybar*ybar

b0=ybar
b1=sxy/sxx
print(c(b0,b1))
# [1] 7.39265432 0.01159718
#
# intercept in original parameterization
print(b0-b1*xbar)
# [1] 6.447484

# agrees with displayed formula near bottom of page 39 of Smith-Young notes

# Estimation of sigma
ssres=syy-sxy*sxy/sxx
s2=ssres/(n-2)
s1=sqrt(s2)
print(c(ssres,s2,s1))

par(mfrow=c(1,1))

# plot the figure 2.1(a)
plot(x,y,xlab='Year',ylab='Mean Temperature',pch=20)
yfit=b0+b1*(x-xbar)
lines(x,yfit)

# compute SEs of parameter estimates
se1=sqrt(s2/sxx)
se0=sqrt(s2*(1/n+xbar*xbar/sxx))
print(c(se0,se1))

# verify this agrees with results of lm command
lm1=lm(y~x)
summary(lm1)

# standardized residuals
sres=lm1$resid/summary(lm1)$sigma

# qq-plot
qqnorm(sres)
abline(a=0,b=1)


# Kolmogorov-Smirnov test
ks.test(sres,y='pnorm')

# Anderson-Darling test
library(nortest)
ad.test(sres)

# Other tests of normality

shapiro.test(sres)
cvm.test(sres)
lillie.test(sres)
pearson.test(sres)
sf.test(sres)


########## charleston data

cha=read.table('C:/Users/rls/aug20/UNC/STOR664/R-Code/charles.txt',header=F)

x=cha[,2]
y=cha[,3]

# calculate initial statistics
n=length(x)
sumx=sum(x)
sumy=sum(y)
sumxx=sum(x*x)
sumxy=sum(x*y)
sumyy=sum(y*y)
xbar=sumx/n
ybar=sumy/n
sxx=sumxx-n*xbar*xbar
sxy=sumxy-n*xbar*ybar
syy=sumyy-n*ybar*ybar

b0=ybar
b1=sxy/sxx
print(c(b0,b1))
# [1] 7.39265432 0.01159718
#
# intercept in original parameterization
print(b0-b1*xbar)
# [1] 6.447484

# agrees with displayed formula near bottom of page 39 of Smith-Young notes

# Estimation of sigma
ssres=syy-sxy*sxy/sxx
s2=ssres/(n-2)
s1=sqrt(s2)
print(c(ssres,s2,s1))


# compute SEs of parameter estimates
se1=sqrt(s2/sxx)
se0=sqrt(s2*(1/n+xbar*xbar/sxx))
print(c(se0,se1))

# verify this agrees with results of lm command
lm1=lm(y~x)
summary(lm1)