##### R-code for Example 10.2. It writes the related data to the file
##### "example102.dat." It also makes "fig103.ps."


x = c(0,0,1,1,2,2,3,3,4,4)
x1 = x-mean(x)
n = length(x)
sxx = var(x)*(n-1)
y = matrix(scan("example102.dat"),ncol=11,byrow=T)

### Fit LS line for each profile

a = rep(0,22)
b = rep(0,22)
mse = rep(0,22)

for(i in 1:22){
   fit = lm(y[i,2:11] ~ x1)
   a[i] = fit$coefficients[1]
   b[i] = fit$coefficients[2]
   mse[i] = sum(fit$residuals^2)/fit$df.residual
}

a0hat = mean(a)
b0hat = mean(b)
sigma02hat = mean(mse)

postscript("fig103.ps",width=7.5,height=7.5,horizontal=F)

par(mfrow=c(2,2),mar=c(5,4,1,1))

j=seq(1,22)

U=a0hat+qt(0.9975,176)*sqrt(21/22/10*sigma02hat)
L=a0hat-qt(0.9975,176)*sqrt(21/22/10*sigma02hat)

plot(j,a,type="o",lty=1,pch=16,xlab="j",xlim=c(0,23),ylab=expression({hat(a)[j]}^"*"),
     ylim=c(190,215),mgp=c(2,1,0),cex=0.9,xaxt="n",yaxt="n")
lines(j,rep(U,22),lty=2,cex=0.9)
lines(j,rep(L,22),lty=2,cex=0.9)
axis(1,at=c(0,5,10,15,20),labels=c("0","5","10","15","20"),cex=.7)
axis(2,at=c(190,200,210,220),labels=c("190","200","210","220"),cex=.7)
title(xlab="(a)",cex=0.9)

U=b0hat+qt(0.9975,176)*sqrt(21/22/sxx*sigma02hat)
L=b0hat-qt(0.9975,176)*sqrt(21/22/sxx*sigma02hat)

plot(j,b,type="o",lty=1,pch=16,xlab="j",xlim=c(0,23),ylab=expression(hat(b)[j]),
     ylim=c(101.4,103.1),mgp=c(2,1,0),cex=0.9,xaxt="n",yaxt="n")
lines(j,rep(U,22),lty=2,cex=0.9)
lines(j,rep(L,22),lty=2,cex=0.9)
axis(1,at=c(0,5,10,15,20),labels=c("0","5","10","15","20"),cex=.7)
axis(2,at=c(101.5,102,102.5,103),labels=c("101.5","102","102.5","103"),cex=.7)
title(xlab="(b)",cex=0.9)

plot(j,mse,type="o",lty=1,pch=16,xlab="j",xlim=c(0,23),ylab=expression({hat(sigma)[j]}^2),
     ylim=c(0,3),mgp=c(2,1,0),cex=0.9,xaxt="n",yaxt="n")
axis(1,at=c(0,5,10,15,20),labels=c("0","5","10","15","20"),cex=.7)
axis(2,at=c(0,1,2,3),labels=c("0","1","2","3"),cex=.7)
title(xlab="(c)",cex=0.9)

f = rep(0,22)
for(i in 1:22){
   f[i] = mse[i]/(sum(mse)-mse[i])*21
}

U=qf(0.9975,8,168)
L=qf(0.0025,8,168)

plot(j,f,type="o",lty=1,pch=16,xlab="j",xlim=c(0,23),ylab=expression(F[j]),
     mgp=c(2,1,0), ylim=c(0,3.2),cex=0.9,xaxt="n",yaxt="n")
lines(j,rep(U,22),lty=2,cex=0.9)
lines(j,rep(L,22),lty=2,cex=0.9)
axis(1,at=c(0,5,10,15,20),labels=c("0","5","10","15","20"),cex=.7)
axis(2,at=c(0,1,2,3),labels=c("0","1","2","3"),cex=.7)
title(xlab="(d)",cex=0.9)

graphics.off()