## **************************************************
## Elephant population size estimation
## **************************************************
rm(list=ls())
gc()

## 1.
vals <- 59:200
likelihoods <- dhyper(x=15, m=27, n=vals, k=74, log=TRUE)
vals[which.max(likelihoods)]
## so, mle number of unmarked individuals is 106, bringing the
## estimate of the total population size is:
106+27

## 2.
chi.2 <- qchisq(0.05, df=1, lower.tail=FALSE)/2
vals.above <- vals[max(likelihoods)-likelihoods<chi.2]
## and thus, the confidence interval for unmarked individuals is
range(vals.above)
## and for the total number is
range(vals.above)+27

## note - what we do below contains some repetition from what we have
## done above - but note that we are not working in log-likelihood
## space here

## 3.
vals <- 59:300

## 4.
ll <- dhyper(x=15, m=27, n=vals, k=74, log=FALSE)

## 5.
vals.prior <- rep(1/length(vals), length(vals))
plot(x=vals, y=vals.prior, xlab='Value', ylab='Prior',
     type='l', las=1)

## 6.
posteriors.non.normalized <- ll*vals.prior
posteriors <- posteriors.non.normalized/sum(posteriors.non.normalized)

plot(x=vals, y=posteriors, xlab='Value', ylab='Posterior',
     type='l', las=1)

## 7.
vals[which.max(posteriors)]
## answer: 106 + 27 = 133

## 8.
##
## First, order the posterior probabilities from highest to lowest
post.ordered <- posteriors[order(posteriors, decreasing=TRUE)]
## Remember to order the corresponding n values the same way
vals.ordered <- vals[order(posteriors, decreasing=TRUE)]
## Obtain the cumulative sum of the posterior probabilities from
## lowest to highest
post.cumsum <- cumsum(post.ordered)
## Figure out which values in the cumsum are needed to exceed 0.95
keep <- 1:which(post.cumsum>0.95)[1]
## Check the sum
sum(post.ordered[keep])
##
## Finally, find n corresponding to a cumulative posterior probability
## of 0.95. 
bci <- range(vals.ordered[keep])
## bci for total:
bci + 27

## could also do this using the uniroot function
f <- function(x) sum((posteriors)[(posteriors-x)>0])-0.95
x <- uniroot(f, interval=c(0,1))$root
## x is the height of the line that interects our posterior at the 95%
## BCI
range(vals[posteriors>x]) + 27

## 9.
##
## Above we calculated the likelihood-based confidence interval and it
## was (105,192).  Here we have calculated the 95% BCI to be
##(103,199).  These are pretty close.

## 10.
vals <- 50:1000
ll <- dhyper(x=15, m=27, n=vals, k=74, log=FALSE)
vals.prior <- rep(1/length(vals), length(vals))
posteriors <- sapply(seq_along(vals), calc.post)
post.ordered <- posteriors[order(posteriors, decreasing=TRUE)]
vals.ordered <- vals[order(posteriors, decreasing=TRUE)]
post.cumsum <- cumsum(post.ordered)
keep <- 1:which(post.cumsum>0.95)[1]
bci <- range(vals.ordered[keep])
bci + 27
## result is the same

## **************************************************
## Biodiversity and ecosystem function
## **************************************************

## 1.
## Read and examine the data
bb <- read.csv('https://www.sfu.ca/~lmgonigl/materials-qm/data/biodiv.csv',
               stringsAsFactors=TRUE,
               strip.white=TRUE)
head(bb)

## 2.
plot(x=bb$species, bb$co2flux, pch=16, las=1,
     xlab='Number of species', ylab=expression(CO[2]~' flux'))

## 3.
out.lm <- lm(co2flux ~ species, data=bb)
summary(out.lm)
curve(coef(out.lm)[1] + coef(out.lm)[2]*x, from=2, to=19, col='red', add=TRUE)

## 4.
out.nls.mm <- nls(co2flux ~ a + b * (species-2)/(c+(species-2)),
                  data=bb,
                  start=list(a=-1000, b=1000, c=1))
summary(out.nls.mm)
coef(out.nls.mm)
curve(coef(out.nls.mm)['a'] + coef(out.nls.mm)['b'] * (x-2) /
      (coef(out.nls.mm)['c']+(x-2)), from=2, to=19, col='blue', add=TRUE)

## 5.
BIC(out.lm)
BIC(out.nls.mm)
## non-linear is a better fit

## 6.
bic.vec <- c(lm=BIC(out.lm), nls=BIC(out.nls.mm))
bic.delta <- bic.vec - min(bic.vec)
L <- exp(-0.5 * bic.delta)
w <- L/sum(L)
w

## 7.
aic.vec <- c(lm=AIC(out.lm), nls=AIC(out.nls.mm))
aic.delta <- aic.vec - min(aic.vec)
L <- exp(-0.5 * aic.delta)
w <- L/sum(L)
w

## 8.
## non-linear has better support under both BIC and AIC

## 9.
##
## null hypothesis testing (if possible) would simply allow you to
## reject (or fail to reject) the null model.  BIC, on the other hand,
## gives you posterior weights for each of the models, which provide a
## measure of how much better one model is compared to the other.

## 10.
##
## I wasn't able to get this working with species on its original
## scale (because 19^c gets very large, very quickly).  Instead, I
## re-scaled species so that it ranges from 0-1 (this is stored in a
## new variable species.sc).  The model runs fine on this new scale.
bb$species.sc <- bb$species-min(bb$species)
bb$species.sc <- bb$species.sc/max(bb$species.sc)
out.nls.pow <- nls(co2flux ~ a + b*(species.sc)^c,
                   data=bb,
                   start=list(a=-1000, b=1500, c=1/5))

plot(x=bb$species.sc, bb$co2flux, pch=16, las=1,
     xlab='Number of species', ylab=expression(CO[2]~' flux'))
curve(coef(out.nls.pow)['a'] +
      coef(out.nls.pow)['b']*x^coef(out.nls.pow)['c'],
      from=0, to=1, col='green4', add=TRUE)

BIC(out.nls.mm)
BIC(out.nls.pow)

## re-run Michaelis-Menten on same scale
out.nls.mm.2 <- nls(co2flux ~ a + b * species.sc/(c+species.sc),
                  data=bb,
                  start=list(a=-1000, b=1000, c=1/5))
summary(out.nls.mm.2)
## same BIC as with non-scaled species (as it should be)
BIC(out.nls.mm)
BIC(out.nls.mm.2)

## add curve to plot
curve(coef(out.nls.mm.2)['a'] + coef(out.nls.mm.2)['b'] * x /
      (coef(out.nls.mm.2)['c']+x), from=0, to=1, col='blue', add=TRUE)

BIC(out.nls.mm)
BIC(out.nls.pow)
bic.vec <- c(mm=BIC(out.nls.mm),
             pow=BIC(out.nls.pow))
bic.delta <- bic.vec - min(bic.vec)
L <- exp(-0.5 * bic.delta)
w <- L/sum(L)
w
## the power function fits better, but not by a lot