Skip to content

waldtest and nested models - poolability (parameter stability)

4 messages · Roberto Patuelli, Achim Zeileis

#
Dear All,

I'm trying to use waldtest to test poolability (parameter stability) between 
two logistic regressions. Because I need to use robust standard errors 
(using sandwich), I cannot use anova. anova has no problems running the 
test, but waldtest does, indipendently of specifying vcov or not. waldtest 
does not appear to see that my models are nested. H0 in my case is the the 
vector of regression parameters beta1 is the same as the vector of 
parameters beta2, where beta1 and beta2 are computed for the two subgroups 
(divided according to a factor).

I was wondering if anyone can help me making waldtest recognize the nesting.

Here's the lines I run:
(BTW, I try to use robust standard errors because what I normally use 
(glm.binomial.disp) to correct for overdispersion does not converge for the 
unpooled model. But this is another story....)

# poolability for leva.fin03.d
# pooled model
inv.log.leva.base = glm(mix.au.bin ~ cat.gap.tot + leva.fin03.d + ... + 
sud0nord1, data = inv.sub.au, family = binomial, maxit = 1000) # I deleted 
almost all variables to make the line more readable

# overdispersed pooled model - NOT THE PROBLEM NOW
inv.log.leva.base.disp = glm.binomial.disp(inv.log.leva.base)

# unpooled model
inv.log.leva = glm(mix.au.bin ~ leva.fin03.d/(cat.gap.tot + ... + 
sud0nord1 - 1), data = inv.sub.au, family = binomial, maxit = 1000) # again 
I deleted most variables for readability

# overdispersed unpooled model - NOT CONVERGING :(
inv.log.leva.disp = glm.binomial.disp(inv.log.leva, maxit = 10000)

# inv.log.leva.disp not converging, so I resort to using waldtest with 
sandwich, BUT IT DOES NOT SEE THE NESTING!
waldtest(inv.log.leva.base, inv.log.leva, test = "Chisq", vcov = sandwich)

# anova would work but its results are not reliable because of the 
overdispersion - basically without correcting for overdispersion almost 
every variable is highly significant
anova(inv.log.leva.base, inv.log.leva, test = "Chisq")


Thanks everyone!
Best regards,
Roberto Patuelli

********************
Roberto Patuelli, Ph.D.
Istituto Ricerche Economiche (IRE) (Institute for Economic Research)
Universit? della Svizzera Italiana (University of Lugano)
#
On Mon, 6 Dec 2010, Roberto Patuelli wrote:

            
Yes. Because waldtest() needs to figure out which contrasts to apply to 
go from the unrestricted model to the restricted model. The current 
implementation can only do so by looking at the names of the coefficients. 
It assumes that unrestricted model has all coefficients from the 
restricted model plus some more (which are set to zero under the null 
hypothesis).

When you use interactions (as you do below), this only works if you use 
the *-coding but not the /-coding.

In pseudo code:

fm0  <- glm(y ~ x, family = binomial)
fm1a <- glm(y ~ a * x, family = binomial)
fm1b <- glm(y ~ a / x, family = binomial)

The restricted model is fm0 and the unrestricted model is fm1a/fm1b. Both 
are equivalent in terms of fitted values. With waldtest() you can compare

   waldtest(fm0, fm1a)

but

   waldtest(fm0, fm1b)

fails because the models do not fulfill the restriction above. So, only 
for the inference with waldtest() you need to compute fm1a as well. If 
significant, you can go on and interpret fm1b.

Hope that helps,
Z
#
Dear Achim,

Thanks a lot for the superquick reply!
Somehow with your suggestion I can get around the problem, but of course I 
run into other problems, such as this:
Error in bread. %*% meat. : non-conformable arguments

Cheers
Roberto

----- Original Message ----- 
From: "Achim Zeileis" <Achim.Zeileis at uibk.ac.at>
To: "Patuelli Roberto" <roberto.patuelli at usi.ch>
Cc: <r-help at r-project.org>
Sent: Monday, December 06, 2010 9:38 PM
Subject: Re: [R] waldtest and nested models - poolability (parameter 
stability)
On Mon, 6 Dec 2010, Roberto Patuelli wrote:

            
Yes. Because waldtest() needs to figure out which contrasts to apply to
go from the unrestricted model to the restricted model. The current
implementation can only do so by looking at the names of the coefficients.
It assumes that unrestricted model has all coefficients from the
restricted model plus some more (which are set to zero under the null
hypothesis).

When you use interactions (as you do below), this only works if you use
the *-coding but not the /-coding.

In pseudo code:

fm0  <- glm(y ~ x, family = binomial)
fm1a <- glm(y ~ a * x, family = binomial)
fm1b <- glm(y ~ a / x, family = binomial)

The restricted model is fm0 and the unrestricted model is fm1a/fm1b. Both
are equivalent in terms of fitted values. With waldtest() you can compare

   waldtest(fm0, fm1a)

but

   waldtest(fm0, fm1b)

fails because the models do not fulfill the restriction above. So, only
for the inference with waldtest() you need to compute fm1a as well. If
significant, you can go on and interpret fm1b.

Hope that helps,
Z
#
On Mon, 6 Dec 2010, Roberto Patuelli wrote:

            
Hmmm, not sure where this comes from. My guess is that it's not 
waldtest() but sandwich().

Try
   sandwich(inv.log.leva.base)
   sandwich(inv.log.leva)
to see which one of the two is the culprit. Is there something else 
non-standard about the fitted model objects that I missed?

If not, try to boil this down to a smaller example that still has the 
problem and try to provide a reproducible version of it.

Best,
Z