[R-sig-eco] WinBUGS linear regression & tree growth
Ben Bolker
bbolker at gmail.com
Wed Nov 24 23:15:03 CET 2010
I tweaked a variety of things. I worked in JAGS rather than WinBUGS
(but I think this code should run in WinBUGS too, if you change the
'jags' call back to a 'bugs' call). I thought WinBUGS might be
case-insensitive (but later googling suggests it is not), so changed
your i index to j ... I switched the priors from lognormal to normal, as
making them lognormal with a big variance led to a spike in the
posterior at zero ...
library(R2jags)
## generate data for sapling growth (y) with respect to irradiance (I)
set.seed(1001)
n = 50 # number of data points to generate
a = 1.7 # specify intercept
b = 0.11 # specify slope
sigma = 2.2 # specify standard deviation
I = 10 + 40*runif(n) # generate independent variable I
y = a + b*I + sigma*rnorm(n) # generate dependent variable y
lm(y~I) ## get best-fit answers for this particular data set
### WinBUGS model (put in separate file, e.g.,"linearmodel.txt" to run
####
jmodel <- function() {
# Priors
sigma ~ dunif(0.01,10)
## a ~ dlnorm(0,0.0001)
## b ~ dlnorm(0,0.0001)
a ~ dnorm(0,0.0001)
b ~ dnorm(0,0.0001)
## Determine means and distribution for dependent variable
## (which determines the likelihood function)
tau2 <- 1/(sigma*sigma)
## Precision parameter for a normal distribution
for (j in 1:n) {
# Mean value
mu[j] <- a + b*I[j]
# Observed value
y[j] ~ dnorm(mu[j], tau2)
}
}
parms <- c("a","b","sigma") ## Give names of parameters
dat = list("n","I","y") ## Set data
startval = list(list(a = 1.7, b = 0.11 , sigma=2.2))
## Set initial parameter estimates
##### run WinBUGS ###
tol_species <- jags(data=dat,
inits = startval,
parameters = parms,
mod=jmodel,
n.chains = 1,
n.iter = 15000,
n.burnin = 5000,
n.thin = 1)
library(coda)
xyplot(as.mcmc(tol_species))
densityplot(as.mcmc(tol_species),layout=c(2,2))
On 10-11-24 12:46 PM, Seth W Bigelow wrote:
> Dear all:
>
> I'm estimating the crossover point irradiance (irradiance at which growth
> rankings reverse) for some shade-tolerant and intolerant trees, in
> Bayesian mode. First I need to estimate growth with respect to light for
> individual species, and I'm using the R interface to WinBUGS. This has
> worked fine with some non-linear models of growth, but I'm getting hung up
> on the simple linear model, and the error-trapping information in WinBUGS
> is really opaque. Attached is some code that stochastically generates data
> for a linear growth model, but which hangs up in WinBUGS. Is anyone able
> to see what I'm doing wrong?
>
>
> ### WinBUGS linear regression
>
> library(R2WinBUGS)
>
> ## generate data for sapling growth (y) with respect to irradiance (I)
>
> n = 50 # number of data
> points to generate
> a = 1.7 # specify
> intercept
> b = 0.11 # specify slope
> sigma = 2.2 # specify standard
> deviation
> I = 10 + 40*runif(n) # generate independent
> variable I
> y = a + b*I + sigma*rnorm(n) # generate dependent
> variable y
>
> ### WinBUGS model (put in separate file, e.g.,"linearmodel.txt" to run
> ####
>
> model {
> # Priors
> sigma ~ dunif(0.01,10)
> log.a ~ dnorm(0,0.0001)
> log.b ~ dnorm(0,0.0001)
> log(a) <- log.a
> log(b) <- log.b
>
> # Determine means and distribution for dependent variable
> # (which determines the likelihood function)
>
> tau2 <- 1/(sigma*sigma) # Precision parameter for a normal
> distribution
>
> for (i in 1:n) {
>
> # Mean value
> mu[i] <- a + b*I[i]
>
> # Observed value
> y[i] ~ dnorm(mu[i], tau2)
> }
> }
>
> ##### specify input data ####
>
> parms <- c("a","b","sigma") # Give
> names of parameters
> dat = list("n","I","y") # Set data
> startval = list(list(a = 1.7, b = 0.11 , sigma=2.2)) #
> Set initial parameter estimates
> directory = (your directory here)
>
> ##### run WinBUGS ###
>
> tol_species <- bugs(
> data=dat,
> inits = startval,
> parameters = parms,
> "linearmodel.txt",
> n.chains = 1,
> n.iter = 15000,
> n.burnin = 5000,
> n.thin = 1,
> codaPkg = FALSE,
> DIC=FALSE,
> working.directory = directory,
> clearWD = FALSE)
>
> ###### END ################################
>
>
>
>
> Dr. Seth W. Bigelow
> Biologist, USDA-FS Pacific Southwest Research Station
> 1731 Research Park Drive, Davis California
> sbigelow at fs.fed.us / ph. 530 759 1718
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
More information about the R-sig-ecology
mailing list