# Germination analysis # Sophie Karrenberg # 180220 # wd, packages #### setwd("~/Documents/Sophie/Projects/MANUSCRIPTS/ms in prep/Reproductive barriers Silene/Supplementary material/Dryad upload") library(lme4) library(multcomp) library(lsmeans) library(nlme) library(boot) # load & adjust germination data#### germ4 <- read.csv("germination.csv") str(germ4) #define dormancy in relatinom to OPT, as 1 - germ/OPT germ4$SD.dorm <- 1-(germ4$SD/germ4$OPT) germ4$SL.dorm <- 1-(germ4$SL/germ4$OPT) # dormancy assuming all hard seeds not germinating are dormant germ4$SD.dorm.tot <- 1-germ4$SD germ4$SL.dorm.tot <- 1-germ4$SL # setting record with dormancy <0 (more seeds germinating under treatment than after cold) sum(germ4$SL.dorm<0) # 20 sum(germ4$SD.dorm<0) # 10 germ4$SL.dorm[germ4$SL.dorm<0] <- 0 germ4$SD.dorm[germ4$SD.dorm<0] <- 0 # plot data dev.off() par(mfrow=c(2,2)) boxplot(germ4$weight*1000 ~ germ4$ct4, lab="Weight (ug)", main="Seed weight") boxplot(germ4$SD.dorm*100~ germ4$ct4, ylab="Dormancy (%)", main="SD conditions") boxplot(germ4$SL.dorm*100~ germ4$ct4,ylab="Dormancy (%)", main="SL conditions") boxplot(germ4$OPT*100~ germ4$ct4,ylab="Germination (%)", main="Dormancy breaking conds.") boxplot(germ4$weight*1000 ~ germ4$ct4, lab="Weight (ug)", main="Seed weight") boxplot(germ4$SD.dorm.tot*100~ germ4$ct4, ylab="Dormancy (%)", main="SD conditions") boxplot(germ4$SL.dorm.tot*100~ germ4$ct4,ylab="Dormancy (%)", main="SL conditions") boxplot(germ4$OPT*100~ germ4$ct4,ylab="Germination (%)", main="Dormancy breaking conds.") # Seed weight#### lme.weight <- lme(log(weight) ~ ct4, random= ~ 1|No, data=germ4) plot(resid(lme.weight), fitted(lme.weight)) qqnorm(resid(lme.weight)) qqline(resid(lme.weight)) anova(lme.weight) #Seed weight, multiple comparisons (seed.weight.171122 <- summary(lsmeans(lme.weight, pairwise~ct4, adjust="fdr", type="r"))$lsmeans) (seed.weight.multi.171122 <- summary(lsmeans(lme.weight, pairwise~ct4, adjust="fdr", type="r"))$contrasts) (cld.seed.weight.180220 <- cld(summary(glht(lme.weight, linfct=mcp(ct4="Tukey")),test=adjusted("fdr")))) # overall models lme #### #combined data frame, for SD & SL conds germ5 <- data.frame(dorm=c(germ4$SD.dorm, germ4$SL.dorm), dorm.tot=c(germ4$SD.dorm.tot, germ4$SL.dorm.tot), ct4=c(germ4$ct4, germ4$ct4), cond=c(rep("SD", times=191), rep("SL", times=191)), weight=rep(germ4$weight[1:191], times=2), No=rep(germ4$No[1:191], times=2)) str(germ5) germ5$ct4 <- factor(germ5$ct4) #IMPORTANT #DORM lme model overall model.all.dorm <- lme(dorm~ct4*cond, ~1|No, data=germ5) plot(resid(model.all.dorm), fitted(model.all.dorm)) qqnorm(resid(model.all.dorm)) qqline(resid(model.all.dorm)) anova(model.all.dorm) plot(model.all.dorm) # DORM.tot lme model overall model.all.dorm.tot <- lme(dorm.tot~ct4*cond, ~1|No, data=germ5) plot(resid(model.all.dorm.tot), fitted(model.all.dorm.tot)) qqnorm(resid(model.all.dorm.tot)) qqline(resid(model.all.dorm.tot)) anova(model.all.dorm.tot) # within SL conditions: DORM#### lme.SL <- lme(SL.dorm~ct4, ~1|No, data=germ4) plot(lme.SL) qqnorm(resid(lme.SL)) qqline(resid(lme.SL)) anova(lme.SL) summary(lme.SL) (SL.DORM.171122 <- summary(lsmeans(lme.SL, pairwise~ct4, adjust="fdr", type="r"))$lsmeans) (SL.DORM.multi.171122 <- lsmeans(lme.SL, pairwise~ct4, adjust="fdr", type="r")$contrasts) (cld.SL.DORM.180220 <- cld(summary(glht(lme.SL, linfct=mcp(ct4="Tukey")),test=adjusted("fdr")))) # within SL conditions: dorm.tot#### lme.SL.tot <- lme(SL.dorm.tot~ct4, ~1|No, data=germ4) plot(lme.SL.tot) qqnorm(resid(lme.SL.tot)) qqline(resid(lme.SL.tot)) anova(lme.SL.tot) (SL.dorm.tot.171122 <- summary(lsmeans(lme.SL.tot, pairwise~ct4, adjust="fdr", type="r")$lsmeans)) (SL.dorm.tot.multi.171122 <- lsmeans(lme.SL.tot, pairwise~ct4, adjust="fdr", type="r")$contrasts) (cld.SL.tot.180220 <- cld(summary(glht(lme.SL.tot, linfct=mcp(ct4="Tukey")),test=adjusted("fdr")))) # within SD conditions: DORM #### lme.SD <- lme(SD.dorm~ct4, ~1|No, data=germ4) plot(lme.SD) qqnorm(resid(lme.SD)) qqline(resid(lme.SD)) anova(lme.SD) (SD.DORM.171122 <- summary(lsmeans(lme.SD, pairwise~ct4, adjust="fdr", type="r"))$lsmeans) (SD.DORM.multi.171122 <- lsmeans(lme.SD, pairwise~ct4, adjust="fdr", type="r")$contrasts) (cld.SD.DORM.180220 <- cld(summary(glht(lme.SD, linfct=mcp(ct4="Tukey")),test=adjusted("fdr")))) # within SD conditions: dorm.tot #### lme.SD.tot <- lme(SD.dorm.tot~ct4, ~1|No, data=germ4) plot(lme.SD.tot) qqnorm(resid(lme.SD.tot)) qqline(resid(lme.SD.tot)) anova(lme.SD.tot) (SD.dorm.tot.171122 <- summary(lsmeans(lme.SD.tot, pairwise~ct4, adjust="fdr", type="r")$lsmeans)) (SD.dorm.tot.multi.171122 <- lsmeans(lme.SD.tot, pairwise~ct4, adjust="fdr", type="r")$contrasts) (cld.SD.tot.180220 <- cld(summary(glht(lme.SD.tot, linfct=mcp(ct4="Tukey")),test=adjusted("fdr")))) # within OPT conditions: OPT (germ)#### lme.OPT <- lme(OPT~ct4, ~1|No, data=germ4) plot(lme.OPT) qqnorm(resid(lme.OPT)) qqline(resid(lme.OPT)) anova(lme.OPT) (OPT.171122 <- summary(lsmeans(lme.OPT, pairwise~ct4, adjust="fdr", type="r"))$lsmeans) (OPT.multi.171122 <- lsmeans(lme.OPT, pairwise~ct4, adjust="fdr", type="r")$contrasts) (cld.OPT.180220 <- cld(summary(glht(lme.OPT, linfct=mcp(ct4="Tukey")),test=adjusted("fdr")))) #Fig. 3 #### dev.off() par(mfrow=c(1,3), xpd=T, mar=c(5, 5.5, 5, 4) + 0.1) # Fig 3a: conditions #### SD<-c(11,10,12, 16, 18, 17, 13, 11) SL<-c(15,13, 16.5, 22, 24, 23, 18, 15) time<- c(0, 6, 9, 13, 16, 19, 22, 24) plot(time, SD, type="l", pch=1, lty=2, col="red", xaxt="n", bty="n", las=1, ylim=c(3, 30), xlab="Time of day (h)", ylab=expression("Temperature ("~degree~"C)"), cex=1.2, cex.lab=1.3, cex.axis=1.3, yaxt="n") lines(time, SL, type="l", pch=1, lty=1, col="black") axis(1, at=c(0, 6, 12, 18, 24), las=1, cex.axis=1.3, lab=c(0, 6, 12, 18, 24)) axis(2, at=c(5, 10, 15, 20,25), las=1, cex.axis=1.3, lab=c(5, 10, 15, 20,25), ylim=c(5,25)) legend(5, 7.5, c("SD cond.", "SL cond."), cex=1.3, col=c("red", "black"), lty=c(2,1), bty="n", horiz=T) polygon(x=c(0, 0,24,24), y=c(29, 30, 30, 29)) polygon(x=c(0, 0,6,6), y=c(29, 30, 30, 29), col="black") polygon(x=c(19, 19,24,24), y=c(29, 30, 30, 29), col="black") text(x=24.5, y=29.5, labels = "Light", adj = 0, cex=1.3) mtext(side=3, line=2, at=-2,text="(a)", cex=1.5) # Fig.3b: SD dorm#### rownames(SD.DORM.171122) <- SD.DORM.171122[,1] mp <- barplot(SD.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"), "lsmean"]*100, beside = T, cex=1, names.arg= c("SD", "SL", expression("F"[1]),expression("BC"[D]),expression("BC"[L]), expression("F"[2]), expression("F"[3])), las=1, ylim=c(0,100), ylab="Dormancy (%)", xlab="Cross type", col=c("red", "white", "pink","red", "black", "pink", "pink"), density=c(NA,NA,NA,15,15,NA,NA), angle=c(0,0,0,45,45,0,0), space=c(0.2, 0.1, 0.6, 0.6, 0.1, 0.1, 0.6), cex.lab=1.3, cex.axis=1.3, cex.names=1.2) arrows(mp, SD.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 - SD.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*100, mp, SD.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 + SD.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*100, angle=90, length=0.1, code=3, lwd=1) text(x=mp, y=SD.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 +SD.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SD","F2", "F3"),"SE"]*100+5, labels=cld.SD.DORM.180220$mcletters$Letters[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3")], cex=1.3) text("SD conditions", y=90, x=7, cex=1.3) mtext(side=3, line=2, at=-0.5,text="(b)", cex=1.5) # Fig. 3c: SL dorm#### rownames(SL.DORM.171122) <- SL.DORM.171122[,1] mp <- barplot(SL.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"), "lsmean"]*100, beside = T, cex=1, names.arg= c("SD", "SL", expression("F"[1]),expression("BC"[D]),expression("BC"[L]), expression("F"[2]), expression("F"[3])), las=1, ylim=c(0,100), ylab="Dormancy (%)", xlab="Cross type", col=c("red", "white", "pink","red", "black", "pink", "pink"), density=c(NA,NA,NA,15,15,NA,NA), angle=c(0,0,0,45,45,0,0), space=c(0.2, 0.1, 0.6, 0.6, 0.1, 0.1, 0.6), cex.lab=1.3, cex.axis=1.3, cex.names=1.3) arrows(mp, SL.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 - SL.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*100, mp, SL.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 + SL.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*100, angle=90, length=0.1, code=3, lwd=1) text(x=mp, y=SL.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 +SL.DORM.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*100+5, labels=cld.SL.DORM.180220$mcletters$Letters[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3")], cex=1.3) mtext(side=3, line=2, at=-0.5,text="(c)", cex=1.5) text("SL conditions", y=90, x=7, cex=1.3) #Figure S2#### par(mfrow=c(1,3), mar=c(6, 7, 4, 4)) # Fig S2 a: climate records#### table(ave.dat$Sp) ave.dat <- read.csv("orig.climate_2.csv") str(ave.dat) ave.SD <- ave.dat[ave.dat$Sp=="SD",] ave.SL <- ave.dat[ave.dat$Sp=="SL",] boxplot(ave.SD[,5],ave.SD[,4], ave.SD[,3],ave.SL[,5],ave.SL[,4], ave.SL[,3], las=1, names=c("min", "mean", "max","min", "mean", "max"), #col=c("lightblue", "purple", "red","lightblue", "purple", "red"), col=c("red", "red", "red", "white", "white", "white"), ylab=expression("Temperature ("~degree~"C)"), ylim=c(3,30), border=T, frame=F, cex.lab=1.3, cex.axis=1.3) text(x=2, y=-3, expression(italic("Silene dioica")), cex=1.3) text(x=5, y=-3, expression(italic("S.latifolia")), cex=1.3) mtext(side=3, line=2, at=-0.5,text="(a)", cex=1.5) # Fig. S2 b: seed weight#### rownames(seed.weight.171122) <- seed.weight.171122[,1] mp <- barplot(seed.weight.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"), "response"]*1000, beside = T, cex=1, names.arg= c("SD", "SL", expression("F"[1]),expression("BC"[D]),expression("BC"[L]), expression("F"[2]), expression("F"[3])), las=1, ylim=c(0,1), ylab="Seed weight (mg)", xlab="Cross type", col=c("red", "white", "pink","red", "black", "pink", "pink"), density=c(NA,NA,NA,15,15,NA,NA), angle=c(0,0,0,45,45,0,0), space=c(0.2, 0.1, 0.6, 0.6, 0.1, 0.1, 0.6), cex.lab=1.3, cex.axis=1.3, cex.names=1.3) arrows(mp, seed.weight.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"response"]*1000 - seed.weight.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*1000, mp, seed.weight.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"response"]*1000 + seed.weight.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*1000, angle=90, length=0.1, code=3, lwd=1) text(x=mp, y=seed.weight.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"response"]*1000 +seed.weight.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*1000+0.05, labels=cld.seed.weight.180220$mcletters$Letters[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3")], cex=1.3) mtext(side=3, line=2, at=-0.5,text="(b)", cex=1.5) # Fig S2 c: OPT dorm#### rownames(OPT.171122) <- OPT.171122[,1] mp <- barplot(OPT.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"), "lsmean"]*100, beside = T, cex=1, names.arg= c("SD", "SL", expression("F"[1]),expression("BC"[D]),expression("BC"[L]), expression("F"[2]), expression("F"[3])), las=1, ylim=c(0,100), ylab="Germination (%)", xlab="Cross type", col=c("red", "white", "pink","red", "black", "pink", "pink"), density=c(NA,NA,NA,15,15,NA,NA), angle=c(0,0,0,45,45,0,0), space=c(0.2, 0.1, 0.6, 0.6, 0.1, 0.1, 0.6), cex.lab=1.3, cex.axis=1.3, cex.names=1.2) arrows(mp, OPT.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 - OPT.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*100, mp, OPT.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 + OPT.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*100, angle=90, length=0.1, code=3, lwd=1) text(x=mp, y=OPT.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 +OPT.171122[c("SD", "SL", "F1", "BC-SD","BC-SD","F2", "F3"),"SE"]*100+5, labels=cld.OPT.180220$mcletters$Letters[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3")], cex=1.3) mtext(side=3, line=2, at=-0.5,text="(c)", cex=1.5) #Figure S5#### par(mfrow=c(1,2)) # Fig. S5 a: SD dorm.tot#### rownames(SD.dorm.tot.171122) <- SD.dorm.tot.171122[,1] mp <- barplot(SD.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"), "lsmean"]*100, beside = T, cex=1, names.arg= c("SD", "SL", expression("F"[1]),expression("BC"[D]),expression("BC"[L]), expression("F"[2]), expression("F"[3])), las=1, ylim=c(0,100), ylab="Dormancy (%)", xlab="Cross type", col=c("red", "white", "pink","red", "black", "pink", "pink"), density=c(NA,NA,NA,15,15,NA,NA), angle=c(0,0,0,45,45,0,0), space=c(0.2, 0.1, 0.6, 0.6, 0.1, 0.1, 0.6), cex.lab=1.3, cex.axis=1.3, cex.names=1.2) arrows(mp, SD.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 - SD.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*100, mp, SD.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 + SD.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*100, angle=90, length=0.1, code=3, lwd=1) text(x=mp, y=SD.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 +SD.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SD","F2", "F3"),"SE"]*100+5, labels=cld.SD.tot.180220$mcletters$Letters[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3")], cex=1.3) text("SD conditions", y=100, x=7, cex=1.3) mtext(side=3, line=2, at=-0.5,text="(a)", cex=1.5) # Fig. S5 b: SL dorm tot #### rownames(SL.dorm.tot.171122) <- SL.dorm.tot.171122[,1] mp <- barplot(SL.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"), "lsmean"]*100, beside = T, cex=1, names.arg= c("SD", "SL", expression("F"[1]),expression("BC"[D]),expression("BC"[L]), expression("F"[2]), expression("F"[3])), las=1, ylim=c(0,100), ylab="Dormancy (%)", xlab="Cross type", col=c("red", "white", "pink","red", "black", "pink", "pink"), density=c(NA,NA,NA,15,15,NA,NA), angle=c(0,0,0,45,45,0,0), space=c(0.2, 0.1, 0.6, 0.6, 0.1, 0.1, 0.6), cex.lab=1.3, cex.axis=1.3, cex.names=1.3) arrows(mp, SL.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 - SL.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*100, mp, SL.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 + SL.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*100, angle=90, length=0.1, code=3, lwd=1) text(x=mp, y=SL.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"lsmean"]*100 +SL.dorm.tot.171122[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3"),"SE"]*100+5, labels=cld.SL.tot.180220$mcletters$Letters[c("SD", "SL", "F1", "BC-SD","BC-SL","F2", "F3")], cex=1.3) mtext(side=3, line=2, at=-0.5,text="(b)", cex=1.5) text("SL conditions", y=100, x=7, cex=1.3) # BOOTstrap CI of RI #### mod<- lmer(SD.dorm~ct4-1+(1|No), data=germ4) summary(lmer.SL)$coef RI.spp.SD.hab <- function(x) { summary(x)$coef -> current 1-2*current[7,1]/(current[7,1]+current[6,1])} RI.F1.SD.hab <- function(x) { summary(x)$coef -> current 1-2*current[3,1]/(current[3,1]+current[6,1])} RI.F2.SD.hab <- function(x) { summary(x)$coef -> current 1-2*current[4,1]/(current[4,1]+current[6,1])} RI.spp.SL.hab <- function(x) { summary(x)$coef -> current 1-2*current[6,1]/(current[6,1]+current[7,1])} RI.F1.SL.hab <- function(x) { summary(x)$coef -> current 1-2*current[3,1]/(current[3,1]+current[7,1])} RI.F2.SL.hab <- function(x) { summary(x)$coef -> current 1-2*current[4,1]/(current[4,1]+current[7,1])} indices.germ <- data.frame(habitat=c("SD/SL"), RI.spp.SD.EST=c(1), RI.spp.SD.lwr=c(1), RI.spp.SD.upr=c(1), RI.F1.SD.EST=c(1), RI.F1.SD.lwr=c(1), RI.F1.SD.upr=c(1), RI.F2.SD.EST=c(1), RI.spp.F2.lwr=c(1), RI.F2.SD.upr=c(1), RI.spp.SL.EST=c(1), RI.spp.SL.lwr=c(1), RI.spp.SL.upr=c(1), RI.F1.SL.EST=c(1), RI.F1.SL.lwr=c(1), RI.F1.SL.upr=c(1), RI.F2.SL.EST=c(1), RI.F2.SL.lwr=c(1), RI.F2.SL.upr=c(1)) str(indices.germ) str(germ4) library(boot) #SD conds mod <- lmer(SD.dorm~ct4-1+(1|No), data=germ4) bt <- bootMer(mod, RI.spp.SD.hab, nsim=1000, type="parametric") ci <- boot::boot.ci(bt, type="perc") ci$t0 -> indices.germ[1, 2] ci$percent[4:5] -> indices.germ[1, 3:4] bt <- bootMer(mod, RI.F1.SD.hab, nsim=1000, type="parametric") ci <- boot::boot.ci(bt, type="perc") ci$t0 -> indices.germ[1, 5] ci$percent[4:5] -> indices.germ[1, 6:7] bt <- bootMer(mod, RI.F2.SD.hab, nsim=1000, type="parametric") ci <- boot::boot.ci(bt, type="perc") ci$t0 -> indices.germ[1, 8] ci$percent[4:5] -> indices.germ[1, 9:10] #SL conds mod <- lmer(SL.dorm~ct4-1+(1|No), data=germ4) bt <- bootMer(mod, RI.spp.SL.hab, nsim=1000, type="parametric") ci <- boot.ci(bt, type="perc") ci$t0 -> indices.germ[1, 11] ci$percent[4:5] -> indices.germ[1, 12:13] bt <- bootMer(mod, RI.F1.SL.hab, nsim=1000, type="parametric") ci <- boot.ci(bt, type="perc") ci$t0 -> indices.germ[1, 14] ci$percent[4:5] -> indices.germ[1, 15:16] bt <- bootMer(mod, RI.F2.SL.hab, nsim=1000, type="parametric") ci <- boot.ci(bt, type="perc") ci$t0 -> indices.germ[1, 17] ci$percent[4:5] -> indices.germ[1, 18:19] indices.germ write.table(indices.germ, "RI.germ.csv")