Breeding seasonality and stacked reproductive investment in female vervet monkeys
Predictions
- As income breeders, vervet births will be concentrated prior to food peaks such that lactation (the 60-120 days you investigated) or weaning (180-240 days) coincides with high food availability as measured by the tree phenology data and climatic factors which may be indicative of phenology we do not record (e.g., insect abundance, crop phenology) .
- If vervets are strict income breeders, they will rely exclusively on external cues of seasonality, such that we predict a strong relationship between climatic factors (rainfall, temperature) and birth seasonality. That is, climate factors will trigger conceptions. If they are relaxed income breeders we expect a less strict relationship between food availability + climate and births. We examine whether the timing of conceptions and births is predicted by the overall availability of naturally occurring vervet plant foods, as well as by the availability of the top five most commonly consumed plant species.
- We further predict that trees, including the top five species, will show intra-annual variation in their fruiting which will be correlated with climatic conditions (i.e., rainfall and temperature).
Load Libraries & Data
## Libraries
library(tidyverse)
library(lmerTest)
library(MuMIn)
library(vegan)
library(ggthemes)
library(corrplot)
library(zoo)
library(readxl)
## Data
concept <- read.csv("data//conceptions.csv")
climate <- read.csv("data//climate.csv")
pheno <- read.csv("data//phenology.csv")
## Functions & Misc
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
source("scripts//functions.r")
source("scripts//compileData.r") ## Conception data = allvars / births data = allvars2
## Drop rows without conception or birth data
allvars <- allvars[12:100,]
allvars2 <- allvars2[17:103,]Pred 1 - Infant vs. Mother strategy
## Create a dataframe with just climate variables
allvarsClimateOnly <- allvars2 %>% select(-births, -freqBirth)
## Push date after birth forward 60 days to compare against other climate variables
motherSurvial <- allvars2 %>% select(date, freqBirth) %>% mutate(date = date + 60) %>% ## add 60 days after birth
mutate(Year = as.numeric(substring(date,1,4)), Month = as.numeric(substring(date,6,7))) %>% ## convert date back to numeric
select(-date) %>% right_join(allvarsClimateOnly) ## drop date column
### Compare the mean after birth
rollingMother<- motherSurvial %>% mutate(Max.temp.roll = rollmean(Max.temp, k=2, fill=NA, align="left"),
Min.temp.roll = rollmean(Min.temp, k=2, fill=NA, align="left"),
Rainfall.roll = rollmean(Rainfall, k=2, fill=NA, align="left"),
UF.roll = rollmean(UF, k=2, fill=NA, align="left"),
RF.roll = rollmean(RF, k=2, fill=NA, align="left"),
)
## Simple models
m1 <- glm(freqBirth ~ Max.temp.roll + Rainfall.roll , data=rollingMother %>% filter(!is.na(Max.temp.roll)), family="binomial")
car::Anova(m1, type=2)## Analysis of Deviance Table (Type II tests)
##
## Response: freqBirth
## LR Chisq Df Pr(>Chisq)
## Max.temp.roll 0.50869 1 0.4757
## Rainfall.roll 0.00977 1 0.9213m2 <- glm(freqBirth ~ UF.roll + RF.roll , data=rollingMother %>% filter(!is.na(Max.temp.roll)), family="binomial")
car::Anova(m2, type=2)## Analysis of Deviance Table (Type II tests)
##
## Response: freqBirth
## LR Chisq Df Pr(>Chisq)
## UF.roll 0.34643 1 0.5561
## RF.roll 0.24999 1 0.6171## Complex models
m3 <- glmer(freqBirth ~ Max.temp.roll + Rainfall.roll + (1|Year), data=rollingMother %>% filter(!is.na(Max.temp.roll)), family="binomial")
car::Anova(m3, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: freqBirth
## Chisq Df Pr(>Chisq)
## Max.temp.roll 0.5643 1 0.4525
## Rainfall.roll 0.4463 1 0.5041## Push date after birth forward 180 days to compare against other climate variables
infantSurvial <- allvars2 %>% select(date, freqBirth) %>% mutate(date = date + 180) %>% ## add 60 days after birth
mutate(Year = as.numeric(substring(date,1,4)), Month = as.numeric(substring(date,6,7))) %>% ## convert date back to numeric
select(-date) %>% right_join(allvarsClimateOnly) ## drop date column
### Compare the mean after birth
rollingInfant<- infantSurvial %>% mutate(Max.temp.roll = rollmean(Max.temp, k=2, fill=NA, align="left"),
Min.temp.roll = rollmean(Min.temp, k=2, fill=NA, align="left"),
Rainfall.roll = rollmean(Rainfall, k=2, fill=NA, align="left"),
UF.roll = rollmean(UF, k=2, fill=NA, align="left"),
RF.roll = rollmean(RF, k=2, fill=NA, align="left"),
)
## Simple models
m1 <- glm(freqBirth ~ Max.temp.roll + Rainfall.roll , data=rollingInfant %>% filter(!is.na(Max.temp.roll)), family="binomial")
car::Anova(m1, type=2)## Analysis of Deviance Table (Type II tests)
##
## Response: freqBirth
## LR Chisq Df Pr(>Chisq)
## Max.temp.roll 0.59111 1 0.4420
## Rainfall.roll 0.03347 1 0.8548m2 <- glm(freqBirth ~ UF.roll + RF.roll , data=rollingInfant %>% filter(!is.na(RF.roll)), family="binomial")
car::Anova(m2, type=2)## Analysis of Deviance Table (Type II tests)
##
## Response: freqBirth
## LR Chisq Df Pr(>Chisq)
## UF.roll 0.0051679 1 0.9427
## RF.roll 0.0000761 1 0.9930## Complex models
m3 <- glmer(freqBirth ~ Max.temp.roll + Rainfall.roll + (1|Year), data=rollingMother %>% filter(!is.na(Max.temp.roll)), family="binomial")
car::Anova(m3, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: freqBirth
## Chisq Df Pr(>Chisq)
## Max.temp.roll 0.5643 1 0.4525
## Rainfall.roll 0.4463 1 0.5041There is no significant relationship with tempearture, rainfall, unripe fruit, or rainfall on the number of infants present 60-120 or 180-240 days after birth.
Pred 2a - Climate & phenology effects on births/conceptions
### Correlations with conceptions
corMat <- allvars %>% select(ML, YL, FW, UF, RF, Max.temp, Mean.temp, Min.temp, Humidity, Rainfall,conceptions,freqConcept) %>% cor(., use = "pairwise.complete.obs", method="pearson")
corrplot(corMat, method="number")
## Try to predict conceptions with climate or food
m1 <- glm(freqConcept ~ UF + RF + Max.temp + Humidity + Rainfall, data= allvars, family="binomial")
car::Anova(m1, type=2)## Analysis of Deviance Table (Type II tests)
##
## Response: freqConcept
## LR Chisq Df Pr(>Chisq)
## UF 0.56156 1 0.4536
## RF 0.15821 1 0.6908
## Max.temp 0.14462 1 0.7037
## Humidity 0.00401 1 0.9495
## Rainfall 0.29423 1 0.5875### Correlations with births
corMat <- allvars2 %>% select(ML, YL, FW, UF, RF, Max.temp, Mean.temp, Min.temp, Humidity, Rainfall,births,freqBirth) %>% cor(., use = "pairwise.complete.obs", method="pearson")
corrplot(corMat, method="number")
m2 <- glm(freqBirth ~ UF + RF + Max.temp + Humidity + Rainfall, data=rollingInfant %>% filter(!is.na(Max.temp.roll)), family="binomial")
car::Anova(m1, type=2)## Analysis of Deviance Table (Type II tests)
##
## Response: freqConcept
## LR Chisq Df Pr(>Chisq)
## UF 0.56156 1 0.4536
## RF 0.15821 1 0.6908
## Max.temp 0.14462 1 0.7037
## Humidity 0.00401 1 0.9495
## Rainfall 0.29423 1 0.5875There does not appear to be a pattern of climate or phenology on conception/births.
Pred 2b - Select food phenology on births/conceptions
## Select foods
## Maesopsis eminii, Pseudospondias microcarpa, Lantana camara, Ficus natalensis, and Pycnanthus angolensis
phenoMonthly <- pheno %>% group_by(Year, Month, Species = SpeciesNameCorrected) %>% summarize_at(vars(ML:RF), mean, na.rm=T)
selectSpecies <- phenoMonthly %>% filter(Species %in% c("Pseudospondias microcarpa", "Lantana camara", "Ficus natalensis", "Pycnanthus angolensis"))
monthlySpecies <- selectSpecies %>% group_by(Year, Month) %>% summarize(UFselect = mean(UF), RFselect = mean(RF))
## merge with original data
allvarsSelect <- merge(allvars, monthlySpecies, by=c("Year","Month"))
allvars2Select <- merge(allvars2, monthlySpecies, by=c("Year","Month"))
### Correlations with conceptions
corMat <- allvarsSelect %>% select(UFselect,RFselect, Max.temp, Mean.temp, Min.temp, Humidity, Rainfall,conceptions,freqConcept) %>% cor(., use = "pairwise.complete.obs", method="pearson")
corrplot(corMat, method="number")
### Correlations with births
corMat <- allvars2Select %>% select(UFselect,RFselect, Max.temp, Mean.temp, Min.temp, Humidity, Rainfall,births,freqBirth) %>% cor(., use = "pairwise.complete.obs", method="pearson")
corrplot(corMat, method="number")
Subseting by key species there is no pattern of phenology on conceptions or births.
Pred 3a - climate effects on vegetation
allvarsNAomit <- allvars[!is.na(allvars$UF),]
m1global <- lm(UF ~ Max.temp + Mean.temp + Min.temp + Rainfall + Humidity + Sun.Hour, data=allvarsNAomit, na.action = "na.fail")
dd1 <- dredge(m1global)## Fixed term is "(Intercept)"plot(dd1)
dd1[1:5,]## Global model call: lm(formula = UF ~ Max.temp + Mean.temp + Min.temp + Rainfall +
## Humidity + Sun.Hour, data = allvarsNAomit, na.action = "na.fail")
## ---
## Model selection table
## (Int) Hmd Max.tmp Men.tmp Min.tmp Rnf Sun.Hor df logLik
## 7 0.42180 -0.1559 0.1958 4 -4.505
## 8 -0.08383 0.003248 -0.1535 0.2055 5 -4.387
## 23 0.43540 -0.1583 0.1985 -8.247e-05 5 -4.471
## 15 0.36790 -0.1492 0.1802 0.01209 5 -4.485
## 39 0.41560 -0.1560 0.1961 1.416e-05 5 -4.504
## AICc delta weight
## 7 17.5 0.00 0.428
## 8 19.6 2.05 0.154
## 23 19.8 2.21 0.141
## 15 19.8 2.24 0.140
## 39 19.8 2.28 0.137
## Models ranked by AICc(x)m2global <- lm(RF ~ Max.temp + Mean.temp + Min.temp + Rainfall + Humidity + Sun.Hour, data=allvarsNAomit, na.action = "na.fail")
dd2 <- dredge(m2global)## Fixed term is "(Intercept)"plot(dd2)
dd2[1:5,]## Global model call: lm(formula = RF ~ Max.temp + Mean.temp + Min.temp + Rainfall +
## Humidity + Sun.Hour, data = allvarsNAomit, na.action = "na.fail")
## ---
## Model selection table
## (Int) Hmd Max.tmp Men.tmp Min.tmp Rnf df logLik AICc delta
## 5 5.493 -0.2122 3 -86.161 178.6 0.00
## 13 5.496 -0.2485 0.04178 4 -86.110 180.8 2.12
## 7 5.257 0.02889 -0.2341 4 -86.127 180.8 2.15
## 6 5.898 -0.002712 -0.2218 4 -86.151 180.8 2.20
## 21 5.502 -0.2134 0.0001263 4 -86.151 180.8 2.20
## weight
## 5 0.425
## 13 0.147
## 7 0.145
## 6 0.141
## 21 0.141
## Models ranked by AICc(x)## Final models
m1 <- lm(UF ~ Max.temp + Mean.temp, data=allvarsNAomit)
summary(m1)##
## Call:
## lm(formula = UF ~ Max.temp + Mean.temp, data = allvarsNAomit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.48245 -0.24829 -0.02345 0.22058 0.55678
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.42177 0.51376 0.821 0.414241
## Max.temp -0.15587 0.04026 -3.872 0.000227 ***
## Mean.temp 0.19582 0.03578 5.472 5.48e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2612 on 76 degrees of freedom
## Multiple R-squared: 0.2975, Adjusted R-squared: 0.279
## F-statistic: 16.09 on 2 and 76 DF, p-value: 1.487e-06m2 <- lm(RF ~Mean.temp, data=allvarsNAomit)
summary(m2)##
## Call:
## lm(formula = RF ~ Mean.temp, data = allvarsNAomit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.9780 -0.6816 -0.1420 0.7049 1.3527
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.49268 1.10241 4.982 3.76e-06 ***
## Mean.temp -0.21224 0.05248 -4.044 0.000124 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7294 on 77 degrees of freedom
## Multiple R-squared: 0.1752, Adjusted R-squared: 0.1645
## F-statistic: 16.35 on 1 and 77 DF, p-value: 0.0001238### Inter-annual variation
yearlyPheno <- allvars %>% group_by(Year) %>% summarize(UFavg = mean(UF, na.rm=T), RFavg = mean(RF, na.rm=T)) %>% ## summarize annually
gather(fruitType, meanAbd, 2:3)
yearlyClim <- allvars %>% group_by(Year) %>% summarize(avgTemp = mean(Mean.temp, na.rm=T), avgPrec = mean(Rainfall, na.rm=T))
## Correlate climate and phenology
corYearly <- yearlyPheno %>% spread(fruitType, meanAbd) %>% left_join(yearlyClim)## Joining, by = "Year"cor(corYearly[-1])## RFavg UFavg avgTemp avgPrec
## RFavg 1.0000000 -0.9172010 -0.4504884 -0.2096165
## UFavg -0.9172010 1.0000000 0.5814415 0.3715938
## avgTemp -0.4504884 0.5814415 1.0000000 0.6461801
## avgPrec -0.2096165 0.3715938 0.6461801 1.0000000## Plot patterns over time
plot1 <- ggplot(yearlyPheno, aes(x=Year, y=meanAbd, color=fruitType)) + geom_line(size=1.5) + theme_minimal() + ylab("Mean Abundance")+ theme(legend.position="top") + ylim(0,3)
plot2 <- ggplot(yearlyClim, aes(x=Year, y=avgTemp)) + geom_line(size=1.5) + theme_minimal() + ylab("Mean Monthly Temperature (°C)")
plot3 <- ggplot(yearlyClim, aes(x=Year, y=avgPrec*12)) + geom_bar(stat="identity") + theme_minimal() + ylab("Yearly Precipitation (mm)")
gridExtra::grid.arrange(plot1, plot2, plot3)
### Intra-annual variation
monthlyPheno <- allvars %>% group_by(Month) %>% summarize(UFavg = mean(UF, na.rm=T), RFavg = mean(RF, na.rm=T)) %>% ## summarize monthly
gather(fruitType, meanAbd, 2:3)
monthlyClim <- allvars %>% group_by(Month) %>% summarize(avgTemp = mean(Mean.temp, na.rm=T), avgPrec = mean(Rainfall, na.rm=T))
## Correlate climate and phenology
corMonthly <- monthlyPheno %>% spread(fruitType, meanAbd) %>% left_join(monthlyClim)## Joining, by = "Month"cor(corMonthly[-1])## RFavg UFavg avgTemp avgPrec
## RFavg 1.00000000 -0.70705758 0.3030119 0.01187883
## UFavg -0.70705758 1.00000000 -0.4663973 -0.03019841
## avgTemp 0.30301190 -0.46639734 1.0000000 -0.34421070
## avgPrec 0.01187883 -0.03019841 -0.3442107 1.00000000## Plot patterns over time
plot1 <- ggplot(monthlyPheno, aes(x=Month, y=meanAbd, color=fruitType)) + geom_line(size=1.5) + theme_minimal() + ylab("Mean Abundance")+ theme(legend.position="top") + ylim(0,3)
plot2 <- ggplot(monthlyClim, aes(x=Month, y=avgTemp)) + geom_line(size=1.5) + theme_classic() + ylab("Mean Monthly Temperature (°C)")
plot3 <- ggplot(monthlyClim, aes(x=Month, y=avgPrec)) + geom_bar(stat="identity") + theme_classic() + ylab("Mean Monthly Precipitation (mm)")
gridExtra::grid.arrange(plot2, plot3)
Climate is seasonal within a year but plant phenology is not. There are differences between years in plant phenology that are determined by mean annual temperature.
Pred 3b - climate effects on vegetation for four select species
## Load data and select species
pheno <- read.csv("data//phenology.csv", stringsAsFactors = F)
phenoMonthly <- pheno %>% group_by(Year, Month, Species = SpeciesNameCorrected) %>% summarize_at(vars(ML:RF), mean, na.rm=T) %>%
filter(Species %in% c("Pseudospondias microcarpa", "Lantana camara", "Ficus natalensis", "Pycnanthus angolensis"))
## Convert month to year
phenoMonthly <- merge(phenoMonthly, data.frame(Month = month.name, Month.num=1:12))
phenoMonthly <- phenoMonthly %>% select(-Month) %>% rename(Month = Month.num) %>%
group_by(Year, Month) %>% summarize(UF = mean(UF), RF = mean(RF, na.rm=T)) ## summarize across species
## Load climate data
monthlyClim <- climate %>% group_by(Year, Month) %>% summarize_at(vars(Max.temp:Sun.Hour), mean) %>% data.frame()
newVars <- merge(phenoMonthly, monthlyClim, by=c("Year","Month"))
### Inter-annual variation
yearlyPheno <- newVars %>% group_by(Year) %>% summarize(UFavg = mean(UF, na.rm=T), RFavg = mean(RF, na.rm=T)) %>% ## summarize annually
gather(fruitType, meanAbd, 2:3)
yearlyClim <- newVars %>% group_by(Year) %>% summarize(avgTemp = mean(Mean.temp, na.rm=T), avgPrec = mean(Rainfall, na.rm=T))
## Correlate climate and phenology
corYearly <- yearlyPheno %>% spread(fruitType, meanAbd) %>% left_join(yearlyClim)## Joining, by = "Year"cor(corYearly[-1])## RFavg UFavg avgTemp avgPrec
## RFavg 1.0000000 -0.9115398 -0.4232939 -0.2174243
## UFavg -0.9115398 1.0000000 0.2182475 0.2766842
## avgTemp -0.4232939 0.2182475 1.0000000 0.3868330
## avgPrec -0.2174243 0.2766842 0.3868330 1.0000000## Plot patterns over time
plot1 <- ggplot(yearlyPheno, aes(x=Year, y=meanAbd, color=fruitType)) + geom_line(size=1.5) + theme_minimal() + ylab("Mean Abundance")+ theme(legend.position="top") + ylim(0,3)
plot2 <- ggplot(yearlyClim, aes(x=Year, y=avgTemp)) + geom_line(size=1.5) + theme_minimal() + ylab("Mean Monthly Temperature (°C)")
plot3 <- ggplot(yearlyClim, aes(x=Year, y=avgPrec*12)) + geom_bar(stat="identity") + theme_minimal() + ylab("Yearly Precipitation (mm)")
gridExtra::grid.arrange(plot1, plot2, plot3)
### Intra-annual variation
monthlyPheno <- newVars %>% group_by(Month) %>% summarize(UFavg = mean(UF, na.rm=T), RFavg = mean(RF, na.rm=T)) %>% ## summarize monthly
gather(fruitType, meanAbd, 2:3)
monthlyClim <- newVars %>% group_by(Month) %>% summarize(avgTemp = mean(Mean.temp, na.rm=T), avgPrec = mean(Rainfall, na.rm=T))
## Correlate climate and phenology
corMonthly <- monthlyPheno %>% spread(fruitType, meanAbd) %>% left_join(monthlyClim)## Joining, by = "Month"cor(corMonthly[-1])## RFavg UFavg avgTemp avgPrec
## RFavg 1.00000000 -0.2300114 0.1628764 0.02411297
## UFavg -0.23001137 1.0000000 -0.1708161 -0.39286735
## avgTemp 0.16287638 -0.1708161 1.0000000 -0.42968292
## avgPrec 0.02411297 -0.3928674 -0.4296829 1.00000000## Plot patterns over time
plot1 <- ggplot(monthlyPheno, aes(x=Month, y=meanAbd, color=fruitType)) + geom_line(size=1.5) + theme_minimal() + ylab("Mean Abundance")+ theme(legend.position="top") + ylim(0,3)
plot2 <- ggplot(monthlyClim, aes(x=Month, y=avgTemp)) + geom_line(size=1.5) + theme_minimal() + ylab("Mean Monthly Temperature (°C)")
plot3 <- ggplot(monthlyClim, aes(x=Month, y=avgPrec)) + geom_bar(stat="identity") + theme_minimal() + ylab("Mean Monthly Precipitation (mm)")
gridExtra::grid.arrange(plot1, plot2, plot3)
Rank model selection
## Load rank data
EloIBI_Mixed_Model <- read_excel("data//EloIBI_Mixed_Model_Input.xlsx", sheet = "IBI") %>% filter(!is.na(IBI))
## Model with everything and mother nested within group
fullNestedIBI <- lmer(IBI ~ RawElo + OrdinalElo + ProportionalElo + Jenks + Group + SexPrevINF + (1|Group/Mother), data = EloIBI_Mixed_Model, na.action = "na.fail")
## Model with just mother
fullMotherIBI <- lmer(IBI ~ RawElo + OrdinalElo + ProportionalElo + Jenks + Group + SexPrevINF + (1|Mother), data = EloIBI_Mixed_Model, na.action = "na.fail")
## Conduct dredge for all model combinations - NESTED
dd1IBI <- dredge(fullNestedIBI, m.lim = c(0, 3), extra = c( R2 = function(x) MuMIn::r.squaredGLMM(x)))
plot(dd1IBI)
dd1IBI[1:10,]## Global model call: lmer(formula = IBI ~ RawElo + OrdinalElo + ProportionalElo +
## Jenks + Group + SexPrevINF + (1 | Group/Mother), data = EloIBI_Mixed_Model,
## na.action = "na.fail")
## ---
## Model selection table
## (Intrc) Group Jenks OrdnE PrprE RawEl SPINF R21 R22 df logLik
## 42 464.2 + -106.70 + 0.07339 0.4213 9 -376.392
## 36 499.9 + -47.82 + 0.07484 0.4099 9 -377.241
## 38 349.9 + 9.288 + 0.08112 0.4113 9 -378.630
## 43 499.4 -28.41 -59.19 + 0.07053 0.3760 8 -381.794
## 12 446.0 + -19.13 -73.72 0.05817 0.4123 8 -382.723
## 50 704.9 + -0.2966 + 0.11040 0.4839 9 -381.425
## 34 385.3 + + 0.03853 0.3684 8 -382.887
## 45 400.8 5.621 -50.74 + 0.06994 0.3738 8 -383.441
## 14 305.8 + 10.660 10.20 0.07134 0.4099 8 -383.831
## 39 424.7 -27.92 5.537 + 0.07608 0.3752 8 -384.377
## AICc delta weight
## 42 774.1 0.00 0.634
## 36 775.8 1.70 0.271
## 38 778.6 4.48 0.068
## 43 782.2 8.09 0.011
## 12 784.1 9.95 0.004
## 50 784.2 10.07 0.004
## 34 784.4 10.28 0.004
## 45 785.5 11.38 0.002
## 14 786.3 12.16 0.001
## 39 787.4 13.25 0.001
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Group/Mother'dd1IBI[1:42,]## Global model call: lmer(formula = IBI ~ RawElo + OrdinalElo + ProportionalElo +
## Jenks + Group + SexPrevINF + (1 | Group/Mother), data = EloIBI_Mixed_Model,
## na.action = "na.fail")
## ---
## Model selection table
## (Intrc) Group Jenks OrdnE PrprE RawEl SPINF R21 R22 df
## 42 464.2 + -106.70 + 0.07339 0.4213 9
## 36 499.9 + -47.8200 + 0.07484 0.4099 9
## 38 349.9 + 9.28800 + 0.08112 0.4113 9
## 43 499.4 -28.4100 -59.19 + 0.07053 0.3760 8
## 12 446.0 + -19.1300 -73.72 0.05817 0.4123 8
## 50 704.9 + -0.2966 + 0.11040 0.4839 9
## 34 385.3 + + 0.03853 0.3684 8
## 45 400.8 5.62100 -50.74 + 0.06994 0.3738 8
## 14 305.8 + 10.66000 10.20 0.07134 0.4099 8
## 39 424.7 -27.9200 5.53700 + 0.07608 0.3752 8
## 8 344.4 + -10.5500 8.69900 0.07288 0.4130 8
## 41 476.2 -111.70 + 0.06515 0.3781 7
## 35 506.0 -50.5900 + 0.06918 0.3697 7
## 57 739.0 20.66 -0.3350 + 0.09897 0.4427 8
## 10 424.5 + -107.40 0.05748 0.4146 7
## 26 746.3 + 51.58 -0.3985 0.10460 0.4900 8
## 51 714.8 -0.1829 -0.2986 + 0.09787 0.4419 8
## 4 461.1 + -47.4200 0.05311 0.3969 7
## 20 710.2 + 17.3000 -0.3696 0.10250 0.4920 8
## 37 346.2 9.06900 + 0.06920 0.3671 7
## 15 362.2 -26.1600 8.65900 29.87 0.06412 0.3652 7
## 6 315.3 + 9.96600 0.07142 0.4242 7
## 53 714.4 0.01933 -0.2987 + 0.09785 0.4408 8
## 22 655.6 + 1.13000 -0.2853 0.09949 0.4802 8
## 11 472.2 -26.6400 -63.36 0.05135 0.3662 6
## 49 715.8 -0.2998 + 0.09941 0.4462 7
## 33 405.7 + 0.02125 0.3117 6
## 27 753.4 -3.1710 43.08 -0.3744 0.08977 0.4427 7
## 18 686.9 + -0.3105 0.10150 0.4884 7
## 2 349.4 + 0.01801 0.3611 6
## 13 340.1 8.70500 -16.73 0.05999 0.3692 6
## 29 680.3 3.69900 65.82 -0.3449 0.09118 0.4383 7
## 7 376.6 -20.1300 7.40500 0.06417 0.3701 6
## 23 683.0 6.8620 1.29400 -0.3074 0.08745 0.4401 7
## 9 450.3 -111.40 0.04756 0.3731 5
## 25 754.9 39.20 -0.3794 0.09113 0.4476 6
## 3 478.4 -50.3200 0.04737 0.3522 5
## 19 718.9 5.7870 -0.3323 0.08896 0.4481 6
## 5 323.1 9.81400 0.06063 0.3682 5
## 21 678.8 1.03800 -0.2892 0.08845 0.4407 6
## 17 709.1 -0.3122 0.08995 0.4475 5
## 1 378.9 0.00000 0.3014 4
## logLik AICc delta weight
## 42 -376.392 774.1 0.00 0.633
## 36 -377.241 775.8 1.70 0.271
## 38 -378.630 778.6 4.48 0.067
## 43 -381.794 782.2 8.09 0.011
## 12 -382.723 784.1 9.95 0.004
## 50 -381.425 784.2 10.07 0.004
## 34 -382.887 784.4 10.28 0.004
## 45 -383.441 785.5 11.38 0.002
## 14 -383.831 786.3 12.16 0.001
## 39 -384.377 787.4 13.25 0.001
## 8 -384.894 788.4 14.29 0.000
## 41 -386.776 789.6 15.43 0.000
## 35 -387.551 791.1 16.98 0.000
## 57 -386.277 791.2 17.06 0.000
## 10 -387.623 791.2 17.13 0.000
## 26 -386.627 791.9 17.75 0.000
## 51 -387.198 793.0 18.90 0.000
## 4 -388.546 793.1 18.97 0.000
## 20 -387.553 793.7 19.61 0.000
## 37 -389.221 794.4 20.32 0.000
## 15 -389.428 794.9 20.74 0.000
## 6 -389.544 795.1 20.97 0.000
## 53 -388.904 796.4 22.31 0.000
## 22 -389.318 797.3 23.14 0.000
## 11 -393.004 799.5 25.36 0.000
## 49 -391.887 799.8 25.66 0.000
## 33 -393.456 800.4 26.27 0.000
## 27 -392.274 800.5 26.43 0.000
## 18 -392.326 800.7 26.53 0.000
## 2 -394.179 801.8 27.71 0.000
## 13 -394.361 802.2 28.08 0.000
## 29 -393.857 803.7 29.60 0.000
## 7 -395.312 804.1 29.98 0.000
## 23 -395.054 806.1 31.99 0.000
## 9 -397.941 806.9 32.80 0.000
## 25 -397.077 807.6 33.51 0.000
## 3 -398.750 808.5 34.42 0.000
## 19 -398.058 809.6 35.47 0.000
## 5 -400.008 811.0 36.93 0.000
## 21 -399.742 813.0 38.84 0.000
## 17 -402.713 816.5 42.34 0.000
## 1 -404.650 818.0 43.86 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Group/Mother'## Conduct dredge for all model combinations - no group
dd2IBI <- dredge(fullMotherIBI, m.lim = c(0, 3), extra = c( R2 = function(x) MuMIn::r.squaredGLMM(x)))
plot(dd2IBI)
dd2IBI[1:10,]## Global model call: lmer(formula = IBI ~ RawElo + OrdinalElo + ProportionalElo +
## Jenks + Group + SexPrevINF + (1 | Mother), data = EloIBI_Mixed_Model,
## na.action = "na.fail")
## ---
## Model selection table
## (Intrc) Group Jenks OrdnE PrprE RawEl SPINF R21 R22 df logLik
## 42 464.2 + -106.70 + 0.07401 0.4164 8 -376.392
## 36 499.9 + -47.82 + 0.07525 0.4067 8 -377.241
## 38 349.9 + 9.288 + 0.08156 0.4079 8 -378.630
## 43 499.4 -28.41 -59.19 + 0.07053 0.3760 7 -381.794
## 12 446.0 + -19.13 -73.72 0.05864 0.4075 7 -382.723
## 50 704.9 + -0.2966 + 0.11100 0.4813 8 -381.425
## 34 385.3 + + 0.03880 0.3639 7 -382.887
## 45 400.8 5.622 -50.72 + 0.06995 0.3736 7 -383.441
## 14 305.8 + 10.660 10.20 0.07174 0.4066 7 -383.831
## 39 424.7 -27.92 5.537 + 0.07608 0.3753 7 -384.377
## AICc delta weight
## 42 771.4 0.00 0.634
## 36 773.1 1.70 0.271
## 38 775.9 4.48 0.068
## 43 779.6 8.18 0.011
## 12 781.4 10.04 0.004
## 50 781.5 10.07 0.004
## 34 781.8 10.37 0.004
## 45 782.9 11.48 0.002
## 14 783.7 12.26 0.001
## 39 784.8 13.35 0.001
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Mother'dd2IBI[1:42,]## Global model call: lmer(formula = IBI ~ RawElo + OrdinalElo + ProportionalElo +
## Jenks + Group + SexPrevINF + (1 | Mother), data = EloIBI_Mixed_Model,
## na.action = "na.fail")
## ---
## Model selection table
## (Intrc) Group Jenks OrdnE PrprE RawEl SPINF R21 R22 df
## 42 464.2 + -106.70 + 0.07401 0.4164 8
## 36 499.9 + -47.820 + 0.07525 0.4067 8
## 38 349.9 + 9.28800 + 0.08156 0.4079 8
## 43 499.4 -28.410 -59.19 + 0.07053 0.3760 7
## 12 446.0 + -19.130 -73.72 0.05864 0.4075 7
## 50 704.9 + -0.2966 + 0.11100 0.4813 8
## 34 385.3 + + 0.03880 0.3639 7
## 45 400.8 5.62200 -50.72 + 0.06995 0.3736 7
## 14 305.8 + 10.66000 10.20 0.07174 0.4066 7
## 39 424.7 -27.920 5.53700 + 0.07608 0.3753 7
## 8 344.4 + -10.550 8.69900 0.07330 0.4096 7
## 41 476.2 -111.70 + 0.06515 0.3781 6
## 57 739.0 20.66 -0.3350 + 0.09897 0.4427 7
## 35 506.0 -50.590 + 0.06918 0.3697 6
## 10 424.5 + -107.40 0.05752 0.4144 6
## 26 746.3 + 51.58 -0.3985 0.10520 0.4870 7
## 51 714.8 -0.181 -0.2986 + 0.09787 0.4420 7
## 4 461.1 + -47.420 0.05340 0.3932 6
## 20 710.2 + 17.300 -0.3696 0.10300 0.4892 7
## 37 346.2 9.06900 + 0.06920 0.3671 6
## 15 362.2 -26.160 8.65900 29.87 0.06412 0.3652 6
## 6 315.3 + 9.96600 0.07284 0.4128 6
## 53 714.4 0.02036 -0.2986 + 0.09785 0.4407 7
## 22 655.6 + 1.13000 -0.2853 0.10020 0.4767 7
## 11 472.2 -26.640 -63.36 0.05135 0.3662 5
## 49 715.7 -0.2997 + 0.09940 0.4456 6
## 33 405.7 + 0.02125 0.3116 5
## 27 753.4 -3.172 43.08 -0.3744 0.08977 0.4427 6
## 18 686.9 + -0.3105 0.10210 0.4855 6
## 2 349.4 + 0.01810 0.3575 5
## 13 340.1 8.70500 -16.73 0.05999 0.3692 5
## 29 680.4 3.69800 65.82 -0.3449 0.09119 0.4385 6
## 7 376.6 -20.130 7.40500 0.06417 0.3701 5
## 23 682.7 6.837 1.29800 -0.3071 0.08744 0.4394 6
## 9 450.3 -111.40 0.04756 0.3731 4
## 25 754.9 39.20 -0.3794 0.09113 0.4476 5
## 3 478.4 -50.320 0.04737 0.3522 4
## 19 718.8 5.768 -0.3322 0.08895 0.4477 5
## 5 323.1 9.81400 0.06063 0.3682 4
## 21 678.5 1.04300 -0.2889 0.08843 0.4400 5
## 17 709.1 -0.3122 0.08995 0.4477 4
## 1 378.9 0.00000 0.3014 3
## logLik AICc delta weight
## 42 -376.392 771.4 0.00 0.633
## 36 -377.241 773.1 1.70 0.271
## 38 -378.630 775.9 4.48 0.068
## 43 -381.794 779.6 8.18 0.011
## 12 -382.723 781.4 10.04 0.004
## 50 -381.425 781.5 10.07 0.004
## 34 -382.887 781.8 10.37 0.004
## 45 -383.441 782.9 11.48 0.002
## 14 -383.831 783.7 12.26 0.001
## 39 -384.377 784.8 13.35 0.001
## 8 -384.894 785.8 14.39 0.000
## 41 -386.776 787.0 15.62 0.000
## 57 -386.277 788.6 17.15 0.000
## 35 -387.551 788.6 17.17 0.000
## 10 -387.623 788.7 17.32 0.000
## 26 -386.627 789.3 17.85 0.000
## 51 -387.198 790.4 18.99 0.000
## 4 -388.546 790.6 19.16 0.000
## 20 -387.553 791.1 19.70 0.000
## 37 -389.221 791.9 20.51 0.000
## 15 -389.428 792.3 20.93 0.000
## 6 -389.544 792.6 21.16 0.000
## 53 -388.904 793.8 22.41 0.000
## 22 -389.318 794.6 23.23 0.000
## 11 -393.004 797.0 25.64 0.000
## 49 -391.887 797.2 25.84 0.000
## 33 -393.456 797.9 26.54 0.000
## 27 -392.274 798.0 26.62 0.000
## 18 -392.326 798.1 26.72 0.000
## 2 -394.179 799.4 27.99 0.000
## 13 -394.361 799.8 28.35 0.000
## 29 -393.857 801.2 29.78 0.000
## 7 -395.312 801.7 30.26 0.000
## 23 -395.054 803.6 32.18 0.000
## 9 -397.941 804.6 33.16 0.000
## 25 -397.077 805.2 33.79 0.000
## 3 -398.750 806.2 34.78 0.000
## 19 -398.058 807.2 35.75 0.000
## 5 -400.008 808.7 37.29 0.000
## 21 -399.742 810.5 39.12 0.000
## 17 -402.713 814.1 42.70 0.000
## 1 -404.650 815.7 44.30 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Mother'## best model
bestModelIBI <- lmer(IBI ~ ProportionalElo + Group + SexPrevINF + (1|Mother), data= EloIBI_Mixed_Model)
r.squaredGLMM(bestModelIBI)## R2m R2c
## [1,] 0.07401252 0.4163747car::Anova(bestModelIBI, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: IBI
## Chisq Df Pr(>Chisq)
## ProportionalElo 2.9712 1 0.08476 .
## Group 0.4584 2 0.79518
## SexPrevINF 1.6315 2 0.44231
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1shapiro.test(residuals(bestModelIBI))##
## Shapiro-Wilk normality test
##
## data: residuals(bestModelIBI)
## W = 0.58476, p-value = 3.563e-12## Plot model
library(ggplot2)
ggplot(EloIBI_Mixed_Model, aes(x= ProportionalElo, y=IBI)) + geom_point()
Income-I environment on weaning date (180 d)
### Compare the mean climate during jestation
rollingJest<- allvars %>% mutate(Max.temp.roll = rollmean(Max.temp, k=6, fill=NA, alight="left"),
Min.temp.roll = rollmean(Min.temp, k=6, fill=NA, alight="left"),
Rainfall.roll = rollmean(Rainfall, k=6, fill=NA, alight="left"),
UF.roll = rollmean(UF, k=6, fill=NA, alight="left"),
RF.roll = rollmean(RF, k=6, fill=NA, alight="left"),
)
m1 <- glm(freqConcept ~ Max.temp.roll +Rainfall.roll , data=rollingJest %>% filter(!is.na(Max.temp.roll)), family="binomial")
anova(m1, test="Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: freqConcept
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 83 12.550
## Max.temp.roll 1 0.168102 82 12.381 0.6818
## Rainfall.roll 1 0.082477 81 12.299 0.7740m2 <- glm(freqConcept ~ UF.roll + RF.roll , data=rollingJest %>% filter(!is.na(Max.temp.roll)), family="binomial")
anova(m2, test="Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: freqConcept
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 51 8.0629
## UF.roll 1 0.043641 50 8.0193 0.8345
## RF.roll 1 0.015977 49 8.0033 0.8994### Compare the mean after birth
rollingBirth<- allvars2 %>% mutate(Max.temp.roll = rollmean(Max.temp, k=6, fill=NA, align="left"),
Min.temp.roll = rollmean(Min.temp, k=6, fill=NA, align="left"),
Rainfall.roll = rollmean(Rainfall, k=6, fill=NA, align="left"),
UF.roll = rollmean(UF, k=6, fill=NA, align="left"),
RF.roll = rollmean(RF, k=6, fill=NA, align="left"),
)
m1 <- glm(freqBirth ~ Max.temp.roll * Rainfall.roll, data=rollingBirth %>% filter(!is.na(Max.temp.roll)), family="binomial")
car::Anova(m1, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: freqBirth
## LR Chisq Df Pr(>Chisq)
## Max.temp.roll 0.22981 1 0.6317
## Rainfall.roll 0.12870 1 0.7198
## Max.temp.roll:Rainfall.roll 0.13121 1 0.7172m2 <- glm(freqBirth ~ UF.roll + RF.roll , data=rollingBirth %>% filter(!is.na(Max.temp.roll)), family="binomial")
anova(m2, test="Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: freqBirth
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 46 9.7440
## UF.roll 1 0.038567 45 9.7055 0.8443
## RF.roll 1 0.084499 44 9.6210 0.7713Income-II environment on mid-lactation date (60 d)
## Create a dataframe with just climate variables
allvarsClimateOnly <- allvars2 %>% select(-births, -freqBirth)
## Push date after birth forward 60 days to compare against other climate variables
motherSurvial <- allvars2 %>% select(date, freqBirth) %>% mutate(date = date + 60) %>% ## add 60 days after birth
mutate(Year = as.numeric(substring(date,1,4)), Month = as.numeric(substring(date,6,7))) %>% ## convert date back to numeric
select(-date) %>% right_join(allvarsClimateOnly) ## drop date column## Joining, by = c("Year", "Month")### Compare the mean after birth
rollingMother<- motherSurvial %>% mutate(Max.temp.roll = rollmean(Max.temp, k=2, fill=NA, align="left"),
Min.temp.roll = rollmean(Min.temp, k=2, fill=NA, align="left"),
Rainfall.roll = rollmean(Rainfall, k=2, fill=NA, align="left"),
UF.roll = rollmean(UF, k=2, fill=NA, align="left"),
RF.roll = rollmean(RF, k=2, fill=NA, align="left"),
)
## Simple models
m1 <- glm(freqBirth ~ Max.temp.roll + Rainfall.roll , data=rollingMother %>% filter(!is.na(Max.temp.roll)), family="binomial")## Warning in eval(family$initialize): non-integer #successes in a binomial glm!summary(m1)##
## Call:
## glm(formula = freqBirth ~ Max.temp.roll + Rainfall.roll, family = "binomial",
## data = rollingMother %>% filter(!is.na(Max.temp.roll)))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5540 -0.4030 -0.2665 0.1202 1.5318
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.2636867 6.7629668 -1.074 0.283
## Max.temp.roll 0.2034986 0.2819022 0.722 0.470
## Rainfall.roll 0.0003921 0.0039457 0.099 0.921
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 15.957 on 84 degrees of freedom
## Residual deviance: 15.389 on 82 degrees of freedom
## (8 observations deleted due to missingness)
## AIC: 26.654
##
## Number of Fisher Scoring iterations: 5car::Anova(m1, type=2)## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Analysis of Deviance Table (Type II tests)
##
## Response: freqBirth
## LR Chisq Df Pr(>Chisq)
## Max.temp.roll 0.50869 1 0.4757
## Rainfall.roll 0.00977 1 0.9213m2 <- glm(freqBirth ~ UF.roll + RF.roll , data=rollingMother %>% filter(!is.na(Max.temp.roll)), family="binomial")## Warning in eval(family$initialize): non-integer #successes in a binomial glm!summary(m2)##
## Call:
## glm(formula = freqBirth ~ UF.roll + RF.roll, family = "binomial",
## data = rollingMother %>% filter(!is.na(Max.temp.roll)))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5799 -0.4216 -0.2863 0.2015 1.4683
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5901 2.9455 -0.200 0.841
## UF.roll -1.5655 2.6416 -0.593 0.553
## RF.roll -0.4960 0.9995 -0.496 0.620
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 12.925 on 65 degrees of freedom
## Residual deviance: 12.579 on 63 degrees of freedom
## (27 observations deleted due to missingness)
## AIC: 22.616
##
## Number of Fisher Scoring iterations: 5car::Anova(m2, type=2)## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Analysis of Deviance Table (Type II tests)
##
## Response: freqBirth
## LR Chisq Df Pr(>Chisq)
## UF.roll 0.34643 1 0.5561
## RF.roll 0.24999 1 0.6171## Complex models
m3 <- glmer(freqBirth ~ Max.temp.roll + Rainfall.roll + (1|Year), data=rollingMother %>% filter(!is.na(Max.temp.roll)), family="binomial")## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!## boundary (singular) fit: see ?isSingularsummary(m3)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: freqBirth ~ Max.temp.roll + Rainfall.roll + (1 | Year)
## Data: rollingMother %>% filter(!is.na(Max.temp.roll))
##
## AIC BIC logLik deviance df.resid
## 17.7 27.5 -4.9 9.7 81
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.2050 -0.0669 0.0871 1.9086 13.1136
##
## Random effects:
## Groups Name Variance Std.Dev.
## Year (Intercept) 0 0
## Number of obs: 85, groups: Year, 8
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -17.26365 19.19742 -0.899 0.369
## Max.temp.roll 0.58132 0.77382 0.751 0.453
## Rainfall.roll -0.01340 0.02005 -0.668 0.504
##
## Correlation of Fixed Effects:
## (Intr) Mx.tm.
## Max.tmp.rll -0.995
## Rainfll.rll 0.000 -0.082
## convergence code: 0
## boundary (singular) fit: see ?isSingularcar::Anova(m3, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: freqBirth
## Chisq Df Pr(>Chisq)
## Max.temp.roll 0.5643 1 0.4525
## Rainfall.roll 0.4463 1 0.5041ggplot(rollingMother, aes(x=RF.roll, y=freqBirth)) + geom_point()## Warning: Removed 28 rows containing missing values (geom_point).
## Push date after birth forward 180 days to compare against other climate variables
infantSurvial <- allvars2 %>% select(date, freqBirth) %>% mutate(date = date + 180) %>% ## add 60 days after birth
mutate(Year = as.numeric(substring(date,1,4)), Month = as.numeric(substring(date,6,7))) %>% ## convert date back to numeric
select(-date) %>% right_join(allvarsClimateOnly) ## drop date column## Joining, by = c("Year", "Month")### Compare the mean after birth
rollingInfant<- infantSurvial %>% mutate(Max.temp.roll = rollmean(Max.temp, k=2, fill=NA, align="left"),
Min.temp.roll = rollmean(Min.temp, k=2, fill=NA, align="left"),
Rainfall.roll = rollmean(Rainfall, k=2, fill=NA, align="left"),
UF.roll = rollmean(UF, k=2, fill=NA, align="left"),
RF.roll = rollmean(RF, k=2, fill=NA, align="left"),
)
## Simple models
m1 <- glm(freqBirth ~ Max.temp.roll + Rainfall.roll , data=rollingInfant %>% filter(!is.na(Max.temp.roll)), family="binomial")## Warning in eval(family$initialize): non-integer #successes in a binomial glm!summary(m1)##
## Call:
## glm(formula = freqBirth ~ Max.temp.roll + Rainfall.roll, family = "binomial",
## data = rollingInfant %>% filter(!is.na(Max.temp.roll)))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5390 -0.4299 -0.1570 0.1393 1.5666
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.1681145 7.4512235 0.425 0.671
## Max.temp.roll -0.2349077 0.3158779 -0.744 0.457
## Rainfall.roll 0.0008226 0.0044350 0.185 0.853
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 15.158 on 80 degrees of freedom
## Residual deviance: 14.562 on 78 degrees of freedom
## (5 observations deleted due to missingness)
## AIC: 25.36
##
## Number of Fisher Scoring iterations: 5car::Anova(m1, type=2)## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Analysis of Deviance Table (Type II tests)
##
## Response: freqBirth
## LR Chisq Df Pr(>Chisq)
## Max.temp.roll 0.59111 1 0.4420
## Rainfall.roll 0.03347 1 0.8548m2 <- glm(freqBirth ~ UF.roll + RF.roll , data=rollingInfant %>% filter(!is.na(RF.roll)), family="binomial")## Warning in eval(family$initialize): non-integer #successes in a binomial glm!summary(m2)##
## Call:
## glm(formula = freqBirth ~ UF.roll + RF.roll, family = "binomial",
## data = rollingInfant %>% filter(!is.na(RF.roll)))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4304 -0.4163 -0.3928 0.1499 1.6234
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.291818 3.415326 -0.671 0.502
## UF.roll -0.215828 2.995268 -0.072 0.943
## RF.roll 0.009644 1.105718 0.009 0.993
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 11.364 on 61 degrees of freedom
## Residual deviance: 11.341 on 59 degrees of freedom
## (5 observations deleted due to missingness)
## AIC: 21.177
##
## Number of Fisher Scoring iterations: 5car::Anova(m2, type=2)## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Analysis of Deviance Table (Type II tests)
##
## Response: freqBirth
## LR Chisq Df Pr(>Chisq)
## UF.roll 0.0051679 1 0.9427
## RF.roll 0.0000761 1 0.9930## Complex models
m3 <- glmer(freqBirth ~ Max.temp.roll + Rainfall.roll + (1|Year), data=rollingInfant %>% filter(!is.na(Max.temp.roll)), family="binomial")## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!## boundary (singular) fit: see ?isSingularsummary(m3)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: freqBirth ~ Max.temp.roll + Rainfall.roll + (1 | Year)
## Data: rollingInfant %>% filter(!is.na(Max.temp.roll))
##
## AIC BIC logLik deviance df.resid
## 16.5 26.1 -4.3 8.5 77
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.3 0.0 0.1 2.5 24765.0
##
## Random effects:
## Groups Name Variance Std.Dev.
## Year (Intercept) 0 0
## Number of obs: 81, groups: Year, 7
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 15.29916 24.76046 0.618 0.537
## Max.temp.roll -0.70980 1.04008 -0.682 0.495
## Rainfall.roll -0.04000 0.04186 -0.956 0.339
##
## Correlation of Fixed Effects:
## (Intr) Mx.tm.
## Max.tmp.rll -0.996
## Rainfll.rll -0.289 0.207
## convergence code: 0
## boundary (singular) fit: see ?isSingularcar::Anova(m3, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: freqBirth
## Chisq Df Pr(>Chisq)
## Max.temp.roll 0.4657 1 0.4950
## Rainfall.roll 0.9132 1 0.3393ggplot(rollingInfant, aes(x=Rainfall.roll, y=freqBirth)) + geom_point()## Warning: Removed 6 rows containing missing values (geom_point).