aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--.gitignore5
-rw-r--r--README.md21
-rw-r--r--analysis_bp_im_prev.R50
-rw-r--r--model_bp_im_prev.R45
-rw-r--r--sample_bp_im_prev.csv11
5 files changed, 131 insertions, 1 deletions
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..63eab7f
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,5 @@
+.Rproj.user
+.Rhistory
+.RData
+bp-intermediate-host.Rproj
+bp_im_prev.csv
diff --git a/README.md b/README.md
index dc12e5c..29b3c1c 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,23 @@
bp-intermediate-host
====================
-Analysis of raccoon roundworm prevalence in suburban small mammal populations
+A mixed-model approach to analysis of larval raccoon roundworm (*Baylisascaris procyonis*) prevalence in suburban small mammal populations, focused on how prevalence is related to human population density. Prevalence of infection (logistic) and infection intensity (Poisson-lognormal) are modeled simultaneously.
+
+The results are published in the following article:
+
+[*Kellner, Kenneth F.; Page, L. Kristen; Downey, Mark; McCord, Sarah E. 2012. Effects of urbanization on prevalence of Baylisascaris procyonis in intermediate host populations. Journal of Wildlife Diseases 48:1083-1087*](http://www.jwildlifedis.org/doi/abs/10.7589/2011-09-267)
+
+Metadata
+====================
+
+The *analysis_bp_im_prev.R* file contains the analysis described above implemented in JAGS (the original publication version used WinBUGS). Additional information is provided in the comments of the file. The model in the BUGS language is contained in the file *model_bp_im_prev.R*.
+
+A sample of the original dataset used for this analysis is contained the file *sample_bp_im_prev.csv*. Each row is an individual mouse captured and sampled for larval raccoon roundworm infection. The columns in the dataset are described below.
+
+1. year07 - Study year; if 1, the mouse was captured in 2007, 0 if captured in 2006.
+2. sex - Sex of the mouse; 1 if male and 0 if female.
+3. mass - Mass in grams of the mouse at time of capture.
+4. dens - Human population density (persons/km2) around capture site.
+5. positive - Infection status; 1 if infected 0 if not.
+6. burden - Number of roundworms collected from mouse after digestion.
+
diff --git a/analysis_bp_im_prev.R b/analysis_bp_im_prev.R
new file mode 100644
index 0000000..9f1a125
--- /dev/null
+++ b/analysis_bp_im_prev.R
@@ -0,0 +1,50 @@
+#Mixed logistic and Poisson-lognormal regression model of
+#Bp prevalence and burden
+
+data <- read.csv('bp_im_prev.csv', header=TRUE)
+
+#Dependent variables
+u <- data$positive
+burden <- data$burden
+
+#Independent variables
+year07 = data$year07
+sex <- data$sex #male = 1
+
+#Standardized mouse mass
+mass.raw <- data$mass
+mass <- (mass.raw - mean(mass.raw,na.rm=TRUE))/sd(mass.raw,na.rm=TRUE)
+
+#Standardized human population density
+dens.raw <- data$dens
+dens <- (dens.raw - mean(dens.raw,na.rm=TRUE))/sd(dens.raw,na.rm=TRUE)
+
+nmice <- length(u)
+
+#Create arguments for JAGS
+data <- list(u=u, burden=burden, year07=year07,
+ sex=sex, mass=mass,
+ nmice=nmice,
+ dens=dens)
+
+parameters <- c("beta0", "p.beta0",
+ "p.beta.year07", "beta.year07",
+ "p.beta.sex", "beta.sex",
+ "p.beta.mass", "beta.mass",
+ "beta.dens", "p.beta.dens",
+ "fit", "fit.new"
+ )
+
+model.file <- "model_bp_im_prev.R"
+
+#Analysis using jagsUI (http://github.com/kenkellner/jagsUI)
+#library(devtools)
+#install_github('kenkellner/jagsUI')
+library(jagsUI)
+
+out <- jags(data, inits=NULL, parameters, model.file,
+ n.thin=2,n.chains=3, n.iter=10000, n.burnin=5000,
+ parallel=TRUE)
+
+#Posterior Predictive Check
+pp.check(out,'fit','fit.new')
diff --git a/model_bp_im_prev.R b/model_bp_im_prev.R
new file mode 100644
index 0000000..a70fd74
--- /dev/null
+++ b/model_bp_im_prev.R
@@ -0,0 +1,45 @@
+model {
+
+ #Priors
+ beta0 ~ dunif(-5,5)
+ beta.year07 ~ dnorm(0,0.001)
+ beta.sex ~ dnorm(0,0.001)
+ beta.mass ~ dnorm(0,0.001)
+ beta.dens ~ dnorm(0,0.001)
+
+ p.beta0 ~ dunif(-5,5)
+ p.beta.year07 ~ dnorm(0,0.001)
+ p.beta.sex ~ dnorm(0,0.001)
+ p.beta.mass ~ dnorm(0,0.001)
+ p.beta.dens ~ dnorm(0,0.001)
+
+ eps.sigma ~ dunif(0,7)
+ eps.prec <- sqrt(1/(eps.sigma*eps.sigma))
+
+ #Likelihood
+ for (i in 1:nmice) {
+ #Prevalence
+ logit(p[i]) <- p.beta0 + p.beta.year07*year07[i] + p.beta.sex*sex[i]
+ + p.beta.dens*dens[i] + p.beta.mass*mass[i]
+ u[i] ~ dbern(p[i])
+
+ #Worm burden
+ eps[i] ~ dnorm(0,eps.prec) #Poisson-lognormal error term
+ log(lambda[i]) <- beta0 + beta.year07*year07[i] + beta.sex*sex[i]
+ + beta.dens*dens[i] + beta.mass*mass[i] + eps[i]
+ c.lambda[i] <- u[i]*lambda[i]
+ burden[i] ~ dpois(c.lambda[i])
+
+ #Posterior Predictive Check Stuff
+ res.raw[i] <- pow((burden[i] - c.lambda[i]), 2) #Pearson residual
+ res[i] <- res.raw[i]/sqrt(c.lambda[i]+0.5)
+ burden.new[i] ~ dpois(c.lambda[i])
+ res.new.raw[i] <- pow((burden.new[i] - c.lambda[i]), 2)
+ res.new[i] <- res.new.raw[i]/sqrt(c.lambda[i]+0.5)
+
+ }
+
+ fit <- sum(res[]) #discrepency of actual dataset
+ fit.new <- sum(res.new[]) #discrepency of simulated dataset
+
+} \ No newline at end of file
diff --git a/sample_bp_im_prev.csv b/sample_bp_im_prev.csv
new file mode 100644
index 0000000..eb89eae
--- /dev/null
+++ b/sample_bp_im_prev.csv
@@ -0,0 +1,11 @@
+"year07","sex","mass","dens","positive","burden"
+0,0,21.2,356.4953,0,0
+0,1,20.9,881.5279,1,10
+0,0,30.2,356.4953,0,0
+0,0,20.5,881.5279,0,0
+1,0,11.7,1590.308,0,0
+1,0,19.5,1791.651,0,0
+1,1,23.6,737.1681,1,3
+0,1,18.5,348.739,1,3
+1,1,20.8,1227.546,0,0
+0,0,26.9,881.5279,1,1