angustifolia_pallidaCrossData.csv <- file.choose()# import data
pacrosses <- read.csv(angustifolia_pallidaCrossData.csv, header=TRUE)####
pacrosses$crossID <- paste(pacrosses$maternalID,pacrosses$paternalID,sep="_")####
pacrosses$styleStatus.corrected <- gsub ("\\?","", pacrosses$styleStatus)#### remove ?'s
pacrosses$styleStatus.corrected2 <- gsub ("-","", pacrosses$styleStatus.corrected)#### remove -'s
pacrosses$cross.length <- nchar(as.character(pacrosses$styleStatus.corrected2))#### count total cross length
vx <- gsub("o","", pacrosses$styleStatus.corrected2)#### just compatible styles
vx####
pacrosses$compatible <- vx####
vo <- gsub("x","", pacrosses$styleStatus.corrected2)#### just incompatible styles
pacrosses$incompatible <- vo####
pacrosses$compatibleLength <- nchar (vx)#### length of compatible styles
pacrosses$incompatibleLength <- nchar (vo)#### length of incompatible styles
pacrosses$cross.length.corrected <- nchar(as.character(pacrosses$styleStatus.corrected2))####
pacrosses$compatibleLength + pacrosses$incompatibleLength == pacrosses$cross.length.corrected#### check on lengths
y <- cbind(pacrosses$compatibleLength,pacrosses$incompatibleLength)####
y####
summary1 <-aggregate(pacrosses[,c("compatibleLength","cross.length.corrected")], by=list(pacrosses$crossID),sum)#### to combine same crosses
summary1$prop <- summary1$compatibleLength / summary1$cross.length.corrected#### proportion compatible
summary1##
pacrosses$crosstype <- "PA"## sets default as pallid angustifolia cross
as.character(pacrosses$maternalID) == as.character(pacrosses$paternalID)##
self <- as.character(pacrosses$maternalID) == as.character(pacrosses$paternalID)##
pacrosses[self, "crosstype"] <- "S"## to sort self crosses
grepl("PAL", pacrosses$maternalID)## to find all id's with PAL
pallida <- grepl("PAL", pacrosses$maternalID)##
pacrosses$maternaltype <- "angustifolia"##
pacrosses[pallida, "maternaltype"] <- "pallida"##
pacrosses$maternaltype##
pacrosses$paternaltype <- "angustifoliap"##
grepl("PAL", pacrosses$paternalID)##
pallida2 <- grepl("PAL", pacrosses$paternalID)##
pacrosses[pallida2, "paternaltype"] <- "pallida"##
pacrosses$paternaltype##
pacrosses$maternaltype == pacrosses$paternaltype##
pallida_pallida<- pacrosses$maternaltype == pacrosses$paternaltype##
pacrosses[pallida_pallida, "crosstype"] <- "PP"## to find pallida pallida crosses
pacrosses$crosstype##
pacrosses$maternaltype2 <- "angustifolia"##
pacrosses[pallida, "maternaltype2"] <- "pallida"##
maternaltype2##
pacrosses$maternaltype2##
pacrosses$paternaltype2 <- "angustifolia"#
pacrosses[pallida2, "paternaltype2"] <- "pallidap"#
pacrosses$paternaltype2##
pacrosses$paternaltype2 == pacrosses$maternaltype2##
angustifolia_angustifolia <- pacrosses$paternaltype2 == pacrosses$maternaltype2##
pacrosses[angustifolia_angustifolia, "crosstype"] <- "AA"## to find angustifolia angustifolia crosses
pacrosses$crosstype##
pacrosses$maternalype3 <- "trueang"#
pacrosses$maternaltype3 <- "trueang"#
pacrosses$paternaltype3 <- "trueang"#
pacrosses[pallida2, "paternaltype3"] <- "angustifolia"#
pacrosses$paternaltype3#
pacrosses$maternaltype2 == pacrosses$paternaltype3#
pacrosses$maternaltype2 == pacrosses$paternaltype3#
angustifolia_pallida <- pacrosses$maternaltype2 == pacrosses$paternaltype3#
pacrosses[angustifolia_pallida, "crosstype"] <- "AP"# to find angustifolia pallida crosses
pacrosses$crosstype#
as.character(pacrosses$maternalID) == as.character(pacrosses$paternalID)#
self <- as.character(pacrosses$maternalID) == as.character(pacrosses$paternalID)#
pacrosses[self, "crosstype"] <- "S"# to find self crosses
pacrosses$crosstype
summary2 <-aggregate(pacrosses[,c("compatibleLength","cross.length.corrected")], by=list(pacrosses$crosstype),sum) # 
summary2 
summary2 <-aggregate(pacrosses[,c("compatibleLength","cross.length.corrected")], by=list(pacrosses$crosstype, pacrosses$crossID),sum) # to combine same crosses sorted by crosstype
summary2
summary2$prop <- summary2$compatibleLength / summary2$cross.length.corrected
summary2
summary2$compatibility <- "unsure"
summary2[summary2$prop>.8 & summary2$cross.length.corrected>5.9, "compatibility"] <-"compatible"
summary2[summary2$prop<.3 & summary2$cross.length.corrected>5.9, "compatibility"] <-"incompatible"
summary2[summary2$compatibleLength==5 & summary2$cross.length.corrected==5 & summary2$prop==1, "compatibility"] <-"compatible"
summary2[summary2$compatibleLength==0 & summary2$cross.length.corrected==5 & summary2$prop==0, "compatibility"] <-"incompatible"
summary2[, c(1,6)]
table(summary2[, c(1,6)])
summary2$compatibility2 <- "unsure"
summary2[summary2$cross.length.corrected>5.9 & summary2$prop<.35, "compatibility2"] <-"incompatible"
summary2[summary2$cross.length.corrected>5.9 & summary2$prop>.65, "compatibility2"] <-"compatible"
summary2[summary2$compatibleLength==5 & summary2$cross.length.corrected==5 & summary2$prop==1, "compatibility2"] <-"compatible"
summary2[summary2$compatibleLength==0 & summary2$cross.length.corrected==5 & summary2$prop==0, "compatibility2"] <-"incompatible"
table(summary2[, c(1,6)])
table(summary2[, c(1,7)])
summary2[summary2$Group.1 =="S" & summary2$prop<.25, "compatibility2"] <-"incompatible"
table(summary2[, c(1,7)])
summary2$compatible <- 0
summary2$incompatible <-0
summary2[summary2$compatibility =="compatible", "compatible"] <- 1
summary2[summary2$compatibility =="unsure", "compatible"] <- ""
summary2[summary2$compatibility =="incompatible", "incompatible"] <- 1
summary2[summary2$compatibility =="unsure", "incompatible"] <- ""
as.character(summary2$compatible)
as.integer(as.character(summary2$compatible))
as.character(summary2$incompatible)
as.integer(as.character(summary2$incompatible))
summary2$compat2 <-as.integer(as.character(summary2$compatible))
summary2$incompat2 <-as.integer(as.character(summary2$incompatible))
y2 <- cbind(summary2$compat2, summary2$incompat2)
summary2$crosstype <- as.factor(summary2$Group.1)
m1 <- glm(y2~crosstype, data=summary2, family=binomial)
m2 <- glm(y2~1, data=summary2, family=binomial)
anova(m2,m1, test= "Chi")
summary(m1)
summary3 <- summary2[summary2$crosstype!="S", ]
y3 <- cbind(summary3$compat2, summary3$incompat2)
m3 <- glm(y3~crosstype, data=summary3, family=binomial)
m4 <- glm(y3~1, data=summary3, family=binomial)
anova(m4,m3, test= "Chi")
summary(m3)
summary3$compatible2 <- 0
summary3$incompatible2 <-0
summary3[summary3$compatibility2 =="compatible", "compatible2"] <- 1
summary3[summary3$compatibility2 =="unsure", "compatible2"] <- ""
summary3[summary3$compatibility2 =="incompatible", "incompatible2"] <- 1
summary3[summary3$compatibility2 =="unsure", "incompatible2"] <- ""
as.character(summary3$compatible2)
as.integer(as.character(summary3$compatible2))
as.character(summary3$incompatible2)
as.integer(as.character(summary3$incompatible2))
summary3$compat3 <-as.integer(as.character(summary3$compatible2))
summary3$incompat3 <-as.integer(as.character(summary3$incompatible2))
y4 <- cbind(summary3$compat3, summary3$incompat3)
m5 <- glm(y4~crosstype, data=summary3, family=binomial)
m6 <- glm(y4~1, data=summary3, family=binomial)
anova(m6,m5, test= "Chi")
summary(m5)



summary4 <- summary3[!is.na(summary3$compat3), ]
y3.1 <- cbind(summary4$compat3, summary4$incompat3)
m3.1 <- glm(y3.1 ~ crosstype, data=summary4, family=binomial)
pred <- predict(m3.1, type="response",se.fit=T)
summary4$fit <- pred$fit


summary4$se <-pred$se.fit
unique(summary4[,c("crosstype","fit","se")])
look <- unique(summary4[,c("crosstype","fit","se")])
bp <- barplot(look$fit, names.arg= c("a", "b", "sdf", ""), ylim = c(0, 1),
ylab = "Mean seed count")
segments(bp, look$fit - look$se, bp, look$fit + look$se)