Hello,
I?m attempting to run zero-inflated and hurdle models using the pscl()
package, but keep encountering problems with running my
fully-parameterized model and I?m not sure why.
My data consist of measurements of seed set (the response, which is
zero-inflated and overdispersed) for two treatments and a control (B,
F, and O) at three elevations (low, mid, and high).? At each elevation
I also measured average temperature (AvgTemp) and abundance of insects
(N).? My goal is to construct a model that best describes seed set
using the four predictors above (and interactions, if necessary).
I find that I have no problems running a partial model with only my
two categorical predictors as follows (just an example using hurdle(),
but results are the same for zeroinfl()):
H1B <- hurdle(AVInt ~ Elev + Treatment ?| ?Elev + Treatment, dist =
"negbin", link = "logit", data = erytseeds)
It?s when I try to add the two continuous predictors that things begin
to get messy;
H4B <- hurdle(AVInt ~ Elev*Treatment*AvgTemp*N | Elev*Treatment*AvgTemp*N)
This results in the error:
Warning message:
In sqrt(diag(object$vcov)) : NaNs produced
Error in optim(fn = loglikfun, gr = gradfun, par = c(start$count, start$zero,? :
non-finite value supplied by optim
Even when I run the model using one continuous predictor and one
categorical predictor, or both continuous predictors, I end up with
NAs in the summary table of the count, but not in the zero-inflated
portion of the model.? For example:
summary(H9B)
Call:
zeroinfl(AVInt ~ AvgTemp+N | AvgTemp+N, data = erytseeds, dist =
"negbin", link = "logit")
Pearson residuals:
Min????? 1Q? Median????? 3Q???? Max
-0.5807 -0.4566 -0.2791 -0.2791? 5.5813
Count model coefficients (negbin with log link):
???????????? Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.770120?? 0.452221? -8.337?? <2e-16 ***
AvgTemp????? 0.404550?? 0.173458?? 2.332?? 0.0197 *
N??????????? 0.001211???????? NA????? NA?????? NA
Log(theta)?? 1.308387???????? NA????? NA?????? NA
Zero-inflation model coefficients (binomial with logit link):
????????????? Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.814e+01? 8.115e+00? -5.932 2.99e-09 ***
AvgTemp????? 4.039e+00? 7.158e-01?? 5.642 1.68e-08 ***
N??????????? 3.010e-03? 4.863e-04?? 6.189 6.04e-10 ***
---
Theta = 3.7002
Number of iterations in BFGS optimization: 31
Log-likelihood: -456.4 on 7 Df
summary(H10B)
Call:
zeroinfl(AVInt ~ Treatment+N | Treatment+N, data = erytseeds, dist =
"negbin", link = "logit")
Pearson residuals:
??? Min????? 1Q? Median????? 3Q???? Max
-0.6608 -0.4684 -0.3565 -0.2956? 6.1730
Count model coefficients (negbin with log link):
???????????? Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.9555134???????? NA????? NA?????? NA
TreatmentF? 0.1407240? 0.1844811?? 0.763??? 0.446
TreatmentO? 0.0430166? 0.1647207?? 0.261??? 0.794
N?????????? 0.0009843???????? NA????? NA?????? NA
Log(theta)? 1.3049021???????? NA????? NA?????? NA
Zero-inflation model coefficients (binomial with logit link):
????????????? Estimate Std. Error z value Pr(>|z|)
(Intercept)? 0.0685048? 1.0056156?? 0.068? 0.94569
TreatmentF? -0.3822091? 0.3631982? -1.052? 0.29264
TreatmentO? -1.1624582? 0.3551596? -3.273? 0.00106 **
N??????????? 0.0009120? 0.0004773?? 1.911? 0.05603 .
---
Theta = 3.6873
Number of iterations in BFGS optimization: 31
Log-likelihood: -458.1 on 9 Df
Reading other posts, I know that my values are all >= 0, and there are
no missing values that could be causing this, but there ARE lots of
zeroes. So my questions are as follows:
1) Why am I getting all these NaNs, and why can?t I run the full
(saturated) model, including continuous as well as categorical
predictors and interactions?? How can I begin the process of model
simplification if I can?t get the full model to run?
2) If this is a result of the many zeroes in my data set, is there a
way to work around this?
Any advice would be greatly appreciated.
Jason Straka
University of Victoria