An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-ecology/attachments/20110408/d4407b73/attachment.pl>
zero-truncated Poisson
6 messages · Ben Bolker, Peter Solymos, Stratford, Jeffrey
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
On 04/08/2011 04:10 PM, Stratford, Jeffrey wrote:
Hi everyone,
I'm attempting to run a zero-truncated Poisson using VGAM and I run into
this error:
Error in if ((temp <- sum(wz[, 1:M, drop = FALSE] < wzepsilon)))
warning(paste(temp, :
argument is not interpretable as logical
Here's the code I'm using
blja <- read.csv("g:\\blja\\blja.csv", header=T)
library(VGAM)
glm1 <- vglm(blja$taken ~ blja$mass + blja$year,family=pospoisson,
data=blja)
Just a *very* quick suggestion: try glm1 <- vglm(taken ~ mass + year,family=pospoisson, data=blja) and see if that helps. Also try summary() on your data (and/or sapply(blja,class)) to see if your data types are what you think they are (hint: R often converts numeric data to factors if it has glitches in it) If neither of those helps, could you use dput() to post your data (or a subset thereof) in an easier-to-use format?
We're looking at the number of acorns taken by jays using year (two years) and how much (mass) they're taking. Our data looks like this: year size distance taken mass one small 30.8735 10 13.8758 one small 121.505 2 2.77516 one small 46.055 12 16.65096 one small 46.055 1 1.38758 one small 46.055 4 5.55032 one small 9.343 1 1.38758 one small 46.055 10 13.8758 one small 46.055 3 4.16274 one small 46.055 5 6.9379 one small 85.271 5 6.9379 one small 85.271 1 1.38758 one small 85.271 5 6.9379 one small 85.271 3 4.16274 ... Any suggestions (and alternatives!) would be appreciated! Thanks a bunch, Jeff ************************************ Dr. Jeffrey A. Stratford Department of Health and Biological Sciences 84 W. South Street Wilkes University, PA 18766 jeffrey.stratford at wilkes.edu 570-408-4761 (office) 570-332-2942 (cell) http://web.wilkes.edu/jeffrey.stratford/ ************************************ [[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
-----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.10 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/ iEYEARECAAYFAk2fbNkACgkQc5UpGjwzenMcAQCfcxwQ51AV3QspEu2goI2ECMnF r1AAn0/PgmgWqvzxcr0ywW1+acUL5q5b =ek5d -----END PGP SIGNATURE-----
Thanks Ben,
That didn't do the trick. I checked all the classes
year size distance taken mass
"factor" "factor" "numeric" "integer" "numeric"
And used the notation that you suggested and I did get the same error.
Just so it's not such a black box for me, what makes vector GLMs
different from "ordinary" GLMs?
Thanks,
Jeff
Here's all the data.
structure(list(year = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label =
c("one",
"two"), class = "factor"), size = structure(c(2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L), .Label = c("large", "small"), class = "factor"), distance =
c(30.8735,
121.505, 46.055, 46.055, 46.055, 9.343, 46.055, 46.055, 46.055,
85.271, 85.271, 85.271, 85.271, 85.271, 85.271, 30.8735, 30.8735,
9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 152.717, 152.717,
152.717, 152.717, 152.717, 152.717, 152.717, 46.055, 46.055,
152.717, 152.717, 20.698, 20.698, 17.217, 17.217, 9.343, 9.343,
17.217, 46.055, 152.717, 46.055, 17.217, 152.717, 17.217, 17.217,
46.055, 30.8735, 30.8735, 30.8735, 5.69, 30.8735, 17.217, 30.8735,
30.8735, 30.8735, 30.8735, 30.8735, 9.343, 9.343, 17.217, 17.217,
17.217, 17.217, 17.217, 17.217, 9.343, 9.343, 17.217, 17.217,
17.217, 20.698, 17.217, 17.217, 17.217, 5.42, 17.217, 17.217,
17.217, 9.343, 30.8735, 30.8735, 30.8735, 9.343, 9.343, 9.343,
9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343,
17.217, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343,
30.8735, 30.8735, 17.217, 17.217, 46.055, 36.239, 17.217, 17.217,
17.217, 46.055, 30.8735, 30.8735, 17.217, 17.217, 17.217, 121.505,
152.717, 152.717, 17.217, 152.717, 121.505, 121.505, 121.505,
121.505, 121.505, 9.343, 121.505, 9.343, 121.505, 121.505, 30.8735,
121.505, 17.217, 17.217, 17.217, 9.343, 30.8735, 85.271, 85.271,
85.271, 85.271, 85.271, 85.271, 9.343, 85.271, 85.271, 85.271,
85.271, 20.698, 9.343, 30.8735, 17.217, 20.698, 30.8735, 17.217,
85.271, 121.505, 121.505, 85.271, 85.271, 85.271, 85.271, 121.505,
121.505, NA, 121.505, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343,
9.343, 30.8735, 17.217, 9.343, 9.343, 9.343, 30.8735, 30.8735,
30.8735, 9.343, 9.343, 30.8735, 17.217, 5.42, 17.217, 85.271,
85.271, 30.8735, 30.8735, 30.8735, 30.8735, 30.8735, 30.8735,
17.217, 85.271, 17.217, 30.8735, 30.8735, 85.271, 30.8735, 30.8735,
85.271, 30.8735, 30.8735, 30.8735, 30.8735, 30.8735, 30.8735,
17.217, 30.8735, 9.343, 30.8735, 30.8735, 30.8735, 30.8735, 9.343,
30.8735, 30.8735, 17.217, 9.343, 9.343, 9.343, 85.271, 30.8735,
30.8735, 46.055, 5.69, 85.271, 85.271, 17.217, 46.055, 85.271,
30.8735, 85.271, 46.055, 121.505, 121.505, 121.505, 17.217, 17.217,
25.508, 9.343, 17.217, 17.217, 9.343, 17.217, 9.343, 17.217,
34.35, 34.35, 46.055, 46.055, 9.343, 20.698, 17.217, 27.25, 20.54,
27.25, 20.698, 17.217, 20.698, 20.698, 15, 15, 8.35, 13.63, 20.34,
17.217, 5.69, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 30,
17.217, 5, 5, 5, 7.93, 7.93, 7.93, 5, 7.71, 5, 17.217, 8.175,
8.175, 6.69, 6.69, 6.69, 10.875, 5.345, 5.345, 5.345, 5.345,
3.54, 10.755, 10.755, 10.755, 19.61, 20.145, 20.145, 20.145,
10.34, 5.35, 6.34, 10.34, 5.35, 10.34, 9.343, 9.343, 9.343, 9.343,
137.111, 17.217, 137.111, 17.217, 137.111, 137.111, 17.217, 46.055,
46.055, 17.217, 17.77, 20.54, 17.217, 17.217, 20.54, 11.75, 11.75,
17.217, 56.89, 20.54, 20.54, 55, 75, 75, 75, 19.61, 19.61, 19.61,
19.61, 19.61, 25.508, 20.698, 5.42, 5.42, 19.61, 19.61, 19.61,
20.698, 5.42, 25.508, 5.42, 17.217, 16.92, 9.343, 9.343, 19.61,
9.343, 19.61, 19.61, 16.92, 16.92, 16.92, 5.42, 5.42, 19.61,
5.42, 9.343, 19.61, 5.42, 5.42, 19.61, 9.343, 46.055, 17.217,
5.42, 19.61, 5.69, 19.61, 17.217, 9.343, 5.42, 19.61, 19.61,
17.217, 35.508, 5.42, 5.42, 5.42, 19.61, 17.217, 5.69, 19.61,
19.61, 8.35, 17.217, 5.69, 19.61, 5.69, 5.69, 17.217, 19.61,
16.92, 19.61, 17.217, 9.343, 5.42, 17.217, 17.217, 9.343, 33.456,
17.217, 46.055, 137.11, 56.89, 54.25, 17.217, 56.89, 55, 17.217,
46.055, 20.698, 46.055, 54.25, 27.25, 53.5, 5.69, 53.5, 20.698,
20.698, 5.69, 17.217, 17.217, 17.217, 17.217, 17.217, 5.69, 17.217,
17.217, 17.217, 11.75, 17.217, 5.69, 11.75, 11.75, 5.42, 5.42,
17.217, 5.69, 20.54, 11.75, 5.69, 11.75, 17.217, 121.505, 121.505,
121.505, 11.75, 5.69, 11.75, 11.75, 20.698, 121.505, 17.217,
11.75, 17.217, 5.42, 17.217, 17.217, 5.42, 5.69, 121.505, 17.217,
5.42, 46.055, 5.69, 20.698, 46.055, 15.86, 15.86, 5.69, 5.69,
11.75, 5.42, 46.055, 5.42, 10.349, 121.505, 15.86, 25.508, 17.217,
17.217, 11.75, 17.217, 17.217, 17.217, 17.217, 17.217, 17.217,
46.055, 17.217, 17.217, 17.217, 17.217, 17.217, 5.42, 17.217,
17.217, 17.217, 9.343, 85.271, 46.055, 17.217, 17.217, 46.055,
25.508, 25.508, 25.508, 20.698, 19.61, 11.75, 11.75, 11.75, 20.698,
11.75, 11.75, 20.45, 11.75, 20.45, 20.45, 7.5, 11.75, 11.75,
20.45), taken = c(10L, 2L, 12L, 1L, 4L, 1L, 10L, 3L, 5L, 5L,
1L, 5L, 3L, 3L, 5L, 2L, 4L, 2L, 1L, 4L, 1L, 1L, 2L, 1L, 3L, 1L,
3L, 2L, 3L, 2L, 2L, 1L, 4L, 1L, 5L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 3L, 1L, 3L, 2L, 2L, 3L, 2L, 3L, 1L, 2L, 3L, 2L,
2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 2L, 1L,
4L, 2L, 3L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L,
2L, 1L, 2L, 3L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, 3L, 4L, 2L,
1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 3L, 2L, 1L, 1L, 5L, 1L, 1L,
5L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 1L, 1L, 1L, 3L, 1L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 2L,
1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L,
1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 1L,
1L, 3L, 1L, 1L, 1L, 2L, 1L, 2L, 3L, 1L, 2L, 1L, 2L, 2L, 2L, 2L,
2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 3L, 2L, 1L, 2L, 2L, 2L, 1L, 3L, 2L,
2L, 3L, 3L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 2L, 4L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 3L, 3L, 2L, 2L, 4L,
1L, 2L, 4L, 3L, 3L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 1L,
1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
3L, 2L, 2L, 2L, 3L, 3L, 2L, 2L, 2L, 3L, 2L, 2L, 2L, 2L, 2L, 2L,
3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 3L, 2L, 3L, 2L, 3L, 1L, 2L,
2L, 1L, 2L, 3L, 2L, 2L, 3L, 2L, 1L, 2L, 1L, 2L), mass = c(13.8758,
2.77516, 16.65096, 1.38758, 5.55032, 1.38758, 13.8758, 4.16274,
6.9379, 6.9379, 1.38758, 6.9379, 4.16274, 4.16274, 6.9379, 2.77516,
5.55032, 2.77516, 1.38758, 5.55032, 1.38758, 1.38758, 2.77516,
1.38758, 4.16274, 1.38758, 4.16274, 2.77516, 4.16274, 2.77516,
2.77516, 1.38758, 5.55032, 1.38758, 6.9379, 2.77516, 2.77516,
5.49268, 5.49268, 5.49268, 2.74634, 2.74634, 2.74634, 2.74634,
2.74634, 2.74634, 8.23902, 1.295, 3.885, 2.59, 2.59, 3.885, 2.59,
3.885, 1.295, 2.59, 3.885, 2.59, 2.59, 2.59, 1.295, 2.59, 2.59,
2.59, 1.295, 1.295, 1.295, 1.295, 1.295, 1.295, 2.59, 3.885,
2.59, 1.295, 5.18, 2.59, 3.885, 2.59, 1.295, 2.59, 1.295, 1.295,
1.295, 2.59, 2.59, 2.59, 2.59, 1.295, 2.59, 2.59, 2.59, 1.295,
2.59, 3.885, 1.295, 1.295, 1.295, 3.885, 1.295, 1.295, 1.295,
1.295, 1.295, 3.885, 5.18, 2.59, 1.295, 1.4429, 2.8858, 1.4429,
1.4429, 2.8858, 1.4429, 1.4429, 2.8858, 4.3287, 2.8858, 1.4429,
1.4429, 7.2145, 1.4429, 1.4429, 7.2145, 1.4429, 1.4429, 1.4429,
1.4429, 1.4429, 2.8858, 4.3287, 1.4429, 1.4429, 1.4429, 4.3287,
1.4429, 2.8858, 2.8858, 2.8858, 1.4429, 2.877, 2.877, 2.877,
2.877, 5.754, 2.877, 2.877, 2.877, 2.877, 2.877, 2.877, 8.631,
2.877, 2.877, 5.754, 2.877, 2.877, 2.877, 2.877, 5.754, 2.877,
2.877, 2.877, 2.877, 5.754, 5.754, 2.877, 2.877, 2.877, 2.877,
2.877, 2.877, 1.3719, 1.3719, 1.3719, 1.3719, 1.3719, 2.7438,
1.3719, 2.7438, 1.3719, 1.3719, 2.7438, 2.7438, 1.3719, 1.3719,
1.3719, 1.3719, 2.7438, 1.3719, 1.3719, 2.7438, 1.3719, 2.7438,
1.3719, 1.3719, 1.3719, 1.3719, 1.3719, 4.1157, 1.3719, 1.3719,
1.3719, 1.3719, 4.1157, 1.3719, 1.3719, 1.3719, 2.7438, 1.3719,
2.7438, 4.1157, 1.3719, 2.7438, 1.3719, 2.7438, 2.7438, 2.7438,
2.7438, 2.7438, 1.3719, 1.3719, 2.7438, 2.7438, 1.3719, 1.3719,
1.3719, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316,
2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316,
2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 1.17028, 2.34056,
1.17028, 1.17028, 1.17028, 1.17028, 1.17028, 1.17028, 2.34056,
2.34056, 3.51084, 2.34056, 1.17028, 2.34056, 2.34056, 2.34056,
1.17028, 3.51084, 2.34056, 2.34056, 3.51084, 3.51084, 2.34056,
1.17028, 2.03456, 4.06912, 2.03456, 2.03456, 2.03456, 4.06912,
2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456,
2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456,
2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456,
2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456,
2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456,
2.03456, 2.03456, 2.03456, 2.03456, 4.06912, 2.03456, 2.03456,
2.03456, 2.03456, 2.03456, 2.03456, 1.189102041, 1.189102041,
1.189102041, 2.378204082, 3.567306123, 2.378204082, 4.756408164,
2.378204082, 2.378204082, 2.378204082, 1.189102041, 1.189102041,
1.189102041, 1.189102041, 2.378204082, 1.189102041, 2.378204082,
1.189102041, 3.567306123, 3.567306123, 2.378204082, 2.378204082,
4.756408164, 1.189102041, 2.378204082, 4.756408164, 3.567306123,
3.567306123, 4.756408164, 2.45232, 2.45232, 2.45232, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 4.90464, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232,
2.45232, 4.90464, 4.90464, 2.45232, 2.45232, 2.45232, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 4.90464, 2.45232, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232,
7.35696, 2.45232, 2.45232, 2.45232, 4.90464, 2.45232, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 4.90464,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 4.90464, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232,
2.64516, 1.76344, 1.76344, 1.76344, 2.64516, 2.64516, 1.76344,
1.76344, 1.76344, 2.64516, 1.76344, 1.76344, 1.76344, 1.76344,
1.76344, 1.76344, 2.64516, 2.64516, 0.88172, 0.88172, 2.46532,
2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532,
2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532,
2.46532, 2.46532, 4.93064, 2.46532, 2.46532, 2.46532, 2.46532,
2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532,
2.46532, 4.93064, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532,
4.93064, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532,
2.46532, 2.46532, 4.93064, 2.46532, 2.46532, 4.93064, 2.46532,
2.46532, 2.46532, 4.93064, 4.93064, 2.46532, 4.93064, 2.46532,
2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.28446939,
2.28446939, 2.28446939, 2.28446939, 2.28446939, 2.28446939, 2.28446939,
2.28446939, 2.28446939, 2.28446939, 2.28446939, 2.28446939, 2.28446939,
2.28446939, 2.28446939, 2.28446939, 2.28446939, 2.28446939, 2.28446939,
2.28446939, 4.56893878, 4.568938778, 6.853408167, 4.56893878,
6.85340817, 4.56893878, 6.85340817, 2.28446939, 4.56893878, 4.568938778,
2.284469389, 4.56893878, 6.85340817, 4.56893878, 4.56893878,
6.85340817, 4.56893878, 2.28446939, 4.56893878, 2.28446939, 4.56893878
)), .Names = c("year", "size", "distance", "taken", "mass"), class =
"data.frame", row.names = c(NA,
-550L))
-----Original Message-----
From: Ben Bolker [mailto:bbolker at gmail.com]
Sent: Friday, April 08, 2011 4:15 PM
To: Stratford, Jeffrey
Cc: r-sig-ecology at r-project.org
Subject: Re: [R-sig-eco] zero-truncated Poisson
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On 04/08/2011 04:10 PM, Stratford, Jeffrey wrote:
Hi everyone, I'm attempting to run a zero-truncated Poisson using VGAM and I run
into
this error:
Error in if ((temp <- sum(wz[, 1:M, drop = FALSE] < wzepsilon)))
warning(paste(temp, :
argument is not interpretable as logical
Here's the code I'm using
blja <- read.csv("g:\\blja\\blja.csv", header=T)
library(VGAM)
glm1 <- vglm(blja$taken ~ blja$mass + blja$year,family=pospoisson,
data=blja)
Just a *very* quick suggestion: try glm1 <- vglm(taken ~ mass + year,family=pospoisson, data=blja) and see if that helps. Also try summary() on your data (and/or sapply(blja,class)) to see if your data types are what you think they are (hint: R often converts numeric data to factors if it has glitches in it) If neither of those helps, could you use dput() to post your data (or a subset thereof) in an easier-to-use format?
We're looking at the number of acorns taken by jays using year (two years) and how much (mass) they're taking. Our data looks like this: year size distance taken mass one small 30.8735 10 13.8758 one small 121.505 2 2.77516 one small 46.055 12 16.65096 one small 46.055 1 1.38758 one small 46.055 4 5.55032 one small 9.343 1 1.38758 one small 46.055 10 13.8758 one small 46.055 3 4.16274 one small 46.055 5 6.9379 one small 85.271 5 6.9379 one small 85.271 1 1.38758 one small 85.271 5 6.9379 one small 85.271 3 4.16274 ... Any suggestions (and alternatives!) would be appreciated! Thanks a bunch, Jeff ************************************ Dr. Jeffrey A. Stratford Department of Health and Biological Sciences 84 W. South Street Wilkes University, PA 18766 jeffrey.stratford at wilkes.edu 570-408-4761 (office) 570-332-2942 (cell) http://web.wilkes.edu/jeffrey.stratford/ ************************************ [[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
-----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.10 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/ iEYEARECAAYFAk2fbNkACgkQc5UpGjwzenMcAQCfcxwQ51AV3QspEu2goI2ECMnF r1AAn0/PgmgWqvzxcr0ywW1+acUL5q5b =ek5d -----END PGP SIGNATURE-----
Jeff,
The zero truncated Poisson distribution can be described as a
probability function conditional on Y>0:
Pr(Y=y|Y>0) = Pr(Y=y) / (1-Pr(Y=0)), y=1,2,3,...
This leads to this function 'nll' for the negative log-likelihood:
nll <- function(theta, y, X) {
lambda <- exp(drop(X %*% theta))
-sum(dpois(y, lambda, log=TRUE) - log(1-exp(-lambda)))
}
where theta is the vector of parameters, y is the observation vector,
X is the design matrix.
A simple simulation is useful to check that the theory works. I used
'optim' with the default Nelder-Mead simplex algorithm, initial values
are based in non-truncated Poisson estimates to speed up convergence:
theta <- c(1, 2)
X <- model.matrix(~rnorm(100))
y <- rpois(100,exp(X %*% theta))
X <- X[y>0,]
y <- y[y>0]
inits <- coef(glm(y ~ X-1, family=poisson))
res <- optim(par=inits, fn=nll, y=y, X=X)
cbind(theta=theta, theta.hat=res$par)
Here I used the data you provided, note that 'hessian=TRUE' for some
subsequent SE calculations:
inits <- coef(glm(taken ~ mass + year, blja, family=poisson))
y <- blja$taken
X <- model.matrix(~ mass + year, blja)
res <- optim(par=inits, fn=nll, y=y, X=X, hessian=TRUE)
Vector of parameters that minimize 'nll' (maximize the log-likelihood)
given the data:
cf <- res$par
Standard error (note that due to minimization, I used the Hessian
instead of -H):
se <- sqrt(diag(solve(res$hessian)))
z and p values and the final output table:
zstat <- cf/se
pval <- 2 * pnorm(-abs(zstat))
out <- cbind(cf, se, zstat, pval)
colnames(out) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)")
printCoefmat(out, digits = 4, signif.legend = FALSE)
For your data set, this gives:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.76392 0.10659 -7.167 7.66e-13 ***
mass 0.24008 0.01262 19.028 < 2e-16 ***
yeartwo -0.29006 0.11974 -2.423 0.0154 *
Cheers,
Peter
P?ter S?lymos
Alberta Biodiversity Monitoring Institute
and Boreal Avian Modelling project
Department of Biological Sciences
CW 405, Biological Sciences Bldg
University of Alberta
Edmonton, Alberta, T6G 2E9, Canada
Phone: 780.492.8534
Fax: 780.492.7635
email <- paste("solymos", "ualberta.ca", sep = "@")
http://www.abmi.ca
http://www.borealbirds.ca
http://sites.google.com/site/psolymos
On Fri, Apr 8, 2011 at 6:40 PM, Stratford, Jeffrey
<jeffrey.stratford at wilkes.edu> wrote:
Thanks Ben,
That didn't do the trick. I checked all the classes
? ?year ? ? ?size ?distance ? ? taken ? ? ?mass
?"factor" ?"factor" "numeric" "integer" "numeric"
And used the notation that you suggested and I did get the same error.
Just so it's not such a black box for me, what makes vector GLMs
different from "ordinary" GLMs?
Thanks,
Jeff
Here's all the data.
structure(list(year = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label =
c("one",
"two"), class = "factor"), size = structure(c(2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L), .Label = c("large", "small"), class = "factor"), distance =
c(30.8735,
121.505, 46.055, 46.055, 46.055, 9.343, 46.055, 46.055, 46.055,
85.271, 85.271, 85.271, 85.271, 85.271, 85.271, 30.8735, 30.8735,
9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 152.717, 152.717,
152.717, 152.717, 152.717, 152.717, 152.717, 46.055, 46.055,
152.717, 152.717, 20.698, 20.698, 17.217, 17.217, 9.343, 9.343,
17.217, 46.055, 152.717, 46.055, 17.217, 152.717, 17.217, 17.217,
46.055, 30.8735, 30.8735, 30.8735, 5.69, 30.8735, 17.217, 30.8735,
30.8735, 30.8735, 30.8735, 30.8735, 9.343, 9.343, 17.217, 17.217,
17.217, 17.217, 17.217, 17.217, 9.343, 9.343, 17.217, 17.217,
17.217, 20.698, 17.217, 17.217, 17.217, 5.42, 17.217, 17.217,
17.217, 9.343, 30.8735, 30.8735, 30.8735, 9.343, 9.343, 9.343,
9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343,
17.217, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343,
30.8735, 30.8735, 17.217, 17.217, 46.055, 36.239, 17.217, 17.217,
17.217, 46.055, 30.8735, 30.8735, 17.217, 17.217, 17.217, 121.505,
152.717, 152.717, 17.217, 152.717, 121.505, 121.505, 121.505,
121.505, 121.505, 9.343, 121.505, 9.343, 121.505, 121.505, 30.8735,
121.505, 17.217, 17.217, 17.217, 9.343, 30.8735, 85.271, 85.271,
85.271, 85.271, 85.271, 85.271, 9.343, 85.271, 85.271, 85.271,
85.271, 20.698, 9.343, 30.8735, 17.217, 20.698, 30.8735, 17.217,
85.271, 121.505, 121.505, 85.271, 85.271, 85.271, 85.271, 121.505,
121.505, NA, 121.505, 9.343, 9.343, 9.343, 9.343, 9.343, 9.343,
9.343, 30.8735, 17.217, 9.343, 9.343, 9.343, 30.8735, 30.8735,
30.8735, 9.343, 9.343, 30.8735, 17.217, 5.42, 17.217, 85.271,
85.271, 30.8735, 30.8735, 30.8735, 30.8735, 30.8735, 30.8735,
17.217, 85.271, 17.217, 30.8735, 30.8735, 85.271, 30.8735, 30.8735,
85.271, 30.8735, 30.8735, 30.8735, 30.8735, 30.8735, 30.8735,
17.217, 30.8735, 9.343, 30.8735, 30.8735, 30.8735, 30.8735, 9.343,
30.8735, 30.8735, 17.217, 9.343, 9.343, 9.343, 85.271, 30.8735,
30.8735, 46.055, 5.69, 85.271, 85.271, 17.217, 46.055, 85.271,
30.8735, 85.271, 46.055, 121.505, 121.505, 121.505, 17.217, 17.217,
25.508, 9.343, 17.217, 17.217, 9.343, 17.217, 9.343, 17.217,
34.35, 34.35, 46.055, 46.055, 9.343, 20.698, 17.217, 27.25, 20.54,
27.25, 20.698, 17.217, 20.698, 20.698, 15, 15, 8.35, 13.63, 20.34,
17.217, 5.69, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 30,
17.217, 5, 5, 5, 7.93, 7.93, 7.93, 5, 7.71, 5, 17.217, 8.175,
8.175, 6.69, 6.69, 6.69, 10.875, 5.345, 5.345, 5.345, 5.345,
3.54, 10.755, 10.755, 10.755, 19.61, 20.145, 20.145, 20.145,
10.34, 5.35, 6.34, 10.34, 5.35, 10.34, 9.343, 9.343, 9.343, 9.343,
137.111, 17.217, 137.111, 17.217, 137.111, 137.111, 17.217, 46.055,
46.055, 17.217, 17.77, 20.54, 17.217, 17.217, 20.54, 11.75, 11.75,
17.217, 56.89, 20.54, 20.54, 55, 75, 75, 75, 19.61, 19.61, 19.61,
19.61, 19.61, 25.508, 20.698, 5.42, 5.42, 19.61, 19.61, 19.61,
20.698, 5.42, 25.508, 5.42, 17.217, 16.92, 9.343, 9.343, 19.61,
9.343, 19.61, 19.61, 16.92, 16.92, 16.92, 5.42, 5.42, 19.61,
5.42, 9.343, 19.61, 5.42, 5.42, 19.61, 9.343, 46.055, 17.217,
5.42, 19.61, 5.69, 19.61, 17.217, 9.343, 5.42, 19.61, 19.61,
17.217, 35.508, 5.42, 5.42, 5.42, 19.61, 17.217, 5.69, 19.61,
19.61, 8.35, 17.217, 5.69, 19.61, 5.69, 5.69, 17.217, 19.61,
16.92, 19.61, 17.217, 9.343, 5.42, 17.217, 17.217, 9.343, 33.456,
17.217, 46.055, 137.11, 56.89, 54.25, 17.217, 56.89, 55, 17.217,
46.055, 20.698, 46.055, 54.25, 27.25, 53.5, 5.69, 53.5, 20.698,
20.698, 5.69, 17.217, 17.217, 17.217, 17.217, 17.217, 5.69, 17.217,
17.217, 17.217, 11.75, 17.217, 5.69, 11.75, 11.75, 5.42, 5.42,
17.217, 5.69, 20.54, 11.75, 5.69, 11.75, 17.217, 121.505, 121.505,
121.505, 11.75, 5.69, 11.75, 11.75, 20.698, 121.505, 17.217,
11.75, 17.217, 5.42, 17.217, 17.217, 5.42, 5.69, 121.505, 17.217,
5.42, 46.055, 5.69, 20.698, 46.055, 15.86, 15.86, 5.69, 5.69,
11.75, 5.42, 46.055, 5.42, 10.349, 121.505, 15.86, 25.508, 17.217,
17.217, 11.75, 17.217, 17.217, 17.217, 17.217, 17.217, 17.217,
46.055, 17.217, 17.217, 17.217, 17.217, 17.217, 5.42, 17.217,
17.217, 17.217, 9.343, 85.271, 46.055, 17.217, 17.217, 46.055,
25.508, 25.508, 25.508, 20.698, 19.61, 11.75, 11.75, 11.75, 20.698,
11.75, 11.75, 20.45, 11.75, 20.45, 20.45, 7.5, 11.75, 11.75,
20.45), taken = c(10L, 2L, 12L, 1L, 4L, 1L, 10L, 3L, 5L, 5L,
1L, 5L, 3L, 3L, 5L, 2L, 4L, 2L, 1L, 4L, 1L, 1L, 2L, 1L, 3L, 1L,
3L, 2L, 3L, 2L, 2L, 1L, 4L, 1L, 5L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 3L, 1L, 3L, 2L, 2L, 3L, 2L, 3L, 1L, 2L, 3L, 2L,
2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 2L, 1L,
4L, 2L, 3L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L,
2L, 1L, 2L, 3L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, 3L, 4L, 2L,
1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 3L, 2L, 1L, 1L, 5L, 1L, 1L,
5L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 1L, 1L, 1L, 3L, 1L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 2L,
1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L,
1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 1L,
1L, 3L, 1L, 1L, 1L, 2L, 1L, 2L, 3L, 1L, 2L, 1L, 2L, 2L, 2L, 2L,
2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 3L, 2L, 1L, 2L, 2L, 2L, 1L, 3L, 2L,
2L, 3L, 3L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 2L, 4L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 3L, 3L, 2L, 2L, 4L,
1L, 2L, 4L, 3L, 3L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 1L,
1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
3L, 2L, 2L, 2L, 3L, 3L, 2L, 2L, 2L, 3L, 2L, 2L, 2L, 2L, 2L, 2L,
3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 3L, 2L, 3L, 2L, 3L, 1L, 2L,
2L, 1L, 2L, 3L, 2L, 2L, 3L, 2L, 1L, 2L, 1L, 2L), mass = c(13.8758,
2.77516, 16.65096, 1.38758, 5.55032, 1.38758, 13.8758, 4.16274,
6.9379, 6.9379, 1.38758, 6.9379, 4.16274, 4.16274, 6.9379, 2.77516,
5.55032, 2.77516, 1.38758, 5.55032, 1.38758, 1.38758, 2.77516,
1.38758, 4.16274, 1.38758, 4.16274, 2.77516, 4.16274, 2.77516,
2.77516, 1.38758, 5.55032, 1.38758, 6.9379, 2.77516, 2.77516,
5.49268, 5.49268, 5.49268, 2.74634, 2.74634, 2.74634, 2.74634,
2.74634, 2.74634, 8.23902, 1.295, 3.885, 2.59, 2.59, 3.885, 2.59,
3.885, 1.295, 2.59, 3.885, 2.59, 2.59, 2.59, 1.295, 2.59, 2.59,
2.59, 1.295, 1.295, 1.295, 1.295, 1.295, 1.295, 2.59, 3.885,
2.59, 1.295, 5.18, 2.59, 3.885, 2.59, 1.295, 2.59, 1.295, 1.295,
1.295, 2.59, 2.59, 2.59, 2.59, 1.295, 2.59, 2.59, 2.59, 1.295,
2.59, 3.885, 1.295, 1.295, 1.295, 3.885, 1.295, 1.295, 1.295,
1.295, 1.295, 3.885, 5.18, 2.59, 1.295, 1.4429, 2.8858, 1.4429,
1.4429, 2.8858, 1.4429, 1.4429, 2.8858, 4.3287, 2.8858, 1.4429,
1.4429, 7.2145, 1.4429, 1.4429, 7.2145, 1.4429, 1.4429, 1.4429,
1.4429, 1.4429, 2.8858, 4.3287, 1.4429, 1.4429, 1.4429, 4.3287,
1.4429, 2.8858, 2.8858, 2.8858, 1.4429, 2.877, 2.877, 2.877,
2.877, 5.754, 2.877, 2.877, 2.877, 2.877, 2.877, 2.877, 8.631,
2.877, 2.877, 5.754, 2.877, 2.877, 2.877, 2.877, 5.754, 2.877,
2.877, 2.877, 2.877, 5.754, 5.754, 2.877, 2.877, 2.877, 2.877,
2.877, 2.877, 1.3719, 1.3719, 1.3719, 1.3719, 1.3719, 2.7438,
1.3719, 2.7438, 1.3719, 1.3719, 2.7438, 2.7438, 1.3719, 1.3719,
1.3719, 1.3719, 2.7438, 1.3719, 1.3719, 2.7438, 1.3719, 2.7438,
1.3719, 1.3719, 1.3719, 1.3719, 1.3719, 4.1157, 1.3719, 1.3719,
1.3719, 1.3719, 4.1157, 1.3719, 1.3719, 1.3719, 2.7438, 1.3719,
2.7438, 4.1157, 1.3719, 2.7438, 1.3719, 2.7438, 2.7438, 2.7438,
2.7438, 2.7438, 1.3719, 1.3719, 2.7438, 2.7438, 1.3719, 1.3719,
1.3719, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316,
2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316,
2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 2.8316, 1.17028, 2.34056,
1.17028, 1.17028, 1.17028, 1.17028, 1.17028, 1.17028, 2.34056,
2.34056, 3.51084, 2.34056, 1.17028, 2.34056, 2.34056, 2.34056,
1.17028, 3.51084, 2.34056, 2.34056, 3.51084, 3.51084, 2.34056,
1.17028, 2.03456, 4.06912, 2.03456, 2.03456, 2.03456, 4.06912,
2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456,
2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456,
2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456,
2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456,
2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456, 2.03456,
2.03456, 2.03456, 2.03456, 2.03456, 4.06912, 2.03456, 2.03456,
2.03456, 2.03456, 2.03456, 2.03456, 1.189102041, 1.189102041,
1.189102041, 2.378204082, 3.567306123, 2.378204082, 4.756408164,
2.378204082, 2.378204082, 2.378204082, 1.189102041, 1.189102041,
1.189102041, 1.189102041, 2.378204082, 1.189102041, 2.378204082,
1.189102041, 3.567306123, 3.567306123, 2.378204082, 2.378204082,
4.756408164, 1.189102041, 2.378204082, 4.756408164, 3.567306123,
3.567306123, 4.756408164, 2.45232, 2.45232, 2.45232, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 4.90464, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232,
2.45232, 4.90464, 4.90464, 2.45232, 2.45232, 2.45232, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 4.90464, 2.45232, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232,
7.35696, 2.45232, 2.45232, 2.45232, 4.90464, 2.45232, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 4.90464,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 4.90464, 2.45232,
2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232, 2.45232,
2.64516, 1.76344, 1.76344, 1.76344, 2.64516, 2.64516, 1.76344,
1.76344, 1.76344, 2.64516, 1.76344, 1.76344, 1.76344, 1.76344,
1.76344, 1.76344, 2.64516, 2.64516, 0.88172, 0.88172, 2.46532,
2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532,
2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532,
2.46532, 2.46532, 4.93064, 2.46532, 2.46532, 2.46532, 2.46532,
2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532,
2.46532, 4.93064, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532,
4.93064, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532,
2.46532, 2.46532, 4.93064, 2.46532, 2.46532, 4.93064, 2.46532,
2.46532, 2.46532, 4.93064, 4.93064, 2.46532, 4.93064, 2.46532,
2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.46532, 2.28446939,
2.28446939, 2.28446939, 2.28446939, 2.28446939, 2.28446939, 2.28446939,
2.28446939, 2.28446939, 2.28446939, 2.28446939, 2.28446939, 2.28446939,
2.28446939, 2.28446939, 2.28446939, 2.28446939, 2.28446939, 2.28446939,
2.28446939, 4.56893878, 4.568938778, 6.853408167, 4.56893878,
6.85340817, 4.56893878, 6.85340817, 2.28446939, 4.56893878, 4.568938778,
2.284469389, 4.56893878, 6.85340817, 4.56893878, 4.56893878,
6.85340817, 4.56893878, 2.28446939, 4.56893878, 2.28446939, 4.56893878
)), .Names = c("year", "size", "distance", "taken", "mass"), class =
"data.frame", row.names = c(NA,
-550L))
-----Original Message-----
From: Ben Bolker [mailto:bbolker at gmail.com]
Sent: Friday, April 08, 2011 4:15 PM
To: Stratford, Jeffrey
Cc: r-sig-ecology at r-project.org
Subject: Re: [R-sig-eco] zero-truncated Poisson
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On 04/08/2011 04:10 PM, Stratford, Jeffrey wrote:
Hi everyone, I'm attempting to run a zero-truncated Poisson using VGAM and I run
into
this error:
Error in if ((temp <- sum(wz[, 1:M, drop = FALSE] < wzepsilon)))
warning(paste(temp, ?:
? argument is not interpretable as logical
Here's the code I'm using
blja <- read.csv("g:\\blja\\blja.csv", header=T)
library(VGAM)
glm1 <- vglm(blja$taken ~ blja$mass + blja$year,family=pospoisson,
data=blja)
?Just a *very* quick suggestion: try glm1 <- vglm(taken ~ mass + year,family=pospoisson, data=blja) ?and see if that helps. ?Also try summary() on your data (and/or sapply(blja,class)) to see if your data types are what you think they are (hint: R often converts numeric data to factors if it has glitches in it) ?If neither of those helps, could you use dput() to post your data (or a subset thereof) in an easier-to-use format?
We're looking at the number of acorns taken by jays using year (two years) and how much (mass) they're taking. Our data looks like this: year ?size ? ?distance ? ? ? ?taken ? mass one ? small ? 30.8735 10 ? ? ?13.8758 one ? small ? 121.505 2 ? ? ? 2.77516 one ? small ? 46.055 ?12 ? ? ?16.65096 one ? small ? 46.055 ?1 ? ? ? 1.38758 one ? small ? 46.055 ?4 ? ? ? 5.55032 one ? small ? 9.343 ? 1 ? ? ? 1.38758 one ? small ? 46.055 ?10 ? ? ?13.8758 one ? small ? 46.055 ?3 ? ? ? 4.16274 one ? small ? 46.055 ?5 ? ? ? 6.9379 one ? small ? 85.271 ?5 ? ? ? 6.9379 one ? small ? 85.271 ?1 ? ? ? 1.38758 one ? small ? 85.271 ?5 ? ? ? 6.9379 one ? small ? 85.271 ?3 ? ? ? 4.16274 ... Any suggestions (and alternatives!) would be appreciated! Thanks a bunch, Jeff ************************************ Dr. Jeffrey A. Stratford Department of Health and Biological Sciences 84 W. South Street Wilkes University, PA 18766 jeffrey.stratford at wilkes.edu 570-408-4761 (office) 570-332-2942 (cell) http://web.wilkes.edu/jeffrey.stratford/ ************************************ ? ? ? [[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
-----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.10 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/ iEYEARECAAYFAk2fbNkACgkQc5UpGjwzenMcAQCfcxwQ51AV3QspEu2goI2ECMnF r1AAn0/PgmgWqvzxcr0ywW1+acUL5q5b =ek5d -----END PGP SIGNATURE-----
_______________________________________________ R-sig-ecology mailing list R-sig-ecology at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
1 day later
On 11-04-09 03:07 AM, Peter Solymos wrote:
Jeff,
The zero truncated Poisson distribution can be described as a
probability function conditional on Y>0:
Pr(Y=y|Y>0) = Pr(Y=y) / (1-Pr(Y=0)), y=1,2,3,...
This leads to this function 'nll' for the negative log-likelihood:
nll <- function(theta, y, X) {
lambda <- exp(drop(X %*% theta))
-sum(dpois(y, lambda, log=TRUE) - log(1-exp(-lambda)))
}
where theta is the vector of parameters, y is the observation vector,
X is the design matrix.
A simple simulation is useful to check that the theory works. I used
'optim' with the default Nelder-Mead simplex algorithm, initial values
are based in non-truncated Poisson estimates to speed up convergence:
theta <- c(1, 2)
X <- model.matrix(~rnorm(100))
y <- rpois(100,exp(X %*% theta))
X <- X[y>0,]
y <- y[y>0]
inits <- coef(glm(y ~ X-1, family=poisson))
res <- optim(par=inits, fn=nll, y=y, X=X)
cbind(theta=theta, theta.hat=res$par)
Here I used the data you provided, note that 'hessian=TRUE' for some
subsequent SE calculations:
inits <- coef(glm(taken ~ mass + year, blja, family=poisson))
y <- blja$taken
X <- model.matrix(~ mass + year, blja)
res <- optim(par=inits, fn=nll, y=y, X=X, hessian=TRUE)
Vector of parameters that minimize 'nll' (maximize the log-likelihood)
given the data:
cf <- res$par
Standard error (note that due to minimization, I used the Hessian
instead of -H):
se <- sqrt(diag(solve(res$hessian)))
z and p values and the final output table:
zstat <- cf/se
pval <- 2 * pnorm(-abs(zstat))
out <- cbind(cf, se, zstat, pval)
colnames(out) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)")
printCoefmat(out, digits = 4, signif.legend = FALSE)
For your data set, this gives:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.76392 0.10659 -7.167 7.66e-13 ***
mass 0.24008 0.01262 19.028 < 2e-16 ***
yeartwo -0.29006 0.11974 -2.423 0.0154 *
Cheers,
Peter
P?ter S?lymos
Alberta Biodiversity Monitoring Institute
and Boreal Avian Modelling project
Department of Biological Sciences
CW 405, Biological Sciences Bldg
University of Alberta
Edmonton, Alberta, T6G 2E9, Canada
Phone: 780.492.8534
Fax: 780.492.7635
email <- paste("solymos", "ualberta.ca", sep = "@")
http://www.abmi.ca
http://www.borealbirds.ca
http://sites.google.com/site/psolymos
I agree. I would just like to point out that, if you define the
truncated Poisson distribution function appropriately, you can then do
this very compactly via the mle2 function in the bbmle package.
dtruncpois <- function(x,lambda,log=FALSE) {
r <- ifelse(x==0,-Inf,dpois(x,lambda,log=TRUE)-
log(1-dpois(0,lambda,log=FALSE)))
if (log) r else exp(r)
}
library(bbmle)
(m1 <- mle2(taken~dtruncpois(exp(logmu)),
parameters=list(logmu~mass+year),
data=blja,
start=list(logmu=1)))
coef(m1)
summary(m1) ## gives a table almost equivalent to Peter's
confint(m1) ## profile confidence intervals
After a bit of digging it seems that VGAM::vglm fails when it hits an
NaN -- I haven't figured out more than that it gets this on the *second*
step, not the first, of the optimization. When it works VGAM should be
faster than the sort of generalized ML estimation that Peter and I are
suggesting, but I have fairly often found it to be not very robust when
I try it.
If you plot the data (with a regular Poisson GLM superimposed):
library(ggplot2)
ggplot(blja,aes(x=mass,y=taken,colour=year))+
stat_sum(alpha=0.5)+geom_smooth(method="glm",family="poisson")+
theme_bw()
You see that there are two observations with very large masses.
It turns out that removing these two observations allows vglm() to work:
glm1 <- vglm(taken ~ mass + year,family=pospoisson,
data=blja,subset=mass<10)
m2 <- mle2(taken~dtruncpois(exp(logmu)),
parameters=list(logmu~mass+year),
data=subset(blja,mass<10),
start=list(logmu=1))
It also changes the estimate of the mass effect considerably:
confint(m1)
confint(m2)
Ben Bolker
Thanks Ben and Peter for your thoughtful replies. This gives me quite a bit to chew on and digest but I think I can analyze these data. -----Original Message----- From: Ben Bolker [mailto:bbolker at gmail.com] Sent: Sunday, April 10, 2011 10:38 AM To: Peter Solymos Cc: Stratford, Jeffrey; r-sig-ecology at r-project.org Subject: Re: [R-sig-eco] zero-truncated Poisson
On 11-04-09 03:07 AM, Peter Solymos wrote:
Jeff,
The zero truncated Poisson distribution can be described as a
probability function conditional on Y>0:
Pr(Y=y|Y>0) = Pr(Y=y) / (1-Pr(Y=0)), y=1,2,3,...
This leads to this function 'nll' for the negative log-likelihood:
nll <- function(theta, y, X) {
lambda <- exp(drop(X %*% theta))
-sum(dpois(y, lambda, log=TRUE) - log(1-exp(-lambda)))
}
where theta is the vector of parameters, y is the observation vector,
X is the design matrix.
A simple simulation is useful to check that the theory works. I used
'optim' with the default Nelder-Mead simplex algorithm, initial values
are based in non-truncated Poisson estimates to speed up convergence:
theta <- c(1, 2)
X <- model.matrix(~rnorm(100))
y <- rpois(100,exp(X %*% theta))
X <- X[y>0,]
y <- y[y>0]
inits <- coef(glm(y ~ X-1, family=poisson))
res <- optim(par=inits, fn=nll, y=y, X=X)
cbind(theta=theta, theta.hat=res$par)
Here I used the data you provided, note that 'hessian=TRUE' for some
subsequent SE calculations:
inits <- coef(glm(taken ~ mass + year, blja, family=poisson))
y <- blja$taken
X <- model.matrix(~ mass + year, blja)
res <- optim(par=inits, fn=nll, y=y, X=X, hessian=TRUE)
Vector of parameters that minimize 'nll' (maximize the log-likelihood)
given the data:
cf <- res$par
Standard error (note that due to minimization, I used the Hessian
instead of -H):
se <- sqrt(diag(solve(res$hessian)))
z and p values and the final output table:
zstat <- cf/se
pval <- 2 * pnorm(-abs(zstat))
out <- cbind(cf, se, zstat, pval)
colnames(out) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)")
printCoefmat(out, digits = 4, signif.legend = FALSE)
For your data set, this gives:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.76392 0.10659 -7.167 7.66e-13 ***
mass 0.24008 0.01262 19.028 < 2e-16 ***
yeartwo -0.29006 0.11974 -2.423 0.0154 *
Cheers,
Peter
P?ter S?lymos
Alberta Biodiversity Monitoring Institute
and Boreal Avian Modelling project
Department of Biological Sciences
CW 405, Biological Sciences Bldg
University of Alberta
Edmonton, Alberta, T6G 2E9, Canada
Phone: 780.492.8534
Fax: 780.492.7635
email <- paste("solymos", "ualberta.ca", sep = "@")
http://www.abmi.ca
http://www.borealbirds.ca
http://sites.google.com/site/psolymos
I agree. I would just like to point out that, if you define the
truncated Poisson distribution function appropriately, you can then do
this very compactly via the mle2 function in the bbmle package.
dtruncpois <- function(x,lambda,log=FALSE) {
r <- ifelse(x==0,-Inf,dpois(x,lambda,log=TRUE)-
log(1-dpois(0,lambda,log=FALSE)))
if (log) r else exp(r)
}
library(bbmle)
(m1 <- mle2(taken~dtruncpois(exp(logmu)),
parameters=list(logmu~mass+year),
data=blja,
start=list(logmu=1)))
coef(m1)
summary(m1) ## gives a table almost equivalent to Peter's
confint(m1) ## profile confidence intervals
After a bit of digging it seems that VGAM::vglm fails when it hits an
NaN -- I haven't figured out more than that it gets this on the *second*
step, not the first, of the optimization. When it works VGAM should be
faster than the sort of generalized ML estimation that Peter and I are
suggesting, but I have fairly often found it to be not very robust when
I try it.
If you plot the data (with a regular Poisson GLM superimposed):
library(ggplot2)
ggplot(blja,aes(x=mass,y=taken,colour=year))+
stat_sum(alpha=0.5)+geom_smooth(method="glm",family="poisson")+
theme_bw()
You see that there are two observations with very large masses.
It turns out that removing these two observations allows vglm() to work:
glm1 <- vglm(taken ~ mass + year,family=pospoisson,
data=blja,subset=mass<10)
m2 <- mle2(taken~dtruncpois(exp(logmu)),
parameters=list(logmu~mass+year),
data=subset(blja,mass<10),
start=list(logmu=1))
It also changes the estimate of the mass effect considerably:
confint(m1)
confint(m2)
Ben Bolker