Skip to content

heteroscedastic model in lme4

7 messages · ONKELINX, Thierry, vito muggeo, Rense Nieuwenhuis +2 more

Original
ONKELINX, Thierry
# 
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:
<B9D1301370916C44B5874AF340C18B9B28AE890D50 at VMAILB.uoa.abdn.ac.uk>
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
structure to allow different variances for each level of my factor
(patch width).
suggestion wouuld be much apprecaited.
Group:sess + Group:VegDensity ,random=~1|Site, data=all,
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



_______________________________________________
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.
vito muggeo
# 
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:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Rense Nieuwenhuis
    

    

      
      #
    

  


  

      
Dear all,

I'm currently working on a similar problem in which I expect the  
within-group variance to change over years. Of course, in an ideal  
world correlation structures would be readily available in lme4, but  
as Thierry already states, this is not the case.

I have been considering a different approach to solve the problem,  
based on statements on heteroscedasticity made by Snijders & Bosker  
(1999, page 119). This only works for heteroscedasticity at the second  
(or higher) level, and if this heteroscedasticity can be related to an  
observed variable.

He argues that the random part of the second level can be defined as:
U0j + U1j*zj
"Thus, strange as it may sound, the level-two variable Z formally gets  
a random slope at level two".

In my understanding, this would lead to something like:

glmer(Count ~ A + B + (B|Group), family=poisson)

assuming that B measures a  'group'-level characteristic.

I'm still trying to get a firm grasp on the matter, but possibly this  
solution can help out in certain situations,

kind regards,

Rense Nieuwenhuis

        
On 15 jan 2009, at 11:13, ONKELINX, Thierry wrote:

        

            

      
      
      
    

        
@book{Snijders:1999,
	Author = {Snijders, Tom A.B. and Bosker, Roel J.},
	Date-Added = {2009-01-15 13:48:53 +0100},
	Date-Modified = {2009-01-15 13:49:35 +0100},
	Publisher = {Sage},
	Title = {Multilevel Analysis, an introduction to basic and advanced  
multilevel modelling},
	Year = {1999}}

  
  


  


      

    

    

      
      

        


  

        

      ONKELINX, Thierry
    

    

      
      #
    

  


  

      
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:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Doran, Harold
    

    

      
      #
    

  


  

      
Arrrg, that whole level 1 level 2 thing eats me up. There are random
effects and everything else is just a covariate. Not quite sure what
your data are, but a random slope can account for the non-constant
variance. Just a random intercept is compound symmetry where random
intercept and slope is a more general covariance structure that allows
for the variance to be different over time.

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Doran, Harold
    

    

      
      #
    

  


  

      
I would think on the transformed data. In a GLMM an offset is applied on
the transformed data and not on the original data, which is what makes
me think the same would be used here.

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Steven J. Pierce
    

    

      
      #
    

  


  

      
Have you looked at the spBayes package? That might offer an alternative
that's designed for doing spatial analyses.  


Steven J. Pierce
E-mail: pierces1 at msu.edu

-----Original Message-----
From: ONKELINX, Thierry [mailto:Thierry.ONKELINX at inbo.be] 
Sent: Thursday, January 15, 2009 8:47 AM
To: vito muggeo
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] heteroscedastic model in lme4

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:

            

      
      
      
    

        

      
      
    

        
--
====================================
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.