Modified from hybridgerm.Rmd (Oct 2017)

Guides

Bolker et al. 2009. Generalized linear mixed models: a practical guide for ecology and evolution. Trends in Ecology & Evolution.
Bolker. GLMM FAQ
Bolker et al. 2011. GLMMs in action: gene-by-environment interaction in total fruit production of wild populations of Arabidopsis thaliana. Part 1, Part 2
Bolker. In: Ecological Statistics: Contemporary Theory and Application
- GLMM:worked examples
- Reporting GLMM results
Bolker et al. 2013. Strategies for fitting nonlinear ecological models in R, AD Model Builder, and BUGS
- Bolker et al. 2012. Owls example: a zero-inflated, generalized linear mixed model for count data
- other examples from Non-Linear Model Working Group
A Practical Guide to Mixed Models in R
Schielzeth and Nakagawa 2012. Nested by design: model fitting and interpretation in a mixed model era. Methods in Ecology and Evolution.
Stroup. Rethinking the Analysis of Non-Normal Data in Plant and Soil Science
Elston et al. 2001. Analysis of aggregation, a worked example: Numbers of ticks on red grouse chicks
Annika Nelson. Ant GLMM Example
R package pscl vignette: Regression Models for Count Data in R

Purpose

Identify reproductive barriers between two sympatric moth-fflinated plant species, Schiedea kaalae and S. hookeri by fitting a generalized linear mixed model (GLMM).

In the experimental design, the following crosstypes were made:

within species, between population (may show outbreeding depression or heterosis)
within species, within populations (may show inbreeding depression)
hybrids between species (indicates species barrier from fflination to seed production)

In this analysis the response variable is the fflen viability of each cross. Other barriers (hybrid survival, flowering) could be analyzed in a similar framework, with appropriate changes to the underlying distribution.

Fixed effects:

crosstype - hybrids, within population, between populations
species - species of the maternal plant that produced the Viability

Potential random effects:

mompop - maternal plant population
mompid - maternal plant, specified by its population and ID
dadpop - paternal plant population

Data Import

ff <- read.table("firstflower6_plusone.csv", header=T, sep="\t")#has one dummy entry that did not flower for each crosstype
ff$firstflower.date[ff$firstflower.date==""] <- NA
ff$firstflower.date[ff$use.firstflower!="yes"] <- NA
ff$firstflower.date <- as.Date(ff$firstflower.date)
ff$firstflower <- as.integer(round(difftime(ff$firstflower.date, "2016-03-01")))
ff <- ff[ff$crossid!=107,] ##FIND OUT WHAT CROSS THIS IS

ff$alive[ff$use.alive.flowered!="yes"] <- NA
ff <- ff[ff$alive=="yes",] #only consider flowering of live individuals
ff$alive <- ff$alive=="yes"
ff$flowered[ff$flowered=="?"] <- NA
ff$flowered <- ff$flowered=="yes"
ff <- ff[!is.na(ff$flowered),]

crosses <- read.table("hybrids.csv", header=T, sep="\t", colClasses=c(mompop="factor", dadpop="factor"))

crosscol <- c("green","blue","orange","red")

#treat populations as factors
ff$mompop <- crosses$mompop[match(ff$crossid, crosses$crossid)]
ff$momid <- crosses$momid[match(ff$crossid, crosses$crossid)]
ff$species <- crosses$momsp[match(ff$crossid, crosses$crossid)]
ff$dadpop <- crosses$dadpop[match(ff$crossid, crosses$crossid)]
ff$dadid <- crosses$dadid[match(ff$crossid, crosses$crossid)]
ff$dadsp <- crosses$dadsp[match(ff$crossid, crosses$crossid)]
ff$crosstype <- crosses$crosstype[match(ff$crossid, crosses$crossid)]
ff$cross <- crosses$cross[match(ff$crossid, crosses$crossid)]

#rename crosstype codes
ff$crosstype <- factor(ff$crosstype, levels=c("between", "within", "hybrid"))
#made "between" the first reference level to facilitate comparison between outcrossing populations and hybridizing species 

ff$mompop <- sapply(ff$mompop, mapvalues, from = c("794","866","899","879","892","904","3587"), to = c("WK","WK","WK","879WKG","892WKG","904WPG","3587WP"))
ff$dadpop <- sapply(ff$dadpop, mapvalues, from = c("794","866","899","879","892","904","3587"), to = c("WK","WK","WK","879WKG","892WKG","904WPG","3587WP"))

