#look for overall effects in card throws
cards=read.table("cards.txt",header=T)
attach(cards)
#plot main effects
par(mfrow=c(2,2))
boxplot(dis~hand,ylab="distance",xlab="hand")
boxplot(dis~rotation,ylab="distance",xlab="rotation")
boxplot(dis~aim,ylab="distance",xlab="aim")
par(mfrow=c(1,1))

#find effects
dis.lm=lm(dis~hand*rotation*aim)
2*coef(dis.lm)
#to assess significance, need to compute se
dis.sp=sqrt(sum(dis.lm$residuals^2)/40)
dis.se=dis.sp*2/sqrt(48)
dis.se*qt(.975,40)

#consider two throwers, with thrower effect
#(assessing significance without replications)
detach()
cards2=cards[c(17:24,33:40),]
#re-code the throwers to be numeric
cards2$thr=-1*(cards2$thrower=="Charae")+1*(cards2$thrower=="Julissa")
attach(cards2)
#main effects
par(mfrow=c(2,2))
boxplot(dis~hand,ylab="distance",xlab="hand")
boxplot(dis~rotation,ylab="distance",xlab="rotation")
boxplot(dis~aim,ylab="distance",xlab="aim")
boxplot(dis~thr,ylab="distance",xlab="thrower")
par(mfrow=c(1,1))
2*coef(lm(dis~hand*rotation*aim*thr))
#how can we tell which are significant?  
#first try assuming all third and higher order intractions are zero
# (these are elements 12-16 of the coefficient vector)
cards2.se=sqrt(mean((2*coef(lm(dis~hand*rotation*aim*thr))[12:16])^2))
cards2.se
cards2.se*qt(.975,5)
#another possible assumption is that most of the effects are zero
qqnorm(2*coef(lm(dis~hand*rotation*aim*thr))[2:16])
qqline(2*coef(lm(dis~hand*rotation*aim*thr))[2:16])
sort(2*coef(lm(dis~hand*rotation*aim*thr))[2:16])

coef(lm(dis~hand*rotation*aim*thr))
coef(lm(dis~hand))
#fitted model is  dis = 51 + 15.19*hand

#consider a different pair of two throwers, with thrower effect
detach()
cards2=cards[c(9:16,33:40),]
#re-code the throwers to be numeric
cards2$thr=-1*(cards2$thrower=="Minh")+1*(cards2$thrower=="Julissa")
attach(cards2)
#main effects
par(mfrow=c(2,2))
boxplot(dis~hand,ylab="distance",xlab="hand")
boxplot(dis~rotation,ylab="distance",xlab="rotation")
boxplot(dis~aim,ylab="distance",xlab="aim")
boxplot(dis~thr,ylab="distance",xlab="thrower")
par(mfrow=c(1,1))
2*coef(lm(dis~hand*rotation*aim*thr))
#but nothing looks significant, as the 3-way terms are as large as the rest
cards2.se=sqrt(mean((2*coef(lm(dis~hand*rotation*aim*thr))[12:16])^2))
cards2.se
cards2.se*qt(.975,5)
qqnorm(2*coef(lm(dis~hand*rotation*aim*thr))[2:16])
qqline(2*coef(lm(dis~hand*rotation*aim*thr))[2:16])
sort(2*coef(lm(dis~hand*rotation*aim*thr))[2:16])

coef(lm(dis~hand*rotation*aim*thr))
#to show what you would do with a bunch of significant terms, suppose that
# thrower, rotation*aim, rotation, hand*rotation*thrower, rotation*aim*thrower,
# and hand*aim are to be included in the model
#then fitted model is  dis = 59.38 + 5.69*hand - 12.5*rotation + 2.06*aim
#  - 17.38*thrower - 7.94*hand*rotation + 14.5*hand*aim - 14.31*rotation*aim
#  + 6.18*hand*thrower + 1.75*rotation*thrower + 4.3*aim*thrower
#  + 12.06*hand*rotation*thrower + 13.94*rotation*aim*thrower

#fitted value for Julissa, right hand, vertical, aim at top of door:
#   59.38 + 5.69 - 12.5 + 2.06 - 17.38 -7.94 + 14.5 - 14.31 + 6.18 + 1.75
#      + 4.3 + 12.06 + 13.94 = 67.73
#fitted value for Minh, right hand, vertical, aim at bottom of window:
#   59.38 + 5.69 - 12.5 - 2.06 + 17.38 -7.94 - 14.5 + 14.31 - 6.18 - 1.75
#      + 4.3 - 12.06 + 13.94 = 58.01
#all fitted values
fitted(lm(dis~hand*rotation*thrower+rotation*aim*thrower+hand*aim))
#longest prediction is for Minh, right hand, horizontal, aim at top of door