Re: [R-sig-eco] Bayesian analysis of canopy cover/percent data

Bob O'Hara Wed, 27 Aug 2014 10:11:41 -0700

On 08/27/2014 05:22 PM, Stratford, Jeffrey wrote:

Hi everyone,


I'd like to estimate the means and sd's of canopy cover for 13 sites (9 are
birds and 4 are potential habitats). I'm using a beta distribution to
account for the fact that the data are bound by 0 and 1. For each sample, I
measured the proportion of points that have some sort of vegetation in the
canopy.

I'm trying to adapt the code that works for the other distributions but I'm
getting the error message:

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
   NA/NaN/Inf in 'y'
In addition: There were 33 warnings (use warnings() to see them)
Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
   NA/NaN/Inf in 'y'

Below are the R code, the model statement, and the data. I'm using JAGS and
R 3.1.1 in Windows Vista. I suspect the problem is with the model statement.

Any help would be greatly appreciated. I feel like the kid that knows all
the rules of baseball but can neither catch or hit!

The beta can't be 0 or 1, and you have a lot of 1's. There are a fewways of dealing with this: the easiest would be to change the 1's into1-e (where e is small, and you can check if the value chosen makes adifference).

BTW, if (as it sounds), you took some points, and counted how many hadvegetation, and how many did not, it might be better to assume that thisis a Bernoulli trial, so that the coverage follows a Binomialdistribution, and go on from there. This then gets around the problembecause a Binomial is allowed to be 0 or N (but if it is, it can besensitive to the priors).

Bob

Thanks,

Jeff



#########  R  code ##########
frag <- read.csv("f:\\brazil\\TIandFRAG.csv", header=T)
library(R2jags)
library(rjags)
setwd("f://brazil")
site <- frag$site
cover20p <- frag$cover20p/100
N <- length(frag$site)

jags.data <- list("site", "cover20p")
jags.params <- c("median", "test100MF","test100MT","test100fc","test100fa",
"test100gv","test100hm","test100mc", "test100ca","test100ct", "test10MF",
"test10MT", "test10fc","test10fa", "test10gv", "test10hm", "test10mc",
"test10ca", "test10ct",
"test1MF", "test1MT", "test1fc",  "test1fa",  "test1gv", "test1hm",
"test1mc", "test1ca", "test1ct", "t1est1_con","t2est10_con","t3est100_con",
"t4est1_100","t5est1_10","t6est10_100")
#inits1 <- list(a=0, sd=0)
#inits2 <- list(a=100, sd=50)
#jags.inits <- list(inits1, inits2)

jags.inits <- function() {
   list(a=c(0,0,0,0,0,0,0,0,0,0,0,0,0), sd=1)}

jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params,
n.iter=1000000, n.burnin=20000, model.file="fragmodelbeta.txt")

my.coda <- as.mcmc(jagsfit)

summary(my.coda, quantiles=c(0.05, 0.25,0.5,0.75, 0.95))

print(jagsfit, digits=3)

######## MODEL ########################################

model{
for (i in 1:227) {

log(mean[i]) <- a[site[i]]
cover20p[i] ~ dbeta(1,0.5)
}
for (i in 1:13){
a[i] ~ dnorm(0, tau)
median[i] <- exp(a[i])
}
sd ~ dunif(0, 10)
tau <- 1 / (sd*sd) # precision

}

######## DATA #################################################

# normally this would be on a flash drive

structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7,
7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8,
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 11, 11,
11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13,
13, 13, 13, 13, 13, 13, 13, 0.95, 0.8, 0.85, 0.9, 0.35, 1, 1,
1, 0.95, 0.55, 0.9, 0.85, 0.7, 0.65, 0.05, 0.6, 1, 1, 0.85, 0.9,
0, 0.45, 1, 0.7, 0.95, 0.5, 0.95, 0.6, 0.65, 0.7, 0.4, 0.85,
0.6, 0.95, 0.75, 0.9, 0.85, 0.75, 0.7, 0.85, 0.3, 0.7, 0.8, 0.7,
0.75, 0.8, 0.75, 0.95, 0.9, 0.05, 0.85, 0.6, 0.65, 0.5, 0.85,
0.95, 0.85, 0.25, 0.75, 1, 0.65, 0.95, 0.8, 0.9, 0.6, 0.8, 1,
0.2, 0.8, 0.4, 1, 0.95, 0.4, 1, 1, 0.95, 0.45, 0.2, 0.7, 0.95,
0.7, 0.8, 0.5, 0.85, 0.55, 0, 0.25, 0.45, 1, 0.95, 1, 0.9, 0.6,
0.35, 0.95, 0.3, 1, 1, 0.5, 0.4, 0.9, 1, 0.7, 1, 0.9, 1, 0.4,
0.55, 0.8, 0.7, 1, 0, 0.8, 0, 0.7, 0.5, 0.8, 0.75, 0, 0.45, 0.1,
0, 0.4, 0.55, 0.4, 1, 0.9, 0.9, 0.15, 0.55, 0.35, 0.9, 0.65,
0.25, 1, 0.85, 1, 0.95, 0.7, 0.5, 0.7, 0.2, 0.95, 1, 1, 0.25,
0.85, 0.5, 0.8, 0.75, 0.85, 0.7, 0.95, 0.05, 0.65, 0.65, 1, 1,
1, 0.65, 0.4, 0.6, 0.9, 0.85, 0.75, 0.5, 0.65, 1, 0.65, 0.55,
0.75, 0.4, 0.9, 0.35, 1, 1, 0.4, 0.5, 0.8, 0.95, 0.95, 0.55,
0.7, 0.85, 0.8, 0.8, 0.65, 1, 0.6, 0.5, 1, 0.8, 1, 0.45, 1, 1,
0.8, 0.85, 1, 1, 1, 1, 0.5, 0.6, 0.15, 0.75, 0.6, 0.1, 0.05,
0, 1, 0.6, 0.1, 0.35, 0.9, 0.9, 0.95, 0.95, 0.9, 0.55, 0.65,
0.9, 0.4, 1, 0.65, 0.5, 0.8), .Dim = c(227L, 2L), .Dimnames = list(
     NULL, c("site2", "cover20p2")))



--
Bob O'Hara

Biodiversity and Climate Research Centre
Senckenberganlage 25
D-60325 Frankfurt am Main,
Germany

Tel: +49 69 7542 1863
Mobile: +49 1515 888 5440
WWW:   http://www.bik-f.de/root/index.php?page_id=219
Blog: http://blogs.nature.com/boboh
Journal of Negative Results - EEB: www.jnr-eeb.org

_______________________________________________
R-sig-ecology mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] Bayesian analysis of canopy cover/percent data

Reply via email to