An update: I did work out the issue shortly after posting the email, and
confirmed my mistake with Jarrod Hadfield, who was very helpful! A scan of
his J. Stat. Soft. paper showed me I was using the wrong link: for
exponential family, the latent variables (l) are mapped onto the data (y)
using y=exp(-l).

David Atkins' point (below: thanks for your help!) is also helpful: for
non-identity links there is a difference between conditional and
marginalised posteriors, and to get the expected marginal value rather than
the conditional mode it's necessary to average over the random effects.
This is something I wasn't aware of! The predict method for MCMCglmm
doesn't support the exponential family yet, but by hacking the code and
augmenting with info from [
http://depts.washington.edu/cshrb/wordpress/wp-content/uploads/2013/04/Tutorials-Tutorial-on-Count-Regression-R-code.txt]
it should be possible... except that I'm not sure on the analytical
solutions for marginalising random effects in exponential models. If anyone
has a solution or suggestion, it'd be very welcome. Failing that I'll
continue to google, or acquiesce to a less-satisfactory polynomial Gaussian
model or gamma lmer.

Thanks all

Iain



On 6 September 2013 21:00, David Atkins <datkins@u.washington.edu> wrote:
>
>>
>> Iain--
>>
>> Just writing back to you as a) I'm not sure about this, and b) it isn't
>> necessarily R specific.
>>
>> One of the things that has bitten me a couple times with mixed models
>> with a non-identity link function is the distinction between conditional
>> and marginal effects.  With normal errors and identity link function,
>> predictions from a mixed model essentially average over the random effects.
>>  With non-identity link functions this is no longer true, and the
>> predictions based on fixed-effects only are conditional on specific levels
>> of the random effects.  It is possible to derive marginal predictions,
>> which do in effect average over the random effects, but this is somewhat
>> more involved and the predictions include random-effects in them.
>>
>> Not sure whether this is what is happening with your data -- since, I'm a
>> psychologist who mostly analyzes psychiatric / drug abuse data, but I would
>> not necessarily assume that your predictions based only on the
>> fixed-effects would give 'good' predictions of your raw data.
>>
>> You might check out an article that colleagues and I wrote on mixed
>> models; in particular, the appendix covers this issue specifically (for
>> Poisson models):
>>
>> http://depts.washington.edu/**cshrb/statistics-resources-**
>> tutorials-tutorial-on-count-**regression/<http://depts.washington.edu/cshrb/statistics-resources-tutorials-tutorial-on-count-regression/>
>>
>> Hope that might help.
>>
>> cheers, Dave
>>
>>
>> Hi R-users,
>>
>> I'm having some problems interpreting mixed effects models using the
>> exponential family. Struggling to find any info on such models on the
>> interwebs.
>>
>> Some background: I'm working with soil carbon data and need robust
>> estimates of concave-up, decreasing depth curves, and confidence intervals
>> on them, which MCMCglmm should be able to provide. The usual way to fit
>> such curves is with a (negative) exponential model, which again is
>> supported by MCMCglmm (my understanding is that the exponential and
>> negative exponential distibutions are the same in everything but name?).
>> My
>> data is pretty much exponentially-distributed (more precisely
>> gamma-distributed with a very strong right skew), therefore is real with
>> range 0<x<Inf.
>>
>> The models: I'm fitting models with one continuous variable (depth), and
>> two categorical variables (soil type / land cover). I've included a single
>> random variable (site), for which I'm fitting either just an intercept or
>> a
>> fully-(co)varying intercept and slope combination using us(). The fixed
>> effects take default priors of N(0, 10e8) and random effects and residuals
>> take uninformative proper priors with V=1, nu=0.002. Below, I've pasted
>> code which will replicate my data and model structure: the range of
>> solutions (intercepts, slopes) I'm getting out of these dummy models are
>> in
>> accordance with what's being computed for my data.
>>
>> The results: Models converge well (albeit after a large number of
>> iterations), but the estimates I'm getting out of them make no sense to me
>> at all. I must be missing something somewhere! Transformed curves do not
>> fit the raw data, and transformed data do not fit the fitted lines. I'm
>> using the canonical inverse link for the exponential distribution,
>> Xbeta=-mu^-1 (and various variations on that theme, to no avail...)
>>
>>
>> I can think of three possible reasons I can't make sense of this:
>>
>> 1. My understanding of the exponential distribution is woefully poor: that
>> it is not equivalent to the negative exponential distribution (in which
>> case Wikipedia would be lying to me) or that I am missing a crucial piece
>> of theory somewhere along the line.
>>
>> 2. Default fixed priors of N(0, 10e8) are completely inadequate for
>> exponential data
>>
>> 3. I've misinterpreted the link function in some way: it's not inverse,
>> the
>> inverse of inverse is for some reason not inverse (!) or I'm using the
>> wrong link altogether.
>>
>>
>> I can't find anything to contradict my thoughts on these points. So, if
>> someone with a better understanding of such things could offer me some
>> words of wisdom, point out to me where I'm making stupid mistakes, it
>> would
>> be hugely appreciated!
>>
>> Thanks
>>
>> Iain
>>
>>
>>
>> CODE:
>>
>>
>> MODELS:
>>
>> resp<-rgamma(406,scale=0.5,**shape=1.4)*75
>> hist(resp,breaks=100)
>>
>> fix1<-rep(c("A","B","C","D","**E","F","G"),58);
>> fix1<-fix1[order(rnorm(406))]; fix1<-as.factor(fix1)
>> fix2<-rep(c("ONE","TWO"),203); fix2<-fix2[order(rnorm(406))];
>> fix2<-as.factor(fix2)
>> fix3<-runif(406)*100
>>
>> ran1<-rep(1:58,7); ran1<-as.factor(ran1)
>>
>> data<-data.frame(resp,fix1,**fix2,fix3,ran1)
>>
>> nitt<-1000000
>> burn<-nitt*0.1
>> thin<-nitt*0.0005
>>
>> prior1.0<-list(R = list(V=1, nu=0.002),
>>                G = list(G1=list(V=1, nu=0.002)))
>>
>>
>> prior1.1<-list(R = list(V=1, nu=0.002),
>>                G = list(G1=list(V=diag(1,2), nu=0.002)))
>>
>> m1.0<-MCMCglmm(resp~fix1+fix2+**fix3+fix1:fix3+fix2:fix3,
>>                random=~ran1, prior=prior1.0, data=data,
>> family="exponential", verbose=F,
>>                nitt=nitt, burnin=burn, thin=thin)
>>
>> m1.1<-MCMCglmm(resp~fix1+fix2+**fix3+fix1:fix3+fix2:fix3,
>>                random=~us(1+fix3):ran1, prior=prior1.1, data=data,
>> family="exponential", verbose=F,
>>                nitt=nitt, burnin=burn, thin=thin)
>>
>> EXAMPLE FIT:
>>
>> sols<-summary(m1.0)$solutions
>>
>> x<-1:100
>> y<-sols["(Intercept)",1]+x***sols["fix3",1]
>> plot(resp~fix3); lines(x,-y^-1)
>> cbind(x,y,-y^-1)
>>
>>
>>
>>
>>
>> - - - - - - - - - - - -
>> Dr. Iain Stott
>> Environment and Sustainability Institute
>> University of Exeter, Penryn Campus
>> Treliever Road, Penryn,
>> Cornwall, TR10 9EZ, UK.
>> http://www.exeter.ac.uk/esi/**people/stott/<http://www.exeter.ac.uk/esi/people/stott/>
>> - - - - - - - - - - - -
>> http://www.exeter.ac.uk/esi/
>> http://biosciences.exeter.ac.**uk/cec/<http://biosciences.exeter.ac.uk/cec/>
>>
>>         [[alternative HTML version deleted]]
>>
>> --
>> Dave Atkins, PhD
>>
>> Research Professor
>> Department of Psychiatry and Behavioral Science
>> University of Washington
>> datkins@u.washington.edu
>> http://depts.washington.edu/**cshrb/david-atkins/#more-48<http://depts.washington.edu/cshrb/david-atkins/#more-48>
>>
>> "The combination of some data and an aching desire for an answer does not
>> ensure that a reasonable answer can be extracted from a given body of
>> data." -- John Tukey
>>
>
>

	[[alternative HTML version deleted]]

