rm(list = ls())
ct=c(104,165,65,100,4,5,13,32,42,142,44,130,3,6,6,23)
s=c(1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2)
y=c(1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2)
r=c(1,1,2,2,3,3,4,4,1,1,2,2,3,3,4,4)

S=factor(s)
Y=factor(y)
R=factor(r)

m7=glm(ct ~ R:Y + R:S + Y:S,family=poisson)
summary(m7)
m5=glm(ct ~ R:Y + Y:S,family=poisson)
summary(m5)
m4=glm(ct ~ R:Y + R:S ,family=poisson)
summary(m4)
m3=glm(ct ~ R:Y + S,family=poisson)
summary(m3)


tab7=c(7,df.residual(m7),
deviance(m7),1-pchisq(deviance(m7),df.residual(m7)))
tab5=c(5,df.residual(m5),
deviance(m5),1-pchisq(deviance(m5),df.residual(m5)))
tab4=c(4,df.residual(m4),
deviance(m4),1-pchisq(deviance(m4),df.residual(m4)))
tab3=c(3,df.residual(m3),
deviance(m3),1-pchisq(deviance(m3),df.residual(m3)))

t(matrix(c(tab7,tab4,tab5,tab3),4,4))

g=c(1,1,2,2,1,1,4,4,1,1,2,2,1,1,4,4)
h=c(1,1,1,1,3,3,1,1,1,1,1,1,3,3,1,1)
G=factor(g)
H=factor(h)
eq2=glm(ct ~ G:H:Y + G:H:S + Y:S,family=poisson)
summary(eq2)
eq3=glm(ct ~ G:H:Y + G:S + Y:S,family=poisson)
summary(eq3)

r=c(1,1,2,2,3,3,4,4,1,1,2,2,3,3,4,4)
p=c(1,1,2,2,2,2,2,2,1,1,2,2,2,2,2,2)
c=c(2,2,1,1,2,2,2,2,2,2,1,1,2,2,2,2)
j=c(2,2,2,2,1,1,2,2,2,2,2,2,1,1,2,2)
o=c(2,2,2,2,2,2,1,1,2,2,2,2,2,2,1,1)
P=factor(p)
C=factor(c)
J=factor(j)
O=factor(o)


eq4=glm(ct ~ P:C:J:O:Y + P:C:J:O:S + Y:S,family=poisson)
summary(eq4)
eq5=glm(ct ~ P:C:J:O:Y + P:S + Y:S,family=poisson)
summary(eq5)
eq6=glm(ct ~ P:C:J:O:Y + C:S + Y:S,family=poisson)
summary(eq6)
eq7=glm(ct ~ P:C:J:O:Y + J:S + Y:S,family=poisson)
summary(eq7)
eq8=glm(ct ~ P:C:J:O:Y + O:S + Y:S,family=poisson)
summary(eq8)
eq9=glm(ct ~ P:C:J:O:Y + Y:S,family=poisson)
summary(eq9)

tab4=c(4,df.residual(eq4),deviance(eq4),
  1-pchisq(deviance(eq4),df.residual(eq4)))
tab5=c(5,df.residual(eq5),deviance(eq5),
  1-pchisq(deviance(eq5),df.residual(eq5)))
tab6=c(6,df.residual(eq6),deviance(eq6),
  1-pchisq(deviance(eq6),df.residual(eq6)))
tab7=c(7,df.residual(eq7),deviance(eq7),
  1-pchisq(deviance(eq7),df.residual(eq7)))
tab8=c(8,df.residual(eq8),deviance(eq8),
  1-pchisq(deviance(eq8),df.residual(eq8)))
tab9=c(9,df.residual(eq9),deviance(eq9),
  1-pchisq(deviance(eq9),df.residual(eq9)))

t(matrix(c(tab4,tab5,tab6,tab7,tab8,tab9),4,6))