Preliminary Results Q1 and Q3
Lara Wootton
2/28/2017
Things to bear in mind
- I haven’t checked whether model assumptions have been violated
- The model fits aren’t always amazing
- Some of the veg variables were causing an error and came up as NA’s in the model output, so I generally just excluded them.
Question 1: What is the species’ preferred habitat?
Part A, Presence/Absence
Presence/Absence of the two species
Modeled using logistical regression. The least significant variables were deleted one by one.
There were two equally good models according to AIC. Both included max and min temperature and NDVI, while one also contained dem30m. The species’ seemed to prefer areas with warm winters and cool summers.
hab6<-glm(Absent ~
dem30m
+ Tmax_jan_mean30m
+ Tmin_jul_mean30m
+ Max_NDVI_2015_Landsat8,
family = binomial, data = count)
summary(hab6)##
## Call:
## glm(formula = Absent ~ dem30m + Tmax_jan_mean30m + Tmin_jul_mean30m +
## Max_NDVI_2015_Landsat8, family = binomial, data = count)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9777 -1.1012 0.5996 1.0321 2.0733
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.691678 5.534581 1.209 0.2266
## dem30m -0.003119 0.002545 -1.225 0.2204
## Tmax_jan_mean30m 0.318583 0.146307 2.178 0.0294 *
## Tmin_jul_mean30m -1.460109 0.360713 -4.048 5.17e-05 ***
## Max_NDVI_2015_Landsat8 -3.867738 1.930400 -2.004 0.0451 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 396.33 on 286 degrees of freedom
## Residual deviance: 359.96 on 282 degrees of freedom
## AIC: 369.96
##
## Number of Fisher Scoring iterations: 4hab7<-glm(count$Absent ~
count$Tmax_jan_mean30m
+ count$Tmin_jul_mean30m
+ count$Max_NDVI_2015_Landsat8,
family = binomial)
summary(hab7)##
## Call:
## glm(formula = count$Absent ~ count$Tmax_jan_mean30m + count$Tmin_jul_mean30m +
## count$Max_NDVI_2015_Landsat8, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.958 -1.100 0.573 1.048 2.017
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7614 2.6241 0.290 0.77169
## count$Tmax_jan_mean30m 0.3834 0.1354 2.832 0.00463 **
## count$Tmin_jul_mean30m -1.1067 0.2085 -5.308 1.11e-07 ***
## count$Max_NDVI_2015_Landsat8 -3.5876 1.9142 -1.874 0.06090 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 396.33 on 286 degrees of freedom
## Residual deviance: 361.46 on 283 degrees of freedom
## AIC: 369.46
##
## Number of Fisher Scoring iterations: 4The two species appear to like cool summers, warm winters and productive areas.
library(visreg)
par(mfrow=c(2,2))
visreg(hab6, "Tmax_jan_mean30m", scale = "response", ylab = "Pr(Absent)")
visreg(hab6, "Tmin_jul_mean30m", scale = "response", ylab = "Pr(Absent)")
visreg(hab6, "Max_NDVI_2015_Landsat8", scale = "response", ylab = "Pr(Absent)")
visreg(hab6, "dem30m", scale = "response", ylab = "Pr(Absent)")
Presence/absence of Mimetes
Max and min temp and solar radiation were the most important variables to predict presence/absence of Mimetes. There was also some indication that aspect and NDVI may play a role. However, there wasn’t that much difference between the AIC of the best 4 models, and overall the models didn’t explain much of the variation.
mim7<-glm(presence_M ~
jan_solar_radiation30m
+ jul_solar_radiation30m
+ Tmax_jan_mean30m
+ Tmin_jul_mean30m,
family = binomial, data = count)
summary(mim7)##
## Call:
## glm(formula = presence_M ~ jan_solar_radiation30m + jul_solar_radiation30m +
## Tmax_jan_mean30m + Tmin_jul_mean30m, family = binomial, data = count)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3041 -0.8101 -0.6451 1.2722 2.1521
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.3655306 5.0497357 0.072 0.9423
## jan_solar_radiation30m 0.0008969 0.0005198 1.726 0.0844 .
## jul_solar_radiation30m -0.0005199 0.0002262 -2.298 0.0215 *
## Tmax_jan_mean30m -0.4552333 0.1434371 -3.174 0.0015 **
## Tmin_jul_mean30m 0.4968661 0.2286484 2.173 0.0298 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 327.63 on 286 degrees of freedom
## Residual deviance: 312.50 on 282 degrees of freedom
## AIC: 322.5
##
## Number of Fisher Scoring iterations: 4What is weird though, is that I would’ve expected the same trend between solar rad and temperature. Eg, high july temps and high july solar rad would increase the probability of M presence. Instead, there’s an inverse relationship between temp and solar rad. This kind of makes sense for summer, where one could imagine that the plants like cooler temps, but still sunny conditions. But why would they prefer less solar rad in winter? #
par(mfrow=c(2,2))
visreg(mim7, "jan_solar_radiation30m", scale = "response", ylab = "Pr(M present)")
visreg(mim7, "jul_solar_radiation30m", scale = "response", ylab = "Pr(M present)")
visreg(mim7, "Tmax_jan_mean30m", scale = "response", ylab = "Pr(M present)")
visreg(mim7, "Tmin_jul_mean30m", scale = "response", ylab = "Pr(M present)")
Presence/absence Leucospermum
The leucs had a different set of variables affecting presence. The most influential variable appeared to be temp in July. Interestingly though, temp in Jan and solar rad had no effect. The other important variables were MaxEnt and NDVI, with the 2nd best model also containing dem30m. These models also seemed to fit quite well, explaining ~20% of the variation, as opposed to ~5% of the Mimetes models.
leu7<-glm(presence_L ~
Tmin_jul_mean30m
+ Max_NDVI_2015_Landsat8
+ A.lig_final_MaxEnt_FrBecker,
family = binomial, data = count)
summary(leu7)##
## Call:
## glm(formula = presence_L ~ Tmin_jul_mean30m + Max_NDVI_2015_Landsat8 +
## A.lig_final_MaxEnt_FrBecker, family = binomial, data = count)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3832 -0.7109 -0.3941 0.1460 2.5262
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.0837 2.2244 -4.983 6.27e-07 ***
## Tmin_jul_mean30m 1.1423 0.2322 4.920 8.67e-07 ***
## Max_NDVI_2015_Landsat8 3.7389 2.1423 1.745 0.080937 .
## A.lig_final_MaxEnt_FrBecker -4.7288 1.2529 -3.774 0.000161 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 323.33 on 286 degrees of freedom
## Residual deviance: 254.28 on 283 degrees of freedom
## AIC: 262.28
##
## Number of Fisher Scoring iterations: 5The probability of finding Leucs sharply increases with warmer winter temps. #
par(mfrow=c(2,2))
visreg(leu7, "Tmin_jul_mean30m", scale = "response", ylab="Pr(L presence)")
visreg(leu7, "Max_NDVI_2015_Landsat8", scale = "response", ylab="Pr(L presence)")
visreg(leu7, "A.lig_final_MaxEnt_FrBecker", scale = "response", ylab="Pr(L presence)")
Question 1: What is the species’ preferred habitat?
Part B, Abundance
Note: The data were extremely overdispersed, so I fitted quasipoisson models. However, that means AICs no longer apply, so it’s a bit harder to determine which are the best models. I’m also not sure whether the quasipoisson has totally corrected the overdispersion problem, as I don’t know how to check that.
Total Abundance
There’s a different set of variables influencing abundance compared to presence/absence. Indexgrid_landsat_30m, NDVI and TPI were the most significant. Unlike the presence/absence models, temp wasn’t influential.
ta6<-glm(Total_count ~
indexgrid_landsat_30m
+ topographicPositionIndex30m
+ veg_national30m
+ Max_NDVI_2015_Landsat8
+ A.lig_final_MaxEnt_FrBecker,
family = quasipoisson, data = count)
summary(ta6)##
## Call:
## glm(formula = Total_count ~ indexgrid_landsat_30m + topographicPositionIndex30m +
## veg_national30m + Max_NDVI_2015_Landsat8 + A.lig_final_MaxEnt_FrBecker,
## family = quasipoisson, data = count)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -8.0754 -3.9024 -2.7139 -0.9241 20.7294
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.848e+00 2.534e+00 -1.519 0.12995
## indexgrid_landsat_30m 1.095e-05 4.159e-06 2.634 0.00891 **
## topographicPositionIndex30m -1.034e-02 5.446e-03 -1.899 0.05853 .
## veg_national30m19 5.863e-01 8.344e-01 0.703 0.48281
## Max_NDVI_2015_Landsat8 2.488e+00 1.710e+00 1.455 0.14674
## A.lig_final_MaxEnt_FrBecker -2.705e+00 9.678e-01 -2.795 0.00554 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 43.9719)
##
## Null deviance: 7891.5 on 286 degrees of freedom
## Residual deviance: 6960.3 on 281 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 7Mimetes Abundance
TPI and indexgrid_landsat_30m were the most significant. None of the temp/solar rad variables were important.
ma6<-glm(No_M ~
Elevation_m
+ indexgrid_landsat_30m
+ slope30m
+ topographicPositionIndex30m,
family = quasipoisson, data = count)
summary(ma6)##
## Call:
## glm(formula = No_M ~ Elevation_m + indexgrid_landsat_30m + slope30m +
## topographicPositionIndex30m, family = quasipoisson, data = count)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6.6399 -3.0624 -1.9613 -0.9584 23.0836
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.208e+01 4.004e+00 -3.017 0.002786 **
## Elevation_m 4.453e-03 2.579e-03 1.727 0.085295 .
## indexgrid_landsat_30m 2.348e-05 6.510e-06 3.607 0.000366 ***
## slope30m -1.743e+00 1.501e+00 -1.161 0.246778
## topographicPositionIndex30m -2.460e-02 9.040e-03 -2.722 0.006899 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 40.41973)
##
## Null deviance: 5857.6 on 286 degrees of freedom
## Residual deviance: 4805.8 on 282 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 7Not sure if visreg works for poisson models - what’s with the different yaxis scales?
par(mfrow=c(2,2))
visreg(ma6, "Elevation_m", scale = "response", ylab = "M abundance")
visreg(ma6, "indexgrid_landsat_30m", scale = "response", ylab = "M abundance")
visreg(ma6, "slope30m", scale = "response", ylab = "M abundance")
visreg(ma6, "topographicPositionIndex30m", scale = "response", ylab = "M abundance")
Leucospermum Abundance
NDVI and MaxEnt were significant. This is similar to the presence/absence, except July temp is no longer important
la7<-glm(No_L ~
aspect30m
+ veg_national30m
+ Max_NDVI_2015_Landsat8
+ A.lig_final_MaxEnt_FrBecker,
family = quasipoisson, data = count)
summary(la7)##
## Call:
## glm(formula = No_L ~ aspect30m + veg_national30m + Max_NDVI_2015_Landsat8 +
## A.lig_final_MaxEnt_FrBecker, family = quasipoisson, data = count)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -8.8516 -2.6854 -1.6454 -0.9901 21.0720
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.88811 0.84862 1.047 0.29621
## aspect30m -0.11058 0.08254 -1.340 0.18140
## veg_national30m19 0.90278 0.63038 1.432 0.15321
## Max_NDVI_2015_Landsat8 5.11252 1.68141 3.041 0.00258 **
## A.lig_final_MaxEnt_FrBecker -7.17619 1.34690 -5.328 2.03e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 25.5887)
##
## Null deviance: 5361.1 on 286 degrees of freedom
## Residual deviance: 4000.0 on 282 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 7par(mfrow = c(1,2))
visreg(la7, "Max_NDVI_2015_Landsat8", scale = "response", ylab = "abundance")
visreg(la7, "A.lig_final_MaxEnt_FrBecker", scale = "response", ylab = "abundance")
Question 1 Summary
- Generally, abundance and presence/absence were affected by different sets of variables, although there was some overlap in the case of the Leucs
- Winter temperatures seemed to have greater effect than summer temperatures.
- Mimetes seems to be more sensitive to temperature and solar radiation than the Leucospermum
- Leucs prefer areas of high NDVI
- Elevation and slope did not have as big an effect as expected from field observations.
- The two species appear to occupy different niches as different sets of variables predicted their presence/absence.
Question 3: Can we predict post-fire mortality and is it higher in “hotter/drier” topoclimates?
The vegetation variables are causing the glm to give this error: glm.fit: fitted rates numerically 0 occurred and the model summaries have NA’s. I googled around, but I can’t find much info on the error.
For the models below, I didn’t differentiate between recently dead and dying, or between species. Ideally, one would eventually differentiate between species. However, I suspect that we don’t have enough observations of each species to fit good models. Plus it’s going to require lots of fiddly subsetting, so I haven’t got around to it yet.
Predicting post-fire mortality
Poisson models. By successively deleting the least significant variables, I chose the best fitting model according to AIC.
mod7<-glm(totR ~
Elevation_m
+ indexgrid_landsat_30m
+ slope30m
+ topographicPositionIndex30m
+ veg_national30m
+ Max_NDVI_2015_Landsat8,
family = poisson, data = count)
summary(mod7)##
## Call:
## glm(formula = totR ~ Elevation_m + indexgrid_landsat_30m + slope30m +
## topographicPositionIndex30m + veg_national30m + Max_NDVI_2015_Landsat8,
## family = poisson, data = count)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4524 -0.5359 -0.3868 -0.2536 4.5686
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.826e+01 3.943e+00 -4.631 3.64e-06 ***
## Elevation_m 3.888e-03 2.117e-03 1.836 0.06629 .
## indexgrid_landsat_30m 2.229e-05 5.556e-06 4.011 6.04e-05 ***
## slope30m 2.609e+00 1.065e+00 2.450 0.01429 *
## topographicPositionIndex30m -2.173e-02 8.142e-03 -2.669 0.00762 **
## veg_national30m19 2.303e+00 9.032e-01 2.550 0.01078 *
## Max_NDVI_2015_Landsat8 5.620e+00 2.168e+00 2.592 0.00953 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 196.73 on 286 degrees of freedom
## Residual deviance: 164.83 on 280 degrees of freedom
## AIC: 238.64
##
## Number of Fisher Scoring iterations: 6Recently deads are more commonly found on steep slopes with low TPI. Interestingly, they are also more likely to be found in areas of higher NDVI. None of the hot/dry variables made it into the final model.
par(mfrow=c(3,2))
visreg(mod7, "Elevation_m", scale = "response")
visreg(mod7, "indexgrid_landsat_30m", scale = "response")
visreg(mod7, "slope30m", scale = "response")
visreg(mod7, "topographicPositionIndex30m", scale = "response")
visreg(mod7, "veg_national30m", scale = "response")
visreg(mod7, "Max_NDVI_2015_Landsat8", scale = "response")
But what if we make a model with just hot/dry variables?
January and July temperatures are significant. However, the AIC for this model is worse than that of the model above (253.8 vs 238.6)
mod10<- glm(totR ~
jan_solar_radiation30m
+ jul_solar_radiation30m
+ Tmax_jan_mean30m
+ Tmin_jul_mean30m,
family = poisson, data = count)
summary(mod10)##
## Call:
## glm(formula = totR ~ jan_solar_radiation30m + jul_solar_radiation30m +
## Tmax_jan_mean30m + Tmin_jul_mean30m, family = poisson, data = count)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0064 -0.5618 -0.4627 -0.3547 5.2269
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.3462191 5.1753903 -0.840 0.40103
## jan_solar_radiation30m 0.0005244 0.0005201 1.008 0.31335
## jul_solar_radiation30m -0.0003627 0.0002364 -1.534 0.12500
## Tmax_jan_mean30m -0.3159015 0.1572617 -2.009 0.04456 *
## Tmin_jul_mean30m 0.8475230 0.3202986 2.646 0.00814 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 196.73 on 286 degrees of freedom
## Residual deviance: 184.02 on 282 degrees of freedom
## AIC: 253.84
##
## Number of Fisher Scoring iterations: 6But now this doesn’t make sense at all. It’s saying that recently deads are more likely to be found in cooler summers and warmer winters, the opposite of our hypothesis. Wait though. That’s the same result we got for presence/absence of Mimetes. Ie this model is saying that where Mimetes are present, you’re more likely to get recently deads, which is actually kinda obvious.
par(mfrow = c(2,2))
visreg(mod10, "jan_solar_radiation30m", scale = "response")
visreg(mod10, "jul_solar_radiation30m", scale = "response")
visreg(mod10, "Tmax_jan_mean30m", scale = "response")
visreg(mod10, "Tmin_jul_mean30m", scale = "response")
Let’s test this with a model with number of M and L as explanatory variables.
model<-glm(totR ~
No_M
+ No_L,
family = poisson, data = count)
summary(model)##
## Call:
## glm(formula = totR ~ No_M + No_L, family = poisson, data = count)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0175 -0.4519 -0.4384 -0.4384 4.1852
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.342543 0.191612 -12.225 < 2e-16 ***
## No_M 0.017648 0.003698 4.772 1.83e-06 ***
## No_L 0.020259 0.004256 4.760 1.94e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 196.73 on 286 degrees of freedom
## Residual deviance: 172.66 on 284 degrees of freedom
## AIC: 238.47
##
## Number of Fisher Scoring iterations: 7Yes, both number of M and L are significant. This suggests that to tease apart what is really causing the deaths, we need to account for abundance.