[R-sig-ME] Extractor Functions
John Maindonald
john.maindonald at anu.edu.au
Sat Aug 13 04:37:15 CEST 2011
Martin, here is a first request!
The following type of plot (cf Figure 10.3 in "Data Analysis and Graphics Using R", 2nd & 3rd editions)
often makes sense for a multilevel model:
library(DAAG)
science <- na.omit(science)
science$class <- factor(science$class)
science$school <- factor(science$school)
science1.lmer <- lmer(like ~ sex + PrivPub + (1 | school:class),
data = science)
ranf <- ranef(obj = science1.lmer, drop=TRUE)[["school:class"]]
ranflev <- science1.lmer at flist[["school:class"]]
privpub <- science[match(names(ranf), ranflev), "PrivPub"]
count <- unclass(table(ranflev))
plot(sqrt(count), ranf, xaxt="n", pch=as.numeric(privpub),
xlab="# in class (square root scale)",
ylab="Estimate of class effect")
Now define the following function:
ranefWithFactors <-
function(obj=science1.lmer, term="school:class", fac="PrivPub",
count=TRUE, data=science){
ranf <- ranef(obj = obj, drop=TRUE)[[term]]
ranflev <- obj at flist[[term]]
## Check that each value (level) corresponds to a unique level of PrivPub
if(count) num <- unclass(table(ranflev)) else num <- NULL
matchfac <- data[match(names(ranf), ranflev), fac]
check <- apply(table(ranflev, data[,fac]),1,function(x)sum(x==0))
if(any(check)!=1)stop(paste("Each value (level) of", term,
"must correspond to a unique level of", fac))
df <- data.frame(effect=ranf, fac=matchfac, count=num)
df
}
library(DAAG)
science <- na.omit(science)
science$class <- factor(science$class)
science$school <- factor(science$school)
science1.lmer <- lmer(like ~ sex + PrivPub + (1 | school:class), data = science)
df <- ranefWithFactors(obj=science1.lmer, term="school:class", fac="PrivPub",
count=TRUE, data=science)
xyplot(effect ~ sqrt(count), xaxt="n", groups=fac, data=df,
auto.key=list(columns=2),
xlab="# in class (square root scale)", ylab="Estimate of class effect")
The function should doubtless allow for 'fac' to be a vector of factor names.
It should be possible to dispense with the argument data -- I've not checked
how to circumvent that.
Thanks for opening up the issue of further extractor functions.
John Maindonald email: john.maindonald at anu.edu.au
phone : +61 2 (6125)3473 fax : +61 2(6125)5549
Centre for Mathematics & Its Applications, Room 1194,
John Dedman Mathematical Sciences Building (Building 27)
Australian National University, Canberra ACT 0200.
http://www.maths.anu.edu.au/~johnm
On 12/08/2011, at 12:11 AM, Martin Maechler wrote:
>>>>>> John Maindonald <john.maindonald at anu.edu.au>
>>>>>> on Thu, 11 Aug 2011 20:54:26 +1000 writes:
>
>> showMethods(class="mer")
>
> yes, indeed, thank you, John.
>
> And...:
> Murray, John, ... <any somewhat experienced *lmer() user> ,
> please do ask for missing "extractor" methods!
> I use quotes because I mean not just extractors in the
> strict sense, but also other sensible (small) utilities.
>
> For *mer objects, I agree that a user should typically never
> have to use str(.) and then access the slots directly...
> (though its clearly necessary now, IIRC).
>
> Notably because current lme4 and future lme4
> result objects will partly have very different slots,
> such that only such "extractor" methods will continue to work
> when you switch from lme4 to "future-lme4"..
>
> Martin Maechler, ETH Zurich
> (digressing from finishing his slides for "useR!" next week...)
>
>
>> John Maindonald email: john.maindonald at anu.edu.au
>> phone : +61 2 (6125)3473 fax : +61 2(6125)5549
>> Centre for Mathematics & Its Applications, Room 1194,
>> John Dedman Mathematical Sciences Building (Building 27)
>> Australian National University, Canberra ACT 0200.
>> http://www.maths.anu.edu.au/~johnm
>
>> On 11/08/2011, at 8:09 PM, Murray Jorgensen wrote:
>
>>> Is there a convenient table or list of available extractor
>>> functions for mer and summary.mer objects in lme4 and lme4a?
>>>
>>> Murray Jorgensen
>>>
>>> On 11/08/11 10:00 PM, r-sig-mixed-models-request at r-project.org
>>> wrote:
>>>> 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. mixed model (Ahmad Rabiee) 2. Re: mixed model (Ben Bolker)
>>>>
>>>>
>>>> ----------------------------------------------------------------------
>>>>
>>>> Message: 1 Date: Thu, 11 Aug 2011 11:53:20 +1000 From: "Ahmad
>>>> Rabiee"<ahmadr at sbscibus.com.au>
>>>> To:<r-sig-mixed-models at r-project.org> Subject: [R-sig-ME]
>>>> mixed model
>>>> Message-ID:<005b01cc57c9$724ca390$56e5eab0$@com.au>
>>>> Content-Type: text/plain
>>>>
>>>> Hi
>>>>
>>>>
>>>>
>>>> I have a binomial dataset (0, 1), and would like to run a
>>>> "mixed model" logistic regression and also a "nested mixed
>>>> model" logistic regression using glmer:
>>>>
>>>>
>>>>
>>>> ket.glm1<- glmer(z_ket_1.4 ~ bcs_pre + bhb_date + lact + (1 |
>>>> herdno) , family = binomial, data = ket)
>>>>
>>>> #-------------------------------
>>>>
>>>>
>>>>
>>>> To account for the overdispersion in the dataset, I used the
>>>> following codes (according to lme4 package), but the output is
>>>> identical to the first model (above= ket.hlm1). Comments
>>>> please?
>>>>
>>>>
>>>>
>>>> # Mixed model accounting for overdispersion
>>>>
>>>> ket$obs<- 1:nrow(ket)
>>>>
>>>> ket.glm2<- glmer(z_ket_1.4 ~ bcs_pre + bhb_date + lact + (1 |
>>>> herdno) + (1|obs), family = binomial, data = ket)
>>>>
>>>> #-------------------------------------------------
>>>>
>>>>
>>>>
>>>> #Nest random effect
>>>>
>>>> When I want to run a nested random effects using "glmer" I get
>>>> an error message (below);
>>>>
>>>>
>>>>
>>>> # herds nested within studies
>>>>
>>>> ket.glm43<- glmer(z_ket_1.4 ~ bcs_pre + bhb_date + lact +
>>>> (1|studyid:herdno) + (1|id), family = binomial, data = ket)
>>>>
>>>>
>>>>
>>>> #Error message (What does this mean?)
>>>>
>>>> Error: length(f1) == length(f2) is not TRUE
>>>>
>>>> In addition: Warning messages:
>>>>
>>>> 1: In study:herdno :
>>>>
>>>> numerical expression has 2695 elements: only the first used
>>>>
>>>> 2: In study:herdno :
>>>>
>>>> numerical expression has 2695 elements: only the first used
>>>>
>>>> #---------------------------------------------------------------------------
>>>> ----------
>>>>
>>>>
>>>>
>>>> #glmmadmb
>>>>
>>>> I believe my dataset (binomial) is zero-inflated- and Ben
>>>> suggested that I should use the "glmmadmb" package to count
>>>> for the zero-inflation (Please correct me if I am wrong). I
>>>> can run this model (below), when I don't have a random effects
>>>> term in the model. But I don't understand the outputs:
>>>>
>>>>
>>>>
>>>> # first model (without random effects term)
>>>>
>>>> ket.glmm1<- glmmadmb(z_ket_1.4 ~ bcs_pre + bhb_date + lact ,
>>>> family = "binomial", data = ket)
>>>>
>>>> summary(ket.glmm2)
>>>>
>>>>
>>>>
>>>> Initial statistics: 10 variables; iteration 0; function
>>>> evaluation 0; phase 1
>>>>
>>>> Function value 1.8680316e+03; maximum gradient component mag
>>>> 1.4283e+01
>>>>
>>>> Var Value Gradient |Var Value Gradient |Var Value Gradient
>>>>
>>>> 1 0.00000 1.42834e+01 | 2 0.00000 -1.25533e-01 | 3 0.00000
>>>> 7.40839e+00
>>>>
>>>> 4 0.00000 9.75790e-01 | 5 0.00000 -1.44553e+00 | 6 0.00000
>>>> -1.77029e+00
>>>>
>>>> 7 0.00000 -1.94537e+00 | 8 0.00000 -1.72752e+00 | 9 0.00000
>>>> -9.14276e-01
>>>>
>>>> 10 0.00000 -1.12861e+00 |
>>>>
>>>>
>>>>
>>>> Intermediate statistics: 10 variables; iteration 10; function
>>>> evaluation 14; phase 1
>>>>
>>>> Function value 1.2444800e+03; maximum gradient component mag
>>>> -4.4890e-02
>>>>
>>>> Var Value Gradient |Var Value Gradient |Var Value Gradient
>>>>
>>>> 1-74.45668 1.75306e-02 | 2 3.38165 1.96037e-02 | 3-39.88270
>>>> -8.55544e-03
>>>>
>>>> 4 -4.77457 -4.48896e-02 | 5 10.47568 -5.10712e-03 | 6 12.37466
>>>> -2.16758e-03
>>>>
>>>> 7 13.04031 3.49481e-03 | 8 10.98531 5.92878e-04 | 9 6.74107
>>>> -4.27903e-03
>>>>
>>>> 10 6.99462 -1.41183e-04 |
>>>>
>>>> 10 variables; iteration 20; function evaluation 24; phase 1
>>>>
>>>> Function value 1.2444697e+03; maximum gradient component mag
>>>> 1.2294e-04
>>>>
>>>> Var Value Gradient |Var Value Gradient |Var Value Gradient
>>>>
>>>> 1-74.57412 6.45335e-05 | 2 3.24452 3.17209e-05 | 3-39.84581
>>>> -2.93776e-05
>>>>
>>>> 4 -4.43904 1.22940e-04 | 5 10.55041 -6.38080e-05 | 6 12.43515
>>>> -6.09376e-05
>>>>
>>>> 7 13.07189 -2.38280e-05 | 8 11.01902 -1.81564e-05 | 9 6.79542
>>>> -2.85983e-05
>>>>
>>>> 10 7.02120 -6.17190e-06 |
>>>>
>>>>
>>>>
>>>> - final statistics:
>>>>
>>>> 10 variables; iteration 21; function evaluation 25
>>>>
>>>> Function value 1.2445e+03; maximum gradient component mag
>>>> 4.5702e-05
>>>>
>>>> Exit code = 1; converg criter 1.0000e-04
>>>>
>>>> Var Value Gradient |Var Value Gradient |Var Value Gradient
>>>>
>>>> 1-74.57453 9.86716e-06 | 2 3.24432 1.00312e-05 | 3-39.84559
>>>> 3.57290e-05
>>>>
>>>> 4 -4.43952 4.57024e-05 | 5 10.55069 -2.73803e-05 | 6 12.43544
>>>> -2.03976e-05
>>>>
>>>> 7 13.07205 -6.19842e-06 | 8 11.01915 -2.28111e-06 | 9 6.79553
>>>> -1.71599e-05
>>>>
>>>> 10 7.02125 -2.24921e-06 |
>>>>
>>>> Estimating row 1 out of 10 for hessian
>>>>
>>>> Estimating row 2 out of 10 for hessian
>>>>
>>>> Estimating row 3 out of 10 for hessian
>>>>
>>>> Estimating row 4 out of 10 for hessian
>>>>
>>>> Estimating row 5 out of 10 for hessian
>>>>
>>>> Estimating row 6 out of 10 for hessian
>>>>
>>>> Estimating row 7 out of 10 for hessian
>>>>
>>>> Estimating row 8 out of 10 for hessian
>>>>
>>>> Estimating row 9 out of 10 for hessian
>>>>
>>>> Estimating row 10 out of 10 for hessian
>>>>
>>>> Estimated covariance matrix may not be positive definite
>>>>
>>>> 4.44173 4.92261 5.06046 5.06419 5.35787 5.45402 5.62318
>>>> 6.84209 8.1491 11.1008
>>>>
>>>> Estimated covariance matrix may not be positive definite
>>>>
>>>> 4.44173 4.92261 5.06046 5.06419 5.35787 5.45402 5.62318
>>>> 6.84209 8.1491 11.1008
>>>>
>>>> #-------------------------------------------
>>>>
>>>>
>>>>
>>>> When I run "glmmadmb" with a random effects term in the model,
>>>> I get an error message. I don't know what I am doing wrong
>>>> here. Any help would be greatly appreciated.
>>>>
>>>>
>>>>
>>>> # Mixed model (herdno is the random effects term)
>>>>
>>>> ket.glmm2<- glmmadmb(z_ket_1.4 ~ bcs_pre + bhb_date + lact +
>>>> (1 | herdno), family = "binomial", data = ket)
>>>>
>>>> summary(ket.glmm2)
>>>>
>>>>
>>>>
>>>> #Error message
>>>>
>>>> Error in process_randformula(formula, random, data = data) :
>>>>
>>>> all grouping variables must be factors
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> Thanks.
>>>>
>>>> Ahmad
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> #---------------------------------------------------------------------------
>>>> ------------
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> "Try not to become a man of success, but rather try to become
>>>> a man of value" Albert Einstein
>>>> <http://www.brainyquote.com/quotes/authors/a/albert_einstein.html>
>>>>
>>>>
>>>>
>>>> Please note my new email address is
>>>> mailto:ahmadr at sbscibus.com.au. Please update your records.
>>>>
>>>>
>>>>
>>>>
>>>> [[alternative HTML version deleted]]
>>>>
>>>>
>>>>
>>>> ------------------------------
>>>>
>>>> Message: 2 Date: Thu, 11 Aug 2011 02:34:52 +0000 (UTC) From:
>>>> Ben Bolker<bbolker at gmail.com> To:
>>>> r-sig-mixed-models at r-project.org Subject: Re: [R-sig-ME] mixed
>>>> model Message-ID:<loom.20110811T042857-460 at post.gmane.org>
>>>> Content-Type: text/plain; charset=us-ascii
>>>>
>>>>> Ahmad Rabiee<ahmadr at ...> writes:
>>>>
>>>> # I have a binomial dataset (0, 1),
>>>>
>>>> this is a key piece of information not stated previously
>>>> (or I missed it ...)
>>>>
>>>>> and would like to run a "mixed model" # logistic regression
>>>> and also a "nested mixed model" logistic regression # using
>>>> glmer:
>>>> #
>>>> # ket.glm1<- glmer(z_ket_1.4 ~ bcs_pre + bhb_date + lact + (1
>>>> | herdno) , # family = binomial, data = ket)
>>>>
>>>> If your data are binomial with values 0/1 (i.e., "binary" or
>>>> "Bernoulli"), it makes sense to incorporate neither
>>>> overdispersion nor zero-inflation.
>>>>
>>>> # To account for the overdispersion in the dataset, I used the
>>>> following codes # (according to lme4 package), but the output
>>>> is identical to the first model # (above= ket.hlm1). Comments
>>>> please?
>>>> #
>>>> # # Mixed model accounting for overdispersion
>>>> #
>>>> # ket$obs<- 1:nrow(ket)
>>>> #
>>>> # ket.glm2<- glmer(z_ket_1.4 ~ bcs_pre + bhb_date + lact + (1
>>>> | herdno) + # (1|obs), family = binomial, data = ket)
>>>>
>>>> As stated above, overdispersion is unidentifiable with binary
>>>> data.
>>>>
>>>> # #Nest random effect # When I want to run a nested random
>>>> effects using "glmer" I get an error # message (below);
>>>> #
>>>> # # herds nested within studies
>>>> #
>>>> # ket.glm43<- glmer(z_ket_1.4 ~ bcs_pre + bhb_date + lact + #
>>>> (1|studyid:herdno) + (1|id), family = binomial, data = ket)
>>>> #
>>>> # #Error message (What does this mean?)
>>>> #
>>>> # Error: length(f1) == length(f2) is not TRUE
>>>> #
>>>> # In addition: Warning messages:
>>>> #
>>>> # 1: In study:herdno :
>>>> #
>>>> # numerical expression has 2695 elements: only the first used
>>>>
>>>> [snip]
>>>>
>>>> It means that you need studyid and herdno to be factors, not
>>>> numeric variables, in order for this to work.
>>>>
>>>> # I believe my dataset (binomial) is zero-inflated- and Ben
>>>> suggested that I # should use the "glmmadmb" package to count
>>>> for the zero-inflation (Please # correct me if I am wrong). I
>>>> can run this model (below), when I don't have a # random
>>>> effects term in the model. But I don't understand the outputs:
>>>>
>>>> When I suggested that, it was before I knew your data were
>>>> binary. Zero-inflation doesn't make sense for binary data.
>>>>
>>>> # # first model (without random effects term)
>>>> #
>>>> # ket.glmm1<- glmmadmb(z_ket_1.4 ~ bcs_pre + bhb_date + lact ,
>>>> family = # "binomial", data = ket)
>>>> #
>>>> # summary(ket.glmm2)
>>>> #
>>>> # Initial statistics: 10 variables; iteration 0; function
>>>> evaluation 0; phase # 1 [snip]
>>>>
>>>> # Estimated covariance matrix may not be positive definite
>>>> #
>>>> # 4.44173 4.92261 5.06046 5.06419 5.35787 5.45402 5.62318
>>>> 6.84209 8.1491 # 11.1008
>>>> #
>>>> # When I run "glmmadmb" with a random effects term in the
>>>> model, I get an # error message. I don't know what I am doing
>>>> wrong here. Any help would be # greatly appreciated.
>>>> #
>>>> # # Mixed model (herdno is the random effects term)
>>>> #
>>>> # ket.glmm2<- glmmadmb(z_ket_1.4 ~ bcs_pre + bhb_date + lact +
>>>> (1 | # herdno), family = "binomial", data = ket)
>>>> #
>>>> # summary(ket.glmm2)
>>>> #
>>>> # #Error message
>>>> #
>>>> # Error in process_randformula(formula, random, data = data) :
>>>> #
>>>> # all grouping variables must be factors
>>>>
>>>> What it says. herdno must be a factor.
>>>>
>>>> Ben Bolker
>>>>
>>>>
>>>>
>>>> ------------------------------
>>>>
>>>> _______________________________________________
>>>> 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 56, Issue 15
>>>> **************************************************
>>>
>>> --
>>> Dr Murray Jorgensen http://www.stats.waikato.ac.nz/Staff/maj.html
>>> Department of Statistics, University of Waikato, Hamilton, New Zealand
>>> Email: maj at waikato.ac.nz majorgensen at ihug.co.nz Fax 7 838 4155
>>> Phone +64 7 838 4773 wk Home +64 7 825 0441 Mobile 021 0200 8350
>>>
>>> _______________________________________________
>>> 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
More information about the R-sig-mixed-models
mailing list