#define interactions
ff <- within(ff, sxc <- interaction(species,crosstype))
ff <- within(ff, sxcxm <- interaction(species,crosstype,mompop,momid))
ff <- within(ff, mompid <- as.factor(paste(mompop,momid,sep=".")))
ff <- within(ff, dadpid <- as.factor(paste(dadpop,dadid,sep=".")))
ff <- within(ff, smompop <- as.factor(paste(species,mompop,sep="")))

#check final structure
str(ff)

   'data.frame':    1368 obs. of  39 variables:
    $ index                  : int  0 0 0 0 0 0 1 2 3 4 ...
    $ crossid                : int  1 4 6 108 115 117 1 1 1 1 ...
    $ plantid                : Factor w/ 33 levels "0","1","10","11",..: 33 33 33 33 33 33 2 15 26 27 ...
    $ crossid.plantid        : Factor w/ 1619 levels "","100-1","100-2",..: 1 1 1 1 1 1 99 155 201 221 ...
    $ death.date             : Factor w/ 133 levels "","100-2: 5/19/16",..: 1 1 1 1 1 1 1 1 1 1 ...
    $ firstflower.day        : Factor w/ 70 levels "","10/11","10/13",..: 1 1 1 1 1 1 3 23 1 26 ...
    $ firstflower.date       : Date, format: NA NA ...
    $ use.alive.flowered     : Factor w/ 3 levels "?","no","yes": 3 3 3 3 3 3 3 3 3 3 ...
    $ alive                  : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
    $ use.firstflower        : Factor w/ 4 levels "missed","never flowered",..: 3 3 3 3 3 3 4 4 1 4 ...
    $ flowered               : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
    $ saved.1                : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 2 1 ...
    $ saved.2                : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 2 1 ...
    $ sampled.VOC            : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
    $ biomass.inflo          : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 2 2 1 2 ...
    $ biomass.firstinflo     : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 2 2 2 2 ...
    $ use.fib                : Factor w/ 4 levels "double","#N/A",..: 3 3 3 3 3 3 4 4 4 4 ...
    $ delay                  : Factor w/ 35 levels "","0","1","10",..: 1 1 1 1 1 1 25 30 1 30 ...
    $ firstinflo.collect.date: Factor w/ 71 levels "","2016-07-12",..: 1 1 1 1 1 1 47 64 21 56 ...
    $ firstinflo.weigh.date  : Factor w/ 49 levels "","2017-07-26",..: 1 1 1 1 1 1 18 21 16 21 ...
    $ firstinflo.biomass.mg  : Factor w/ 1094 levels "","100","100.1",..: 1 1 1 1 1 1 509 630 1003 462 ...
    $ comments.fib           : Factor w/ 31 levels "","10/13/17?",..: 1 1 1 1 1 1 1 1 1 1 ...
    $ biomass.veg            : logi  NA NA NA NA NA NA ...
    $ comments.SS            : Factor w/ 44 levels "","?","\"1 terminlal infl\"",..: 1 1 1 1 1 1 1 1 30 1 ...
    $ comments.JP.SGW        : Factor w/ 78 levels "","105-9: 7/12/16 spray",..: 1 1 1 1 1 1 1 1 1 1 ...
    $ firstflower            : int  NA NA NA NA NA NA 226 271 NA 250 ...
    $ mompop                 : Factor w/ 5 levels "3587WP","WK",..: 3 3 3 4 4 4 3 3 3 3 ...
    $ momid                  : Factor w/ 17 levels "1","10","10-1",..: 7 7 7 1 1 1 7 7 7 7 ...
    $ species                : Factor w/ 2 levels "hook","kaal": 1 1 1 2 2 2 1 1 1 1 ...
    $ dadpop                 : Factor w/ 5 levels "3587WP","WK",..: 3 2 4 3 4 5 3 3 3 3 ...
    $ dadid                  : Factor w/ 23 levels "1","10","10-1",..: 19 7 1 8 14 7 19 19 19 19 ...
    $ dadsp                  : Factor w/ 2 levels "hook","kaal": 1 1 2 1 2 2 1 1 1 1 ...
    $ crosstype              : Factor w/ 3 levels "between","within",..: 2 1 3 3 2 1 2 2 2 2 ...
    $ cross                  : Factor w/ 4 levels "HH","HK","KH",..: 1 1 2 3 4 4 1 1 1 1 ...
    $ sxc                    : Factor w/ 6 levels "hook.between",..: 3 1 5 6 4 2 3 3 3 3 ...
    $ sxcxm                  : Factor w/ 510 levels "hook.between.3587WP.1",..: 195 193 197 24 22 20 195 195 195 195 ...
    $ mompid                 : Factor w/ 23 levels "3587WP.10","3587WP.14",..: 8 8 8 12 12 12 8 8 8 8 ...
    $ dadpid                 : Factor w/ 24 levels "3587WP.10","3587WP.14",..: 9 20 12 8 15 17 9 9 9 9 ...
    $ smompop                : Factor w/ 5 levels "hook879WKG","hookWK",..: 1 1 1 4 4 4 1 1 1 1 ...

