[R-sig-ME] heteroscedastic model in lme4

Thu Jan 15 14:47:22 CET 2009

Dear Vito,

I aggree that the glmer() and nlme() examples assume different distribution. That's why I called the nlme() version an apprioximation. If the steps described in P&B are only valid for linear models but not in the generalised models, then I have a dilemma. With glmer() I can use the appropriate distribution but a wrong correlation structure. A structure of which I'm certain that it is there (spatially clustered points). nlme() allows me to model the correlation structure but only unther the gaussion distribution. The latter is in my opinion a better alternative given that with enough data the residuals will behave approximately gaussian. Please do correct me if that is an incorrect statement.

Why nlme() and not lme() with a log-transformation? Well: the zero's in counts. Using a log(x + 1) transformation complicates the interpretation of the model. And what transformation would you suggest with binomial data? nlme() handles zero's with a log-link as well as true-false data with a logit-link.

Thierry

----------------------------------------------------------------------------
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium 
tel. + 32 54/436 185
Thierry.Onkelinx at inbo.be 
www.inbo.be 

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

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

-----Oorspronkelijk bericht-----
Van: vito muggeo [mailto:vmuggeo at dssm.unipa.it] 
Verzonden: donderdag 15 januari 2009 13:21
Aan: ONKELINX, Thierry
CC: Doran, Harold; Alan Cobo-Lewis; r-sig-mixed-models at r-project.org
Onderwerp: Re: [R-sig-ME] heteroscedastic model in lme4

dear Thierry,
I am adding a simple comment only on your second point.

If I am not wrong, I think that the two alternatives underlie different 
models

1)glmer(.., family = poisson) assumes a real Poisson distribution for 
your response y (conditioned to random effects), i.e. y=rpois(n,exp(mu)).

2) nlme(..) assumes a gaussian distribution for your response with a 
nonlinear mean model, i.e. y=rnorm(n,exp(mu))

Another (different) approach would be lmer() with log-transformed data, 
i.e. y=exp(rnorm(n,mu))

Probably, in a pure likelihood framework the first approach should be 
preferred if you have real count data..

Hope this helps,

vito

