y1 = c(8.4,12.8,9.6,9.8,8.4,8.6,8.9,7.9)
y2 = c(9.4,15.2,9.1,8.8,8.2,9.9,9.0,8.1)
y3 = c(9.8,12.9,11.2,9.9,8.5,9.8,9.2,8.2)
y4 = c(12.2,14.4,9.8,12.0,8.5,10.9,10.4,10)
y=c(y1,y2,y3,y4)

c0=rep(1,32)
e1=e2=e3=e4=e5=e6=e7=rep(0,8)
e1[1]=e2[2]=e3[3]=e4[4]=e5[5]=e6[6]=e7[7]=1
c1=c2=c3=rep(0,32)
c1[1:8]=1
c2[9:16]=1
c3[17:24]=1
b1=rep(e1,4) 
b2=rep(e2,4)
b3=rep(e3,4)
b4=rep(e4,4)
b5=rep(e5,4)
b6=rep(e6,4)
b7=rep(e7,4)
X = data.frame(c1,c2,c3,b1,b2,b3,b4,b5,b6,b7)

X=X[-15,]
y=y[-15]

attach(X)
lm.fit = lm(y~c1+c2+c3+b1+b2+b3+b4+b5+b6+b7, data=X)
anova(lm.fit)

detach(X)
 vv=vcov(lm.fit)
 V=vv[2:4,2:4]
 V= solve(V)
 bb=lm.fit$coef[2:4]
 xx = t(bb)%*%V%*%bb/3
 xx

#############################################################

Correct SS
        SS          F-ratio
tmt   12.8421       6.2218
sub   78.4220
res   13.7604 


y1 = c(8.4,12.8,9.6,9.8,8.4,8.6,8.9,7.9)
y2 = c(9.4,15.2,9.1,8.8,8.2,9.9,9.0,8.1)
y3 = c(9.8,12.9,11.2,9.9,8.5,9.8,9.2,8.2)
y4 = c(12.2,14.4,9.8,12.0,8.5,10.9,10.4,10)
y=c(y1,y2,y3,y4)

tmt=factor(c(rep(1,8),rep(2,8),rep(3,8),rep(4,8)))
blk=factor(rep(1:8,4))
X=data.frame(blk,tmt)


X=data.frame(X[-15,])
y=y[-15]

options(contrasts=c("contr.treatment","contr.poly"))
attach(X)


lm.fit= lm(y~tmt+blk,data=X) 

anova(lm.fit)

 vv=vcov(lm.fit)
 V=vv[2:4,2:4]
 V= solve(V)
 bb=lm.fit$coef[2:4]
 xx = t(bb)%*%V%*%bb/3
 xx



lm.fit= lm(y~blk+tmt,data=X)
anova(lm.fit)

 vv=vcov(lm.fit)
 V=vv[9:11,9:11]
 V= solve(V)
 bb=lm.fit$coef[9:11]
 xx = t(bb)%*%V%*%bb/3
 xx