Rooks

R Markdown

library(rjags)
library(tidyverse)
library(rjags)
library(ggmcmc)
library(polspline)
library(propagate)

d<-read.csv("/home/aqm/course/data/rooks.csv")
d$sward<-as.factor(as.numeric(as.factor(d$sward))-1) ## Code as zeros and ones

mod<-lm(data=d,rate~density*sward)
anova(mod)

## Analysis of Variance Table
## 
## Response: rate
##               Df  Sum Sq Mean Sq  F value Pr(>F)    
## density        1 1.11596 1.11596 193.7391 <2e-16 ***
## sward          1 0.01638 0.01638   2.8442 0.0999 .  
## density:sward  1 0.00017 0.00017   0.0301 0.8632    
## Residuals     38 0.21888 0.00576                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

g0<-ggplot(d,aes(x=density,y=rate,col=sward)) + geom_point()
g0

HDmod<-nls(rate~a*density/(1+a*density*H),data = d,start = list(a =1,H=2)) 

summary(HDmod)

## 
## Formula: rate ~ a * density/(1 + a * density * H)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a  0.22643    0.03401   6.657 5.65e-08 ***
## H  0.46670    0.28002   1.667    0.103    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07498 on 40 degrees of freedom
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 2.355e-06

confint(HDmod)

##         2.5%     97.5%
## a  0.1699599 0.3067778
## H -0.1215280 0.9841328

g0<-ggplot(d,aes(x=density,y=rate,col=sward)) +geom_point()
fig2<-g0+geom_smooth(method="nls",formula=y~a*x/(1+a*x*H),method.args=list(start = c(a = 1,H=2)), se=FALSE) 
fig2

hol_mod<-"
model{

#Likelihood
for( i in 1:n)
{
  y[i] ~ dgamma( ((mu[i]^2)/sig), (mu[i]/sig)) ## Reparameterised gamma
  mu[i]<- a*x[i]/(1+a*x[i]*H)  ## The non linear model to be fit

}

#priors
# The user must provide their own prior estimates of both the best guess expected value for a (a_mu) and H (H_mu) and an estimated standard deviation representing uncertainty regarding these values (a_sig and H_sig.

a~dgamma( ((a_mu^2)/a_sig), (a_mu/a_sig)) 
H~dgamma( ((H_mu^2)/H_sig), (H_mu/H_sig))


## Hyper parameter for the gamma on the likelihood.

sig~dunif(0,10) 

}"

a_mu=0.2
a_sig= 0.3
H_mu =1
H_sig= 0.002

prior_dat<-list(
  n=1,
  x=d$density,
  y=d$rate,
  a_mu=a_mu,
  a_sig=a_sig,
  H_mu=H_mu,
  H_sig=H_sig)

dat<-list(
  n=length(d$density),
  x=d$density,
  y=d$rate,
  a_mu=a_mu,
  a_sig=a_sig,
  H_mu=H_mu,
  H_sig=H_sig)

prior_mod <- jags.model(textConnection(hol_mod), data = prior_dat, n.chains = 1)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 3
##    Total graph size: 110
## 
## Initializing model

update(prior_mod, 10000)
prior_mod_sim <- coda.samples(model = prior_mod, variable.names = c("a","H"), n.iter = 10000)
mod <- jags.model(textConnection(hol_mod), data = dat, n.chains = 1)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 42
##    Unobserved stochastic nodes: 3
##    Total graph size: 234
## 
## Initializing model

update(mod, 10000)
    
mod_sim <- coda.samples(model = mod, variable.names = c("a","H"), n.iter = 10000)

ms <-ggs(mod_sim) 
ms %>% spread(Parameter, value) -> ms
    
x<-seq(0,10,length=100)

pred<-function(x=x,mt=ms){
      pred<-mt$a*x/(1+mt$a*x*mt$H)
      return(list(quantile(pred,c(0.025,0.5,0.975))))
}
    
pred_mod<-data.frame(x=x,do.call("rbind",sapply(x,pred)))
names(pred_mod)<-c("x","lwr","mid","upr")
    

fig8<-g0+geom_point()+ geom_line(data=pred_mod,aes(x=x,y=lwr),col="blue") + geom_line(data=pred_mod,aes(x=x,y=upr),col="red") + geom_line(data=pred_mod,aes(x=x,y=mid),col="black")

fig8

summary(prior_mod_sim)

## 
## Iterations = 11001:21000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##     Mean     SD Naive SE Time-series SE
## H 0.9997 0.0445 0.000445      0.0005908
## a 1.6962 1.3550 0.013550      0.0291205
## 
## 2. Quantiles for each variable:
## 
##     2.5%    25%   50%   75% 97.5%
## H 0.9136 0.9694 0.999 1.029 1.088
## a 0.1800 0.7215 1.324 2.296 5.247

summary(mod_sim)

## 
## Iterations = 11001:21000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##     Mean      SD  Naive SE Time-series SE
## H 0.9988 0.04437 0.0004437      0.0006634
## a 0.3119 0.02296 0.0002296      0.0004008
## 
## 2. Quantiles for each variable:
## 
##     2.5%    25%    50%    75%  97.5%
## H 0.9142 0.9685 0.9977 1.0280 1.0886
## a 0.2721 0.2960 0.3102 0.3255 0.3623

Subsetting

a_mu=0.3
a_sig= 0.3
H_mu = 1
H_sig= 0.02

