Hello all,

 

Thanks to everyone who responded to my question regarding the behavior of deviances when terms are added to a hierarchical logistic regression model.  I received a number of very helpful comments and suggestions, and I have included all of them at the end of the email.

 

Some main points most relevant to my problem include:

-       In hierarchical models, deviances don’t necessarily decline when adding fixed effects.

-       The random effect at the data level (my attempt to add extra-binomial variance) is inappropriate, and the L2 random effect is suspect.

-       This problem occurs when there are few data points per group (the data level effect is an extreme example, with 1 point per group).

-       The deviance calculations of BUGS may be problematic for discrete data distributions.

 

Based on these responses, I’m guessing that the deviances may in fact be perfectly reasonable estimates, but the effective number of parameters is varying greatly from model to model as I add fixed effects.  The problem is exacerbated by my model formulation and by the fact that the species I’m modeling are rare.  So, a model with no fixed effects may still have a “good” deviance score, because the variance in the random effects is high, and so the effective number of parameters is high.  As I add fixed effects this variance (and therefore the number of parameters) declines, but the deviance stays the same (or under some conditions gets slightly better or worse).  DIC should address this in cases where it works, but I’m afraid it may not be appropriate for my situation.

 

To go out on a limb into the realm of speculation, it seems to me that someone clever could develop an estimated criterion based on the change in variance in the random effects as well as the deviance and number of fixed parameters.  Alternatively, it seems plausible to me that one could first fit a model fully saturated with random effects but no fixed effects, and compare different fixed effects models based on the decline in random effect variance and the number of fixed parameters.  But this has probably already been fully explored and I’m just ignorant of the literature.

 

All responses follow.  The original question is pasted at the end.  Thanks again.

 

Seth Wenger

 

___

 

If I get you right, this is a problem commonly present both in Bayesian and maximum likelihood estimations when observations and/or cell frequencies within each cluster are few. Try fitting a prior distribution (well justified) to the cluster that have the problem. Note than when comparing

models the data and priors should be the same among models.

 

Polychronis P. Kostoulas

___

 

The model description you gave isn't complete, but it looks like your model is overfitted, as you suggested.  You need  at least two binary outcomes at each site in order for the site-level random effects to be identifiable, but you only have one.

 

Martyn Plummer

 

I asked Dr. Plummer a follow-up question as to whether the L2 effect should be restricted to only those groups with more than one sample.  His reply:

 

I don't think you should limit random effects to reaches with multiple sites. What I said wasn't quite correct, in fact. To be precise, the *variance* of the random effects is unidentifiable if there are no replicates. As long as you have some replicates, the variance will - formally speaking - be identifiable and the model will work correctly, even for reaches with only 1 site.

 

I say "formally speaking" because it might still be highly sensitive to the prior... but that's another can of worms.

 

Martyn Plummer

___

 

I don't have any hard scientific evidence for this so, if you send this email out, please delete my name!), but I have always been suspect of BUGS calculations of deviances.  This is certainly true for models with data distributions that are discrete (I think this is a known problem), but I never use the automated deviance calculations anymore.

 

What I do to make sure my likelihoods are correct is to construct a logical node in my model that computes the log likelihood of the data. For a logistic regression model, you know the vector of probabilities, so you can take the log of the element of the probability vector that corresponds to the observed data point.  Sum them up over observations and you have a log likelihood for that particular MCMC draw.

 

To compare models, I use the Newton and Raftery method of computing log marginal likelihoods (this involves taking the harmonic mean of the log likelihood draws).  I typically save the draws and then compute the LML in Matlab.  Congdon's book describes this method in detail.

 

Others may have different opinions of this (for example, those that support DIC instead of Bayes factors and LML), but this approach gives me more confidence in the numbers.

 

(name withheld as requested)

___

 

I think you may be encountering a phenomenon similar to one I encountered with F-statistics and dominant markers (see the attached

paper). When you add predictors to a hierarchical model, it constrains the relationships among parameters in a way that may reduce the fit of

the model to data. In my case it occurred because assuming a common inbreeding coefficient across all populations means that the within