ONKELINX, Thierry ha scritto:
> Dear all,
> 
> I would like to analyse some spatial data with mixed model. As I'm
> dealing with presence/absence data or counts I should use the bionomial
> or poisson family. These families are implemented in lme4 but
> correlation structures are not. I'm wondering if the steps from section
> 5 in Pinheiro and Bates can be applied in case of a GLMM. If one can do
> that, should one apply the transformation on the response in the
> original scale or the transformed (logit / log) scale?
> 
> Another, more approximate, solution might be to code the GLMM as a NLMM.
> E.g. glmer(Count ~ A + B + (1|Group), family = poisson) versus
> nlme(model = Count ~ exp(mu), fixed = mu ~ A + B, random = mu ~ Group)
> Any ideas on that?
> 
> Thierry
> 
> ------------------------------------------------------------------------
> ----
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
> and Forest
> Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
> methodology and quality assurance
> Gaverstraat 4
> 9500 Geraardsbergen
> Belgium 
> tel. + 32 54/436 185
> Thierry.Onkelinx at inbo.be 
> www.inbo.be 
> 
> To call in the statistician after the experiment is done may be no more
> than asking him to perform a post-mortem examination: he may be able to
> say what the experiment died of.
> ~ Sir Ronald Aylmer Fisher
> 
> The plural of anecdote is not data.
> ~ Roger Brinner
> 
> 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
> 
> -----Oorspronkelijk bericht-----
> Van: r-sig-mixed-models-bounces at r-project.org
> [mailto:r-sig-mixed-models-bounces at r-project.org] Namens Doran, Harold
> Verzonden: vrijdag 19 december 2008 20:52
> Aan: Alan Cobo-Lewis; r-sig-mixed-models at r-project.org
> Onderwerp: Re: [R-sig-ME] heteroscedastic model in lme4
> 
> This isn't an entirely accurate statement. nlme has built-in functions
> that implement the methods for correlational and variance structures as
> described in section 5 of Pinhiero and Bates. lme4 doesn't have these
> functions built in as does nlme, but those same methods can be
> implemented by the user and then the data can be analyzed using
> functions in lme4. So, functions in lme4 can "handle" the same issues as
> nlme, it just requires the user to perform the steps described in PB
> section 5 et seq on their own. 
> 
> 
> 
> 
> -----Original Message-----
> From: r-sig-mixed-models-bounces at r-project.org on behalf of Alan
> Cobo-Lewis
> Sent: Fri 12/19/2008 11:19 AM
> To: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] heteroscedastic model in lme4
> 
> 
> Anna,
> 
> lme4 cannot handle certain kinds of heteroscedasticity, but I believe it
> can handle the kind you have in mind. Search the r-sig-mixed-models
> archive for a discussion involving me and David Afshartous, especially
> the summary message titled
> "[R-sig-ME] random effect variance per treatment group in lmer" that
> David posted 07/13/2007 04:18:08 PM
> 
> I can't be certain that the suggestion below would work without knowing
> more about your design, but if width were a factor with three levels
> then you might try setting up indicator variables Wind1, Wind2, and
> Wind3 (that each take on the value 1
> when a site is at the indicator's target width and 0 otherwise) and then
> fit the model with something like
> mrem <- lmer( log(Nhat+1)~Group + GreenPerc + sess + crop + VegDensity +
> Group:sess + Group:VegDensity + (0+Wind1|site) + (0+Wind2|site) +
> (0+Wind3|site), data=all, method="REML" )
> 
> alan
> 
> 
> r-sig-mixed-models at r-project.org on Friday, December 19, 2008 at 6:00 AM
> -0500 wrote:
>> Message: 1
>> Date: Thu, 18 Dec 2008 11:23:46 +0000
>> From: "Renwick, A. R." <a.renwick at abdn.ac.uk>
>> Subject: [R-sig-ME] heteroscedastic model in lme4
>> To: "'r-sig-mixed-models at r-project.org'"
>> 	<r-sig-mixed-models at r-project.org>
>> Message-ID:
>>
> <B9D1301370916C44B5874AF340C18B9B28AE890D50 at VMAILB.uoa.abdn.ac.uk>
>> Content-Type: text/plain; charset="us-ascii"
>>
>> I have been using the nlme package to run some LMM's, however I would
> like to try rerunning them using the lme4 package so that I can use mcmc
> sampling.  The data I am using shows some heteroscesdasticity of the
> within error group and so I have
>> been using the 'weights' argument and the varIdent variance function
> structure to allow different variances for each level of my factor
> (patch width).
>> My problem is how to code for a heteroscedastic model in lme4 and any
> suggestion wouuld be much apprecaited.
>> The code I used in the nlme package:
>>
>> # model fit using "REML"
>> mrem<-lme(log(Nhat+1)~Group + GreenPerc + sess + crop + VegDensity +
> Group:sess + Group:VegDensity ,random=~1|Site, data=all,
>>       method="REML",correlation=NULL,weights=varIdent(form=~1|width))
>>
>>
>> Many thanks,
>> Anna
>>
>> Anna Renwick
>> Institute of Biological & Environment Sciences
>> University of Aberdeen
>> Zoology Building
>> Tillydrone Avenue
>> Aberdeen
>> AB24 2TZ
>>
>>
>> The University of Aberdeen is a charity registered in Scotland, No
> SC013683.
> 
> 
> --
> Alan B. Cobo-Lewis, Ph.D.		(207) 581-3840 tel
> Department of Psychology		(207) 581-6128 fax
> University of Maine
> Orono, ME 04469-5742     		alanc at maine.edu
> 
> http://www.umaine.edu/visualperception
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer 
> en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
> door een geldig ondertekend document. The views expressed in  this message 
> and any annex are purely those of the writer and may not be regarded as stating 
> an official position of INBO, as long as the message is not confirmed by a duly 
> signed document.
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> 

-- 
====================================
Vito M.R. Muggeo
Dip.to Sc Statist e Matem `Vianelli'
Università di Palermo
viale delle Scienze, edificio 13
90128 Palermo - ITALY
tel: 091 6626240
fax: 091 485726/485612
http://dssm.unipa.it/vmuggeo
====================================

Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer 
en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document. The views expressed in  this message 
and any annex are purely those of the writer and may not be regarded as stating 
an official position of INBO, as long as the message is not confirmed by a duly 
signed document.