Q1.

The mallee ningaui is a small marsupial that inhabits mallee vegetation. The aim of this assessment is to build models of ningaui occurrence as a function of habitat variables. Firstly, write out a hypothesis for the relationship between ningaui occurrence and the cover of triodia. Remember to think about the direction and shape of the expected relationship. (one point)

### My hypothesis is that ningaui occurrence and the cover of triodia will have a weak positive relationship. It will be more strongly correlated initially and as the triodia cover increases there will be a weaker correlation.

Q2.

Of course, we don’t need to be restricted to just one hypothesis. Now, write out two more hypotheses: One that includes an additional habitat variable, and one that includes an interaction. (one point)

### My second hypothesis that the shrub cover will have a weak negative relationship with ningaui occurence.As shrub cover increases it is predicted that ningaui occurrences will be less likely

### My third hypothesis is that time since fire and occurrence of ningaui will have a strong positive correlation. As the years without fire increases as will the likelihood of ningaui occurrence.

Q3.

As for the example model described in this practical, we should start by doing some data exploration. Produce summary statistics and plots that you think will be useful for data exploration. Consider summaries of the response variable (ningaui_pa), the predictor variables, their relationships, and interactions. (three points)

### Shrub Exploration


model_shrub <- glm(ningaui_pa ~ shrub_cov, family=binomial(link='logit'),data=mallee)
summary(model_shrub)
## 
## Call:
## glm(formula = ningaui_pa ~ shrub_cov, family = binomial(link = "logit"), 
##     data = mallee)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4025  -1.3361   0.9679   1.0068   1.1702  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.51507    0.23102   2.230   0.0258 *
## shrub_cov   -0.02490    0.02715  -0.917   0.3590  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 189.22  on 139  degrees of freedom
## Residual deviance: 188.38  on 138  degrees of freedom
## AIC: 192.38
## 
## Number of Fisher Scoring iterations: 4
100*((189.22-183.66)/189.22)
## [1] 2.938379
new_tsfdata_shrub <- data.frame(shrub_cov = seq(from = 0, to = 30, by = 1))
new_ningauidata_shrub <- predict(model_shrub, newdata = new_tsfdata_shrub, se.fit = TRUE, type =  "response")
new_tsfdata_shrub$ningaui_pa <- new_ningauidata_shrub$fit
new_tsfdata_shrub$ningaui_pa_SE <- new_ningauidata_shrub$se.fit

pred_ningaui_shrub <- ggplot(data = new_tsfdata_shrub, aes(x = shrub_cov , y = ningaui_pa )) + 
  geom_line() + 
  labs(x ='Shrub cover (%)', y = 'Ningaui occurrence', title= 'Response of N. yvonneae to shrub cover')
pred_ningaui_shrub

pred_ningaui_shrub2 <- pred_ningaui_shrub + geom_point(data = mallee, aes(x = shrub_cov, y= ningaui_pa))
pred_ningaui_shrub3 <- pred_ningaui_shrub2 + 
  geom_line(data = new_tsfdata_shrub, aes(x = shrub_cov, y = (ningaui_pa+1.96*ningaui_pa_SE)), linetype = 'dashed') + 
  geom_line(data = new_tsfdata_shrub, aes(x = shrub_cov, y = (ningaui_pa-1.96*ningaui_pa_SE)), linetype = 'dashed') + 
  theme_bw()
pred_ningaui_shrub3 

## Time since fire data exploration

mallee_clean <- subset(mallee, time_since_fire < 150) 
plot(mallee_clean$time_since_fire, mallee_clean$ningaui_pa)

model4 <- glm(ningaui_pa ~ time_since_fire, family=binomial(link='logit'),data=mallee_clean)
summary(model4)
## 
## Call:
## glm(formula = ningaui_pa ~ time_since_fire, family = binomial(link = "logit"), 
##     data = mallee_clean)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8749  -1.1610   0.9208   1.0057   1.1604  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)     -0.086009   0.303132  -0.284   0.7766  
## time_since_fire  0.015757   0.008385   1.879   0.0602 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 187.42  on 138  degrees of freedom
## Residual deviance: 183.66  on 137  degrees of freedom
## AIC: 187.66
## 
## Number of Fisher Scoring iterations: 4
100*((187.42-183.66)/187.42)
## [1] 2.006189
##2.00
new_tsfdata2 <- data.frame(time_since_fire = seq(from = 0, to = 100, by = 1))
new_ningauibdata <- predict(model4, newdata = new_tsfdata2, se.fit = TRUE, type =  "response")
new_tsfdata2$ningaui_pa <- new_ningauibdata$fit
new_tsfdata2$ningaui_pa_SE <- new_ningauibdata$se.fit
pred_ningauib <- ggplot(data = new_tsfdata2, aes(x = time_since_fire, y = ningaui_pa)) + 
  geom_line() + 
  labs(x ='Time since fire (years)', y = 'Ningaui occurrence', title= 'Response of N. yvonneae to time since fire')