Data Inspection

Replication

The sample sizes are unbalanced at all levels, including maternal population:

reptab <- with(ff, table(smompop,crosstype))
mosaic(reptab, pop=F)
labeling_cells(text = reptab, margin = 0)(reptab)

Replication is low for some within-population crosses. The replication is even lower for each maternal plant, so we need to be wary of estimates when subsetting at this level:

with(ff, kable(table(mompid,crosstype)))

	between	within	hybrid
3587WP.10	0	1	0
3587WP.14	17	0	8
3587WP.15	7	3	0
3587WP.7	25	16	16
3587WP.A	5	2	0
3587WP.C	16	1	0
879WKG.10-1	55	7	43
879WKG.2-2	41	54	108
879WKG.G-2	20	11	35
879WKG.H-2	0	1	8
879WKG.N-5	10	10	25
892WKG.1	29	6	11
892WKG.10	1	0	0
892WKG.2	3	0	0
892WKG.3	4	1	3
892WKG.4	1	0	1
892WKG.5	13	7	6
904WPG.2	17	10	17
904WPG.3	22	32	5
904WPG.5	82	36	28
WK.2	171	20	132
WK.2E- 1	7	0	34
WK.4	70	14	40

Overall data distribution

To identify the best-fitting distribution, we make quantile-quantile plots of the raw data against various distributions. The more points within the confidence interval envelopes, the better the fit. Later, we present quantile-quantile plots of the model residuals to assess model fit.

Fixed effects

Effects and interactions in these plots are simply given by the mean, which may be unduly influenced by high values.

intplot <- ggplot(ff,aes(fill=factor(flowered))) + geom_bar(position="fill")  + labs(y="Proportion") + guides(fill=guide_legend(title="Flowered")) + facet_grid(~species)
intplot + aes(x=crosstype)

Random effects

Maternal population

intplot + aes(x=mompop)

Maternal plant

intplot + aes(x=mompid)

Paternal population

intplot + aes(x=dadpop)

Models

We constructed the following models with the package glmmADMB. They all have the same fixed effects, species x crosstype, and response variable, flowered

X = standard GL(M)M

distribution, Random Effects:	None	Maternal plant	Maternal population
normal (norm)	X	X	X

sc.bin          <- glm(flowered~species*crosstype, data=ff, family="binomial")
sc.mix.mompid.bin         <- glmer(flowered~species*crosstype + (1|mompid), data=ff, family="binomial")
sc.mix.mompop.bin    <- glmer(flowered~species*crosstype + (1|mompop), data=ff, family="binomial")
sc.mix.momdadpid.bin         <- glmer(flowered~species*crosstype + (1|mompid) + (1|dadpid), data=ff, family="binomial")

Model comparison

AIC

We will use the Aikake Information Criterion to pick the model the best fits the data, penalized by the number of parameters. Differences of 2 units are significant.

#AICtab(sc.b, sc.mix.mompop.b, sc.mix.mompid.b)
sc.names <- c("sc.bin", "sc.mix.mompid.bin", "sc.mix.mompop.bin","sc.mix.momdadpid.bin")
sc.list <- sapply(sc.names, get, USE.NAMES=T)
sc.AIC <- ICtab(sc.list,mnames=sc.names,type="AIC", base=T, delta=F) # for AICc, nobs=nobs(sc.list[[1]])
class(sc.AIC)<-"data.frame"
all.names <- c(sc.names)
all.list <- sapply(all.names, get, USE.NAMES=T)
all.AIC <- dfun(rbind(sc.AIC))
all.AIC <- all.AIC[order(all.AIC$dAIC),]
kable(all.AIC, format.arg=list(digits=3))

	dAIC	df
sc.bin	0.00	6
sc.mix.mompid.bin	2.00	7
sc.mix.mompop.bin	2.00	7
sc.mix.momdadpid.bin	3.46	8

The best-fiting model is a model with the following components:

- response: flowered
- fixed effects: species, crosstype, species x crosstype
- random effect:

Coefficients

The coefficients estimated for each model agree qualitatively.

