[R-sig-ME] fitting mixed effects logistic regression with weights

Thu Jun 9 10:46:41 CEST 2011

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 11-06-10 05:42 PM, Gebregziabher, Mulugeta wrote:
> Dear Ben, Let me state my questions if I could get more definitive
> answers on these problems that I have been encountering recently:
> 
> 1. is it currently possible to fit mixed effects logit model for
> survey sampled data in R as in Proc GLIMMIX?
> 

  PROC GLIMMIX is equivalent to glmmPQL in the MASS package in R:
however, as far as I know glmmPQL does not support sampling weights.
You should also be aware that there are issues with penalized
quasi/marginal likelihood for binary data (see e.g. Breslow 2003),
although there is some debate as to whether these problems have been
overemphasized (see Murray et al 2004 ref, below).

   In older versions, (g)lmer allowed a PQL option, but Doug Bates (the
author) found that it was sufficiently unreliable for the kinds of
problems he was working on that he removed the option.

> 2. Does MCMCglmm support inclusion of sampling weights for binomial
> data? The problem with lmer is that the responses become non integer
> with the weights and it does not seem to like it.

  No, MCMCglmm does not (again as far as I know).
  Does lmer actually give an error, or does it just give a warning about
"non-integer number of successes in binomial ..."?  It may be OK, if you
know what you're doing ...  lmer is also likely to be your best bet for
doing mixed model analyses of very large data.
> 
> 3. How do we specify the parametric bootstrap for my model given
> below? I am thinking that this will provide an alternative SE/CI
> estimates to the REML/RSPL. For my large data (n>800K) this may not
> be useful, actually.

  I think for n>800K you are going to have a hard/slow time with almost
any approach ... n>800K seems to be the realm of "big data"
algorithms/approaches like biglm, bigglm.

  If your data set is big enough it may be that the Wald approximation
