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/
http://biosciences.exeter.ac.uk/cec/

	[[alternative HTML version deleted]]

