[R-sig-ME] R-sig-mixed-models Digest, Vol 72, Issue 1
chico3 at sapo.pt
chico3 at sapo.pt
Sat Dec 1 18:33:53 CET 2012
Quoting r-sig-mixed-models-request at r-project.org:
> Send R-sig-mixed-models mailing list submissions to
> r-sig-mixed-models at r-project.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> or, via email, send a message with subject or body 'help' to
> r-sig-mixed-models-request at r-project.org
>
> You can reach the person managing the list at
> r-sig-mixed-models-owner at r-project.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of R-sig-mixed-models digest..."
>
>
> Today's Topics:
>
> 1. mer_optimize (Ivar Herfindal)
> 2. glmer, split-plot design simple question (chico3 at sapo.pt)
> 3. Re: glmer, split-plot design simple question (Ben Bolker)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Fri, 30 Nov 2012 13:57:09 +0000
> From: Ivar Herfindal <ivar.herfindal at ntnu.no>
> To: "r-sig-mixed-models at r-project.org"
> <r-sig-mixed-models at r-project.org>
> Subject: [R-sig-ME] mer_optimize
> Message-ID:
> <D6C55C60CCF1AC44AF5805862F86F997022BA378 at WAREHOUSE04.win.ntnu.no>
> Content-Type: text/plain
>
> Hi all
>
> Some time ago, I raised a question on this mailing list about using
> .Call("mer_optimize", mymodel, PACKAGE="lme4") for a model that
> failed to converge, see:
> https://stat.ethz.ch/pipermail/r-sig-mixed-models/2008q3/001354.html
>
> As a "horrible hack" solution to my problem, I added a few lines
> with .Call("mer_optimize", mymodel, PACKAGE="lme4") in my script,
> and was satisfied that the model seemed to converge after a few of
> these runs. Then I forgot it all, until I recently ran the same
> script but on a different model.
>
> This new model did not have problem with convergence (at least not
> any problem that are reported), it is a simple linear mixed model
> with one random factor. I was therefore surprised that the
> log-likelihood and parameter estimates change after I run
> .Call("mer_optimize"...) on the model. Since I use log-likelihoods
> to compare models (likelihood-ratio tests or AICc), it is quite
> important to have the most "correct" log-likelihood. The change is
> quite big (change of 2 - 3 in AICc) and differs between candidate
> models, such that it can have impact on the ranking of models. Also,
> the parameter estimates and their SE change, which affect the
> corresponding t-values. It seems that more complex models have the
> largest change after calling mer_optimize.
>
> My question is then: Should I continue to use
> .Call("mer_optimize"...) on my models until log-Likelihood not
> change for a model, or can I trust the first run of lmer as long as
> I do not get any warnings about convergence?
>
> I attach below example for a model summary before and after running
> mer_optimize. Unfortunately, I have not been able to make a small,
> reproducible example.
>
> Best regards
>
> Ivar
>
> ##########################################################################
>> mod <- lmer(LRS ~ Apr*ALR + I(ALR^2) + (1|Lok), data=LRStab, REML=FALSE)
>> mod
>> mod
> Linear mixed model fit by maximum likelihood
> Formula: LRS ~ Apr * ALR + I(ALR^2) + (1 | Lok)
> Data: LRStab
> AIC BIC logLik deviance REMLdev
> 1724 1751 -855 1710 1730
> Random effects:
> Groups Name Variance Std.Dev.
> Lok (Intercept) 4.1578e-09 6.4481e-05
> Residual 6.6801e+00 2.5846e+00
> Number of obs: 361, groups: Lok, 99
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) -1.611300 1.359076 -1.186
> Apr 0.075193 0.205390 0.366
> ALR 2.677924 0.365207 7.333
> I(ALR^2) 0.004942 0.013571 0.364
> Apr:ALR -0.118061 0.049655 -2.378
>
>> .Call("mer_optimize", mod, PACKAGE="lme4")
> NULL
>> mod
> Linear mixed model fit by maximum likelihood
> Formula: LRS ~ Apr * ALR + I(ALR^2) + (1 | Lok)
> Data: LRStab
> AIC BIC logLik deviance REMLdev
> 1723 1750 -854.3 1709 1728
> Random effects:
> Groups Name Variance Std.Dev.
> Lok (Intercept) 0.33923 0.58243
> Residual 6.34722 2.51937
> Number of obs: 361, groups: Lok, 99
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) -1.620101 1.355985 -1.195
> Apr 0.074193 0.204893 0.362
> ALR 2.714618 0.364657 7.444
> I(ALR^2) 0.004308 0.013566 0.318
> Apr:ALR -0.121868 0.049614 -2.456
>
>
>> sessionInfo()
> R version 2.14.2 (2012-02-29)
> Platform: i386-pc-mingw32/i386 (32-bit)
>
> locale:
> [1] LC_COLLATE=Norwegian (Bokmål)_Norway.1252 LC_CTYPE=Norwegian
> (Bokmål)_Norway.1252
> [3] LC_MONETARY=Norwegian (Bokmål)_Norway.1252 LC_NUMERIC=C
> [5] LC_TIME=Norwegian (Bokmål)_Norway.1252
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] arm_1.5-03 foreign_0.8-49 abind_1.4-0
> R2WinBUGS_2.1-18 coda_0.14-6
> [6] car_2.0-12 nnet_7.3-1 MASS_7.3-17
> languageR_1.4 lme4_0.999375-42
> [11] Matrix_1.0-4 lattice_0.20-0
>
> loaded via a namespace (and not attached):
> [1] grid_2.14.2 nlme_3.1-103 stats4_2.14.2
>>
> ########################################################
>
> --------------------------------------------------
> Ivar Herfindal (PhD, Researcher)
> Centre for Conservation Biology
> Department of Biology
> Norwegian University for Science and Technology
> N-7491 Trondheim, Norway
>
> email: Ivar.Herfindal at bio.ntnu.no
> web: http://www.ntnu.no/employees/ivar.herfindal
> phone: +47 73596253
> Cellphone: +47 91625333
> fax: +47 73596090
> --------------------------------------------------
>
> [[alternative HTML version deleted]]
>
>
>
> ------------------------------
>
> Message: 2
> Date: Fri, 30 Nov 2012 19:33:41 +0000
> From: chico3 at sapo.pt
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] glmer, split-plot design simple question
> Message-ID:
> <20121130193341.Horde.p_32ZuanEAJQuQoVyXuT4aA at mail.sapo.pt>
> Content-Type: text/plain; charset=UTF-8; format=flowed; DelSp=Yes
>
>
> Dear experts,
>
> I have a question concerning the implementation of a mixed model in my
> specific experiment. I have search the archives and found no answer to
> my problem.
>
> I have four luminaires randomly distributed. Two luminaires have light
> on and two luminaires have no light. Below the luminaires I have two
> treatments. Each treatment as four replicates that were randomly
> distributed below the luminaires. Moreover, I have made all the
> possible combinations of presence absence of light and treatments.
>
> I want to check if there is a effect of each treatment and the
> interactions between all the combinations.The response variable is a
> proportion (number of specific specie/total number of species)
>
> Treatment1<-as.factor(?Yes?,?No?)
> Treatment2<- as.factor(?Yes?,?No?)
> Light<- as.factor(?Yes?,?No?)
> Response<-cbind(number_of_specific_species, total_species)
>
> I have made this model,
>
> model<-glmer(Response ~ Treatmen1*Treatment2 + (1|Light), family=binomial)
>
> However this doesn?t allow to see the effect of light or the
> interaction between light and treatments. Moreover, I don?t know how
> to include overdispersion, since glmer doesn?t allow quasi families
> such as glm.
>
> I know this is a simple question, but I greatly would appreciate some
> clues on how to proceed.
>
>
>
> ------------------------------
>
> Message: 3
> Date: Sat, 1 Dec 2012 02:30:07 +0000 (UTC)
> From: Ben Bolker <bbolker at gmail.com>
> To: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] glmer, split-plot design simple question
> Message-ID: <loom.20121201T032348-561 at post.gmane.org>
> Content-Type: text/plain; charset=utf-8
>
> <chico3 at ...> writes:
>
>>
>
> [snip]
>
>> I have four luminaires randomly distributed. Two luminaires have light
>> on and two luminaires have no light. Below the luminaires I have two
>> treatments. Each treatment has four replicates that were randomly
>> distributed below the luminaires. Moreover, I have made all the
>> possible combinations of presence absence of light and treatments.
>>
>> I want to check if there is a effect of each treatment and the
>> interactions between all the combinations.The response variable is a
>> proportion (number of specific specie/total number of species)
>>
>> Treatment1<-as.factor(?Yes?,?No?)
>> Treatment2<- as.factor(?Yes?,?No?)
>> Light<- as.factor(?Yes?,?No?)
>> Response<-cbind(number_of_specific_species, total_species)
>>
>> I have made this model,
>>
>> model<-glmer(Response ~ Treatment1*Treatment2 + (1|Light),
>> family=binomial)
>>
>> However this doesn?t allow to see the effect of light or the
>> interaction between light and treatments. Moreover, I don?t know how
>> to include overdispersion, since glmer doesn?t allow quasi families
>> such as glm.
>>
>> I know this is a simple question, but I greatly would appreciate some
>> clues on how to proceed.
>
> I may be missing something here, but it doesn't make sense to me
> to treat Light as a random effect (one criterion -- the levels "Light"
> and "No light" aren't *exchangeable*, i.e. you couldn't switch
> the factor labels without changing the meaning of the experiment).
> I would say you are looking for something more like
>
> glmer(Response ~ Treatment1*Treatment2*Light + (1|luminaire),
> family=binomial)
>
> where the luminaires are labeled uniquely (e.g.
> "Light1","Light2","Nolight1","Nolight2")
>
> It sounds like this is a randomized block design (i.e. every level of
> Treat1*Treat2 is present at each luminaire), so you could test the
> interaction between Treat1*Treat2 and luminaire via
> (Treat1*Treat2|luminaire) instead of (1|luminaire) (see e.g.
> Schielzeth 2009), although you may find that this leads to estimates
> of zero variance (overfitting the data).
>
>
>
> ------------------------------
>
> _______________________________________________
> R-sig-mixed-models mailing list
> R-sig-mixed-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
> End of R-sig-mixed-models Digest, Vol 72, Issue 1
> *************************************************
More information about the R-sig-mixed-models
mailing list