dd<-d
d0<-filter(d,sward==0)
dat0<-list(
  n=length(d0$density),
  x=d0$density,
  y=d0$rate,
  a_mu=a_mu,
  a_sig=a_sig,
  H_mu=H_mu,
  H_sig=H_sig)

mod0 <- jags.model(textConnection(hol_mod), data = dat0, n.chains = 1)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 21
##    Unobserved stochastic nodes: 3
##    Total graph size: 150
## 
## Initializing model

update(mod0, 10000)
mod_sim0 <- coda.samples(model = mod0, variable.names = c("a","H"), n.iter = 10000)

ms <-ggs(mod_sim0) 
ms %>% spread(Parameter, value) -> ms
x<-seq(0,3,length=100)

pred0<-function(x=x,mt=ms){
      pred<-mt$a*x/(1+mt$a*x*mt$H)
      return(list(quantile(pred,c(0.025,0.5,0.975))))
}
    
pred_mod0<-data.frame(x=x,do.call("rbind",sapply(x,pred)))
names(pred_mod0)<-c("x","lwr","mid","upr")

d1<-filter(d,sward==1)
dat1<-list(
  n=length(d1$density),
  x=d1$density,
  y=d1$rate,
  a_mu=a_mu,
  a_sig=a_sig,
  H_mu=H_mu,
  H_sig=H_sig)

mod1 <- jags.model(textConnection(hol_mod), data = dat1, n.chains = 1)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 21
##    Unobserved stochastic nodes: 3
##    Total graph size: 150
## 
## Initializing model

update(mod1, 10000)
mod_sim1 <- coda.samples(model = mod1, variable.names = c("a","H"), n.iter = 10000)

ms <-ggs(mod_sim1) 
ms %>% spread(Parameter, value) -> ms
x<-seq(0,3,length=100)

pred1<-function(x=x,mt=ms){
      pred<-mt$a*x/(1+mt$a*x*mt$H)
      return(list(quantile(pred,c(0.025,0.5,0.975))))
}
    
pred_mod1<-data.frame(x=x,do.call("rbind",sapply(x,pred)))
names(pred_mod1)<-c("x","lwr","mid","upr")

g0+geom_point() + 
 geom_line(data=pred_mod0,aes(x=x,y=mid),col="blue") +
 geom_line(data=pred_mod1,aes(x=x,y=upr),col="red") 

summary(mod_sim0)

## 
## Iterations = 11001:21000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##     Mean      SD  Naive SE Time-series SE
## H 0.9796 0.12327 0.0012327       0.002879
## a 0.3402 0.04096 0.0004096       0.001047
## 
## 2. Quantiles for each variable:
## 
##     2.5%    25%    50%    75%  97.5%
## H 0.7528 0.8939 0.9735 1.0597 1.2404
## a 0.2755 0.3118 0.3351 0.3627 0.4353

summary(mod_sim1)

## 
## Iterations = 11001:21000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##     Mean      SD  Naive SE Time-series SE
## H 1.0116 0.13443 0.0013443       0.002634
## a 0.2938 0.04534 0.0004534       0.001635
## 
## 2. Quantiles for each variable:
## 
##     2.5%    25%    50%   75%  97.5%
## H 0.7654 0.9189 1.0065 1.098 1.2886
## a 0.2272 0.2640 0.2875 0.316 0.3905

Using a common prior for H and a sward specific a

hol_mod2<-"
model{

#Likelihood
for( i in 1:n)
{
  y[i] ~ dgamma( ((mu[i]^2)/sig), (mu[i]/sig)) ## Reparameterised gamma
  mu[i]<- a[sward[i]]*x[i]/(1+a[sward[i]]*x[i]*H)  ## The non linear model to be fit
  
}

#priors
# 

a[1]~dgamma( ((a_mu^2)/a_sig), (a_mu/a_sig)) 
a[2]~dgamma( ((a_mu^2)/a_sig), (a_mu/a_sig)) 
difa<-a[1]-a[2]
H~dgamma( ((H_mu^2)/H_sig), (H_mu/H_sig))


## Hyper parameter for the gamma on the likelihood.

sig~dunif(0,10) 

}"

a_mu=0.3
a_sig= 0.3
H_mu = 1
H_sig= 0.1

dat<-list(
  n=length(d$density),
  x=d$density,
  sward=as.numeric(d$sward),
  y=d$rate,
  a_mu=a_mu,
  a_sig=a_sig,
  H_mu=H_mu,
  H_sig=H_sig)

mod <- jags.model(textConnection(hol_mod2), data = dat, n.chains = 1)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 42
##    Unobserved stochastic nodes: 4
##    Total graph size: 328
## 
## Initializing model

update(mod, 10000)
mod_sim <- coda.samples(model = mod, variable.names = c("a","H"), n.iter = 10000)
summary(mod_sim)

## 
## Iterations = 11001:21000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##        Mean      SD  Naive SE Time-series SE
## H    0.9736 0.18462 0.0018462       0.006885
## a[1] 0.3472 0.05335 0.0005335       0.001922
## a[2] 0.2739 0.03801 0.0003801       0.001323
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%    75%  97.5%
## H    0.6226 0.8456 0.9680 1.0950 1.3434
## a[1] 0.2618 0.3093 0.3410 0.3780 0.4711
## a[2] 0.2124 0.2474 0.2687 0.2956 0.3629

ms <-ggs(mod_sim) 
mt<-filter(ms,Parameter!="H")
ggs_caterpillar(mt)