[R-sig-ME] glmmPQL simplification
Ben Bolker
bbolker at gmail.com
Sun Jan 23 21:04:33 CET 2011
On 11-01-23 02:12 PM, Iker Vaquero Alba wrote:
>
> Â Â Hello, Toby:
>
> Â Â Thank you very much for your reply and useful comments.
  Are
> you sure t-values from "summary()" can be used to test model
> significance? I thought you couldn't use those values as they only
> give you information about wether the linear slope with each variable
> is significantly different from 0, or at least that's the case in lme
> (and of course, in lmer, which doesn't even return any p-values when
> calling "summary()"). That's why you use model simplification and
> anova between the simplified and the complete model to test the
> significance of each term. Actually, I was about to reply to your 20
> January post with this information (you say there that you can't use
> lmer as you need to get p-values for your fits and you won't get
> those from lmer). Â Â So, do you know of any equivalent way in
> glmmPQLÂ to this model simplification you can carry out in lme and
> lmer, or will I have too look at the p-values from the t-tests
> (although I think this is not correct)?
There are a few distinctions that you might be confusing here:
The main distinction is between marginal tests of parameters
(typically Wald tests), as they show up in summary() vs. nested model
comparisons (either F or likelihood ratio ('test="Chisq"') as done by
anova() or drop1()
For linear models of continuous covariates, or for factors with only
two levels, in linear models (i.e. from lm()), at the top level of
marginality (i.e. without higher-level interactions present), these are
identical. In most other cases (glm(), lmer(), they are not): when you
have a choice, the model comparison approach is usually better (although
less convenient), for two reasons: (a) the approximations involved in
the test may be more accurate (Wald vs LRT) and (b) it is easier to be
precise about the actual hypothesis you're testing (continuing confusion
about marginality, "type III SS", etc.). The main advantage of the
marginal tests is convenience, and applicability in some cases where the
other (generally preferred) statistical tests such as LRT cannot be
applied).
For factors (categorical predictors) with more than two levels, some
sort of model comparison approach is generally necessary, because the
hypothesis you want to test (stated in terms of parameter values) is
that several parameters are jointly equal to zero. anova() and drop1()
can do this, although the wald.test() function in the aod package can do
a multi-variable equivalent of the Wald tests typically encoded in
summary().
The other issue is the dreaded degrees-of-freedom problem. The general
question is whether one can do a finite-sample correction of the test
statistic, and if so whether there exists a sensible framework and a
sensible procedure for estimating the effective amount of information in
the data set (and in some cases also the effective difference in amount
of information used by two different models).
For glmmPQL, summary() shows Wald t-tests based on the degrees of
freedom estimated in lme() via a "between-within" method. It's not clear
how coherent this is (John Maindonald or Brian Ripley might comment),
although it should be a *reasonable* order-of-magnitude estimate of the
effects of different parameters.
Likelihood ratio tests are not possible because there is no
likelihood defined for a model fitted via quasi-likeihood, although *in
principle* one might be able to come up with a "quasi-" test analogous
to that used by anova() for models fitted with quasi-likelihood in glm().
Because it depends on the likelihood, AIC in the strict sense is also
not available for glmmPQL. For the simpler case of GLMs fitted with
quasi-likelihood, researchers have proposed a quasi-AIC (QAIC) that may
(?) have similar properties, but it hasn't been implemented for glmmPQL
(and would take some careful thinking to see how it extended to this
slightly different case).
I personally think it would be a reasonable philosophical stance to
say of GLMM fits via PQL, "this is an approximation to a well-defined
quantity (likelihood), here is my estimate in this case", but that's not
what has been done.
Bottom line: I think you are stuck with marginal tests in this case.
Please don't disregard all the very important warnings given by others
about whether stepwise regression is appropriate, and whether you can
reasonably expect to do the analysis you are attempting (estimate 15
parameters from a data set that has a total of 48 observations, and
based on the reported df from your fit, may have far fewer 'effective
degrees of freedom').
good luck,
Ben Bolker
>
> Â Â Thank you very much.
>
> Â Â Iker Vaquero
>
> --- El dom, 23/1/11, Toby Marthews <toby.marthews at ouce.ox.ac.uk>
> escribió:
>
> De: Toby Marthews <toby.marthews at ouce.ox.ac.uk> Asunto: RE:
> [R-sig-ME] glmmPQL simplification Para: "Iker Vaquero Alba"
> <karraspito at yahoo.es> CC: "r-sig-mixed-models at r-project.org"
> <r-sig-mixed-models at r-project.org> Fecha: domingo, 23 de enero, 2011
> 02:57
>
> Hi Iker, Not sure you're approaching things the right way here: (i)
> Brian Ripley stated once “I have no idea where you got the idea
> that anova could be used with glmmPQL†(on
> http://www.biostat.wustl.edu/archives/html/s-news/2002-06/msg00055.html)
> which would be a short answer to this question. Also see
> http://www.mail-archive.com/r-help@stat.math.ethz.ch/msg46894.html .
> (ii) However, people do get around this in various ways, e.g. by
> changing the class of the glmmPQL object returned, though note that
> this is probably inadvisable:
> http://markmail.org/message/jddj6aq66wdidrog#query:how%20to%20use%20anova%20with%20glmmPQL+page:1+mid:jddj6aq66wdidrog+state:results
>
>
(iii) Lastly, if you're trying to do model simplification, the p-values
are already there in the summary table you posted. Why not remove the
least significant term (starting with sexM:weatherpc1 not
weatherpc1:tlength) and look at the new summary table rather than that
anova() test?
> Toby
>
> ________________________________________ From:
> r-sig-mixed-models-bounces at r-project.org
> [r-sig-mixed-models-bounces@ Sent: 22 January 2011 17:50 To:
> r-sig-mixed-models at r-project.org Subject: [R-sig-ME] glmmPQL
> simplification
>
> Â Â Â Hello:
>
> Â Â Â I am trying to fit a mixed effects model with the glmmPQL
> function, as I get an error message trying to fit the same model with
> lmer.
>
> Â Â Â When I try to simplify the model, removing some of the
> interactions, I get this message:
>
> Â Â "Error in anova.glmmPQL(fledgecoltailmodel1,
> s.fledgecoltailmodel1.1) : Â 'anova' is not available for PQL fits"
>
> Â Â Â So, I don't know what to do or what function to use to reliably
> simplify the model and drop some of the interaction terms. Can you
> give me any indication? Any little help would be greatly
> appreciated.
>
> Â Â Â >
> fledgecoltailmodel1<-glmmPQL(nfledge~sex*briventral*inslarge*weatherpc1*tlength-sex:briventral:inslarge:weatherpc1:tlength-sex:briventral:inslarge:weatherpc1-sex:briventral:inslarge:tlength-sex:briventral:weatherpc1:tlength-sex:inslarge:weatherpc1:tlength-briventral:inslarge:weatherpc1:tlength-sex:briventral:inslarge-sex:briventral:weatherpc1-sex:briventral:tlength-sex:inslarge:weatherpc1-sex:inslarge:tlength-sex:weatherpc1:tlength-briventral:inslarge:weatherpc1-briventral:inslarge:tlength-briventral:weatherpc1:tlength-inslarge:weatherpc1:tlength,random=~1|site/pair,family=poisson)
>
>
iteration 1
> iteration 2 iteration 3
>> summary(fledgecoltailmodel1)
> Linear mixed-effects model fit by maximum likelihood Data: NULL Â
> AIC BIC logLik    NA NA    NA
>
> Random effects: Formula: ~1 | site       (Intercept) StdDev:
> 3.425823e-07
>
> Formula: ~1 | pair %in% site     (Intercept)   Â
> Residual StdDev:Â Â Â 0.2358495 2.329028e-07
>
> Variance function: Structure: fixed weights Formula: ~invwt Fixed
> effects: nfledge ~ sex * briventral * inslarge * weatherpc1 * tlength
> -Â Â Â sex:briventral:inslarge:weatherpc1:tlength -
> sex:briventral:inslarge:weatherpc1 -Â Â Â
> sex:briventral:inslarge:tlength - sex:briventral:weatherpc1:tlength
> -Â Â Â sex:inslarge:weatherpc1:tlength -
> briventral:inslarge:weatherpc1:tlength -Â Â Â
> sex:briventral:inslarge - sex:briventral:weatherpc1 -
> sex:briventral:tlength -Â Â Â sex:inslarge:weatherpc1 -
> sex:inslarge:tlength - sex:weatherpc1:tlength -Â Â Â
> briventral:inslarge:weatherpc1 - briventral:inslarge:tlength -Â Â Â
> briventral:weatherpc1:tlength - inslarge:weatherpc1:tlength   Â
>             Value Std.Error DF   t-value
> p-value (Intercept)Â Â Â Â Â Â 1.3502758 0.29547025 15Â
> 4.569921Â 0.0004 sexMÂ Â Â Â Â Â Â Â Â Â Â 0.0000013
> 0.00000137 11 0.913264 0.3807 briventral     Â
> -0.0000001 0.00000002 11 -2.154738 0.0542 inslarge     Â
> Â Â Â 0.0179169 0.16643724Â 6Â 0.107649Â 0.9178 weatherpc1Â Â Â
>      0.0192141 0.21848655 6 0.087942 0.9328 tlength Â
> Â Â Â Â Â Â 0.0000001 0.00000004 11Â 1.901240Â 0.0838
> sexM:briventral    0.0000000 0.00000000 11 -2.486618Â
> 0.0302 sexM:inslarge     0.0000009 0.00000029 11Â
> 3.284681 0.0073 briventral:inslarge  0.0000000 0.00000001 11Â
> 2.072806Â 0.0625 sexM:weatherpc1Â Â Â Â 0.0000001 0.00000010 11Â
> 0.639981Â 0.5353 briventral:weatherpc1Â 0.0000000 0.00000000 11Â
> 2.295970Â 0.0423 inslarge:weatherpc1Â Â Â -0.0589347 0.13022509Â 6
> -0.452560 0.6668 sexM:tlength       0.0000000 0.00000001
> 11 -1.594436 0.1391 briventral:tlength    0.0000000 0.00000000
> 11 1.142253 0.2776 inslarge:tlength     0.0000000
> 0.00000001 11 -3.174971 0.0088 weatherpc1:tlength    0.0000000
> 0.00000001 11 -0.653040Â 0.5271 Correlation: Â Â Â Â Â Â Â Â
> Â Â Â (Intr) sexMÂ Â Â brvntr inslrg wthrp1 tlngth sxM:br sxM:ns
> sexM           0.000 briventral      Â
>  0.000 0.224 inslarge       -0.980 0.000 0.000
> weatherpc1Â Â Â Â Â Â Â Â 0.811Â 0.000Â 0.000 -0.817 tlengthÂ
> Â Â Â Â Â Â Â 0.000 -0.088Â 0.637Â 0.000Â 0.000
> sexM:briventral    0.000 -0.230 0.748 0.000 0.000Â
> 0.588 sexM:inslarge     0.000 -0.432 -0.054 0.000Â
> 0.000 0.644 -0.001 briventral:inslarge  0.000 -0.443 -0.406Â
> 0.000Â 0.000Â 0.133Â 0.026Â 0.317 sexM:weatherpc1Â Â Â Â
> 0.000Â 0.426Â 0.239Â 0.000Â 0.000 -0.295Â 0.149 -0.706
> briventral:weatherpc1Â 0.000Â 0.203Â 0.021Â 0.000Â 0.000 -0.298
> -0.183 -0.253 inslarge:weatherpc1Â Â Â -0.815Â 0.000Â 0.000Â 0.826
> -0.979 0.000 0.000 0.000 sexM:tlength       0.000
> -0.895 -0.469Â 0.000Â 0.000 -0.341 -0.062Â 0.107
> briventral:tlength    0.000 0.049 -0.889 0.000 0.000
> -0.766 -0.874 -0.117 inslarge:tlength     0.000 0.403
> -0.039 0.000 0.000 -0.747 -0.117 -0.952 weatherpc1:tlength  Â
> Â 0.000 -0.414 -0.350Â 0.000Â 0.000Â 0.039 -0.270Â 0.518 Â Â Â
> Â Â Â Â Â Â Â Â brvntrl:n sxM:w1 brvn:1 insl:1 sxM:tl
> brvntrl:t inslr: sexM briventral inslarge weatherpc1 tlength
> sexM:briventral sexM:inslarge briventral:inslarge sexM:weatherpc1Â Â
> Â Â Â -0.265 briventral:weatherpc1 -0.453Â Â Â Â 0.282
> inslarge:weatherpc1  0.000    0.000 0.000 sexM:tlength Â
>      0.373  -0.242 -0.060 0.000 briventral:tlength Â
> -0.039  -0.097 0.198 0.000 0.268 inslarge:tlength  Â
> -0.306Â Â Â Â 0.615Â 0.305Â 0.000 -0.054Â 0.213
> weatherpc1:tlength    0.310  -0.831 -0.017 0.000 0.352Â
> 0.219Â Â -0.369
>
> Standardized Within-Group Residuals:      Min    Â
>  Q1       Med      Q3       Max
> -1.980407e+00 -3.989934e-01 -6.236178e-07Â 3.989821e-01Â
> 1.979438e+00
>
> Number of Observations: 48 Number of Groups: Â Â Â Â Â site pair
> %in% site       10        25
>> s.fledgecoltailmodel1.1<-update(fledgecoltailmodel1,~.-weatherpc1:tlength)
>
>>
iteration 1
> iteration 2 iteration 3
>> anova(fledgecoltailmodel1,s.fledgecoltailmodel1.1)
> Error in anova.glmmPQL(fledgecoltailmodel1, s.fledgecoltailmodel1.1)
> : Â 'anova' is not available for PQL fits
>
> [[elided Yahoo spam]]
>
>
>
> Â Â Â Iker Vaquero-Alba
>
> Â Â Â Centre for Ecology and Conservation
>
> Â Â Â Daphne du Maurier Building
>
> Â Â Â University of Exeter, Cornwall Campus
>
> Â Â Â Treliever Road
>
> Â Â Â TR10 9EZ Penryn
>
> Â Â Â U.K.
>
>
>
>
>
>
>
> Â Â Â Â [[alternative HTML version deleted]]
>
>
>
>
>
> [[alternative HTML version deleted]]
>
>
>
>
> _______________________________________________
> 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