%load_ext rmagic

%%R
# changes some default display options
options(show.signif.stars=F)
options(digits=3)

%%R
require(arm)
require(car)
require(lattice)

%%R -nd dll.data
dll.data <- read.csv("/Users/olsonran/Documents/dll.csv", header=TRUE)

%%R
dll.data$temp <- factor(dll.data$temp)
dll.data$replicate <- factor(dll.data$replicate)

%%R
dll.data$genotype <- relevel(dll.data$genotype, "wt")  
 # Setting the wild-type (wt) as reference

%%R
dll.data <- na.omit(dll.data) 

%%R
# The plot
plot(dll.data$SCT ~ dll.data$tarsus, 
      col=c("red", "blue")[dll.data$genotype], ylim=c(5,25), 
      pch=c(1,2)[dll.data$temp],cex=0.8, ylab="# of SCT", 
      xlab="Tarsus length", 
      main="scatterplot of Sex comb teeth and basi tarsus length")

# Add a legend    
legend(0.23,25, legend=c("Dll 25","Dll 30", "wt 25", "wt 30"),
      col=c("red", "red", "blue", "blue"), pch=c(1,2,1,2))

# lowess curves
lines(lowess(dll.data$tarsus[dll.data$genotype=="Dll" & dll.data$temp=="25"],
     dll.data$SCT[dll.data$genotype=="Dll" & dll.data$temp=="25"]), 
     lwd=2, type="o", col="red")

lines(lowess(dll.data$tarsus[dll.data$genotype=="Dll" & dll.data$temp=="30"],
     dll.data$SCT[dll.data$genotype=="Dll" & dll.data$temp=="30"]), 
     lwd=2, type="o", col="red", pch=2) 

lines(lowess(dll.data$tarsus[dll.data$genotype=="wt" & dll.data$temp=="25"],
     dll.data$SCT[dll.data$genotype=="wt" & dll.data$temp=="25"]), 
     lwd=2, type="o", col="blue")

lines(lowess(dll.data$tarsus[dll.data$genotype=="wt" & dll.data$temp=="30"],
     dll.data$SCT[dll.data$genotype=="wt" & dll.data$temp=="30"]), 
     lwd=2, type="o", col="blue", pch=2) 

%%R
regression.1 <- lm(SCT ~ tarsus, data=dll.data) # The intercept is implicit
regression.1.Intercept <- lm(SCT ~ 1 + tarsus, data=dll.data) 
  # With explicit intercept, but the same as above!

%%R
print(head(model.matrix(regression.1)))

%%R
print(model.matrix(regression.1)[500:508,])

%%R
print(coef(regression.1))

%%R
print(confint(regression.1))

%%R
print(summary(regression.1))

%%R
print(vcov(regression.1))

%%R
print(cov2cor(vcov(regression.1)))

%%R
print(coef(regression.1))

%%R
plot(SCT ~ tarsus, data=dll.data, xlim=c(0, 0.258))
abline(regression.1)

# We can look at the confidence region (confidence bands ) as well to see how this influences our uncertainty.

# To make confidence bands, we generally need to generate a set of "new" values for our explanatory variable (tarsus), and then use these to predict fitted values (and Confidence intervals on fitted values)

# We make the range go from 0, since we are fitting the intercept back to zero.  Clearly a crazy idea!

# make a new data frame with our "new" values of tarsus length
new_tarsus <- data.frame(tarsus = seq(0, range(dll.data$tarsus)[2], by=0.005) )   # Makes a new data frame

# "predict" values for our new values of tarsus
predicted.model.1  <- predict(regression.1, new_tarsus, interval="confidence")

#Then we add these lines back on
lines(x = new_tarsus[,1], y = predicted.model.1[,1], lwd=4) # this would be the same as abline(model.1)
lines(x = new_tarsus[,1], y = predicted.model.1[,3], lwd=3, lty=2, col="grey")
lines(x = new_tarsus[,1], y = predicted.model.1[,2], lwd=3, lty=2, col = "grey")

%%R
confidenceEllipse(regression.1)  
    # The estimates are highly correlated with one another.

%%R
dll.data$cent.tarsus <- dll.data$tarsus - mean(dll.data$tarsus)

%%R
par(mfrow=c(2,1))
hist(dll.data$tarsus, breaks=15)
hist(dll.data$cent.tarsus, breaks=15) #just changes the mean

%%R
print(mean(dll.data$tarsus))

%%R
print(mean(dll.data$cent.tarsus))

%%R
print(sd(dll.data$tarsus))

%%R
print(sd(dll.data$cent.tarsus))

%%R
regression.2 <-lm(SCT ~ cent.tarsus, data=dll.data)

par(mfrow=c(2,1))
plot(SCT~cent.tarsus,data=dll.data, cex=0.5, pch=16)
abline(regression.2)

new_tarsus.2 <- data.frame(cent.tarsus = seq(range(dll.data$cent.tarsus)[1], range(dll.data$cent.tarsus)[2], by=0.005) )   # Makes a new data frame

# "predict" values for our new values of tarsus
predicted.model.2 <- predict(regression.2, new_tarsus.2, interval="confidence")

#Then we add these lines back on
lines(x = new_tarsus.2[,1], y = predicted.model.2[,1], lwd=3) # this would be the same as abline(model.1)
lines(x = new_tarsus.2[,1], y = predicted.model.2[,3], lwd=2, lty=2)
lines(x = new_tarsus.2[,1], y = predicted.model.2[,2], lwd=2, lty=2)

confidenceEllipse(regression.2)

%%R
print(cov2cor(vcov(regression.2)))

%%R
print(summary(regression.2))

%%R
coefplot(regression.2, intercept=T, vertical=F, ylim=c(0, 35)) 

%%R
print(quantile(dll.data$tarsus))

%%R
print(range(dll.data$tarsus))

%%R
print(max(dll.data$tarsus) - min(dll.data$tarsus))

%%R
print(sd(dll.data$SCT))

%%R
print(confint(regression.2))

%%R
dll.data$cent.tarsus.microM <- 1000*dll.data$cent.tarsus

regression.3 <- lm(SCT ~ cent.tarsus.microM, data=dll.data )

print(coef(regression.3))

%%R
dll.data$tarsus.std <- as.numeric(scale(dll.data$tarsus))
regression.4 <- lm(SCT ~ tarsus.std, data=dll.data )
print(summary(regression.4)$coef[,1:2])

%%R
print(range(dll.data$tarsus.std))

%%R
print(sd(dll.data$tarsus))

%%R
plot(jitter(SCT) ~ tarsus.std , data=dll.data )
abline(regression.4)

new_tarsus.4 <- data.frame(tarsus.std = seq(-4.5, 4, by=0.01) )   # Makes a new data frame

# "predict" values for our new values of tarsus
predicted.model.4  <-predict(regression.4, new_tarsus.4 , interval="confidence")

#Then we add these lines back on
lines(x = new_tarsus.4[,1], y = predicted.model.4[,1], lwd=3) # this would be the same as abline(model.1)
lines(x = new_tarsus.4[,1], y = predicted.model.4[,3], lwd=2, lty=2)
lines(x = new_tarsus.4[,1], y = predicted.model.4[,2], lwd=2, lty=2)