###############  Handout #2 for ST440/540  ###############

### Simple Linear Regression####################

########## Example 1 Page 15  #############
##################################

## Input data
x <- c(20,55,30)
y <- c(5,12,10)
#par(mfrow=c(2,2)) ##to print 2 by 2 graphs in the same plot
## Plot the data.
plot(x,y, main="Scatterplot of the Data")  


## Fit a Least Squares Line
myfit <- lm(y~x)
summary(myfit)

## print myfit
myfit

##calculate by formulas

coef<- function(a,b){ 
abar<-sum(a)/3
bbar<-sum(b)/3

deno<-sum((a-abar)^2)
nume<-(a-abar)%*%(b-bbar)
b1<-nume/deno
b0<-bbar-b1*abar
list(b0,b1)
}

coef(x,y)

## print just the estimated b_0 and b_1
myfit$coef

## Assign the values to b.0 and b.1
b.0 <- myfit$coef[1]
b.1 <- myfit$coef[2] 


## Plot the data with the fitted LS line
plot(x,y, main="Fitted Line Plot")  
abline(b.0, b.1)

##predicted value
pv<- b.0 + b.1*x
pv
  
 
## Assign the residuals to res
res <- myfit$residuals
res

##use formula to calculate residuals
res1<-y-pv
res1

## Plot the residuals by x
plot(x, res, main="Residuals Plot")  

## Sum of residuals equals 0
sum(res)

## Estimate sigma^2
sse <-  sum(res^2)
n <- length(res)
mse <- sse/(n-2)
mse


############## Example 2, P 36, 1.27 ################
## Import Data from local computer ##

## Read in data
ex.data <- read.table(file="W:/teaching/stat440540/data/CH1/CH01PR27.txt")
ex.data[1:10,]
## How many observations
n <- nrow(ex.data)
n

##scatterplot
names(ex.data)[1]<-"mass"
names(ex.data)[2]<-"age"
ex.data
plot(ex.data$age, ex.data$mass, main="Scatterplot of Muscle Mass")
plot(age, mass, main="Scatterplot of Muscle Mass") ##this doesn't work, since age and mass are not global name

## Define the variables globally
age <- ex.data$age
mass <- ex.data$mass
plot(age, mass, main="Scatterplot of Muscle Mass") 


## Fit a Least Square Line
myfit2 <- lm(mass~age)
myfit2
summary(myfit2)

## Plot the data with the fitted LS line
plot(age, mass, main="Fitted Line Plot")  
abline(myfit2$coef[1], myfit2$coef[2])

## Plot the residuals by x
plot(age, myfit2$residuals, main="Residuals Plot")