# importing data
sp <- read.csv("http://blog.lib.umn.edu/wage0005/echinacea/hhelianthoides.csv")


# the data
#sp$headId = identification number assigned to flower head
#sp$row = which row of styles (1=bottom row)
#sp$revRow = which row of styles (1=top)
#sp$treatment = pollination treatment applied
#sp$fristDayFlw = first day flower head put out anthers
#sp$styleDate = date when styles emerged
#sp$dateTreat = date when treatment was first applied
#sp$styleShriv = whether styles were shriveled on first visit after treatment (0=no, 1-yes)
#sp$dateShriv = date when styles were observed to have shriveled
#sp$daysPerstLong = large estimate of how long styles persisted
#sp$daysPerstAvg = middle estimate of how long styles persisted
#sp$daysPerstShort = small estimate of how long styles persisted
#sp$lastDayFlw = last day flower head put out anthers
#sp$endShrivLong = large estimate of how many days styles persisted after last day of flowering (0=shriveled on last day)
#sp$endShrivAvg = middle estimate of how many days styles persisted after last day of flowering
#sp$endShrivShort = small estimate of how many days styles persisted after last day of flowering
#sp$shrivNotes = notes about shriveling
#sp$repollinated = date when pollination treatment was reapplied
#sp$notes = other notes


# checking the data file
str(sp) #structure
sp$headId
levels(sp$treatment)


# how many rows received each treatment
table(sp$treatment)
barplot(table(sp$treatment))


# how many rows do heads have
x <- tapply(sp$row, sp$headId, max)
quartz()
hist(x)


# phenology
#This is not a representative sample of the population. For my experiment, I selected heads that flowered later than most, but early than some heads that started flowering at the beginning of August.
fs <- subset(sp, row == 1)
table(fs$firstDayFlw) #how many heads started flowering each day
quartz()
barplot(table(fs$firstDayFlw))
table(fs$lastDayFlw) #how many heads finished flowering each day
quartz()
barplot(table(fs$lastDayFlw))


# subsets for each treatment
sc <- subset(sp, treatment == "control")
so <- subset(sp, treatment == "outcross")
ss <- subset(sp, treatment == "self")


# style persistence after treatment
summary(sp$styleShriv)
summary(sc$styleShriv)
summary(so$styleShriv)
summary(ss$styleShriv)

tapply(sp$styleShriv, sp$treatment, summary) #another way of getting the summaries above

hsp <- subset(sp, styleShriv == 1)
table(hsp$treatment)
prop.test(table(hsp$treatment), table(sp$treatment))


# length of style persistence
summary(sp$daysPerstAvg) #overall
hist(sp$daysPerstAvg, breaks= 0:18)
summary(sc$daysPerstAvg) #control
quartz()
hist(sc$daysPerstAvg, breaks= 0:18)
summary(so$daysPerstAvg) #outcross
quartz()
hist(so$daysPerstAvg, breaks= 0:18)
summary(ss$daysPerstAvg) #self
quartz()
hist(ss$daysPerstAvg, breaks= 0:18)

m1 <- glm(sp$daysPerstShort~sp$treatment, family = poisson)
m0 <- glm(sp$daysPerstShort~1, family = poisson)
anova(m0, m1, test = "Chisq")
summary(m1)
par(mfcol = c(2,2))
plot(m1)
par(mfcol = c(1,1))
mean(sp$daysPerstShort, na.rm = TRUE); var(sp$daysPerstShort, na.rm = TRUE)

m3 <- glm(sp$daysPerstLong~sp$treatment, family = poisson)
m4 <- glm(sp$daysPerstLong~1, family = poisson)
anova(m3, m4, test = "Chisq")
summary(m3)
par(mfcol = c(2,2))
plot(m3)
par(mfcol = c(1,1))
mean(sp$daysPerstLong, na.rm = TRUE); var(sp$daysPerstLong, na.rm = TRUE)


# length of style persistence after last day of flowering
summary(sp$endShrivAvg) #overall
quartz()
hist(sp$endShrivAvg, breaks= -12:12)
summary(sc$endShrivAvg) #control
quartz()
hist(sc$endShrivAvg, breaks= -12:12)
summary(so$endShrivAvg) #outcross
quartz()
hist(so$endShrivAvg, breaks= -12:12)
summary(ss$endShrivAvg) #self
quartz()
hist(ss$endShrivAvg, breaks= -12:12)
summary(aov(sp$endShrivAvg~sp$treatment))


# analysis of first four rows only


# subsets of first four rows only
ep <- subset(sp, row < 5)
ec <- subset(sc, row < 5)
eo <- subset(so, row < 5)
es <- subset(ss, row < 5)


# style persistence after treatment
summary(ep$styleShriv)
summary(ec$styleShriv)
summary(eo$styleShriv)
summary(es$styleShriv)

hep <- subset(ep, styleShriv == 1)
table(hep$treatment)
prop.test(table(hep$treatment), table(ep$treatment))


# length of style persistence
summary(ep$daysPerstAvg) #overall
quartz()
hist(ep$daysPerstAvg, breaks= 0:15)
summary(ec$daysPerstAvg) #control
quartz()
hist(ec$daysPerstAvg, breaks= 0:15)
summary(eo$daysPerstAvg) #outcross
quartz()
hist(eo$daysPerstAvg, breaks= 0:15)
summary(es$daysPerstAvg) #self
quartz()
hist(es$daysPerstAvg, breaks= 0:15)

