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

[Bob O'Hara]

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 few 
ways of dealing with this: the easiest would be to change the 1's into 
1-e (where e is small, and you can check if the value chosen makes a 
difference).

BTW, if (as it sounds), you took some points, and counted how many had 
vegetation, and how many did not, it might be better to assume that this 
is a Bernoulli trial, so that the coverage follows a Binomial 
distribution, and go on from there. This then gets around the problem 
because a Binomial is allowed to be 0 or N (but if it is, it can be 
sensitive 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