Score tests for marginal regression models fit using GEE

Key ideas: Simulation study, score test, Gaussian GEE

This notebook uses simulation to check the performance of the GEE score test. The score test can be used to test null hypotheses that take the form L*params = R, where L is a q x p matrix, and R is a q-vector. Score tests are especially useful in GEE analyses since in GEE there is no likelihood ratio test.

Here are the import statements:

In [1]:
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy.stats.distributions import norm

Testing two linear constraints

Score tests can be used to test models that are fit under an arbitrary number of linear constraints on the parameter vector against the unconstrained parent model. First we evaluate the performance of score testing when the mean structure is constrained using two linear constraints.

This is the number of replicates in the simulation study. Increase it to get more precise results (but the simulation will take a longer time to run).

In [2]:
nrep = 400

These are the matrices that define the two-dimensional constraint params[1] == params[2] and params[3] == params[4].

In [3]:
L = np.array([[0, 1, -1, 0, 0], [0, 0, 0, 1, -1]])
R = np.array([0, 0])

We will simulate data under the null hypothesis and under the alternative ypothesis. The p-values for the null hypothesis will be placed in pvalues[0] and the p-values for the alternative hypothesis will be placed in pvalues[1].

In [4]:
pvalues = [[], []]

Here is a function to simulate data suitable for GEE analysis.

In [5]:
def gendat(params):
    
    exog = np.random.normal(size=(500,5))
    groups = np.kron(np.arange(100), [1,1,1,1,1])

    endog = np.dot(exog, params) + 2*np.random.normal(size=500)
    
    # The cluster effects
    endog += np.random.normal(size=100)[groups]

    return endog, exog, groups

Here is the main simulation loop that constructs the data sets, fits marginal regression models using GEE, and collects the p-values for the score test.

In [6]:
for j in 0,1:
    for iter in range(nrep):

        # Null hypothesis, the constraints are satisfied
        if j == 0:
            params = np.r_[0, 1, 1, -2, -2]
            
        # Alternative hypothesis, the first constraint is not satisfied
        else:
            params = np.r_[0, 1, 0.8, -2, -2]

        endog, exog, groups = gendat(params)
        
        ex = sm.cov_struct.Exchangeable()
        ga = sm.families.Gaussian()

        model = sm.GEE(endog, exog, groups, family=ga, cov_struct=ex,
                 constraint=(L,R))
        result = model.fit()
        pvalues[j].append(model.score_test_results["p-value"])

These are the results of one run in the simulation study.

In [7]:
print result.summary()
print model.score_test_results
                               GEE Regression Results                              
===================================================================================
Dep. Variable:                           y   No. Observations:                  500
Model:                                 GEE   No. clusters:                      100
Method:                        Generalized   Min. cluster size:                   5
                      Estimating Equations   Max. cluster size:                   5
Family:                           Gaussian   Mean cluster size:                 5.0
Dependence structure:         Exchangeable   Num. iterations:                     5
Date:                     Sat, 18 Oct 2014   Scale:                           5.193
Covariance type:                    robust   Time:                         23:45:40
==============================================================================
                 coef    std err          z      P>|z|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
x1             0.1763      0.095      1.863      0.062        -0.009     0.362
x2             0.9634      0.062     15.520      0.000         0.842     1.085
x3             0.9634      0.062     15.520      0.000         0.842     1.085
x4            -2.1134      0.063    -33.614      0.000        -2.237    -1.990
x5            -2.1134      0.063    -33.614      0.000        -2.237    -1.990
==============================================================================
Skew:                          0.0282   Kurtosis:                      -0.0719
Centered skew:                 0.0668   Centered kurtosis:             -0.2998
==============================================================================
{'df': 2, 'p-value': 0.18965154372410886, 'statistic': 3.3251337420443026}

Our assessment will be based on the order statistics of the p-values, so we sort the p-values here.

In [8]:
pvalues = [np.sort(np.asarray(x)) for x in pvalues]

Next we plot the order statistics. The p-values for the null model (orange) should follow the grey line (the p-values should be uniform, so the order statistics should follow a linear pattern), and the p-values for the alternative model (purple) should fall well below the grey line (the smaller the p-values, the greater the power).

In [9]:
import matplotlib.pyplot as plt
plt.clf()
a1, = plt.plot(pvalues[0], lw=5, color='orange')
a2, = plt.plot(pvalues[1], lw=5, color='purple')
plt.plot([0, len(pvalues[0])], [0, 1], '-', color='grey')
leg = plt.figlegend((a1, a2), ("Null", "Alternative"), "upper center",
                    numpoints=1, ncol=2)
leg.draw_frame(False)
plt.xlabel("Position", size=17)
plt.ylabel("Order statistic", size=17)
Out[9]:
<matplotlib.text.Text at 0x7fdd911558d0>

Comparing score tests for one constraint with Wald tests

When the score test is testing only one linear constraint, it is easy to compare the score test to a Wald test.

In [10]:
# Constraint matrices for score test, testing that params[1] == params[2]
L = np.array([[0., 1, -1, 0, 0]])
R = np.r_[0.]

wald_tests = []
score_tests = []

nrep = 200
for k in range(nrep):
    
    params = np.r_[0, 1, 1 + np.random.normal(), 0, 2]
    endog, exog, groups = gendat(params)
    
    # Unconstrained model
    ex = sm.cov_struct.Exchangeable()
    ga = sm.families.Gaussian()
    model0 = sm.GEE(endog, exog, groups, family=ga, cov_struct=ex)
    result0 = model0.fit()

    # Wald test for the null hypothesis params[1] == params[2]
    se = np.sqrt(np.dot(L[0,:], np.dot(result0.cov_params(), L[0,:])))
    wald_test = np.dot(L[0, :], result0.params) / se
    pv = 2*norm.cdf(-np.abs(wald_test))
    wald_tests.append(pv)
    
    # Score test for the hypothesis params[1] == params[2]
    ex = sm.cov_struct.Exchangeable()
    ga = sm.families.Gaussian()
    model = sm.GEE(endog, exog, groups, family=ga, cov_struct=ex,
              constraint=(L,R))
    result = model.fit()
    score_tests.append(model.score_test_results["p-value"])

We plot the p-values from the two tests against each other to confirm that they are operating similarly.

In [11]:
plt.clf()
plt.plot(wald_tests, score_tests, 'o', color='orange', alpha=0.6)
plt.grid(True)
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xlabel("Wald test p-value", size=14)
plt.ylabel("Score test p-value", size=14)
Out[11]:
<matplotlib.text.Text at 0x7fdd910dca90>