BUG: GEE score #1839

Merged
merged 4 commits into from Jul 20, 2014

Projects

None yet

3 participants

@kshedden
Contributor

Fixes a bug in GEE score test that crept in somewhere during refactoring. There was previously no test coverage for GEE score testing. I added a simple regression test here.

kshedden added some commits Jul 19, 2014
@kshedden kshedden Fixed bug in score test 99b030b
@kshedden kshedden Simplified test code
c66deba
@kshedden kshedden Remove unwanted print statement
a4a7e36
@coveralls

Coverage Status

Coverage increased (+0.01%) when pulling a4a7e36 on kshedden:gee_score into b2c1259 on statsmodels:master.

@coveralls

Coverage Status

Coverage increased (+0.01%) when pulling a4a7e36 on kshedden:gee_score into b2c1259 on statsmodels:master.

@josef-pkt josef-pkt commented on the diff Jul 19, 2014
statsmodels/genmod/tests/test_gee.py
+ n = 200 # Must be divisible by 4
+ exog = np.random.normal(size=(n, 4))
+ endog = exog[:, 0] + exog[:, 1] + exog[:, 2]
+ endog += 3*np.random.normal(size=n)
+ group = np.kron(np.arange(n/4), np.ones(4))
+
+ # Test under the null.
+ L = np.array([[1., -1, 0, 0]])
+ R = np.array([0.,])
+ family = Gaussian()
+ va = Independence()
+ mod1 = GEE(endog, exog, group, family=family,
+ cov_struct=va, constraint=(L, R))
+ rslt1 = mod1.fit()
+ assert_almost_equal(mod1.score_test_results["statistic"],
+ 1.08126334)
@josef-pkt
josef-pkt Jul 19, 2014 Member

add assert also for p-values?

Are these close to a Wald test, so we can have a coarse cross-check?

@kshedden kshedden Added test comparing score test to Wald test
cfef60e
@coveralls

Coverage Status

Coverage increased (+0.02%) when pulling cfef60e on kshedden:gee_score into b2c1259 on statsmodels:master.

@josef-pkt josef-pkt commented on the diff Jul 20, 2014
statsmodels/genmod/tests/test_gee.py
+ np.random.normal(size=n)
+ family = Gaussian()
+ va = Independence()
+ mod0 = GEE(endog, exog, group, family=family,
+ cov_struct=va)
+ rslt0 = mod0.fit()
+ family = Gaussian()
+ va = Independence()
+ mod1 = GEE(endog, exog, group, family=family,
+ cov_struct=va, constraint=(L, R))
+ rslt1 = mod1.fit()
+ se = np.sqrt(np.dot(f, np.dot(rslt0.cov_params(), f)))
+ wald_z = np.dot(f, rslt0.params) / se
+ wald_p = 2*norm.cdf(-np.abs(wald_z))
+ score_p = mod1.score_test_results["p-value"]
+ assert(np.abs(wald_p - score_p) < 0.02)
@josef-pkt
josef-pkt Jul 20, 2014 Member

same but better is assert_allclose(wald_p, score_p, atol=0.02)
(more informative failure message, and assert get optimized away in byte compiling python code with option -o)

@josef-pkt
Member

notebook looks very good.

merging in spite of assert

@josef-pkt josef-pkt merged commit 3f49c10 into statsmodels:master Jul 20, 2014

2 checks passed

continuous-integration/appveyor AppVeyor build succeeded
Details
continuous-integration/travis-ci The Travis CI build passed
Details
@josef-pkt josef-pkt added the PR label Aug 11, 2014
@josef-pkt josef-pkt added this to the 0.6 milestone Aug 24, 2014
@kshedden kshedden deleted the kshedden:gee_score branch Sep 27, 2014
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment