#Program to analyze the data in Table2.1 for the continuation of Example 2.1 in Chapter 3
yd=read.table("C:\\Research\\Data\\table21.txt")
id  = yd[[1]]
age = yd[[2]]
sbp = yd[[3]]
con = rep(1,length(age))
dma = cbind(con, age, sbp)
cprod = t(dma)%*%dma
cprod  #cross products X'X, X'Y, Y'Y with the intercept
av.age=mean(age)
av.sbp=mean(sbp)
sd.age=sqrt(var(age))
sd.sbp=sqrt(var(sbp))
nage = age - av.age
nsbp = sbp - av.sbp
ndma = cbind(nage, nsbp)
ncprod=t(ndma)%*%ndma
ncprod  #corrected cross products
desc.age=summary(age)
desc.sbp=summary(sbp)
rbind(desc.age, desc.sbp)  #Descriptives statistics for age and sbp
rbind(sd.age, sd.sbp)
#Regression analysis with sbp as response and age as predictor
simple = sbp ~ age
lreg = lm(simple)
lreg
summary.lm(lreg)
anova.lm(lreg)
confint(lreg, parm=c(1, 2), level=0.95)
confint(lreg, parm=c(1, 2), level=0.99)
fitted=predict.lm(lreg)
c.int = predict(lreg, interval="confidence", level=0.95, type="response")
c.int = c.int[,-1]
p.int = predict(lreg, interval="prediction", level=0.95, type="response")
p.int = p.int[,-1]
co1 = c("conf.lwr","conf.upr")
co2 = c("pred.lwr","pred.upr")
colnames(c.int) = co1
colnames(p.int) = co2
residual = residuals(lreg)
c1.int = predict(lreg, se.fit=TRUE, interval="confidence", level=0.95, type="response")$se.fit
p1.int = predict(lreg, se.fit=TRUE, interval="prediction", level=0.95, type="response")$se.fit
cbind(age,sbp,fitted,residual,c1.int,c.int,p1.int,p.int)