population band frequencies for dominant markers can't match as closely as if allele frequencies were estimated independently in every population.

 

Fundamentally, that does mean that using DIC as a model choice criterion is a bit dicey. I don't mention it in the attached paper (because I hadn't explored it yet), but it appears that a posterior predictive approach (like the one Gelfand & Ghosh proposed) may be more suitable. It is, however, an open question as to which criterion is the "best" and one that Bayesians are still struggling with. Alan Gelfand isn't convinced there is a "best" criterion. He thinks different criteria tell you different things -- and he's probably right.

 

Kent E. Holsinger

___

 

I agree that deviance scores should not increase as you add more species.  My question is whether or not you have normalized the scores of each species to each other's mean.  Otherwise, I would expect an increase in deviance score, if the median scores differ. 

 

Daniel M. Byrd III, Ph.D., D.A.B.T.

___

 

I'm not sure why this would occur, although I do thing that the ordinary theory of 'deviance goes down as number of parameters goes up' doesn't necessarily hold for random effects, hierarchical models.

 

When you say 'deviance scores' though, do you mean the mean or mode of the posterior distribution of these? How are the distributions shaped under the different models?

 

Finally, I can't help but wonder how this problem would behave if fit using pymc.  This should be easy to do and would require little reformatting of the data. For instance your first piece could be like so (with appropriate initiatlizaiton of self.a and self.tau)

 

like=self.normal_like(self.a,0,self.tau)

for i in range(no_sites):

   like +=self.bernoulli_like(A[i,51],invlogit(self.a[i]))

 

Dr. Michael J. Conroy

___

You are fitting a random effect at the data level with a Bernoulli model, which cannnot ever be appropriate and I think means that none of

the deviances are interpretable.  You need to remove err1 from the model.  I think the second level random effect is dubious too as there

appear to be hardly any replicates, and I would imagine tau2 is unestimable without strong prior inofrmation.

 

David Spiegelhalter

___

As to me, with Poisson multilevel models, I have found a similar problem, in that the deviance (and DIC) was nearly constant between non-nested models, although some of the models included clearly non-zero variables and others not, as you...

The problem seems to come from a trade-off between the different levels of variations: in models with good fixed effects variables, more variation is in the residual, whereas in poorer fixed effect models, more variation goes into the random effects... I have not managed to solve this problem so far.

Frederic Gosselin

___

 

Hello all,

 

We are trying to compare various multilevel logistic regression models to

explain presence/absence of different fish species based on predictors at

multiple scales.  The problem is, with some models the deviance scores do

not appear credible.  Specifically, when adding certain predictor

variables the deviance increases, rather than decreases.  This happens

even though examination of the credible interval for the parameter in

question shows that the predictor is strong and clearly non-zero.  The

problem disappears in a non-hierarchical formulation of the model, when

the random effects are removed.  I should also note that the problem only

occurs for species with very few positive occurrences.

 

Does anyone else have examples of such a problem?  Is this somehow a side

effect of overfitting the data?  Without credible deviances we are unable

to use an information criterion, which makes model selection trickier.

 

I have tried reformulating the models in various ways but nothing helps. 

Here is one version-- the full version also incorporates incomplete

detectability, but the problem is the same with or without that added

complication.

 

# Level 1. Site

for (i in 1:232) {

 A[i,51] ~ dbern(P[i])

 logit(P[i]) <- a[1]*A[i,10] + a[2]*A[i,32] + L2[A[i,8]] + err1[i]  err1[i] ~ dnorm(0, tau[1])  }

 

# Level 2. Reach-level. 

 for (j in 1 : 204) {

    L2[j] ~ dnorm(mu1[j], tau[2])

    mu1[j] <- L3[B[j,2]]  }

 

# Level 3. Trib system.

 for ( k in 1 : 24) {

    L3[k] ~ dnorm(mu2[k], tau[3])

    mu2[k]<- c0 + c1*C[k,2]  }

 

Thanks for any suggestions.

 

Seth Wenger

Institute of Ecology

University of Georgia

Athens, GA USA