# Read in the data with read.table function
data=read.table("http://www.stat.unm.edu/~ghuerta/stat574/carbo.txt",header=T)
data
data$carbohydrate
data$age
y=data$carbohydrate
x1=data$age
x2=data$weight
x3=data$protein
pairs(data)
plot(x1,y)
plot(x2,y)
plot(x3,y)
n=length(y)
# Lets do some calculations explicitly
# 
X=cbind(rep(1,n),x1,x2,x3)
XtX=t(X)%*%X
XtX
Xty=t(X)%*%y
Xty
inXtX=solve(XtX)
b=inXtX%*%Xty
# For sum of squares of residuals
ssr1=sum(y*y)-t(b)%*%t(X)%*%y
# If we wished for an estimate of sigma
sigma2=ssr1/(n-4)
#  R square
Shat0=(n-1)*var(y)
Rsq=(Shat0-ssr1)/Shat0
# To test the hypothesis that the response
# does not depend on age
X=cbind(rep(1,n),x2,x3)
b=solve(t(X)%*%X)%*%t(X)%*%y
ssr0=sum(y*y)-t(b)%*%t(X)%*%y
# now lets compute the F-statistic
F=((ssr0-ssr1)/1)/(ssr1/(n-4))
pvalue=1-pf(F,1,16)
table6.3=data
res=lm(carbohydrate ~ age + weight + protein,data=table6.3)
summary(res)
anova(res)
names(res)
plot(res$residuals)
qqnorm(res$residuals)
qqline(res$residuals)
residuals(res)
lm.influence(res)
dffits(res)
d=cooks.distance(res)
plot(d)
plot(res$fitted,res$residuals)
abline(h=0)
res1=lm(carbohydrate ~ age + protein,data=table6.3)
res0=lm(carbohydrate ~ protein,data=table6.3)
anova(res1)
anova(res0)
res2=lm(carbohydrate ~ weight + protein,data=table6.3)
res3=lm(carbohydrate ~ weight + protein + weight*protein, data=table6.3)
#  Here is something more on logistic regression
#  Beetle mortality data on page 127
y=c(6,13,18,28,52,53,61,60)
n=c(59,60,62,56,63,59,62,60)
x=c(1.6907,1.7242,1.7552,1.7842,1.8113,1.8369,1.8610,1.8839)
n_y=n-y
beetle.mat=cbind(y,n_y)
# Figure 7.3
plot(x,y/n,xlab="Dose",ylab="Proportion Killed")