read.csv("halverson.data.091.csv")

hh <- read.csv("halverson.data.091.csv")

levels(hh$tt.color)

badData <- hh[is.na(hh$sort.order.1),]
badData

hh <- hh[!is.na(hh$sort.order.1),]
str(hh)

hh$corrected.shrivel.48

levels(hh$shrivel.48.count)

hh$corrected.shrivel.48.length <- nchar(as.character(hh$corrected.shrivel.48))

head(hh)

vx <- gsub("o","",hh$corrected.shrivel.48) # vx is 48 hr corrected shrivel data with only x's (no o's)
vx

hh$corrected.shrivel.48.lengthx <- nchar(vx) # length of vx

head(hh)

vo <- gsub("x","",hh$corrected.shrivel.48)
vo

hh$corrected.shrivel.48.lengtho <- nchar(vo)
head(hh)

hh$corrected.shrivel.48.lengthx + hh$corrected.shrivel.48.lengtho == hh$corrected.shrivel.48.length

y <- cbind(hh$corrected.shrivel.48.lengthx,hh$corrected.shrivel.48.lengtho)
y

summary1 <- aggregate(hh[,c("corrected.shrivel.48.lengthx","corrected.shrivel.48.length")],by=list(hh$paint.color),sum)
summary1$prop <- summary1$corrected.shrivel.48.lengthx / summary1$corrected.shrivel.48.length
summary1


summary2 <- aggregate(hh[,c("corrected.shrivel.48.lengthx","corrected.shrivel.48.length")],by=list(hh$cross.id),sum)
summary2$prop <- summary2$corrected.shrivel.48.lengthx / summary2$corrected.shrivel.48.length
summary2

summary3 <- aggregate(summary2$prop,by=list(hh$paint.color),na.rm=TRUE,mean)
cbind(summary1,summary3)


m1 <- glm(y~paint.color,family=binomial,data=hh)
summary(m1)
newd <- data.frame(paint.color = unique(hh$paint.color))
p1 <- predict(m1,newd,type = "response", se.fit = TRUE)
newd$fit <- p1$fit
newd$se.fit <- p1$se.fit
newd

barplot(newd$fit, names.arg=newd$paint.color)
# start of 24 hr corrected shrivel analysis, ask stuart about changing the definition of y?

hh$corrected.shrivel.24.length <- nchar(as.character(hh$corrected.shrivel.24))
hh$corrected.shrivel.24.length

sx <- gsub("o","",hh$corrected.shrivel.24) # sx is 24 hr corrected shrivel data with only x's (no o's)
sx
hh$corrected.shrivel.24.lengthx <- nchar(sx)

so <- gsub("x","",hh$corrected.shrivel.24)
hh$corrected.shrivel.24.lengtho <- nchar(so)
hh$corrected.shrivel.24.lengthx + hh$corrected.shrivel.24.lengtho == hh$corrected.shrivel.24.length

y <- cbind(hh$corrected.shrivel.24.lengthx,hh$corrected.shrivel.24.lengtho)
y

summary4 <- aggregate(hh[,c("corrected.shrivel.24.lengthx","corrected.shrivel.24.length")],by=list(hh$paint.color),sum)
summary4$prop <- summary4$corrected.shrivel.24.lengthx / summary4$corrected.shrivel.24.length
summary4


summary5 <- aggregate(hh[,c("corrected.shrivel.24.lengthx","corrected.shrivel.24.length")],by=list(hh$cross.id),sum)
summary5$prop <- summary5$corrected.shrivel.24.lengthx / summary5$corrected.shrivel.24.length
summary5

summary6 <- aggregate(summary5$prop,by=list(hh$paint.color),na.rm=TRUE,mean)
cbind(summary4,summary6)

m2 <- glm(y~paint.color,family=binomial,data=hh)
summary(m2)
newdd <- data.frame(paint.color = unique(hh$paint.color))
p2 <- predict(m2,newdd,type = "response", se.fit = TRUE)
newdd$fit <- p2$fit
newdd$se.fit <- p2$se.fit
newdd
barplot(newdd$fit, names.arg=newd$paint.color)




head(hh)