Hybrids First Inflo Biomass GLMM
John Powers
Feb 4, 2017
Modified from hybridInflo Biomass.Rmd (April 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
- 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-pollinated 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 pollination to seed production)
In this analysis the response variable is the biomass of the ofibspring produced by 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 Inflo Biomass
Potential random effects:
- mompop - maternal plant population
- mompid - maternal plant, specified by its population and ID
- dadpop - paternal plant population
Data Import
fib <- read.table("firstinflobiomass.csv", header=T, sep="\t",
colClasses=c(firstflower.date="Date", firstinflo.collect.date="Date", firstinflo.weigh.date="Date"), na.strings=c("#N/A", "#VALUE!"))
fib <- fib[fib$crossid!=107,] ##FIND OUT WHAT CROSS THIS IS
fib$firstflower.date[fib$use.firstflower!="yes"] <- NA
fib$firstinflo.biomass.mg[fib$use.fib!="yes"] <- NA
fib <- fib[!is.na(fib$firstinflo.biomass.mg),]
fib$alive[fib$use.alive.flowered!="yes"] <- NA
fib <- fib[!is.na(fib$alive),]
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
fib$mompop <- crosses$mompop[match(fib$crossid, crosses$crossid)]
fib$momid <- crosses$momid[match(fib$crossid, crosses$crossid)]
fib$species <- crosses$momsp[match(fib$crossid, crosses$crossid)]
fib$dadpop <- crosses$dadpop[match(fib$crossid, crosses$crossid)]
fib$dadid <- crosses$dadid[match(fib$crossid, crosses$crossid)]
fib$dadsp <- crosses$dadsp[match(fib$crossid, crosses$crossid)]
fib$crosstype <- crosses$crosstype[match(fib$crossid, crosses$crossid)]
fib$cross <- crosses$cross[match(fib$crossid, crosses$crossid)]
#rename crosstype codes
fib$crosstype <- factor(fib$crosstype, levels=c("between", "within", "hybrid"))
#made "between" the first reference level to facilitate comparison between outcrossing populations and hybridizing species
fib$mompop <- sapply(fib$mompop, mapvalues, from = c("794","866","899","879","892","904","3587"), to = c("WK","WK","WK","879WKG","892WKG","904WPG","3587WP"))
fib$dadpop <- sapply(fib$dadpop, mapvalues, from = c("794","866","899","879","892","904","3587"), to = c("WK","WK","WK","879WKG","892WKG","904WPG","3587WP"))
#define interactions
fib <- within(fib, sxc <- interaction(species,crosstype))
fib <- within(fib, sxcxm <- interaction(species,crosstype,mompop,momid))
fib <- within(fib, mompid <- as.factor(paste(mompop,momid,sep=".")))
fib <- within(fib, dadpid <- as.factor(paste(dadpop,dadid,sep=".")))
fib <- within(fib, smompop <- as.factor(paste(species,mompop,sep="")))
fib$firstflower <- as.integer(round(difftime(fib$firstflower.date, "2016-03-10")))
#check final structure
fib$mass <- fib$firstinflo.biomass.mg/1000 #convert mg to g
str(fib) 'data.frame': 1213 obs. of 40 variables:
$ index : int 1 2 3 4 5 7 11 12 13 14 ...
$ crossid : int 1 1 1 1 1 1 1 1 1 1 ...
$ plantid : Factor w/ 32 levels "0","1","10","11",..: 2 15 26 27 28 30 4 5 7 8 ...
$ crossid.plantid : Factor w/ 1618 levels "100-1","100-2",..: 98 154 200 220 270 372 101 104 107 110 ...
$ 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",..: 3 23 1 26 1 1 1 43 36 1 ...
$ firstflower.date : Date, format: "2016-10-13" "2016-11-27" ...
$ use.alive.flowered : Factor w/ 3 levels "?","no","yes": 3 3 3 3 3 3 3 3 3 3 ...
$ alive : Factor w/ 3 levels "?","no","yes": 3 3 3 3 3 3 3 3 3 3 ...
$ use.firstflower : Factor w/ 4 levels "missed","never flowered",..: 4 4 1 4 1 3 1 4 4 1 ...
$ flowered : Factor w/ 3 levels "?","no","yes": 3 3 3 3 3 3 3 3 3 3 ...
$ saved.1 : Factor w/ 2 levels "no","yes": 1 1 2 1 2 1 1 1 2 2 ...
$ saved.2 : Factor w/ 2 levels "no","yes": 1 1 2 1 2 1 1 1 2 2 ...
$ sampled.VOC : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 2 1 ...
$ biomass.inflo : Factor w/ 2 levels "no","yes": 2 2 1 2 1 2 2 2 1 1 ...
$ biomass.firstinflo : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
$ use.fib : Factor w/ 4 levels "?","double","no",..: 4 4 4 4 4 4 4 4 4 4 ...
$ delay : int 5 7 NA 7 NA NA NA 2 6 NA ...
$ firstinflo.collect.date: Date, format: "2016-10-18" "2016-12-04" ...
$ firstinflo.weigh.date : Date, format: "2017-08-18" "2017-08-24" ...
$ firstinflo.biomass.mg : num 34.6 45.9 84.5 31.1 39.9 ...
$ comments.fib : Factor w/ 28 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 30 1 30 35 30 1 1 30 ...
$ comments.JP.SGW : Factor w/ 76 levels "","105-9: 7/12/16 spray",..: 1 1 1 1 1 1 1 1 1 1 ...
$ mompop : Factor w/ 5 levels "3587WP","WK",..: 3 3 3 3 3 3 3 3 3 3 ...
$ momid : Factor w/ 17 levels "1","10","10-1",..: 7 7 7 7 7 7 7 7 7 7 ...
$ species : Factor w/ 2 levels "hook","kaal": 1 1 1 1 1 1 1 1 1 1 ...
$ dadpop : Factor w/ 5 levels "3587WP","WK",..: 3 3 3 3 3 3 3 3 3 3 ...
$ dadid : Factor w/ 23 levels "1","10","10-1",..: 19 19 19 19 19 19 19 19 19 19 ...
$ dadsp : Factor w/ 2 levels "hook","kaal": 1 1 1 1 1 1 1 1 1 1 ...
$ crosstype : Factor w/ 3 levels "between","within",..: 2 2 2 2 2 2 2 2 2 2 ...
$ cross : Factor w/ 4 levels "HH","HK","KH",..: 1 1 1 1 1 1 1 1 1 1 ...
$ sxc : Factor w/ 6 levels "hook.between",..: 3 3 3 3 3 3 3 3 3 3 ...
$ sxcxm : Factor w/ 510 levels "hook.between.3587WP.1",..: 195 195 195 195 195 195 195 195 195 195 ...
$ mompid : Factor w/ 22 levels "3587WP.10","3587WP.14",..: 8 8 8 8 8 8 8 8 8 8 ...
$ dadpid : Factor w/ 24 levels "3587WP.10","3587WP.14",..: 9 9 9 9 9 9 9 9 9 9 ...
$ smompop : Factor w/ 5 levels "hook879WKG","hookWK",..: 1 1 1 1 1 1 1 1 1 1 ...
$ firstflower : int 217 262 NA 241 NA NA NA 159 138 NA ...
$ mass : num 0.0346 0.0459 0.0845 0.0311 0.0399 ...Data Inspection
ggplot(fib, aes(x=firstflower, y=log10(mass), color=cross)) + geom_point(cex=0.5) + geom_smooth(se=F) + scale_color_manual(values=crosscol) + xlab("Days to first flower") + ylab("log First infloresence biomass (g)") `geom_smooth()` using method = 'loess' and formula 'y ~ x'
ggplot(fib, aes(x=log10(mass), fill=cross, color=cross)) + geom_density(alpha=0.1) + scale_fill_manual(values=crosscol) +scale_color_manual(values=crosscol) +xlab("log First infloresence biomass (g)")
ggplot(fib, aes(x=firstflower, fill=cross, color=cross)) + geom_density(alpha=0.1) + scale_fill_manual(values=crosscol) +scale_color_manual(values=crosscol) + xlab("Days to first flower")
ggplot(fib[fib$delay>0 & fib$delay <25,], aes(x=delay, y=log10(mass), color=cross)) + geom_point(alpha=0.8) + geom_smooth(se=F) + scale_color_manual(values=crosscol) `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Replication
The sample sizes are unbalanced at all levels, including maternal population:
reptab <- with(fib, 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(fib, kable(table(mompid,crosstype)))| between | within | hybrid | |
|---|---|---|---|
| 3587WP.10 | 0 | 1 | 0 |
| 3587WP.14 | 15 | 0 | 5 |
| 3587WP.15 | 7 | 3 | 0 |
| 3587WP.7 | 22 | 13 | 14 |
| 3587WP.A | 5 | 1 | 0 |
| 3587WP.C | 15 | 1 | 2 |
| 879WKG.10-1 | 49 | 7 | 34 |
| 879WKG.2-2 | 38 | 44 | 91 |
| 879WKG.G-2 | 20 | 10 | 31 |
| 879WKG.H-2 | 0 | 0 | 5 |
| 879WKG.N-5 | 7 | 8 | 17 |
| 892WKG.1 | 27 | 5 | 7 |
| 892WKG.10 | 1 | 0 | 0 |
| 892WKG.2 | 3 | 0 | 0 |
| 892WKG.3 | 4 | 1 | 2 |
| 892WKG.5 | 13 | 7 | 5 |
| 904WPG.2 | 16 | 10 | 16 |
| 904WPG.3 | 22 | 29 | 4 |
| 904WPG.5 | 75 | 35 | 26 |
| WK.2 | 161 | 20 | 115 |
| WK.2E- 1 | 6 | 0 | 31 |
| WK.4 | 66 | 14 | 27 |
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.
#QQ plots against various distributions
set.seed(1)
par(mfrow=c(1,3))
normal <- fitdistr(log10(fib$mass+1), "normal")
qqp(log10(fib$mass+1), "norm", main="Normal")
lognormal <- fitdistr(fib$mass+1, "lognormal")
qqp(fib$mass+1, "lnorm", main="Log Normal")
#pois <- fitdistr(fib$mass+1, "Poisson")
#qqp(fib$mass, "pois", pois$estimate, main="Poisson")
#nbinom <- fitdistr(fib$mass+1, "Negative Binomial")
#qqp(fib$mass+1, "nbinom", size = nbinom$estimate[[1]], mu=nbinom$estimate[[2]], main="Negative Binomial")
gamma <- fitdistr(fib$mass+1, "gamma")
qqp(fib$mass+1, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]], main="Gamma")
Distributions by fixed factors
ggplot(fib, aes(x = log10(mass), fill=species)) +
geom_histogram(data=subset(fib,species == "hook"), aes(y=-..density..),binwidth=0.05)+
geom_histogram(data=subset(fib,species == "kaal"), aes(y= ..density..),binwidth=0.05)+
coord_flip() + facet_grid(~crosstype) + labs(y="Histogram", x="Log Inflo Biomass")
Distributions by random factors
ggplot(aes(y=mass, x=mompid, color=crosstype), data=fib) + geom_count(alpha=0.8) + coord_flip() + labs(x="Maternal plant", y="Mass")
Homogeneity of variances across subsets
Our mixed model uses one parameter to capture random efibect variance, which is assumed to be homogeneous. Plotting on a log scale should uncouple variances from means to assess this visually. Subsets are species * crosstype * maternal plant.
Subset variances are not homogeneous:
ggplot(aes(y=log10(mass+1), x=sxcxm, color=crosstype), data=fib) + geom_boxplot() + coord_flip() + labs(y="ln(Inflo Biomass + 1)",x="Subsets")
Subset mean-variance relationship
Various distributions make difiberent assumptions about the mean-variance (µ-Var) ratio.
grpVars <- with(fib, tapply(mass, list(sxcxm), var))
grpMeans <- with(fib, tapply(mass, list(sxcxm), mean))
grpCounts <- with(fib, tapply(mass, list(sxcxm), length))
#set weight=grpCounts to weight loess by sample sizes
ggplot(na.omit(data.frame(grpMeans,grpVars,grpCounts)),
aes(x=grpMeans,y=grpVars, weight=1))+geom_point(aes(size=grpCounts))+
guides(colour=guide_legend(title="Fit"),size=guide_legend(title="Sample size")) + labs(x="Subset Mean", y="Subset Variance") + labs(subtitle="Subset: species*crosstype*mompid")
Fixed effects
effects and interactions in these plots are simply given by the mean, which may be unduly influenced by high values.
intplot <- ggplot(fib,aes(x=crosstype,y=mass))+
geom_count(aes(size = ..prop.., group=sxc),alpha=0.5)+
stat_summary(aes(x=as.numeric(crosstype)),fun.y=mean,geom="line")+ facet_grid(~species)
intplot + aes(group=species, color=species)
Random effects
Maternal population
intplot + aes(group=mompop, color=mompop)
Maternal plant
intplot + aes(group=mompid, color=mompop)
Paternal population
intplot + aes(group=dadpop, color=dadpop)
Run models on subsets
Run many generalized linear models on subsets of the data defined by crosstype | mompid to see if effects estimates are consistent within maternal plants.
Most maternal plant subsets agree, but some are problematic outliers. These plants can be picked out visually from the random effects interaction plot above, the estimated parameters of each subset model, and the QQ plot of the estimated parameters:
#had to get rid of species or mompid since mompid is nested inside species. dadpop also works
glm.lis <- lmList(log10(mass)~crosstype|mompid,data=fib, family="gaussian")
plot.lmList(glm.lis,scale=list(x=list(relation="free"))) Loading required package: reshape
Attaching package: 'reshape' The following objects are masked from 'package:plyr':
rename, round_any The following object is masked from 'package:Matrix':
expand Using grp as id variables
qqmath.lmList(glm.lis)# Using as id variables
Models
We constructed the following models with the package glmmADMB. They all have the same fixed effects, species x crosstype, and response variable, log10(mass)
- X = standard GL(M)M
| Distribution, Random effects: | None | Maternal plant | Maternal population |
|---|---|---|---|
| normal (norm) | X | X | X |
#Normal (Gaussian) distribution, identity link
sc.norm.l <- lm(log10(mass)~species*crosstype, data=fib)
sc.mix.mompid.l <- lmer(log10(mass)~species*crosstype + (1|mompid), data=fib)
sc.mix.mompop.l <- lmer(log10(mass)~species*crosstype + (1|mompop), data=fib)
sc.mix.momdadpid.l <- lmer(log10(mass)~species*crosstype + (1|mompid) + (1|dadpid), data=fib)
sc.mix.momdadpop.l <- lmer(log10(mass)~species*crosstype + (1|mompop) + (1|dadpop), data=fib)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. Difiberences of 2 units are significant.
sc.names <- c("sc.norm.l","sc.mix.mompid.l","sc.mix.mompop.l","sc.mix.momdadpid.l","sc.mix.momdadpop.l")
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.mix.momdadpid.l | 0.000 | 9 |
| sc.norm.l | 0.337 | 7 |
| sc.mix.mompid.l | 20.139 | 8 |
| sc.mix.momdadpop.l | 37.235 | 9 |
| sc.mix.mompop.l | 39.294 | 8 |
The best-fiting model is a mixed model with the following components:
- response: log10(mass)
- fixed effects: species, crosstype, species x crosstype
- random efibect: mompid
Overdispersion
Looking at the normal, fixed effects model, we see that the residuals are not normal:
shapiro.test(sc.norm.l$residuals)#raw residuals!
Shapiro-Wilk normality test
data: sc.norm.l$residuals
W = 0.98563, p-value = 1.42e-09Coefficients
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=F)
Inference
We chose the model with nearly the best (lowest) AIC, to carry out inference tests and parameter estimation.
Description
mod <- sc.mix.momdadpid.l
print(mod) Linear mixed model fit by REML ['lmerMod']
Formula: log10(mass) ~ species * crosstype + (1 | mompid) + (1 | dadpid)
Data: fib
REML criterion at convergence: -53.1832
Random effects:
Groups Name Std.Dev.
dadpid (Intercept) 0.05054
mompid (Intercept) 0.08122
Residual 0.22873
Number of obs: 1213, groups: dadpid, 24; mompid, 22
Fixed Effects:
(Intercept) specieskaal
-1.05562 0.75427
crosstypewithin crosstypehybrid
-0.02738 0.25946
specieskaal:crosstypewithin specieskaal:crosstypehybrid
0.05251 -0.66987Test significance of random effects
Using a likelihood ratio test, with a null hypothesis of zero variance, the random efibect (maternal plant) is significant for both model parts:
anova(sc.norm.l, sc.mix.momdadpid.l) #double this p-value. or simulate null by permuting data.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 significant for the count model but not the binary model:
sxc.chisq <- drop1(mod, test="Chisq") #load from file
dfun(sxc.chisq) Single term deletions
Model:
log10(mass) ~ species * crosstype + (1 | mompid) + (1 | dadpid)
Df dAIC LRT Pr(Chi)
<none> 0.000
species:crosstype 2 39.716 43.716 3.215e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model summary
The model estimated the following parameters, with individual parameter significance determined by the Wald z-test, and fixed efibect significance determined by analysis of deviance Wald test.
summary(mod) Linear mixed model fit by REML ['lmerMod']
Formula: log10(mass) ~ species * crosstype + (1 | mompid) + (1 | dadpid)
Data: fib
REML criterion at convergence: -53.2
Scaled residuals:
Min 1Q Median 3Q Max
-4.0638 -0.6008 0.1069 0.6429 3.0984
Random effects:
Groups Name Variance Std.Dev.
dadpid (Intercept) 0.002554 0.05054
mompid (Intercept) 0.006597 0.08122
Residual 0.052318 0.22873
Number of obs: 1213, groups: dadpid, 24; mompid, 22
Fixed effects:
Estimate Std. Error t value
(Intercept) -1.05562 0.03836 -27.519
specieskaal 0.75427 0.05075 14.863
crosstypewithin -0.02738 0.02721 -1.006
crosstypehybrid 0.25946 0.03050 8.507
specieskaal:crosstypewithin 0.05251 0.04059 1.294
specieskaal:crosstypehybrid -0.66987 0.05979 -11.203
Correlation of Fixed Effects:
(Intr) spcskl crsstypw crsstyph spcskl:crsstypw
specieskaal -0.756
crsstypwthn -0.196 0.149
crsstyphybr -0.470 0.506 0.242
spcskl:crsstypw 0.136 -0.210 -0.677 -0.170
spcskl:crsstyph 0.386 -0.516 -0.124 -0.816 0.191Anova(mod, type=3) Analysis of Deviance Table (Type III Wald chisquare tests)
Response: log10(mass)
Chisq Df Pr(>Chisq)
(Intercept) 757.304 1 < 2.2e-16 ***
species 220.899 1 < 2.2e-16 ***
crosstype 82.384 2 < 2.2e-16 ***
species:crosstype 137.781 2 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Predicted random effects
These are box and QQ (to check normality) plots of the estimated random efibect 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 Inflo Biomass than either crosses between or within S. hookeri populations. The other difiberences are not significant, but remember that the fixed efibect of hybrid (vs. between) was significant (model summary).
#Count
rg <- ref.grid(mod) Loading required namespace: lmerTest#summary(rg)
sxc.lsm <- lsmeans(rg, ~ crosstype*species)
plot(sxc.lsm)
options(digits=4)
cld.mod <- cld(sxc.lsm, Letters=letters) #tukey letterings
cld.mod$response <- 10 ^ cld.mod$lsmean
cld.mod$uSE <- 10 ^ (cld.mod$lsmean+cld.mod$SE)
cld.mod$lSE <- 10 ^ (cld.mod$lsmean-cld.mod$SE)
cld.mod[rev(order(cld.mod$species, cld.mod$crosstype)),] crosstype species lsmean SE df lower.CL upper.CL .group response
hybrid kaal -0.7118 0.04169 24.92 -0.7976 -0.6259 b 0.19419
within kaal -0.2762 0.03948 20.91 -0.3583 -0.1941 c 0.52940
between kaal -0.3013 0.03324 11.78 -0.3739 -0.2288 c 0.49963
hybrid hook -0.7962 0.03608 7.17 -0.8811 -0.7113 b 0.15990
within hook -1.0830 0.04246 13.18 -1.1746 -0.9914 a 0.08260
between hook -1.0556 0.03836 8.70 -1.1429 -0.9684 a 0.08798
uSE lSE
0.21376 0.17642
0.57978 0.48340
0.53937 0.46282
0.17375 0.14715
0.09109 0.07491
0.09610 0.08054
Degrees-of-freedom method: satterthwaite
Results are given on the log10 (not the response) scale.
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 6 estimates
significance level used: alpha = 0.05H.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
with(fib, wilcox.test(fib[species=="kaal" & crosstype=="hybrid","mass"], mu=intermed))
Wilcoxon signed rank test with continuity correction
data: fib[species == "kaal" & crosstype == "hybrid", "mass"]
V = 820, p-value = 8e-05
alternative hypothesis: true location is not equal to 0.2938with(fib, wilcox.test(fib[species=="hook" & crosstype=="hybrid","mass"], mu=intermed))
Wilcoxon signed rank test with continuity correction
data: fib[species == "hook" & crosstype == "hybrid", "mass"]
V = 7900, p-value <2e-16
alternative hypothesis: true location is not equal to 0.2938round(c(H.wb,K.wb,K.H,HK.int,KH.int),2) [1] -0.06 0.06 4.68 0.17 0.26ggplot(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="First inflorescence biomass (g)",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=-1) +
scale_y_continuous(expand = expand_scale(add=c(0,0)), breaks = scales::pretty_breaks(n = 5)) +
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))