is pretty good (this may be what you're saying).
> 
> 4. How do we account for the sampling weights in the computation of
> the HPD intervals?

   Hmm.  Not sure here.

==============

@article{murray_design_2004,
	title = {Design and Analysis of {Group-Randomized} Trials: A Review of
Recent Methodological Developments},
	volume = {94},
	shorttitle = {Design and Analysis of {Group-Randomized} Trials},
	url = {http://ajph.aphapublications.org/cgi/content/abstract/94/3/423},
	doi = {<p>10.2105/AJPH.94.3.423</p>},
	abstract = {We review recent developments in the design and analysis of
group-randomized trials {(GRTs).} Regarding design, we summarize
developments in estimates of intraclass correlation, power analysis,
matched designs, designs involving one group per condition, and designs
in which individuals are randomized to receive treatments in groups.
Regarding analysis, we summarize developments in marginal and
conditional models, the sandwich estimator, model-based estimators,
binary data, survival analysis, randomization tests, survey methods,
latent variable methods and nonlinear mixed models, time series methods,
global tests for multiple endpoints, mediation effects, missing data,
trial reporting, and software. We encourage investigators who conduct
{GRTs} to become familiar with these developments and to collaborate
with methodologists who can strengthen the design and analysis of their
trials.},
	number = {3},
	journal = {Am J Public Health},
	author = {David M. Murray and Sherri P. Varnell and Jonathan L. Blitstein},
	month = mar,
	year = {2004},
	pages = {423--432}
}
> 
> Thank you!
> 
> ________________________________________ From: Ben Bolker
> [bbolker at gmail.com] Sent: Friday, June 10, 2011 4:45 PM To:
> Gebregziabher, Mulugeta Cc: r-sig-mixed-models at r-project.org Subject:
> Re: [R-sig-ME] fitting mixed effects logistic regression with
> weights
> 
> On 06/10/2011 04:44 PM, Gebregziabher, Mulugeta wrote:
>> Dear Ben,
> 
>> Thank you for your prompt response.  Does MCMCglmm provide the
>> usual REML/RSPL estimates or Bayesian estimates for the fixed
>> effects? I mean are the coefficient and SE estimates REML/RSPL or
>> Bayesian? I want to be sure what I am getting.
> 
> Bayesian.
> 
> 
>> Mulugeta ________________________________________ From: 
>> r-sig-mixed-models-bounces at r-project.org 
>> [r-sig-mixed-models-bounces at r-project.org] On Behalf Of Ben Bolker 
>> [bbolker at gmail.com] Sent: Friday, June 10, 2011 4:33 PM To: 
>> r-sig-mixed-models at r-project.org Subject: Re: [R-sig-ME] fitting 
>> mixed effects logistic regression with  weights
> 
>> On 06/10/2011 04:24 PM, Gebregziabher, Mulugeta wrote:
>>> Dear all,
>>> 
>>> I have been trying to fit a logistic regression model with
>>> random intercept for survey sampled data using the lmer in the
>>> lme4 package. I also want to get HPD intervals. I tried
>>> family=binomial as well as family=quasibinomial. I am
>>> encountering the following errors. "Eror in .local(object, n,
>>> verbose, ...) : Update not yet written" and "Error in HPD
>>> interval(mcmcsamp(fit, 1000)) :   error in evaluating the
>>> argument 'object' in selecting a method for function
>>> 'HPDinterval'.
> 
>> As the error message suggests (obliquely!), mcmcsamp() does not
>> yet work for GLMMs (i.e. [g]lmer with 'family' specified).
>> Furthermore, quasi-likelihood fitting no longer works with the
>> latest releases of lme4 (because the author decided that he didn't
>> really understand what it was doing, hence safer to omit it).
> 
>> Have you considered MCMCglmm, or parametric bootstrapping 
>> (help("simulate-mer")) ?
> 
> 
>> See more below:
>>> 
>>> Does anyone know how to go around this problem? Thanks.
>>> 
>>>> library(lme4)
>>> Loading required package: Matrix Loading required package:
>>> lattice Attaching package: 'Matrix' The following object(s) are
>>> masked from 'package:base': det Attaching package: 'lme4' The
>>> following object(s) are masked from 'package:stats': AIC
>>>> attach(cohort)
>>>> 
>>>> f  <- a1cge8 ~  (1|id)   +  time  +  nhb  + hispanic  + other
>>>> + male  + mstat  +svcpct   +urban   +  comor1   +comor2  +
>>>> comor3
>>>> 
>>>> 
>>>> est        <- round(slot(summary(fit), "coefs")[,1], digits=2)
>>>> t <- slot(summary(fit), "coefs")[,3] df      <- c(1, 1, 3, 3,
>>>> 3, 1, 1, 1, 1, 3, 3, 3) p.value    <- ifelse(t<0, round(pt(t, 
>>>> df=df), digits=3), round(1-pt(t, df), digits=3) ) se
>>>> <- round(slot(summary(fit), "coefs")[,2], digits=4) ci
>>>> <- HPDinterval(mcmcsamp(fit, 1000))
>>> Error in .local(object, n, verbose, ...) : Update not yet
>>> written Error in HPDinterval(mcmcsamp(fit, 1000)) : error in
>>> evaluating the argument 'object' in selecting a method for
>>> function 'HPDinterval'
>>>> lower <- round(ci$fixef[,1], digits=2)
>>> Error: object 'ci' not found
>>>> upper <- round(ci$fixef[,2], digits=2)
>>> Error: object 'ci' not found 
>>> __________________________________________ Mulugeta
>>> Gebregziabher, PhD Assistant Professor Division of Biostatistics
>>> and Epidemiology 135 cannon St. Suite 303 Charleston, SC 29425
>>> Tel: 843-876-1112; Fax: 843-876-1126 E-mail: gebregz at musc.edu 
>>> _______________________________________________ 
>>> 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
> 
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iEYEARECAAYFAk3wiHEACgkQc5UpGjwzenO3fwCeIoFoopWJWrIfebM+EmwGOec0
ezIAoI1wBcy+28cBRYOLGQzOKV/ExaGc
=yGex
-----END PGP SIGNATURE-----