pred_ningauib

pred_ningauib2 <- pred_ningauib + geom_point(data = mallee_clean, aes(x = time_since_fire, y= ningaui_pa))
pred_ningauib3 <- pred_ningauib2 + 
  geom_line(data = new_tsfdata2, aes(x = time_since_fire, y = (ningaui_pa+1.96*ningaui_pa_SE)), linetype = 'dashed') + 
  geom_line(data = new_tsfdata2, aes(x = time_since_fire, y = (ningaui_pa-1.96*ningaui_pa_SE)), linetype = 'dashed') + 
  theme_bw()
pred_ningauib3

ggplot(data=mallee, aes(x=shrub_cov, y=triodia_cov)) + 
  geom_point() +
  xlab("shrub cover (%)") +
  ylab("Triodia cover (%)")

cor.test(mallee$shrub_cov, mallee$triodia_cov, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  mallee$shrub_cov and mallee$triodia_cov
## t = -2.6957, df = 138, p-value = 0.007898
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.37562583 -0.05998298
## sample estimates:
##        cor 
## -0.2236605
model_veg <- glm(triodia_cov ~ shrub_cov, family = gaussian, data = mallee)
summary(model_veg)
## 
## Call:
## glm(formula = triodia_cov ~ shrub_cov, family = gaussian, data = mallee)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -15.548  -10.953   -2.034   10.452   33.805  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  15.5485     1.4269  10.897   <2e-16 ***
## shrub_cov    -0.4596     0.1705  -2.696   0.0079 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 161.3305)
## 
##     Null deviance: 23436  on 139  degrees of freedom
## Residual deviance: 22264  on 138  degrees of freedom
## AIC: 1113
## 
## Number of Fisher Scoring iterations: 2
100*((23436-22264 )/23436)
## [1] 5.000853

Q4.

Build a logistic regression model of ningaui occurrence as a function of triodia cover. Summarise it. (one point)

model_ninguai_triodia <- glm(ningaui_pa ~ triodia_cov, family=binomial(link='logit'),data=mallee)
summary(model_ninguai_triodia)
## 
## Call:
## glm(formula = ningaui_pa ~ triodia_cov, family = binomial(link = "logit"), 
##     data = mallee)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5539  -0.7724   0.2111   0.6316   1.6462  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.05669    0.27955  -3.780 0.000157 ***
## triodia_cov  0.14263    0.02459   5.801 6.61e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 189.22  on 139  degrees of freedom
## Residual deviance: 129.01  on 138  degrees of freedom
## AIC: 133.01
## 
## Number of Fisher Scoring iterations: 5
new_tsfdata <- data.frame(triodia_cov = seq(from = 0, to = 100, by = 1))
new_ningauidata <- predict(model_ninguai_triodia, newdata = new_tsfdata, se.fit = TRUE, type =  "response")
new_tsfdata$ningaui_pa <- new_ningauidata$fit
new_tsfdata$ningaui_pa_SE <- new_ningauidata$se.fit

pred_ningaui <- ggplot(data = new_tsfdata, aes(x = triodia_cov , y = ningaui_pa )) + 
  geom_line() + 
  labs(x ='Triodia cover (%)', y = 'Ningaui occurrence', title= 'Response of N. yvonneae to time since fire')
pred_ningaui

pred_ningaui2 <- pred_ningaui + geom_point(data = mallee, aes(x = triodia_cov, y= ningaui_pa))
pred_ningaui3 <- pred_ningaui2 + 
  geom_line(data = new_tsfdata, aes(x = triodia_cov, y = (ningaui_pa+1.96*ningaui_pa_SE)), linetype = 'dashed') + 
  geom_line(data = new_tsfdata, aes(x = triodia_cov, y = (ningaui_pa-1.96*ningaui_pa_SE)), linetype = 'dashed') + 
  theme_bw()
pred_ningaui3

100*((189.22-129.01)/189.22)
## [1] 31.8201
### Using this modelit indicates that 31.82% of variation in ningaui occurrence can be explained by triodia cover. This means that as triodia cover increases it is likely that ningaui occurrence will become more probable. This model has an AIC value of 133.01 which compared to the AIC values for time since fire (187.66) and shrub cover (192.38) suggests that this model is more parsimonious and a better fit

Q5.

Now build on this model. First, by adding two more predictor variables (- and then summarise that model (but don’t use time since fire as a variable, stick to habitat measures). Second, by adding an interaction term. Summarise both of those models. Which model is the best ‘fit’ to the data? Use the deviance explained and AIC to inform your answer. Include a brief description of how you are interpreting the model outputs in the context of ningaui ecology and mallee ecosystems.(four points)

### 

model_triodia_shrub_litter <- glm(ningaui_pa ~ triodia_cov + shrub_cov + litter_cov, family=binomial(link='logit'),data=mallee)
summary(model_triodia_shrub_litter)
## 
## Call:
## glm(formula = ningaui_pa ~ triodia_cov + shrub_cov + litter_cov, 
##     family = binomial(link = "logit"), data = mallee)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4722  -0.7590   0.2203   0.5746   1.8034  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.98629    0.57112  -3.478 0.000505 ***
## triodia_cov  0.14750    0.02600   5.672 1.41e-08 ***
## shrub_cov    0.03661    0.03442   1.064 0.287415    
## litter_cov   0.03161    0.01793   1.763 0.077853 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 189.22  on 139  degrees of freedom
## Residual deviance: 124.82  on 136  degrees of freedom
## AIC: 132.82
## 
## Number of Fisher Scoring iterations: 5
100*((189.22-124.82)/189.22)
## [1] 34.03446
##34.03446
## AIC 133.82


model_triodia_shrub_litter_reserve <- glm(ningaui_pa ~ triodia_cov + shrub_cov + litter_cov + reserve, family=binomial(link='logit'),data=mallee)
summary(model_triodia_shrub_litter_reserve)
## 
## Call:
## glm(formula = ningaui_pa ~ triodia_cov + shrub_cov + litter_cov + 
##     reserve, family = binomial(link = "logit"), data = mallee)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4616  -0.7713   0.2303   0.5835   1.7935  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -1.67832    0.63234  -2.654  0.00795 ** 
## triodia_cov          0.14013    0.02577   5.437 5.41e-08 ***
## shrub_cov            0.01869    0.03728   0.501  0.61614    
## litter_cov           0.03982    0.01980   2.011  0.04429 *  
## reserveeasternblock -0.37630    0.53218  -0.707  0.47950    
## reservewesternblock -0.82874    0.67072  -1.236  0.21661    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 189.22  on 139  degrees of freedom
## Residual deviance: 123.24  on 134  degrees of freedom
## AIC: 135.24
## 
## Number of Fisher Scoring iterations: 5
100*((189.22-123.24)/189.22)
## [1] 34.86946
##34.58937
## AIC 135.24
### The model using triodia cover, litter cover and shrub cover explains 34.03% of the variation in ningaui occurrence. The model using triodia cover, litter cover, shrub cover and reserve explains 34.87% of variation in ningaui occurrence. Of all the models explored in this prac this model using 4 variables explains the most variation. While it explains more of the variation it is not as good as a fit to to the data indicated by its AIC value of 135.24 being higher than the triodia cover, litter cover and shrub cover model which has a AIC value of 132.82. The triodia cover, litter cover and shrub cover model has a lower AIC value as it is a simpler model and therefore more parsimonious. 

## Both these models indicate that triodia cover is the most significant coeffecient for explaining ningaui occurrence with z values of 5.672 and 5.437 respectively. This would suggest that triodia is an important habiat requirement for ningaui. Shrub cover however is not a significant coefficent in influencing the occurrence of ningaui. This means they are mostly just as likely to occur in areas with shrubs as they are in areas without shrubs. Litter cover, while not as significant as triodia cover, influences the occurrence of ningaui occurrence. This is likely due to litter cover providing habitat and food for invertebrates that the ningaui feeds on.