# HW 1 Solutions
#
# DG 6.1
library(ISwR)
data(zelazo)
zelazo
attach(zelazo)
?zelazo
Y=c(active,passive,none,ctr.8w);Y
fg=factor(c(rep(1,6),rep(2,6),rep(3,6),rep(4,5)),labels=c("active","passive","none","ctr.8w"));fg
# note the following results
c(rep(1,6),rep(2,6),rep(3,6),rep(4,5))
c(rep("active",6),rep("passive",6),rep("none",6),rep("ctr.8w",5))
as.factor(c(rep("active",6),rep("passive",6),rep("none",6),rep("ctr.8w",5)))
as.numeric(as.factor(c(rep("active",6),rep("passive",6),rep("none",6),rep("ctr.8w",5))))
as.numeric(fg)

anova(lm0<-lm(Y~fg))
lm0
summary(lm0)
print(lm0)

plot(Y~fg,ylab="Age at Walking (months)",xlab="treatment")
# which is identical to this
plot(fg,Y)
# but slightly different here, since a different plot method is now applied.
plot(as.numeric(fg),Y)

t.test(active,ctr.8w)
t.test(active,c(none,ctr.8w)) # lost some significance but this pooling may have been most logical prior to seeing the data
plot(lm0) # normality assumptions OK

pairwise.t.test(Y,fg,p.adj = "bonf")
pairwise.t.test(Y,fg,p.adj = "fdr")  # neither of these pick up anything


#
# 10.3
#
a=gl(2,2,8);a
b=gl(2,4,8);b
gl(2,4,9)
x=1:8;x
y=c(1:4,8:5);y
z=rnorm(8);z

tx=tapply(x,list(a,b),mean);tx # note that x yields an additive model situations, +4 on columns and +2 on rows
ty=tapply(y,list(a,b),mean);ty # and that y yields an interaction model situation 
tz=tapply(z,list(a,b),mean);tz # z is noise, so the parameters should all have non-significant P values

a
b
tapply(x,a,mean);
tapply(x,b,mean);


model.matrix(~a*b)
model.matrix(~a+b+a:b) # this is equivalent to the one above it
model.matrix(~a:b)     # this design matrix is singular, the first column is the sum of the others
model.matrix(~a:b-1)   # this design matrix removes the first column (intercept), so the interaction parameters are the cell means
# e.g. 
summary(lm(x~a:b-1))

# let's first look at the non-degenerate design matrix results
summary(lm(x~a*b))
summary(lm(y~a*b)) # note that the upper left corner is the reference cell or intercept for x and y
summary(lm(z~a*b))

summary(lm(x~a:b))
tx
summary(lm(y~a:b))
ty
summary(lm(z~a:b))
tz
# The degeneracy of the design matrix manifests itself as a non-estimable 5th parameter (4 means => at most 4 estimable parameters ).
# Note that the sum of the first two parameters is the (1,1) cell mean, the sum of the first and third the (2,1) cell mean, the sum
# of the first and fourth the (1,2) cell mean, and that the first term equals the (2,2) cell mean, i.e. that NA=0 in these results.