m5 <- glm(ep$daysPerstShort~ep$treatment, family = poisson)
m6 <- glm(ep$daysPerstShort~1, family = poisson)
anova(m5, m6, test = "Chisq")
summary(m5)
par(mfcol = c(2,2))
plot(m5)
par(mfcol = c(1,1))
mean(ep$daysPerstShort, na.rm = TRUE); var(ep$daysPerstShort, na.rm = TRUE)

m7 <- glm(ep$daysPerstLong~ep$treatment, family = poisson)
m8 <- glm(ep$daysPerstLong~1, family = poisson)
anova(m7, m8, test = "Chisq")
summary(m7)
par(mfcol = c(2,2))
plot(m7)
par(mfcol = c(1,1))
mean(ep$daysPerstLong, na.rm = TRUE); var(ep$daysPerstLong, na.rm = TRUE)


# length of style persistence after last day of flowering
summary(ep$endShrivAvg) #overall
quartz()
hist(ep$endShrivAvg, breaks= -12:5)
summary(ec$endShrivAvg) #control
quartz()
hist(ec$endShrivAvg, breaks= -12:5)
summary(eo$endShrivAvg) #outcross
quartz()
hist(eo$endShrivAvg, breaks= -12:5)
summary(es$endShrivAvg) #self
quartz()
hist(es$endShrivAvg, breaks= -12:5)
summary(aov(ep$endShrivAvg~ep$treatment))


# analysis of last four rows only


# subsets of last four rows only
ep <- subset(sp, revRow < 5)
ec <- subset(sc, revRow < 5)
eo <- subset(so, revRow < 5)
es <- subset(ss, revRow < 5)


# style persistence after treatment
summary(ep$styleShriv)
summary(ec$styleShriv)
summary(eo$styleShriv)
summary(es$styleShriv)

hep <- subset(ep, styleShriv == 1)
table(hep$treatment)
prop.test(table(hep$treatment), table(ep$treatment))


# length of style persistence
summary(ep$daysPerstAvg) #overall
quartz()
hist(ep$daysPerstAvg, breaks= 0:15)
summary(ec$daysPerstAvg) #control
quartz()
hist(ec$daysPerstAvg, breaks= 0:15)
summary(eo$daysPerstAvg) #outcross
quartz()
hist(eo$daysPerstAvg, breaks= 0:15)
summary(es$daysPerstAvg) #self
quartz()
hist(es$daysPerstAvg, breaks= 0:15)

m5 <- glm(ep$daysPerstShort~ep$treatment, family = poisson)
m6 <- glm(ep$daysPerstShort~1, family = poisson)
anova(m5, m6, test = "Chisq")
summary(m5)
par(mfcol = c(2,2))
plot(m5)
par(mfcol = c(1,1))
mean(ep$daysPerstShort, na.rm = TRUE); var(ep$daysPerstShort, na.rm = TRUE)

m7 <- glm(ep$daysPerstLong~ep$treatment, family = poisson)
m8 <- glm(ep$daysPerstLong~1, family = poisson)
anova(m7, m8, test = "Chisq")
summary(m7)
par(mfcol = c(2,2))
plot(m7)
par(mfcol = c(1,1))
mean(ep$daysPerstLong, na.rm = TRUE); var(ep$daysPerstLong, na.rm = TRUE)


# length of style persistence after last day of flowering
summary(ep$endShrivAvg) #overall
quartz()
hist(ep$endShrivAvg, breaks= -1:12)
summary(ec$endShrivAvg) #control
quartz()
hist(ec$endShrivAvg, breaks= -1:12)
summary(eo$endShrivAvg) #outcross
quartz()
hist(eo$endShrivAvg, breaks= -1:12)
summary(es$endShrivAvg) #self
quartz()
hist(es$endShrivAvg, breaks= -1:12)
summary(aov(ep$endShrivAvg~ep$treatment))


# Do the top rows shrivel faster than the bottom rows?
plot(sp$row, sp$daysPerstAvg)
plot(sc$row, sc$daysPerstAvg)
plot(so$row, so$daysPerstAvg)
plot(ss$row, ss$daysPerstAvg)

m4 <- lm(daysPerstAvg ~ treatment + row + treatment:row, data = sp)
m3 <- lm(daysPerstAvg ~ treatment + row, data = sp)
m2 <- lm(daysPerstAvg ~ treatment      , data = sp)
m1 <- lm(daysPerstAvg ~             row, data = sp)
m0 <- lm(daysPerstAvg ~              1 , data = sp)

anova(m3, m4)

summary(m4)
cc <- coef(m4)

plot(daysPerstAvg ~ row, data = sp)
points(daysPerstAvg ~ row, data = sp[sp$treatment == "control",], pch = 19, col = "red")
points(daysPerstAvg ~ row, data = sp[sp$treatment == "outcross",], pch = 19, col = "blue")
points(daysPerstAvg ~ row, data = sp[sp$treatment == "self",], pch = 19, col = "green")
abline(cc[1], cc[4], col = "red", lwd = 2)
abline(cc[1] + cc[2], cc[4] + cc[5], col = "blue", lwd = 2)
abline(cc[1] + cc[3], cc[4] + cc[6], col = "green", lwd = 2)