Jags non linear model
Duncan Golicher
12/27/2018
Buntings
library(rjags)## Loading required package: coda## Linked to JAGS 4.2.0## Loaded modules: basemod,bugslibrary(tidyverse)## ── Attaching packages ───────────────────────────────────────────────────────── tidyverse 1.2.1 ──## ✔ ggplot2 3.0.0 ✔ purrr 0.2.5
## ✔ tibble 1.4.2 ✔ dplyr 0.7.6
## ✔ tidyr 0.8.1 ✔ stringr 1.3.1
## ✔ readr 1.1.1 ✔ forcats 0.3.0## ── Conflicts ──────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(rjags)
library(ggmcmc)
library(polspline)
library(propagate)## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: tmvtnorm## Loading required package: mvtnorm## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following object is masked from 'package:tidyr':
##
## expand## Loading required package: stats4## Loading required package: gmm## Loading required package: sandwich## Loading required package: Rcpp## Loading required package: ff## Loading required package: bit## Attaching package bit## package:bit (c) 2008-2012 Jens Oehlschlaegel (GPL-2)## creators: bit bitwhich## coercion: as.logical as.integer as.bit as.bitwhich which## operator: ! & | xor != ==## querying: print length any all min max range sum summary## bit access: length<- [ [<- [[ [[<-## for more help type ?bit##
## Attaching package: 'bit'## The following object is masked from 'package:base':
##
## xor## Attaching package ff## - getOption("fftempdir")=="/tmp/Rtmp4oCmZo"## - getOption("ffextension")=="ff"## - getOption("ffdrop")==TRUE## - getOption("fffinonexit")==TRUE## - getOption("ffpagesize")==65536## - getOption("ffcaching")=="mmnoflush" -- consider "ffeachflush" if your system stalls on large writes## - getOption("ffbatchbytes")==16777216 -- consider a different value for tuning your system## - getOption("ffmaxbytes")==536870912 -- consider a different value for tuning your system##
## Attaching package: 'ff'## The following objects are masked from 'package:bit':
##
## clone, clone.default, clone.list## The following objects are masked from 'package:utils':
##
## write.csv, write.csv2## The following objects are masked from 'package:base':
##
## is.factor, is.ordered## Loading required package: minpack.lmd<-read.csv("/home/aqm/course/data/buntings.csv")
plot(rate~density,data=d,pch=21,bg=2)
Convert density to cm2
d$densitycm2<-d$density/10000
plot(rate~densitycm2,data=d,pch=21,bg=2)
Model fit using least squares
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.976978655HDmod<-nls(rate~a*densitycm2/(1+a*densitycm2*H),data = d,start = list(a =1,H=2))
summary(HDmod)##
## Formula: rate ~ a * densitycm2/(1 + a * densitycm2 * H)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 39.6636 8.6245 4.599 9.22e-06 ***
## H 1.8419 0.0638 28.868 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1393 on 144 degrees of freedom
##
## Number of iterations to convergence: 10
## Achieved convergence tolerance: 2.286e-06confint(HDmod)## Waiting for profiling to be done...## 2.5% 97.5%
## a 25.930863 66.949394
## H 1.713496 1.976979require(propagate)
newdata<-data.frame(densitycm2=seq(0,0.4,length=100))
pred_model <- predictNLS(HDmod, newdata=newdata,nsim = 10000)## predictNLS: Propagating predictor value #1...## predictNLS: Propagating predictor value #2...## predictNLS: Propagating predictor value #3...## predictNLS: Propagating predictor value #4...## predictNLS: Propagating predictor value #5...## predictNLS: Propagating predictor value #6...## predictNLS: Propagating predictor value #7...## predictNLS: Propagating predictor value #8...## predictNLS: Propagating predictor value #9...## predictNLS: Propagating predictor value #10...## predictNLS: Propagating predictor value #11...## predictNLS: Propagating predictor value #12...## predictNLS: Propagating predictor value #13...## predictNLS: Propagating predictor value #14...## predictNLS: Propagating predictor value #15...## predictNLS: Propagating predictor value #16...## predictNLS: Propagating predictor value #17...## predictNLS: Propagating predictor value #18...## predictNLS: Propagating predictor value #19...## predictNLS: Propagating predictor value #20...## predictNLS: Propagating predictor value #21...## predictNLS: Propagating predictor value #22...## predictNLS: Propagating predictor value #23...## predictNLS: Propagating predictor value #24...## predictNLS: Propagating predictor value #25...## predictNLS: Propagating predictor value #26...## predictNLS: Propagating predictor value #27...## predictNLS: Propagating predictor value #28...## predictNLS: Propagating predictor value #29...## predictNLS: Propagating predictor value #30...## predictNLS: Propagating predictor value #31...## predictNLS: Propagating predictor value #32...## predictNLS: Propagating predictor value #33...## predictNLS: Propagating predictor value #34...## predictNLS: Propagating predictor value #35...## predictNLS: Propagating predictor value #36...## predictNLS: Propagating predictor value #37...## predictNLS: Propagating predictor value #38...## predictNLS: Propagating predictor value #39...## predictNLS: Propagating predictor value #40...## predictNLS: Propagating predictor value #41...## predictNLS: Propagating predictor value #42...## predictNLS: Propagating predictor value #43...## predictNLS: Propagating predictor value #44...## predictNLS: Propagating predictor value #45...## predictNLS: Propagating predictor value #46...## predictNLS: Propagating predictor value #47...## predictNLS: Propagating predictor value #48...## predictNLS: Propagating predictor value #49...## predictNLS: Propagating predictor value #50...## predictNLS: Propagating predictor value #51...## predictNLS: Propagating predictor value #52...## predictNLS: Propagating predictor value #53...## predictNLS: Propagating predictor value #54...## predictNLS: Propagating predictor value #55...## predictNLS: Propagating predictor value #56...## predictNLS: Propagating predictor value #57...## predictNLS: Propagating predictor value #58...## predictNLS: Propagating predictor value #59...## predictNLS: Propagating predictor value #60...## predictNLS: Propagating predictor value #61...## predictNLS: Propagating predictor value #62...## predictNLS: Propagating predictor value #63...## predictNLS: Propagating predictor value #64...## predictNLS: Propagating predictor value #65...## predictNLS: Propagating predictor value #66...## predictNLS: Propagating predictor value #67...## predictNLS: Propagating predictor value #68...## predictNLS: Propagating predictor value #69...## predictNLS: Propagating predictor value #70...## predictNLS: Propagating predictor value #71...## predictNLS: Propagating predictor value #72...## predictNLS: Propagating predictor value #73...## predictNLS: Propagating predictor value #74...## predictNLS: Propagating predictor value #75...## predictNLS: Propagating predictor value #76...## predictNLS: Propagating predictor value #77...## predictNLS: Propagating predictor value #78...## predictNLS: Propagating predictor value #79...## predictNLS: Propagating predictor value #80...## predictNLS: Propagating predictor value #81...## predictNLS: Propagating predictor value #82...## predictNLS: Propagating predictor value #83...## predictNLS: Propagating predictor value #84...## predictNLS: Propagating predictor value #85...## predictNLS: Propagating predictor value #86...## predictNLS: Propagating predictor value #87...## predictNLS: Propagating predictor value #88...## predictNLS: Propagating predictor value #89...## predictNLS: Propagating predictor value #90...## predictNLS: Propagating predictor value #91...## predictNLS: Propagating predictor value #92...## predictNLS: Propagating predictor value #93...## predictNLS: Propagating predictor value #94...## predictNLS: Propagating predictor value #95...## predictNLS: Propagating predictor value #96...## predictNLS: Propagating predictor value #97...## predictNLS: Propagating predictor value #98...## predictNLS: Propagating predictor value #99...## predictNLS: Propagating predictor value #100...conf_model <- pred_model$summary
newdata<-data.frame(newdata,conf_model)g0<-ggplot(data=d,aes(x=densitycm2,y=rate))
g1<-g0+geom_point()
g1+geom_smooth(method="nls",formula=y~a*x/(1+a*x*H),method.args=list(start = c(a = 1,H=2)), se=FALSE,col="black") + geom_line(data=newdata,aes(x=densitycm2,y=Prop.2.5.),col="blue") + geom_line(data=newdata,aes(x=densitycm2,y=Prop.97.5.),col="blue")
x<-seq(0,0.4,length=100)
fit<-function(x,a=25,H=1.7) a*x/(1+a*x*H)
dd<-data.frame(x=x,y1=fit(x),y2=fit(x,a=67,H=1.97))
library(ggplot2)
g0<-ggplot(data=d,aes(x=densitycm2,y=rate))
g1<-g0+geom_point()
g1+geom_smooth(method="nls",formula=y~a*x/(1+a*x*H),method.args=list(start = c(a = 1,H=2)), se=FALSE) +geom_line(data=dd,aes(x=x,y=y1),col="red") +geom_line(data=dd,aes(x=x,y=y2),col="blue")
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.
Priors
priors<-"
model{
#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)
}"
dat<-list(a_mu=40,a_sig= 4,H_mu=2.64950,H_sig=0.28)
mod1 <- jags.model(textConnection(priors), 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: 17
##
## Initializing modelupdate(mod1, 10000)
priors_mod <- coda.samples(model = mod1, variable.names = c("a","H"), n.iter = 50000)
plot(priors_mod)
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$densitycm2,y=d$rate,a_mu=100,a_sig=10,H_mu=2.6,H_sig=0.4)
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 modelupdate(mod1, 10000)
mod1_sim <- coda.samples(model = mod1, variable.names = c("a","H"), n.iter = 50000)
ms <-ggs(mod1_sim) post_samples = filter(ms, Parameter == "a")$value
fit.posterior <- logspline(post_samples, ubound=200, lbound=0)## Warning in logspline(post_samples, ubound = 200, lbound = 0): too much data
## close together## Warning in logspline(post_samples, ubound = 200, lbound = 0): re-ran with
## oldlogsplineposterior_w <- dlogspline(100, fit.posterior)
posterior_w## [1] 0.1264911summary(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.599 0.6336 0.001267 0.001267
## a 99.997 3.1632 0.006326 0.006318
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## H 1.51 2.147 2.547 2.993 3.982
## a 93.89 97.847 99.963 102.108 106.299plot(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