data coleman;
input verbal ses @@;
label verbal='Verbal test score' ses='Socioeconomic status';

Headbreadth <- scan(what=list(breadth=0,group=""))
141 English 148 English 132 English 138 English 154 English 142 English 150 English
146 English 155 English 158 English 150 English 140 English 147 English 148 English
144 English 150 English 149 English 145 English 133 Celt 138 Celt 130 Celt
138 Celt 134 Celt 127 Celt 128 Celt 138 Celt 136 Celt 131 Celt
126 Celt 120 Celt 124 Celt 132 Celt 132 Celt 125 Celt

Celt <- data.frame(Headbreadth$breadth,Headbreadth$group)
names(Celt)=c("Skull_Breadth","Skull_Origin")

#Plot the data
boxplot(Skull_Breadth ~ Skull_Origin,data=Celt,xlab="Skull Origin",ylab="Skull Breadth (mm)")
Celt.lm <- lm(Skull_Breadth ~ 0 + Skull_Origin,data=Celt)
summary(Celt.lm)

par(mfrow=c(2,2))
plot(Celt.lm)

par(mfrow=c(1,1))
Data  <- scan(what=list(verbal=0,ses=0))
37.01 7.20 23.30 -12.86 26.51 -11.71 35.20 0.92
36.51 12.32 34.90 4.77 40.70 14.28 33.10 -0.96
37.10 6.31 22.70 -16.04 33.90 6.16 39.70 10.62
41.80 12.70 31.80 2.66 33.40 -0.17 31.70 -10.99
41.01 9.85 43.10 15.03 37.20 -0.05 41.01 12.77

Coleman <- data.frame(Data$verbal,Data$ses)
names(Coleman) <- c("Verbal_Score","Socioeconomic_Status")
with(Coleman,plot(Verbal_Score~Socioeconomic_Status,xlab="Socioeconomic Status",ylab="Verbal Test Score"))
 
Coleman.lm <- lm(Verbal_Score~Socioeconomic_Status,data=Coleman)
summary(Coleman.lm)

#Add regression line using linear models object
abline(Coleman.lm,col="red")

#Page through diagnostic plots
plot(Coleman.lm)

#Prediction and confidence intervals
newdata=data.frame(Socioeconomic_Status=seq(min(Coleman$Socioeconomic_Status),max(Coleman$Socioeconomic_Status)))
Coleman.conf <- predict(Coleman.lm, newdata, interval="confidence",level=0.95,type="response")
Coleman.predict <- predict(Coleman.lm, newdata, interval="predict",level=0.95,type="response")
with(Coleman,plot(Socioeconomic_Status,Verbal_Score,xlab="Socioeconomic Status",ylab="Verbal Test Score"))
abline(Coleman.lm,col="black")
lines(newdata[,1],Coleman.conf[,2],col="red",lty=2)
lines(newdata[,1],Coleman.conf[,3],col="red",lty=2)
lines(newdata[,1],Coleman.predict[,2],col="blue",lty=2)
lines(newdata[,1],Coleman.predict[,3],col="blue",lty=2)