3
$\begingroup$

There are some related posts on this issue, but no answers actually demonstrate the mechanics of how to accomplish the task that I could find. I want to compare two parameter estimates in a binomial GLM (but I expect the answer to this question will work for any GLM).

My model is of the form:

y ~ b0 + b1x1 + b2x2 + e(binom)

My question: How do you test for a difference between b1 and b2?

Running this model using glm() in R provides parameter estimates and standard errors for all b. So comparing them should be easy, but I'm just having a hard time figuring out exactly what to do.

Here are some example data to work with:

set.seed(1987)
y <- rbinom(n = 100, size = 1, prob = 0.5) # binomial response variable
set.seed(1988)
x1a <- rnorm(n = 100, mean = 2, sd = 3)
set.seed(1988)
x1b <- rnorm(n = 100, mean = 0, sd = 3)
x1 <- ifelse(y == 0, x1a, x1b) # negative relationship between y and x1
set.seed(1990)
x2a <- rnorm(n = 100, mean = 2, sd = 5)
set.seed(1990)
x2b <- rnorm(n = 100, mean = 15, sd = 5)
x2 <- ifelse(y == 0, x2a, x2b) # a strong, positive relationship between y and x2

Run this model and return relevant output:

an1 <- glm(y ~ x1 + x2, family=binomial)
s_an1 <- summary(an1)

s_an1$coefficients 
    paste("residual df = ", an1$df.residual)
paste("null.df = ", an1$df.null)

And get these estimates:

              Estimate Std. Error   z value     Pr(>|z|)
(Intercept) -4.2711210  1.0451565 -4.086585 4.377686e-05
x1          -0.1638442  0.1433595 -1.142891 2.530841e-01
x2           0.6130596  0.1378431  4.447519 8.686793e-06

residual df = 97
null df = 99

Notice that x2 is a significant positive predictor of y, and x1 is non-significant, negative predictor of y. As we all know, however, this is not good evidence that the parameter estimates are different.

(How) Can I use the standard errors to compare the estimates?

I'll point out that there appears to be a lot of information on the web about how to calculate confidence intervals, which are nice, but I want:

1) a parameter estimate for the difference (yes, I know -- it's just the difference between them),

2) a standard error (probably sqrt(se.x1^2 + se.x2^2)),

3) a t-value (or equivalent distribution statistic), and

4) a p-value.

In my real-world scenario, x1 and x2 are measured on the same scale, but don't have the same mean or standard deviations (as in the example data above -- see sd arguments to random number generators). Is standardization of some kind necessary in this case?

Help much appreciated, cheers!

$\endgroup$
2
$\begingroup$

You can simulate parameter estimates using the mvtnorm package with mean vector coef(an1) and covariance matrix vcov(an1) and then summarise them. Or you could bootstrap.

However, it's probably easier just to use the multcomp package to examine the contrast. In your model that would be:

library(multcomp)
cont <- glht(an1, linfct="x1 - x2 = 0")
summary(cont) ## estimate, standard error, z-statistic and p-value
confint(cont) ## confidence interval
$\endgroup$
  • $\begingroup$ That's great thanks. What are the advantages of the multcomp vs simulation vs. bootstrap methods? Can you briefly describe the mechanics of the multcomp method? Is it more or less like comparing two means with a given standard error (e.g., t-test)? $\endgroup$ – tim.farkas Feb 14 '15 at 22:43
  • 1
    $\begingroup$ The theory for testing linear functions of parameters is fairly standard. The package vignette: vignette('generalsiminf', package='multcomp') summarises it and provides references. Your problem of testing whether b1-b2=0 is a special case. The package also takes care of the multiple comparisons issues that can arise in more complex situations - hence its title. The actual mechanics vary depending on the model, but they all basically depend on making some kind of parameteric distributional assumptions and using the model's variance covariance matrix. $\endgroup$ – conjugateprior Feb 15 '15 at 16:57
  • $\begingroup$ Is it more or less like a t-test? Yes. That test is a special case of the general framework. $\endgroup$ – conjugateprior Feb 15 '15 at 16:58
  • $\begingroup$ Bootstrapping may give you more accurate estimate of uncertainty about the contrast than the calculation used in the standard linear contrasts testing approach. Or it may not. This is a large subject that can't usefully be summarised here. Practically speaking, if you want to compare to some bootstrap approaches on your data the rms package may be helpful. $\endgroup$ – conjugateprior Feb 15 '15 at 17:04
0
$\begingroup$

In essence, you are asking to compare the (sampling) distribution of two estimates. Then I think bootstrapping is what you can do.

Basically once you've generated your simulation data set, you use that as your population, do draw random samples with replacement, and fit the GLM, obtain the estimates of the coefficients. By that way, you can obtain bootstrap sampling distribution of the estimates and do the common inference.

Here it seems to me that you imposed an implicit relationship between X's and y. I mean the true model is not well specified. So I am not sure what kind of sampling mechanism is better. Maybe you can try resampling residuals or block resampling.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.