Hi Wayne,

Thanks for the suggestion, I will check out the cplm package.

Do you (or anyone else from the list) have any advice for me about using
the slope parameters from the mixed-effects models? My original email was
probably far to long for anyone to read, so I will paraphrase my primary
question here very briefly:

I am interested in change over time, so I would like to use slope
parameters for a number of different species (abundance data for 154
locations of 42 different species, measured in 5 equally-spaced years, the
same locations are measured in each year, thus the need for a random effect
of location).  I would like to model each species separately and use the
slope parameter from each model as an estimate of the direction and
magnitude of change over time for that species.  However, the distribution
of the data for each species are quite different - some are highly
zero-inflated, others are nearly normal.

My question is this: should I use the same model specifications for every
species (e.g., use negative binomial/zero-inflated distributions even for
those species that appear more normally distributed), or can I use the
"best" model for each species (based on the distribution of the data for
that species) and still compare the slope values to each other?  The
important thing is that the slope values are comparable, because I am
interested in how much each species has changed relative to the other
species.

I don't often see slope parameters from mixed models used in this way, so I
am a little hesitant about the feasibility of what I'm trying to do. Any
insight would be hugely appreciated.

Thanks very much!
Anne

Ph.D. Candidate, University of British Columbia
Vancouver, Canada

On Tue, Mar 6, 2012 at 7:37 PM, Zhang,Yanwei <Yanwei.Zhang@cna.com> wrote:

