# fruitfly code - part 2

data <- read.table("http://www.amstat.org/publications/jse/datasets/fruitfly.dat")
names(data) <- c("id","partners","type","longevity","thorax","sleep")
attach(data)

# define vectors for the colors (col) and plotting characters (pch)
col.partn <- 1*(partners==0)+2*(partners==1)+3*(partners==8)  
	# 1=black, 2=red, 3=green 
pch.type <- 1*(type==0) + 2*(type==1)+3*(type==9)
  	# 1='circle', 2='triangle', 3='plus'
	# type=0: pregnant, type=1: virgin, type=9: NA


############
# Method 1 #
############

# create dummy variables for the different categories:
# (type=0: pregnant, type=1: virgin, type=9: NA)
dummy.1.p <- (partners==1)*(type==0)*1  # 1 pregant female
dummy.1.v <- (partners==1)*(type==1)*1  # 1 virgin female
dummy.8.p <- (partners==8)*(type==0)*1  # 8 pregnant females
dummy.8.v <- (partners==8)*(type==1)*1  # 8 virging females
dummy.0   <- (partners==0)              # 0 females

# we fit a model with a common slope (effect of thorax length) and 
# different intercepts for the 5 groups, taking the group with 
# zero females as the baseline:
# y = alpha + beta*thorax + gamma.1.p*dummy.1.p + 
#     gamma.1.v*dummy.1.v + gamma.8.p*dummy.8.p +
#     gamma.8.v*dummy.8.v + epsilon
m6 <- lm(longevity ~ thorax + dummy.1.p + dummy.1.v + 
                     dummy.8.p + dummy.8.v)
summary(m6)

# vector with estimated coefficients:
m6$coef

# show fitted regression lines:

# 0 females:
par(mfrow=c(2,2))
plot(thorax[partners==0],longevity[partners==0],pch=3,col=1,
     main="0 partners",ylim=c(16,97),xlim=c(.64,.94), 
     xlab="thorax length", ylab="longevity")
abline(m6$coef[1], m6$coef[2])

# 1 female:
plot(thorax[partners==1],longevity[partners==1],pch=pch.type[partners==1],
     col=2,main="1 partners",ylim=c(16,97),xlim=c(.64,.94), 
     xlab="thorax length", ylab="longevity")
legend(0.635,97,c("pregnant","virgin"),pch=c(1,2),col=2)
abline(m6$coef[1]+m6$coef[3], m6$coef[2], col="red") # 1 pregnant female
abline(m6$coef[1]+m6$coef[4], m6$coef[2], col="red") # 1 virgin female

# 8 females:
plot(thorax[partners==8],longevity[partners==8],pch=pch.type[partners==8],
     col=3,main="8 partners",ylim=c(16,97),xlim=c(.64,.94), 
     xlab="thorax length", ylab="longevity")
legend(0.635,97,c("pregnant","virgin"),pch=c(1,2),col=3)
abline(m6$coef[1]+m6$coef[5], m6$coef[2], col="green") # 8 pregnant females
abline(m6$coef[1]+m6$coef[6], m6$coef[2], col="green") # 8 virgin females

# Note that all lines are parallel:
# Does this model allow for interaction between thorax and type of female?
# Does this model allow for interaction between thorax and number of females?
# Does this model allow for interaction between type of female and number of females?

# If there is no interaction, then the effect of type of female does not 
# depend on the number of females:
# This implies: gamma.1.p - gamma.1.v = gamma.8.p - gamma.8.v
# or equivalently: gamma.1.p = gamma.1.v + gamma.8.p - gamma.8.v

# Fitting this reduced linear model (see derivation on the board):
m7 <- lm(longevity ~ thorax + I(dummy.1.v + dummy.1.p) + 
                     I(dummy.8.p + dummy.1.p) + I(dummy.8.v - dummy.1.p))
anova(m7,m6)


############
# Method 2 #
############

dummy.1 <- (partners==1)*1
dummy.8 <- (partners==8)*1
dummy.v <- (type==1)*1

# we again fit a model with a common slope (effect of thorax length) and 
# different intercepts for the 5 groups, but we now use a 
# different parametrization:
# y = alpha + beta*thorax + gamma.1*dummy.1 + 
#     gamma.8*dummy.8 + gamma.v*dummy.v + gamma.8.v*dummy.8.v

m7 <- lm(longevity ~ thorax + dummy.1 + dummy.8 + dummy.v + dummy.8*dummy.v)
summary(m7)

# show fitted regression lines:

# 0 females:
par(mfrow=c(2,2))
plot(thorax[partners==0],longevity[partners==0],pch=3,col=1,
     main="0 partners",ylim=c(16,97),xlim=c(.64,.94), 
     xlab="thorax length", ylab="longevity")
abline(m7$coef[1], m7$coef[2])

# 1 female:
plot(thorax[partners==1],longevity[partners==1],pch=pch.type[partners==1],
     col=2,main="1 partners",ylim=c(16,97),xlim=c(.64,.94), 
     xlab="thorax length", ylab="longevity")
legend(0.635,97,c("pregnant","virgin"),pch=c(1,2),col=2)
abline(m7$coef[1]+m7$coef[3], m7$coef[2], col="red") # 1 pregnant female
abline(m7$coef[1]+m7$coef[3]+m7$coef[5], m7$coef[2], col="red") # 1 virgin female

# 8 females:
plot(thorax[partners==8],longevity[partners==8],pch=pch.type[partners==8],
     col=3,main="8 partners",ylim=c(16,97),xlim=c(.64,.94), 
     xlab="thorax length", ylab="longevity")
legend(0.635,97,c("pregnant","virgin"),pch=c(1,2),col=3)
abline(m7$coef[1]+m7$coef[4], m7$coef[2], col="green") # 8 pregnant females
abline(m7$coef[1]+m7$coef[4]+m7$coef[5]+m7$coef[6], m7$coef[2], col="green") # 8 virgin females

# We can test for interaction by looking at the p-value 
# in the summary table, or by using an incremental F-test:
m8 <- lm(longevity ~ thorax + dummy.1 + dummy.8 + dummy.v)
anova(m8,m7)