oyster <- read.table('oyster.txt') names(oyster) <- c('TRT','REP','INIT','FINAL') oyster$TRT <- factor(oyster$TRT) m.lm <- lm(FINAL~TRT+INIT, oyster) summary(m.lm) anova(m.lm) drop1(m.lm) X<-cbind(rep(c(0,1,0),c(0,4,16)), rep(c(0,1,0),c(4,4,12)), rep(c(0,1,0),c(8,4,8)), rep(c(0,1,0),c(12,4,4))) X summary(lm(oyster$FINAL~oyster$INIT+X)) XX<-cbind(X,rep(c(0,1),c(16,4))) tmpX<-cbind(XX,oyster$INIT) colnames(tmpX) <- c(paste0('TRT', c(1,2,3,4,5)), 'INIT') summary(lm(oyster$FINAL~tmpX-1)) SSE <- sum(m.lm$residuals^2) MSE <- SSE/14 B <- solve(t(tmpX)%*%tmpX)%*%t(tmpX)%*%oyster$FINAL COV <- MSE*solve(t(tmpX)%*%(tmpX)) sqrt(COV[1,1]+COV[5,5]-2*COV[1,5]) TRT1.coef = c(4, -1, -1, -1, -1) / 5 E.TRT1 <- sum(TRT1.coef * B[1:5]) S.TRT1 <- sqrt(sum(TRT1.coef %*% t(TRT1.coef) * COV[1:5,1:5])) E.TRT1 S.TRT1 E.TRT1 / S.TRT1 TRT1.X <- XX[,2:5] TRT1.X[1:4,]=1/4 TRT1.X SSE.1 <- sum(lm(oyster$FINAL~oyster$INIT+TRT1.X-1)$residuals^2) SSE.1 - SSE LVL.X <- cbind(rep(c(0,1,0),c(8,8,4)),rep(c(1,0,1),c(8,8,4))) SSE.2 <- sum(lm(oyster$FINAL~oyster$INIT+LVL.X-1)$residuals^2) SSE.2 - SSE m2.lm <- lm(FINAL~TRT+INIT+TRT*INIT, oyster) summary(m2.lm)