### Hargreaves, Langston, Johnson 2019 SAJB Aloe kraussii analyses ##### library(reshape) #need for reshaping data from wide to long (eg 'melt') library(lme4) #for mixed models library(car) #for Anova - to see results in Anova style table with type III estimates library(lsmeans) #for post-hoc contrasts of factor means (least squared means - package recently updated to emmeans) library(MASS) #for negative binomial mixed model (glm.nb) library(MuMIn) #for R squared for GLMM R.Version() setwd("~/Documents/Hargreaves etal 2019 SAJB Aloe kraussii data & Rcode") # DESIGN ############ #trial = ID of sunbird-feeding-on-plants trial (each trial lasted until each inflorescence visited) #trials 1-6 had opened & control flowers on each infloresence, trial 7 had only control flowers #plants randomly assigned to be male (stigmas removed) or female (anthers removed) #on each plant some flowers were experimentally widened/opened ('open') and some were left in their natural shape ('control') #trial - 1-6 had opened & control flowers on each infloresence, trial 7 had only control flowers. don't have data from trial 7. #1) WIDTH ############ #read in flower widths #1st columns = male plants, next group = paired female plants width1 <- read.csv("kraussii aviary flower widths.csv", skip=6) head(width1) #originally kept track of fertility (# anthers open on & stima state) in case needed to control for them later. Did try and keep these constant among trials but didn't use in analyses #anthers = # anthers that were freshly open #sitgma state = code describing age/state of stigma #cut into male vs female plants then restitch data back together #male widthM <- width1[,c('Date','Trial','plant','open.mm','control.mm')]; head(widthM) colnames(widthM) <- c('date','trial','plant','width_open','width_ct'); head(widthM) widthM$sex <- 'male' #female widthF <- width1[,c('Date','Trial','plant.1','open.mm.1','control.mm.1')]; head(widthF) colnames(widthF) <- c('date','trial','plant','width_open','width_ct'); head(widthF) widthF$sex <- 'female' widthF$plant <- as.factor(as.character(widthF$plant)) #paste back to 1 long data frame summary(widthM) summary(widthF) width2 <- rbind(widthF, widthM) head(width2) width2$plantID <- as.factor(paste(width2$sex, width2$plant, sep='')) #get widths into single column width3 <- melt(width2, id.vars=c('trial','sex','plantID'), measure.vars =c('width_open','width_ct')) summary(width3) width3$treat <- ifelse(width3$variable=='width_open', 'open', ifelse(width3$variable=='width_ct', 'ct', NA)) head(width3) width3$width <- width3$value width <- width3[,c('trial','sex','plantID','treat','width')] width$sex <- as.factor(width$sex) width$treat <- as.factor(width$treat) summary(width) dim(width) #244 flowers levels(width$plantID) #analysis: compare widths hist(width$width, col='grey') #data look ~normal width two modes (as expected) wlm <- lmer(width ~ treat*sex + (1|plantID), data=width) Anova(wlm, type='III') #int NS wlm.noX <- lmer(width ~ treat+sex + (1|plantID), data=width) anova(wlm, wlm.noX, test='Chisq') wlm.nosex <- lmer(width ~ treat + (1|plantID), data=width) anova(wlm.nosex, wlm.noX, test='Chisq') #sex NS wlm.notr <- lmer(width ~ sex + (1|plantID), data=width) anova(wlm.notr, wlm.noX, test='Chisq') #treat signif lsm.flwr <- data.frame(summary(lsmeans(wlm.noX, ~ treat))) lsm.flwr #1.8 vs 6.4 #2) POLLEN ON ANTHERS ############ #read in pollen on anthers (male plants only) #data generated from Elzone particle counter #FileName - name of file generated by particle counter #poll_treat - CT = control, OP = experimentally opened/widened, new = fresh anthers, newCT & newOP = plant prepared for experiment but not used anth <- read.csv("kraussii aviary pollen in anthers.csv", header=T) head(anth) summary(anth) #trial numbers - only 3 letters (lower case & upper case) but had 6 trials (+ trial 7 but no anthers collected) #many CT and OP flower with no trial number - trial# rubbed off vials during transport levels(anth$FileName) #AK = aloe kraussii, CT/OP/NEW = poll treatment, # = plant, letter=trial head(anth[,1:8], n=10L) anth$trial2 <- as.factor(ifelse(anth$trial=='b'|anth$trial=='B','b', ifelse(anth$trial=='c'|anth$trial=='C','c', ifelse(anth$trial=='d'|anth$trial=='D','d', '.')))) table(anth$trial, anth$poll.treat) table(anth$trial2, anth$poll.treat) anth$plant <- as.factor(anth$plant) levels(anth$plant) summary(anth) hist(anth$polperanther[anth$poll.treat=='OP'], col='grey') #poisson? almost normal hist(anth$polperanther[anth$poll.treat=='CT'], col='grey') #normal? #create flower-level random factor to deal with overdispersion anth$flwrID <- as.factor(seq.int(nrow(anth))) head(anth) #how are data distributed? panth.pf <- glm(round(polperanther) ~ poll.treat, family=poisson, data=anth, subset=(poll.treat=='CT'|poll.treat=='OP')) summary(panth.pf) #overdispersed #analysis 1A): compare pollen remaining on anthers from ct vs widened flowers panth.nb <- glm.nb(round(polperanther) ~ poll.treat, data=anth, subset=(poll.treat=='CT'|poll.treat=='OP')) panth.nbnull <- glm.nb(round(polperanther) ~ 1, data=anth, subset=(poll.treat=='CT'|poll.treat=='OP')) anova(panth.nb, panth.nbnull, test='Chisq') lsm.anth.nb <- data.frame(summary(lsmeans(panth.nb, ~ poll.treat, transform='response'))) lsm.anth.nb #CI dont overlap lsm.anth.nbSE <- data.frame(summary(lsmeans(panth.nb, ~poll.treat))); lsm.anth.nbSE lsm.anth.nbSE$btLSE <- exp(lsm.anth.nbSE$lsmean - lsm.anth.nbSE$SE) lsm.anth.nbSE$btUSE <- exp(lsm.anth.nbSE$lsmean + lsm.anth.nbSE$SE) lsm.anth.nbSE #analysis 1B): compare pollen remaining on ct flowers when half flowers widened vs none widened #cant do (no anthers from trial 7) #analysis 1C): does widening flowers knock off a bunch of pollen? levels(anth$poll.treat) table(anth$plant[!grepl('new', anth$poll.treat)], anth$poll.treat[!grepl('new', anth$poll.treat)]) table(anth$plant[anth$poll.treat=='CT'], anth$trial[anth$poll.treat=='CT']) panth2.p <- glmer(round(polperanther) ~ poll.treat + (1|plant) + (1|flwrID), family=poisson, data=anth, subset=(poll.treat=='newCT'|poll.treat=='newOP')) panth2.pnull <- glmer(round(polperanther) ~ (1|plant) + (1|flwrID), family=poisson, data=anth, subset=(poll.treat=='newCT'|poll.treat=='newOP')) anova(panth2.p, panth2.pnull, test='Chisq') #NS lsm.anth2 <- data.frame(summary(lsmeans(panth2.p, ~ poll.treat, transform='response'))) lsm.anth2 #n table(anth$poll.treat, anth$plant) table(anth$poll.treat) #3) POLLEN ON STIGMAS ############ #read in pollen on stigmas (female plants only) #trial 1-6 had opened (treat=opened) and control (treat=control) flowers on each inflorescens #trial 7 had only control flowers stig <- read.csv("kraussii aviary pollen on stigma.csv", header=T) head(stig) summary(stig) #make column to denote the type of experiment (other than trial 7 each plant had some opened and some control flowers) stig$exp <- as.factor(ifelse(stig$trial==7, 'none.widened', 'half.widened')) table(stig$exp, stig$trial) #how are data distributed? hist(stig$pollen[stig$treat=='opened'], col='grey') #poisson / negbin hist(stig$pollen[stig$treat=='control'], col='grey') #create flower-level random factor to deal with overdispersion stig$flwrID <- as.factor(seq.int(nrow(stig))) head(stig) #analysis 1A): compare pollen loads of ct vs widened flowers #if residual deviance > residual degrees of freedom, this is an indication of overdispersion #summary doest give this info for glmer - look at fixed model pstig.pf <- glm(pollen ~ treat, family=poisson, data=stig, subset=exp=='half.widened') summary(pstig.pf) #overdispersed pstig.nb <- glmer.nb(pollen ~ treat + (1|plant) + (1|trial), data=stig, subset=exp=='half.widened') pstig.nbnull <- glmer.nb(pollen ~ (1|plant) + (1|trial), data=stig, subset=exp=='half.widened') anova(pstig.nb, pstig.nbnull, test='Chisq') lsm.stig.nb <- data.frame(summary(lsmeans(pstig.nb, ~treat, transform='response'))) lsm.stig.nb #72 v 31 CI overlap (more) lsm.stig.nbSE <- data.frame(summary(lsmeans(pstig.nb, ~treat))); lsm.stig.nbSE lsm.stig.nbSE$btLSE <- exp(lsm.stig.nbSE$lsmean - lsm.stig.nbSE$SE) lsm.stig.nbSE$btUSE <- exp(lsm.stig.nbSE$lsmean + lsm.stig.nbSE$SE) lsm.stig.nbSE #analysis 1B): compare pollen loads on ct flowers when half flowers widened (trials 1-6) vs none widened (trial 7) pstigCT.nb <- glmer.nb(pollen ~ exp + (1|plant) + (1|trial), data=stig, subset=treat=='control') pstigCT.nbnull <- glmer.nb(pollen ~ (1|plant) + (1|trial), data=stig, subset=treat=='control') anova(pstigCT.nb, pstigCT.nbnull, test='Chisq') #NS lsm.stigCT.nb <- data.frame(summary(lsmeans(pstigCT.nb, ~ exp, transform='response'))) lsm.stigCT.nb #61 v 76 r.squaredGLMM(pstigCT.nb) ## PLOTTING ######## lsm.anth.nb lsm.stig.nb quartz(width=4.1, height=7) #3 plots (flower width, anthers, stigmas) par(mfrow=c(3,1), oma=c(4,5.4,1,1), mar=c(1,2,0,0), bty='l'); lylab=4.3; cxyax=1.3 colp='grey'; cxp=2.4; Ew=.06 #Ew=width of arrows at top of error bars xlim=c(.5,2.5); x=c(1,2); #flower width lsm.flwr ylimF=c(0,7); yF=lsm.flwr$lsmean; yU=lsm.flwr$upper.UCL; yL=lsm.flwr$lower.LCL #normal w CI yU= yF+lsm.flwr$SE; yL= yF-lsm.flwr$SE; yU; yL plot(yF ~ x, pch=20, cex=cxp, xlim=xlim, ylim=ylimF, xaxt='n', ylab=NA, las=1, cex.axis=cxyax) arrows(x0=c(1,2),x1=c(1,2), y0=yF, y1=yU, length=Ew, angle=90) arrows(x0=c(1,2),x1=c(1,2), y0=yF, y1=yL, length=Ew, angle=90) points(y=yF, x=x, pch=20, cex=cxp-.3, col=colp) axis(side=1, at=x, labels=NA) mtext(side=2, line=lylab, text='Width of corolla opening\n (mm)') #pollen remaining #ylimA=c(2000,11000); y=lsm.anth.nb$response; yU=lsm.anth.nb$asymp.UCL; yL=lsm.anth.nb$asymp.LCL #negbin w CI ylimA=c(4000,12000); y=lsm.anth.nb$response; yU=lsm.anth.nbSE$btUSE; yL=lsm.anth.nbSE$btLSE; #negbin w SE plot(y ~ x, pch=20, cex=cxp, xlim=xlim, ylim=ylimA, xaxt='n', ylab=NA, las=1, cex.axis=cxyax) arrows(x0=c(1,2),x1=c(1,2), y0=y, y1=yU, length=Ew, angle=90) arrows(x0=c(1,2),x1=c(1,2), y0=y, y1=yL, length=Ew, angle=90) points(y=y, x=x, pch=20, cex=cxp-.3, col=colp) axis(side=1, at=x, labels=NA) mtext(side=2, line=lylab, text='Pollen remaining\n (grains per anther)') #pollen deposited #ylimS=c(0,105);y=lsm.stig.nb$response; yU=lsm.stig.nb$asymp.UCL; yL=lsm.stig.nb$asymp.LCL #negbin w CI ylimS=c(0,105);y=lsm.stig.nb$response; yU=lsm.stig.nbSE$btUSE; yL=lsm.stig.nbSE$btLSE; #negbin w SE plot(y ~ x, pch=20, cex=cxp, xlim=xlim, ylim=ylimS, xaxt='n', ylab=NA, las=1, cex.axis=cxyax) arrows(x0=c(1,2),x1=c(1,2), y0=y, y1=yU, length=Ew, angle=90) arrows(x0=c(1,2),x1=c(1,2), y0=y, y1=yL, length=Ew, angle=90) points(y=y, x=x, pch=20, cex=cxp-.3, col=colp) axis(side=1, at=x, labels=NA) mtext(side=2, line=lylab, text='Pollen deposited\n (grains per stigma)') mtext(side=1, line=1, at=c(1,2), text=c('unmanipulated','widened')) mtext(side=1, line=2.8, text='Corolla treatment', cex=1.2)