############################################################ # Statistics and figures for # Grossenbacher et al. 2018, Madrono # Self-fertilization and herbivory in a rare alpine plant # in California, Claytonia megarhiza (Montiaceae) ############################################################ library(MASS) #neg. binomial models library(multcomp) #pre planned contrasts library(sciplot) #for Fig 2 barcharts library(RgoogleMaps) #for Fig 1 map df <- read.table("data.csv", sep=",", header=TRUE) ############################################################ ###### stats set fruits ################################### ############################################################ m1=glm(as.integer(set.fruit)~ treatment + block + treatment:block, family="poisson",data=df) #poisson distribution for count data summary(m1) #deviance residuals skewed (median not 0), suggesting not normally distributed with(m1, cbind(res.deviance = deviance, df = df.residual, p = pchisq(deviance, df.residual, lower.tail=FALSE))) #data do not fit model well, possibly due to overdispersion # compare neg binomial vs poisson (which is nested in the neg binomial model) summary(m1 <- glm.nb(as.integer(set.fruit) ~ treatment + block + treatment:block, data = df)) #negative binomial m3 <- glm(as.integer(set.fruit) ~ treatment + block + treatment:block, data = df, family = "poisson") # poisson pchisq(2 * (logLik(m1) - logLik(m3)), df = 1, lower.tail=FALSE) #neg binomial model much better, possibly because it copes with over-dispersion summary(aov(m1)) # neg binomial model main effects and interaction terms #preplanned contrasts with neg. binomial model df$trt.block <- interaction(df$treatment, df$block) levels(df$trt.block) aovRes <- glm.nb(as.integer(set.fruit) ~ trt.block, data = df) c1 <- c(1,-.25,-.25,-.25,-.25,1,-.25,-.25,-.25,-.25,1,-.25,-.25,-.25,-.25) #study question 1 c2 <- c(0,3,-1,-1,-1,0,3,-1,-1,-1,0,3,-1,-1,-1) #study question 2 c3 <- c(0,0,0,.5,-.5,0,0,0,.5,-.5,0,0,0,.5,-.5) #effect of flashing c4 <- c(0,0,1,-.5,-.5,0,0,1,-.5,-.5,0,0,1,-.5,-.5) #effect of netting cntrMat <- rbind(c1,c2,c3,c4) summary(glht(aovRes, linfct=mcp(trt.block=cntrMat), alternative="two.sided"), test=adjusted("none")) ############################################################ ###### stats percent leaf damage ######################### ############################################################ m1=glm(as.integer(percent.leaves.damaged*100)~ treatment + block + treatment:block, family="poisson",data=df) #poisson distribution for count data summary(m1) #deviance residuals skewed (median not 0), suggesting not normally distributed with(m1, cbind(res.deviance = deviance, df = df.residual, p = pchisq(deviance, df.residual, lower.tail=FALSE))) #data do not fit model well, possibly due to overdispersion # compare neg binomial vs poisson (which is nested in the neg binomial model) summary(m1 <- glm.nb(as.integer(percent.leaves.damaged*100) ~ treatment + block + treatment:block, data = df)) #negative binomial m3 <- glm(as.integer(percent.leaves.damaged*100) ~ treatment + block + treatment:block, data = df, family = "poisson") # poisson pchisq(2 * (logLik(m1) - logLik(m3)), df = 1, lower.tail=FALSE) #neg binomial model much better, possibly because it copes with over-dispersion summary(aov(m1)) # neg binomial model main effects and interaction terms #preplanned contrasts with neg. binomial model df$trt.block <- interaction(df$treatment, df$block) levels(df$trt.block) aovRes <- glm.nb(as.integer(percent.leaves.damaged*100) ~ trt.block, data = df) c1 <- c(1,-.25,-.25,-.25,-.25,1,-.25,-.25,-.25,-.25,1,-.25,-.25,-.25,-.25) #study question 1 c2 <- c(0,3,-1,-1,-1,0,3,-1,-1,-1,0,3,-1,-1,-1) #study question 2 c3 <- c(0,0,0,.5,-.5,0,0,0,.5,-.5,0,0,0,.5,-.5) #effect of flashing c4 <- c(0,0,1,-.5,-.5,0,0,1,-.5,-.5,0,0,1,-.5,-.5) #effect of netting cntrMat <- rbind(c1,c2,c3,c4) summary(glht(aovRes, linfct=mcp(trt.block=cntrMat), alternative="two.sided"), test=adjusted("none")) ################################################# ###### Fig 1 map ############################### ################################################# site=c(-119.31603, 37.73978) clme = read.table(file = "occurrences.csv", header=TRUE, sep=",") #get map extent lon.low=-2+min(clme$Longitude) lon.high=2+max(clme$Longitude) lat.low=-2+min(clme$Latitude) lat.high=2+max(clme$Latitude) pdf(file="Fig1.pdf", width=4,height=4) lat <- c(lat.low, lat.high) long <- c(lon.low,lon.high) center = c(mean(c(lat.low, lat.high)),mean(c(lon.low,lon.high))) zoom <-4 terrmap <- GetMap(center=center,zoom=zoom,maptype="terrain", destfile = "terrain.png") PlotOnStaticMap(terrmap) # Translate original data coords <- LatLon2XY.centered(terrmap, clme$Latitude, clme$Longitude, zoom) coords <- data.frame(coords) points(coords$newX, coords$newY, pch=21, cex=1, col="black", bg="white", lwd=0.8) # Translate original data site=c(-119.31603, 37.73978) site.coords <- LatLon2XY.centered(terrmap, site[2], site[1], zoom) site.coords <- data.frame(site.coords) points(site.coords$newX, site.coords$newY, pch=23, cex=1.2, col="black", bg="grey50", lwd=0.8) dev.off() ############################################ ###### Fig 2 ############################## ############################################ pdf(file="Fig2.pdf", width=8, height= 4) par(mfrow=c(1,2), cex.main=.8) #reorder factors trt<-factor(df$treatment,levels=c("No treatment","All animals excluded","Crawling insects excluded", "Control A", "Control B")) bargraph.CI( df$treatment, #categorical factor for the x-axis df$set.fruit, #numerical DV for the y-axis legend=F, x.leg=19, ylab="fruit per plant", xlab="", main="", space=0.5, cex.lab=1) bargraph.CI( df$treatment, #categorical factor for the x-axis df$percent.leaves.damaged*100, #numerical DV for the y-axis legend=F, x.leg=19, ylab="% leaves damaged per plant", xlab="", main="", space=0.5, cex.lab=1) dev.off()