glmmPQL simplification
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 at 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