sc.log.names <- sc.names
sc.log <- sapply(sc.log.names, get, USE.NAMES=T)

coefplot2(sc.log, legend.x="topright",legend=T,legend.args=list(cex=0.8, xpd=T, inset=c(-0.1,0)), col.pts=sample(gg_color_hue(length(sc.log.names))), spacing=0.05, lwd.2=2, lwd.1=4, intercept=T)

Inference

We chose the model with nearly the best (lowest) AIC, to carry out inference tests and parameter estimation.

Description

mod <- sc.bin
print(mod)

   
   Call:  glm(formula = flowered ~ species * crosstype, family = "binomial", 
       data = ff)
   
   Coefficients:
                   (Intercept)                  specieskaal  
                        4.8176                       0.6672  
               crosstypewithin              crosstypehybrid  
                       -1.4765                      -2.0450  
   specieskaal:crosstypewithin  specieskaal:crosstypehybrid  
                        0.7279                      -0.5494  
   
   Degrees of Freedom: 1367 Total (i.e. Null);  1362 Residual
   Null Deviance:       354.4 
   Residual Deviance: 323.6     AIC: 335.6

Test significance of interaction

By dropping it from the model and performing a likelihood-ratio test, we see that the species x crosstype interaction is not significant.

sxc.chisq <- drop1(mod, test="Chisq") #load from file
dfun(sxc.chisq)

   Single term deletions
   
   Model:
   flowered ~ species * crosstype
                     Df Deviance   dAIC    LRT Pr(>Chi)
   <none>                 323.59 2.7283                
   species:crosstype  2   324.86 0.0000 1.2717   0.5295

Model summary

The model estimated the following parameters, with individual parameter significance determined by the Wald z-test, and fixed effect significance determined by analysis of deviance Wald test.

summary(mod)

   
   Call:
   glm(formula = flowered ~ species * crosstype, family = "binomial", 
       data = ff)
   
   Deviance Residuals: 
       Min       1Q   Median       3Q      Max  
   -3.3133   0.1269   0.1322   0.3482   0.3482  
   
   Coefficients:
                               Estimate Std. Error z value Pr(>|z|)    
   (Intercept)                   4.8176     0.5797   8.311  < 2e-16 ***
   specieskaal                   0.6672     1.1577   0.576 0.564384    
   crosstypewithin              -1.4765     0.7713  -1.914 0.055577 .  
   crosstypehybrid              -2.0450     0.6152  -3.324 0.000888 ***
   specieskaal:crosstypewithin   0.7279     1.6149   0.451 0.652171    
   specieskaal:crosstypehybrid  -0.5494     1.2625  -0.435 0.663414    
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   
   (Dispersion parameter for binomial family taken to be 1)
   
       Null deviance: 354.37  on 1367  degrees of freedom
   Residual deviance: 323.59  on 1362  degrees of freedom
   AIC: 335.59
   
   Number of Fisher Scoring iterations: 8

Anova(mod, type=3)

   Analysis of Deviance Table (Type III tests)
   
   Response: flowered
                     LR Chisq Df Pr(>Chisq)    
   species             0.3661  1  0.5451265    
   crosstype          17.5840  2  0.0001519 ***
   species:crosstype   1.2717  2  0.5294877    
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Predicted random effects

These are box and QQ (to check normality) plots of the estimated random effect of each maternal plant.

predre <- setNames(data.frame(RE=ranef(mod)$mompid,SD=ranef(mod, sd=T)$`1`),c("RE","SD"))
ggplot(predre, aes(x = rownames(predre),y=RE)) +
  geom_point(size = 2) + coord_flip()+
  geom_errorbar(aes(ymin = RE-SD, ymax = RE+SD)) + labs(x="Maternal plants", y="Predicted random effects")

#Count
reStack <- ldply(ranef(mod))
print( qqmath( ~`(Intercept)`|.id, data=reStack, scales=list(relation="free"),
                 prepanel = prepanel.qqmathline,
                 panel = function(x, ...) {
                   panel.qqmathline(x, ...)
                   panel.qqmath(x, ...)
                 },
                 layout=c(1,1)))

Least square means

The least square means procedure can generate predictor estimates of each type, and give their significance groupings with a post-hoc Tukey test. S. hookeri-produced hybrids produce less Viability than either crosses between or within S. hookeri populations. The other differences are not significant, but remember that the fixed effect of hybrid (vs. between) was significant (model summary).

#Count
rg <- ref.grid(mod)
#summary(rg)
sxc.lsm <- lsmeans(rg, ~ crosstype*species)
plot(sxc.lsm)

