[R-sig-eco] WinBUGS linear regression & tree growth

[Ben Bolker]

  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