
   Dear Ben:

   Thank you so, so much for your quick, extense and certainly clarifying reply. I think I have some concepts definitely clearer now, although I'm glad to see that I wasn't so mistaken in the fact that it's more correct to use a model comparison approach whenever it's possible. 
   I already knew about the difficulty of the analysis I'm attempting with the data set I have, but the problem is that, according to the type of hypothesis I am trying to test, almost all the possible interactions between the variables are quite important and can have an important biological significance. I suspected the problem with the df's, as before realizing about the convenience of removing the 3 and 4-way interactions (for the impossibility of getting a rational biological explanation for them) I used to get an error message which, according to the local colleagues I asked to, was due to the fact that the model "ran out of df's" trying to estimate so many parameters with such a limited number of observations for a nested random effects model.
   Nevertheless, I guess I will have to manage to find out which of the interactions are really biologically important to reply to my question, and remove the rest of them, so that the df problem with the parameter estimation disappears. Anyway, I think I can do that using lmer instead.
   So, if I got it right, and correct me if I'm wrong, I can use the Wald test p-values encoded in summary() from the glmmPQL function, but being *very* careful about the fact that that's not the very best possible test for that type of data, although it can still be a good approximation. And of course, after trying to write a simpler model to avoid the issue with the "df-observations ratio".  

   I wouldn't like to impose on you at all, but there's another issue I would like to find an answer to, if it was possible. When trying to fit the same model with lmer, I get a warning like this:
"Warning message: In mer_finalize(ans) : false convergence (8)"
   and when calling "summary()", this error message:

   "in asMethod(object) : matrix is not symmetric [1,2]"

   I have posted it before and no one seemed to know the reason very well. So, sorry for the insistence, but actually, after trying all the possible options in terms of removing interactions, I realized the term which is causing the problem is "briventral:tlength", this is, the interaction between the two possible sexual signals I'm testing, tail streamer length and ventral colouration, an important interaction for testing my hypothesis. When not including it in the model, I don't get the warning or the error message at all, and I cannot include it in the model until I only have 6 terms left. I can have an idea of the reason for the warning thanks to this post:
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q1/003306.html
   but I still don't know what that error message about the matrix not being symmetric means. 

   Anyway, thank you VERY much again for your helpful and kind reply. 

   Iker Vaquero

  

--- El dom, 23/1/11, Ben Bolker <bbolker@gmail.com> escribió:

De: Ben Bolker <bbolker@gmail.com>
Asunto: Re: [R-sig-ME] glmmPQL simplification
Para: "Iker Vaquero Alba" <karraspito@yahoo.es>
CC: "Toby Marthews" <toby.marthews@ouce.ox.ac.uk>, "r-sig-mixed-models@r-project.org" <r-sig-mixed-models@r-project.org>
Fecha: domingo, 23 de enero, 2011 21:04

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@ouce.ox.ac.uk>
> escribiÃ³:
> 
> De: Toby Marthews <toby.marthews@ouce.ox.ac.uk> Asunto: RE:
> [R-sig-ME] glmmPQL simplification Para: "Iker Vaquero Alba"
> <karraspito@yahoo.es> CC: "r-sig-mixed-models@r-project.org"
> <r-sig-mixed-models@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@r-project.org
> [r-sig-mixed-models-bounces@ Sent: 22 January 2011 17:50 To:
> r-sig-mixed-models@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@r-project.org mailing list 
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models




      
	[[alternative HTML version deleted]]