cld.mod <- cld(sxc.lsm, Letters=letters) #tukey letterings
library(boot)
cld.mod$response <- inv.logit(cld.mod$lsmean)
#cld.mod$SE[c(5,6)] <- 0
cld.mod$uSE <- inv.logit(cld.mod$lsmean+cld.mod$SE)
cld.mod$lSE <- inv.logit(cld.mod$lsmean-cld.mod$SE)
options(digits=4)
cld.mod[rev(order(cld.mod$species, cld.mod$crosstype)),]

    crosstype species lsmean     SE df asymp.LCL asymp.UCL .group response
    hybrid    kaal     2.890 0.4595 NA     1.990     3.791  ab      0.9474
    within    kaal     4.736 1.0044 NA     2.768     6.705  ab      0.9913
    between   kaal     5.485 1.0021 NA     3.521     7.449  ab      0.9959
    hybrid    hook     2.773 0.2062 NA     2.369     3.177  a       0.9412
    within    hook     3.341 0.5088 NA     2.344     4.338  ab      0.9658
    between   hook     4.818 0.5797 NA     3.681     5.954   b      0.9920
       uSE    lSE
    0.9661 0.9192
    0.9968 0.9766
    0.9985 0.9888
    0.9516 0.9287
    0.9792 0.9444
    0.9955 0.9858
   
   Results are given on the logit (not the response) scale. 
   Confidence level used: 0.95 
   Results are given on the log odds ratio (not the response) scale. 
   P value adjustment: tukey method for comparing a family of 6 estimates 
   significance level used: alpha = 0.05

H.wb <-  with(cld.mod[cld.mod$species=="hook",], response[crosstype=="within"]/response[crosstype=="between"] - 1)
K.wb <- with(cld.mod[cld.mod$species=="kaal",], response[crosstype=="within"]/response[crosstype=="between"] - 1)
K.H <- with(cld.mod[cld.mod$crosstype=="between",], response[species=="kaal"]/response[species=="hook"] - 1)
maxsp <- ifelse(K.H>0, "kaal","hook")
minsp <- ifelse(K.H<0, "kaal","hook")
maxresp <- with(cld.mod, response[species==maxsp & crosstype=="between"])
minresp <- with(cld.mod, response[species==minsp & crosstype=="between"])
HK.resp <-  with(cld.mod, response[species=="hook" & crosstype=="hybrid"])
KH.resp <-  with(cld.mod, response[species=="kaal" & crosstype=="hybrid"])
HK.int <-   with(cld.mod, ifelse(HK.resp > minresp & HK.resp < maxresp, (HK.resp-minresp)/(maxresp-minresp), 
                                 ifelse(HK.resp < minresp, HK.resp/minresp-1, HK.resp/maxresp-1)))
KH.int <-   with(cld.mod, ifelse(KH.resp > minresp & KH.resp < maxresp, (KH.resp-minresp)/(maxresp-minresp), 
                                 ifelse(KH.resp < minresp, KH.resp/minresp-1, KH.resp/maxresp-1)))

intermed <- (minresp + maxresp) / 2

round(c(H.wb,K.wb,K.H,HK.int,KH.int),2)

   [1] -0.03  0.00  0.00 -0.05 -0.04

ggplot(as.data.frame(cld.mod), aes(y=response, x=relevel(crosstype, "within"), fill=species)) +
  geom_col(position=position_dodge2()) +
  geom_linerange(aes(ymin=lSE, ymax=uSE), position=position_dodge(0.9)) +
  labs(x="", y="Proportion flowered",fill="Maternal species") +
  scale_fill_manual(labels = c("S. hookeri  ", "S. kaalae  "), values=brewer.pal(name="Set1", n=3)[c(3,2)]) +
  scale_x_discrete(labels = c("Intrapopulation", "Interpopulation", "Hybrid")) +
  geom_text(aes(label=.group), position=position_dodge(0.9), hjust=0, vjust=-0.2) +
  scale_y_continuous(expand = expand_scale(mult=c(0,0.04)), breaks = scales::pretty_breaks(n = 5), limits=c(0,1)) +
  theme_classic() + theme(legend.text=element_text(face="italic", size=rel(1)), legend.position="bottom", axis.text = element_text(colour="black", size=rel(1)), text=element_text(size=14), axis.ticks.x = element_blank()) + geom_segment(aes(x=2.5, y=intermed, xend=3.5, yend=intermed))

Hybrids Flowering GLMM No Add

John Powers

Jan 18, 2018