[R-sig-ME] Dealing with heteroscedasticity in repeated measure models

Sun Sep 26 16:14:45 CEST 2010

Hi Manuel,

Your first problem to overcome is that your response data are not from
a normal (gaussian) distribution. Far from it.

See this:

ipaDat <- read.table("ipaloc.txt", header=T, sep="\t")

library(lattice)
histogram(~ipa|tiempo, data=ipaDat)

Using an appropriate transformation or error distribution will likely
solve the heteroscedasticity too.

The quick and dirty method would be to do a log transformation of your
response variable. You have some zero values so you would have to
either exclude some data or apply a kludge like adding 1 to all values
(or perhaps 1/2 the minimum non-zero value).

How is ipa calculated? You say it is an incidence. Did you start with
a raw count?

It's better if you reply to the list rather than entering private
email chats with respondents.

andydolman at gmail.com

2010/9/26 Manuel Spínola <mspinola10 at gmail.com>:
>  Thank you very much Andy.
>
> I am attaching the data.  "tempo" in my model is "tiempo" centered.
>
> I am trying to see is ipa shows a declining trend for the 4 years period.
>  Ipa refer to incidence of malaria in these counties (localidades).
> May be I need to see the difference in IPA between year 1 and 4 and model
> that?
> IPA usually decrease in time, for several counties.
> I tried lmer, but I found that lmer do not have the capability to deal with
> covariance structures and that lme has much more options, for example the
> function "intervals" to get the CI.
>
> If I use a poisson, should I get the count and an offset variable (the way
> that the incidence rate was calculated), right?
>
> Best,
>
> Manuel
>
>
> On 26/09/2010 06:51 a.m., Andrew Dolman wrote:
>>>
>>> I found an increase of variance when increasing time, is that a
>>> heteroscedasticity problem?
>>
>> Yes, but it could be due to several things. For example, does your
>> dependent variable, ipa, increase with time? For many measurements the
>> variance increases with the mean so if the mean increases with time so
>> will the variance. To deal with this you would need to look at using
>> either a transformation, e.g. log, or using a generalized linear model
>> with an appropriate error distributions, e.g. poisson. It could also
>> be that you have divergent time series due to autocorrelation, in
>> which case your third model might be appropriate.
>>
>> I'm having to guess because you haven't given enough information about
>> the data you are modeling. Ideally this would include a sample of the
>> data or dummy data that reproduces the problem.
>>
>>
>>> I fitted the model in the way you suggested me but I got an error (lack
>>> of
>>> convergence).  Is there any way to get around this problem?
>>
>> This is fairly likely because you only have 4 data points per
>> localidad. You may have more luck using lmer from the lme4 package.
>>
>>
>> Andy.
>>
>>
>>
>>>> modelo2 = lme(ipa ~ tempo, random = ~ tempo | localidad, data=ipa) #
>>>> intercepto y pendiente diferentes para cada localidad
>>>
>>> Error en lme.formula(ipa ~ tempo, random = ~tempo | localidad, data =
>>> ipa) :
>>>  nlminb problem, convergence error code = 1
>>>  message = iteration limit reached without convergence (9)
>>>
>>>
>>>
>>> On 26/09/2010 05:58 a.m., Andrew Dolman wrote:
>>>>
>>>> Hi Manuel,
>>>>
>>>> First of all a quick correction to your basic model specifications.
>>>>
>>>> First model is correct:
>>>> model1 = lme(ipa ~ tempo, data=ipa, random = ~1 | localidad) # random
>>>> intercept
>>>>
>>>> Second model is wrong
>>>> model2 = lme(ipa ~ tempo, random = ~1 | tempo/localidad, data=ipa) #
>>>> random intercept and slope
>>>>
>>>> should be
>>>> model2 = lme(ipa ~ tempo, random = ~tempo|localidad, data=ipa) #
>>>> random intercept and slope
>>>>
>>>> Your 3rd model fits an autoregressive model for temporally correlated
>>>> data. Whether this deals with heteroscedasticity depends on the type
>>>> of heteroscedasticity .
>>>>
>>>> model3 = lme(ipa ~ tempo, random = ~1 | localidad, data=ipa,
>>>> correlation=corAR1(form=~ tempo))
>>>>
>>>>
>>>>
>>>> andydolman at gmail.com
>>>>
>>>>
>>>>
>>>> On 26 September 2010 13:24, Manuel Spínola<mspinola10 at gmail.com>
>>>>  wrote:
>>>>>
>>>>>  Dear list members,
>>>>>
>>>>> I am fitting a repeated measure model using lme.
>>>>> I have 4 measurements of a rate (variable called ipa) measured each
>>>>> year
>>>>> (variable called tempo, which was centered) on 48 counties (all the
>>>>> counties
>>>>> from a province, variable called localidad).
>>>>> I am considering county like a random factor.
>>>>> My models are:
>>>>>
>>>>> model1 = lme(ipa ~ tempo, data=ipa, random = ~1 | localidad) # random
>>>>> intercept
>>>>>
>>>>> model2 = lme(ipa ~ tempo, random = ~1 | tempo/localidad, data=ipa) #
>>>>> random
>>>>> intercept and slope
>>>>>
>>>>> model3 = lme(ipa ~ tempo, random = ~1 | localidad, data=ipa,
>>>>> correlation=corAR1(form=~ tempo))
>>>>>
>>>>> I have heteroscedasticity.
>>>>> Is my last model dealing with heteroscedasticity?
>>>>> Thank you very much in advance.
>>>>> Best,
>>>>>
>>>>> Manuel
>>>>>
>>>>> --
>>>>> Manuel Spínola, Ph.D.
>>>>> Instituto Internacional en Conservación y Manejo de Vida Silvestre
>>>>> Universidad Nacional
>>>>> Apartado 1350-3000
>>>>> Heredia
>>>>> COSTA RICA
>>>>> mspinola at una.ac.cr
>>>>> mspinola10 at gmail.com
>>>>> Teléfono: (506) 2277-3598
>>>>> Fax: (506) 2237-7036
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-mixed-models at r-project.org mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>>
>>>
>>> --
>>> Manuel Spínola, Ph.D.
>>> Instituto Internacional en Conservación y Manejo de Vida Silvestre
>>> Universidad Nacional
>>> Apartado 1350-3000
>>> Heredia
>>> COSTA RICA
>>> mspinola at una.ac.cr
>>> mspinola10 at gmail.com
>>> Teléfono: (506) 2277-3598
>>> Fax: (506) 2237-7036
>>>
>>>
>
>
> --
> Manuel Spínola, Ph.D.
> Instituto Internacional en Conservación y Manejo de Vida Silvestre
> Universidad Nacional
> Apartado 1350-3000
> Heredia
> COSTA RICA
> mspinola at una.ac.cr
> mspinola10 at gmail.com
> Teléfono: (506) 2277-3598
> Fax: (506) 2237-7036
>
>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: ipaloc.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20100926/380ec7cc/attachment.txt>