> Hi,
>
> The cplm package may be a good place to look at for the first problem. It
> handles zero-inflated continuous data using the compound Poisson
> distribution, and the mixed-model version (cpglmm) is developed based on
> glmer.
>
> Regards,
> Wayne Zhang
>
>
> -----Original Message-----
> From: r-sig-mixed-models-bounces@r-project.org [mailto:
> r-sig-mixed-models-bounces@r-project.org] On Behalf Of Anne Bjorkman
> Sent: Tuesday, March 06, 2012 7:31 PM
> To: r-sig-mixed-models@r-project.org
> Subject: [R-sig-ME] zero-inflation with mixed-effects, glmmADMB, and
> cloglog
>
> Dear Mixed Modelers,
>
>
> I have a few questions about various aspects of a dataset I am working on,
> and I was hoping to get some feedback from the experts.  I apologize for
> not being able to provide my data, but if anyone is interested I would be
> happy to share it off-line.
>
>
> My overall question is about changes in crop species abundances and
> presence/absence over a 40 year period in 154 countries in the world, with
> the measurements occurring at 10 year intervals (so, 1967, 1977, 1987,
> 1997, and 2007).  I would like to analyze these data in terms of both
> abundance AND presence/absence, as they answer slightly different questions
> (the first about abundance, the second about spread/distribution).  A
> subset of my data look like this:
>
>
>            Country             Item_general Total    Percent_of_Total Year
> pres.abs
>
> 28465 Albania                Apples  6.45       0.002903365 1967        1
>
> 28466 Albania   Bananas_and_plantains  0.00       0.000000000 1967        0
>
> 28467 Albania                Barley 17.09      0.007692792 1967        1
>
> 28468 Albania                 Beans 41.80      0.018815610 1967        1
>
> 28469 Albania   Beverages_Alcoholic  5.60       0.002520751 1967        1
>
> 28470 Albania   Beverages_Fermented  0.00       0.000000000 1967        0
>
>
> The abundance metric ("Total") here is number of Kcalories consumed of that
> crop in that country in that year.  I have already looked at multivariate
> changes in composition, but now I want to know about the change that occurs
> for each crop species (i.e., which crop species have increased or decreased
> the most over the past 40 years?).  My goal is to model each crop species
> separately and extract slope parameters and standard errors for the slopes,
> which would then tell me about the magnitude and direction of change for
> each crop.  However, there are three complications:
>
>
> 1) The data for most crops are zero-inflated.  The zero-inflation is not
> equal throughout the years, because overall countries have been increasing
> in the number of crops they use, so there are more 0's in 1967 than in 2007
> for a given crop.  This is true of most crops, but there are some that have
> been decreasing.
>
>
> 2) I have repeated measurements on individual countries over the 40 year (5
> sampling times) period. Thus I have been trying to account for this by
> using "Country" as a random effect in a mixed model.
>
>
> 3) The data for each crop are not distributed in the same way.  Some crops
> (such as corn) approach a nearly-normal distribution, while other crops
> (e.g., Tea) are extremely zero-inflated.
>
>
> I have attempted to use glmmADMB to model these data, as I understand this
> is one of the few ways you can have random effects and account for
> zero-inflation, using a negative binomial distribution and
> zeroInflation=TRUE.
>
>
> admb.model<-glmmadmb(Total_Rnd~I(Year-1967),random=~1|Country,family=
> "nbinom",zeroInflation=TRUE,data=mydata)
>
>
> where Total_Rnd is the "Total" column (total abundance) rounded to the
> nearest integer using "ceiling".  This usually works on a subset of the
> data that includes only one crop species, but if I try to model the entire
> dataset (all crops together) I get the error:
>
>
> Error in glmmadmb(Rel.Total.Rnd ~ I(Year - 1967), random = ~1 | Country,  :
>
>  The function maximizer failed (couldn't find STD file)
>
> In addition: Warning message:
>
> running command './glmmadmb -maxfn 500 -maxph 5 -noinit -shess' had status
> 1
>
> Need to increase the maximum number of separable calls allowed to at least
> 20001
>
> Current value is 20000
>
> Use the -ndi N command line option
>
>
> The output from just one crop species (Peas in this case) looks like this:
>
>
> Call:
>
> glmmadmb(formula = Total_Rnd ~ I(Year - 1967), data = subset(kcal.sub.test,
>
>    Item_general == "Pulses_Other"), family = "nbinom", random = ~1 |
>
>    Country, zeroInflation = TRUE)
>
>
>
> Coefficients:
>
>               Estimate Std. Error z value Pr(>|z|)
>
> (Intercept)     2.64530    0.13035   20.29   <2e-16 ***
>
> I(Year - 1967) -0.00107    0.00136   -0.78     0.43
>
> ---
>
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
>
> Number of observations: total=745, Country=154
>
> Random effect variance(s):
>
> $Country
>
>            (Intercept)
>
> (Intercept)      2.3298
>
>
> Negative binomial dispersion parameter: 5.8899 (std. err.: 0.48002)
>
> Zero-inflation: 0.013045  (std. err.:  0.0050869 )
>
>
> Log-likelihood: -2788
>
>
> This seems to match up somewhat with the raw data:
>
>
> > mean(kcal.sub.test$Total_Rnd[kcal.sub.test$Item_general=="Pulses_Other" &
> kcal.sub.test$Year==1967])
>
> [1] 31.81208
>
> > mean(kcal.sub.test$Total_Rnd[kcal.sub.test$Item_general=="Pulses_Other" &
> kcal.sub.test$Year==2007])
>
> [1] 29.69799
>
>
> (so there has been a decrease in mean abundance from 1967 to 2007)
>
>
> In addition, when I plot fitted values vs. residuals for a given crop to
> try to check the assumptions for these models, they do not look good (a
> very serious fanning-out from left to right):
>
>
> admb.model.fit<-fitted(admb.model)
>
> admb.model.res<-resid(admb.model)
>
> plot(admb.model.fit,admb.model.res)
>
>
> I have read on the mixed-models help list and elsewhere that plotting the
> fitted v. residual values is not as straightforward with mixed effects
> models, so I'm not entirely sure this is the right code to use, or what is
> a better way to check how well this model fits (other than comparing AIC
> values of a bunch of different models with the same data).
>
>
> Of course, as I mentioned above, the data for each crop are distributed
> differently, and I would be enormously appreciative if someone could
> comment on whether I could use different models (in other words, some of
> them with zeroInflation=TRUE and some with just a normal family="nbinom"
> without zero inflation) for the different crops and still compare their
> slopes.  Is my priority to find the best model for each crop, or is it to
> use the same model for each crop so that the slopes are directly
> comparable?  I have the same question about transforming the data.  Some of
> the crops (not the zero-inflated ones, of course) would be better fit by
> transforming to a square root of abundance, but I'm unsure whether I could
> then compare the slope parameters with non-square-root-transformed crops. I
> remember reading somewhere (but can't remember where) that
> back-transforming the parameters is not a simple matter with mixed models,
> thus making differently-transformed slope values comparable to each other
> might not be possible.
>
>
> To make matters more complicated, I would really like to model this as *
> relative* abundance (or, "Percent of Total") to account for the fact that
> people have been eating more over the years (total Kcalories has increased
> over time). I would still be dealing with zero-inflation here, and my
> understanding is that I could still use the same model described above,
> since negative binomial and poisson distributions are for counts, and
> proportions can be modeled like counts, even though they have an upper
> limit?
>
>
> The second part of the data, the presence-absence data, are a bit more
> straightforward.  I transformed the data to presence-absence to try to get
> away from the zero-inflation issue, and have been performing logistic
> regression using lmer:
>
>
> pres.abs.mod<-lmer(pres.abs~I(Year-1967)+(1|Country),data=mydata,family=
> binomial)
>
>
> For one crop, Peas, the output looks like this:
>
>
> Generalized linear mixed model fit by the Laplace approximation
>
> Formula: pres.abs ~ I(Year - 1967) + (1 | Country)
>
>   Data: subset(pres.abs.kcal, Item_general == "Pulses_Other")
>
>   AIC   BIC logLik deviance
>
>  255.2 269.1 -124.6    249.2
>
> Random effects:
>
>  Groups  Name        Variance Std.Dev.
>
>  Country (Intercept) 60.785 7.7965
>
> Number of obs: 760, groups: Country, 152
>
>
> Fixed effects:
>
>               Estimate Std. Error z value Pr(>|z|)
>
> (Intercept)     7.64468    1.47113   5.196 2.03e-07 ***
>
> I(Year - 1967)  0.06714    0.01747   3.843 0.000122 ***
>
> ---
>
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
>
> Correlation of Fixed Effects:
>
>            (Intr)
>
> I(Yer-1967) -0.196
>
>
> When I try to transform the Intercept term using plogis(7.64468) I get a
> value of 0.9995216 (which, if I understand correctly, is like the expected
> proportion of countries that have that crop) which I know from the data is
> much too high.
>
>
> My data do NOT have equal numbers of 1's and 0's, so if I try
> family=binomial(link="cloglog") I get a much more reasonable estimate:
>
>
> Generalized linear mixed model fit by the Laplace approximation
>
> Formula: pres.abs ~ I(Year - 1967) + (1 | Country)
>
>   Data: subset(pres.abs.kcal, Item_general == "Pulses_Other")
>
>   AIC   BIC logLik deviance
>
>  229.8 243.7 -111.9    223.8
>
> Random effects:
>
>  Groups  Name        Variance Std.Dev.
>
>  Country (Intercept) 10.482   3.2376
>
> Number of obs: 760, groups: Country, 152
>
>
> Fixed effects:
>
>               Estimate Std. Error z value Pr(>|z|)
>
> (Intercept)    2.458405   0.696478    3.53 0.000416 ***
>
> I(Year - 1967) 0.028910   0.009901    2.92 0.003503 **
>
> ---
>
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
>
> Correlation of Fixed Effects:
>
>            (Intr)
>
> I(Yer-1967) -0.292
>
>
> which matches up almost perfectly with the mean of the raw data (which is
> what I think the Intercept term estimates...):
>
>
> > plogis(2.458405)
>
> [1] 0.9211739
>
>
> > mean(pres.abs.kcal$pres.abs[pres.abs.kcal$Item_general=="Pulses_Other"])
>
> [1] 0.9276316
>
>
>
> the problem is, I get an error that says
>
>
> Warning message:
>
> In mer_finalize(ans) : false convergence (8)
>
>
> Can I still use the output of the link="cloglog" model, since it seems to
> fit my model better (incidentally, the AIC is lower and the plot of fitted
> v. residuals looks more even - even though of course for logistic models
> the plots of fitted v. residuals aren't exactly straightforward to
> interpret!).  Or can someone advise me as to what I am doing wrong to get
> this error?
>
>
> Finally, similarly to my question posed above, can I still compare the
> slope parameters from my logistic regression models if I use link="logit"
> for some and link="cloglog" for others, depending on which one provides a
> better model fit, or should I use the same model for all crop species, even
> if some are fit better than others.
>
>
> I sincerely apologize for the novel-like length of this email, I have been
> battling these questions for over a month now, searching R-list archives
> and reading every (layman's) book I can get my hands on, and I've gotten to
> the point now where, for every new thing I learn, I forget something I read
> yesterday.
>
>
> Enormous thanks in advance to anyone who feels willing to tackle my
> statistical novel!
>
>
> Best,
>
> Anne
>
>         [[alternative HTML version deleted]]
>
>
> NOTICE:  This e-mail message, including any attachments and appended
> messages, is for the sole use of the intended recipients and may contain
> confidential and legally privileged information.
> If you are not the intended recipient, any review, dissemination,
> distribution, copying, storage or other use of all or any portion of this
> message is strictly prohibited.
> If you received this message in error, please immediately notify the
> sender by reply e-mail and delete this message in its entirety.
>

	[[alternative HTML version deleted]]


