set.seed(2332)
nk<-4
m.A<-3
m.B<-5
n<-nk*m.A*m.B
A<-rep(gl(m.A,m.B,labels=c(0:(m.A-1))),nk)
B<-rep(gl(m.B,1,length=m.A*m.B,labels=c(0:(m.B-1))),nk)
table(A,B)

beta.A<-c(0,2,-1)
beta.B<-c(1,1,2,3,1)
beta.AB<-matrix(0,nrow=m.A,ncol=m.B)
beta.AB[2,2:5]<-0.1
beta.AB[3,2:5]<--0.1

Y<-beta.A[as.numeric(A)]+beta.B[as.numeric(B)]+diag(beta.AB[as.numeric(A),as.numeric(B)])+rnorm(n)*2

table(A,B)



fit1<-lm(Y~A)
summary(fit1)
n*mean(Y)^2

fit2<-lm(Y~-1+A)
summary(fit2)


fit3<-lm(Y~A+B)
summary(fit3)

fit4<-lm(Y~-1+A+B)
summary(fit4)

xtabs(Y~as.numeric(A)+as.numeric(B))

tapply(Y,INDEX=list(A=A,B=B),mean)
tapply(fitted(fit3),INDEX=list(A=A,B=B),mean)

anova(fit1)
anova(fit2)


fit.bal1<-lm(Y~A*B)
summary(fit.bal1)
anova(lm(Y~A*B))
anova(lm(Y~B*A))

drop1(fit.bal1,test='F')

fit.bal2<-lm(Y~-1+A*B)
anova(fit.bal1)
anova(fit.bal2)
summary(fit.bal2)

summary(lm(Y~A))
summary(lm(Y~-1+A))

anova(lm(Y~A))

anova(lm(Y~-1+A))

summary(lm(Y~A+B))
summary(lm(Y~-1+A+B))


fit5<-lm(Y~A*B)
summary(fit5)

fit6<-lm(Y~-1+A*B)
summary(fit6)

tapply(Y,INDEX=list(A=A,B=B),mean)
tapply(fitted(fit6),INDEX=list(A=A,B=B),mean)