d<-read.csv("/home/aqm/course/data/Hollings.csv")
plot(d,pch=21,bg=2)

Jags linear

library(rjags)
## Loading required package: coda
## Linked to JAGS 4.2.0
## Loaded modules: basemod,bugs
reg_mod<-"
model{

#Likelihood
  for( i in 1:n)
    {
      y[i]~dnorm(mu[i], tau)
      mu[i]<-b0+b1*x[i]
    }

#priors
b0~dnorm(0,.01)
b1~dnorm(0,.01)
tau<-pow(sd, -2)
sd~dunif(0,100)


}"
dat<-list(n=length(d$Resource),x=d$Resource,y=d$Consumption)

mod1 <- jags.model(textConnection(reg_mod), data = dat, n.chains = 3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 30
##    Unobserved stochastic nodes: 3
##    Total graph size: 133
## 
## Initializing model
update(mod1, 1000)

mod1_sim <- coda.samples(model = mod1, variable.names = c("b0","b1"), n.iter = 5000)
summary(mod1_sim)
## 
## Iterations = 2001:7000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##        Mean       SD  Naive SE Time-series SE
## b0 13.83641 0.383798 3.134e-03      0.0100769
## b1  0.01091 0.001183 9.659e-06      0.0000311
## 
## 2. Quantiles for each variable:
## 
##         2.5%      25%     50%      75%    97.5%
## b0 13.085589 13.58563 13.8369 14.08744 14.60060
## b1  0.008574  0.01013  0.0109  0.01168  0.01321
confint(lm(Consumption~Resource,data=d))
##                    2.5 %      97.5 %
## (Intercept) 13.102428065 14.59200222
## Resource     0.008568358  0.01318626
hyp_mod<-"
model{

#Likelihood
  for( i in 1:n)
    {
      y[i]~dnorm(mu[i], tau)
      mu[i]<-b0*x[i]/(b1+x[i])
    }

#priors
b0~dnorm(0,.01)
b1~dnorm(0,.01)
tau<-pow(sd, -2)
sd~dunif(0,100)


}"
mod1 <- jags.model(textConnection(hyp_mod), data = dat, n.chains = 3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 30
##    Unobserved stochastic nodes: 3
##    Total graph size: 163
## 
## Initializing model
update(mod1, 1000)

mod1_sim <- coda.samples(model = mod1, variable.names = c("b0","b1"), n.iter = 5000)
summary(mod1_sim)
## 
## Iterations = 2001:7000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##     Mean     SD Naive SE Time-series SE
## b0 19.80 0.1884 0.001538       0.004803
## b1 36.54 2.4159 0.019725       0.065050
## 
## 2. Quantiles for each variable:
## 
##     2.5%   25%   50%   75% 97.5%
## b0 19.42 19.68 19.81 19.93 20.16
## b1 31.69 34.96 36.60 38.17 41.15
nlmod<-nls(Consumption~s*Resource/(F+Resource),data = d,start = list( F = 20,s=20))
confint(nlmod)
## Waiting for profiling to be done...
##       2.5%    97.5%
## F 33.97300 43.66697
## s 19.58576 20.32930

Buntings

d<-read.csv("/home/aqm/course/data/buntings.csv")
plot(rate~density,data=d,pch=21,bg=2)

Trying to write the Bayesian model.

Use Gamma errors to avoid negative rates. Set gammas on both paprameters as priors. Use the rep-paramaterisation of gamma in terms of mean and standard deviation instead of shape and scale for clarity.

hol_mod<-"
model{

#Likelihood
  for( i in 1:n)
    {
      
    y[i] ~dgamma( ((mu[i]^2)/sig), (mu[i]/sig))
     mu[i]<- a*x[i]/(1+a*x[i]*H)
     
    }

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


sig~dunif(0.0001,1)


}"

The model still does not really work, even with informative priors derived from the quantile regression.

dat<-list(n=length(d$Resource),x=d$density,y=d$rate,a_mu=0.00117,a_sig= 0.00053,H_mu=2.64950,H_sig=0.28727)
mod1 <- jags.model(textConnection(hol_mod), data = dat, n.chains = 5)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 0
##    Unobserved stochastic nodes: 3
##    Total graph size: 310
## 
## Initializing model
update(mod1, 10000)

mod1_sim <- coda.samples(model = mod1, variable.names = c("a","H"), n.iter = 50000)
summary(mod1_sim)
## 
## Iterations = 10001:60000
## Thinning interval = 1 
## Number of chains = 5 
## Sample size per chain = 50000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD  Naive SE Time-series SE
## H 2.648940 0.53585 1.072e-03      1.068e-03
## a 0.001209 0.02251 4.502e-05      4.514e-05
## 
## 2. Quantiles for each variable:
## 
##    2.5%        25%        50%       75%    97.5%
## H 1.703  2.272e+00  2.613e+00 2.988e+00 3.80e+00
## a 0.000 1.308e-234 4.067e-118 7.514e-50 1.51e-05
plot(mod1_sim)

HDmod<-nls(rate~a*density/(1+a*density*H),data = d,start = list(a =0.001,H=2)) 
confint(HDmod)
## Waiting for profiling to be done...
##          2.5%       97.5%
## a 0.002593086 0.006694939
## H 1.713495694 1.976978655