10.2.3 - Joint Independence

With three variables, there are three ways to satisfy joint independence. The assumption is that one variable is independent of the other two, but the latter two can have an arbitrary association. For example, if \(A\) is jointly independent of \(B\) and \(C\), which we denote by \((A, BC)\), then \(A\) is independent of each of the others, and we can factor the three-way joint distribution into the product of the marginal distribution for \(A\) and the two-way joint distribution of \((B,C)\).

Main assumptions

(these are stated for the case \((A,BC)\) but hold in a similar way for \((B, AC)\) and \((C, AB)\) as well)

  • The \(N = IJK\) counts in the cells are assumed to be independent observations of a Poisson random variable, and
  • there are no partial interactions involving \(A\): \(\lambda_{ij}^{AB} =\lambda_{ik}^{AC} =0\), for all \(i, j, k\), and
  • there is no three-way interaction: \(\lambda_{ijk}^{ABC}=0\) for all \(i, j, k\).

Note the constraints above are in addition to the usual set-to-zero or sum-to-zero constraints (present even in the saturated model) imposed to avoid overparameterization.

Model Structure

\(\log(\mu_{ijk})=\lambda+\lambda_i^A+\lambda_j^B+\lambda_k^C +\lambda_{jk}^{BC}\)

In the example below, we consider the model where admission status is jointly independent of department and sex, which we denote by \((A, DS)\).

In SAS, the model can be fitted like this:

proc genmod data=berkeley order=data;
class D S A;
model count = D S A D*S / dist=poisson link=log lrci type3 obstats;
run;

This model implies that the association between D and S does NOT depend on the level of the variable A. That is, the association between department and sex is independent of the rejection/acceptance decision.

\(\theta_{DS(A=reject)}=\theta_{DS(A=accept)}=\dfrac{\text{exp}(\lambda_{ij}^{DS})\text{exp}(\lambda_{i'j'}^{DS})}{\text{exp}(\lambda_{i'j}^{DS})\text{exp}(\lambda_{ij'}^{DS})}\)

Since we are assuming that the (DS) distribution is independent of A, then we are assuming that the conditional distribution of DS, given A, is the same as the unconditional distribution of DS. And equivalently, the conditional distribution of A, given DS, is the same as the unconditional distribution of A. In other words, if this model fits well, neither department nor sex tells us anything significant about admission status.

The first estimated coefficient 1.9436 for the DS associations...

D*S        DeptA   Male     1    1.9436    0.1268    1.6950    2.1921   234.84

implies that the estimated odds ratio between sex and department (specifically, A versus F) is \(\exp(1.9436) = 6.98\) with 95% CI

\((\exp(1.695), \exp(2.192))= (5.45, 8.96)\)

But, we should really check the overall fit of the model first, to determine if these estimates are meaningful.

Model Fit

The goodness-of-fit statistics indicate that the model does not fit since the "Value/df" is much larger than 1.

 
Criteria For Assessing Goodness Of Fit
Criterion DF Value Value/DF
Deviance 11 877.0564 79.7324
Scaled Deviance 11 877.0564 79.7324
Pearson Chi-Square 11 797.7041 72.5186
Scaled Pearson X2 11 797.7041 72.5186
Log Likelihood   20074.6774  

How did we get these degrees of freedom? As usual, we're comparing this model to the saturated model, and the df is the difference in the numbers of parameters involved:

\(DF = (IJK-1) - [(I-1)+(J-1)+(K-1)+(J-1)(K-1)] = (I-1)(JK-1) \)

With \(I=2\), \(J=2\), and \(K=6\) corresponding to the levels of A, S, and D, respectively, this works out to be 11.

As before, we can look at the residuals for more information about why this model fits poorly. Recall that adjusted residuals have approximately a N(0, 1) distribution (i.e., "Std Pearson Residual"). In general, we have a lack of fit if (1) we have a large number of cells and adjusted residuals are greater than 3, or (2) we have a small number of cells and adjusted residuals are greater than 2. Here is only part of the output. Notice that the absolute value of the standardized residual for the first six cells are all large (e.g., in cell 1, the value is \(-15.1793\)). Many other such residuals are rather large as well, indicating that this model fits poorly.

Evaluate the residuals

 
Observation Statistics
Observation count D S A Predicted Value Linear Predictor Standard Error of the Linear Predictor HessWgt Lower Upper Raw Residual Pearson Residual Deviance Residual Std Deviance Residual Std Pearson Residual Likelihood Residual Leverage CookD DFBETA_Intercept DFBETA_DDeptA DFBETA_DDeptB DFBETA_DDeptC DFBETA_DDeptD DFBETA_DDeptE DFBETA_SMale DFBETA_AReject DFBETA_DDeptASMale DFBETA_DDeptBSMale DFBETA_DDeptCSMale DFBETA_DDeptDSMale DFBETA_DDeptESMale DFBETAS_Intercept DFBETAS_DDeptA DFBETAS_DDeptB DFBETAS_DDeptC DFBETAS_DDeptD DFBETAS_DDeptE DFBETAS_SMale DFBETAS_AReject DFBETAS_DDeptASMale DFBETAS_DDeptBSMale DFBETAS_DDeptCSMale DFBETAS_DDeptDSMale DFBETAS_DDeptESMale
1 313 DeptA Male Reject 505.09831 6.2247531 0.0367703 505.09831 469.97739 542.84378 -192.0983 -8.547431 -9.199152 -16.3367 -15.17931 -15.55562 0.6829213 38.173694 0.1338572 1.051E-15 0 0 -2.63E-16 0 0 -0.218635 -0.734349 5.255E-16 0 0 0 2.3367482 9.518E-15 0 0 -3.51E-15 0 0 -7.166703 -5.790199 2.414E-15 0 0 0
2 512 DeptA Male Accept 319.90169 5.7680137 0.0395092 319.90169 296.06444 345.65816 192.09831 10.740272 9.8692317 13.948262 15.179309 14.575998 0.4993589 17.678562 0.1338572 -1.33E-15 -4.99E-16 -4.99E-16 -4.99E-16 -4.99E-16 -6.66E-16 -0.218635 0.4650965 3.328E-16 6.656E-16 6.656E-16 6.656E-16 2.3367481 -1.21E-14 -2.41E-15 -7.35E-15 -6.67E-15 -6.75E-15 -8.88E-15 -7.166703 3.6671961 1.529E-15 6.534E-15 6.441E-15 5.751E-15
3 19 DeptA Female Reject 66.122038 4.1915021 0.0969494 66.122038 54.679277 79.959431 -47.12204 -5.794967 -6.845121 -11.12613 -9.419201 -10.09928 0.6214932 11.205917 0.0275065 -1.152726 -7.77E-16 -7.77E-16 -7.86E-16 -7.77E-16 -1.13E-15 -0.044928 1.1527259 1.062E-15 1.074E-15 1.075E-15 1.128E-15 0.4801819 -10.4398 -3.75E-15 -1.14E-14 -1.05E-14 -1.05E-14 -1.51E-14 -1.472697 9.0890214 4.878E-15 1.054E-14 1.04E-14 9.745E-15
4 89 DeptA Female Accept 41.878088 3.7347627 0.0980209 41.878088 34.558213 50.748407 47.121912 7.2816446 6.3202597 8.1755761 9.4191761 8.6973682 0.402369 4.5948769 0.0275065 0.7300717 3.761E-16 3.761E-16 3.761E-16 3.761E-16 5.813E-16 -0.044928 -0.730072 -5.47E-16 -5.47E-16 -5.47E-16 -5.81E-16 0.4801807 6.6119817 1.815E-15 5.535E-15 5.027E-15 5.083E-15 7.759E-15 -1.472693 -5.756474 -2.51E-15 -5.37E-15 -5.29E-15 -5.02E-15
5 207 DeptB Male Reject 342.85461 5.8373065 0.0438822 342.85461 314.59903 373.64795 -135.8546 -7.337015 -7.925271 -13.59608 -12.58691 -12.93864 0.6602177 23.679966 0.0883403 -6.94E-16 0 -6.94E-16 -6.94E-16 -3.47E-16 -6.94E-16 -0.14429 1.04E-15 -0.713979 6.936E-16 1.04E-15 6.936E-16 1.5421589 -6.28E-15 0 -1.02E-14 -9.27E-15 -4.69E-15 -9.26E-15 -4.729733 8.203E-15 -3.279442 6.809E-15 1.007E-14 5.993E-15
6 353 DeptB Male Accept 217.14539 5.3805671 0.0462014 217.14539 198.34621 237.72635 135.85461 9.2193238 8.4461135 11.531262 12.586906 12.032087 0.4635118 10.529199 0.0883403 -1.1E-16 0 1.098E-16 0 -1.1E-16 0 -0.14429 -2.2E-16 0.4521955 0 -2.2E-16 0 1.5421589 -9.95E-16 0 1.616E-15 0 -1.48E-15 0 -4.729733 -1.73E-15 2.0770196 0 -2.13E-15 0
7 8 DeptB Female Reject 15.30601 2.7282455 0.2003495 15.30601 10.335325 22.667301 -7.30601 -1.867451 -2.056977 -3.312462 -3.007258 -3.128479 0.6143822 1.108357 0.0041861 1.682E-16 -0.75785 1.939E-16 1.926E-16 1.939E-16 1.981E-16 -0.006837 -1.95E-16 0.7578499 -2.47E-16 -2.31E-16 -2.39E-16 0.0730767 1.523E-15 -3.657546 2.853E-15 2.574E-15 2.62E-15 2.644E-15 -0.224123 -1.54E-15 3.4809481 -2.43E-15 -2.23E-15 -2.07E-15
8 17 DeptB Female Accept 9.6939907 2.2715062 0.2008702 9.6939907 6.5391536 14.37089 7.3060093 2.3465452 2.1180239 2.7143924 3.0072581 2.8325522 0.3911414 0.4469052 0.0041861 -1.3E-16 0.4799807 -1.41E-16 -1.41E-16 -1.41E-16 -1.46E-16 -0.006837 1.457E-16 -0.479981 1.769E-16 1.665E-16 1.717E-16 0.0730767 -1.18E-15 2.3164901 -2.07E-15 -1.88E-15 -1.9E-15 -1.94E-15 -0.224123 1.149E-15 -2.204642 1.737E-15 1.612E-15 1.484E-15
9 205 DeptC Male Reject 198.97812 5.2931949 0.0567174 198.97812 178.04404 222.3736 6.0218769 0.426903 0.4247764 0.7080436 0.7115884 0.7103146 0.6400844 0.0692709 -0.003697 2.177E-17 2.902E-17 0 2.177E-17 1.451E-17 2.902E-17 0.006038 -2.9E-17 -4.35E-17 0.0514811 -2.9E-17 -2.9E-17 -0.064534 1.971E-16 1.401E-16 0 2.909E-16 1.961E-16 3.874E-16 0.197922 -2.29E-16 -2E-16 0.5053783 -2.81E-16 -2.51E-16
10 120 DeptC Male Accept 126.02188 4.8364555 0.0585301 126.02188 112.36343 141.34059 -6.021879 -0.536425 -0.540785 -0.717372 -0.711589 -0.714881 0.431723 0.029591 -0.003697 9.191E-18 -4.6E-18 1.838E-17 4.596E-18 4.596E-18 0 0.006038 0 9.191E-18 -0.032605 0 0 -0.064534 8.324E-17 -2.22E-17 2.705E-16 6.142E-17 6.21E-17 0 0.1979221 0 4.222E-17 -0.320079 0 0
11 391 DeptC Female Reject 363.05854 5.8945641 0.0427349 363.05854 333.88785 394.77779 27.941455 1.4664278 1.4481968 2.4948336 2.5262405 2.5157016 0.6630449 0.9659997 -0.018322 1.899E-17 -9.36E-17 0.1398371 -1.58E-17 1.433E-17 -3.99E-18 0.0299254 -1E-17 1.916E-16 -0.139837 3.936E-17 3.992E-18 -0.31984 1.719E-16 -4.52E-16 2.0575648 -2.11E-16 1.936E-16 -5.33E-17 0.9809346 -7.92E-17 8.802E-16 -1.372749 3.809E-16 3.449E-17
12 202 DeptC Female Accept 229.94146 5.4378247 0.0451131 229.94146 210.48294 251.19886 -27.94146 -1.84264 -1.881985 -2.580183 -2.526241 -2.555081 0.4679758 0.4318154 -0.018322 9.111E-17 1.367E-16 -0.088565 9.111E-17 6.833E-17 9.111E-17 0.0299254 -9.11E-17 -2.05E-16 0.0885652 -1.14E-16 -9.11E-17 -0.31984 8.251E-16 6.595E-16 -1.303149 1.218E-15 9.233E-16 1.216E-15 0.9809347 -7.18E-16 -9.42E-16 0.8694243 -1.1E-15 -7.87E-16
13 279 DeptD Male Reject 255.30424 5.5424559 0.0503787 255.30424 231.29997 281.79967 23.695762 1.4830018 1.4609054 2.4622378 2.4994795 2.4864328 0.6479663 0.8845534 -0.014872 2.919E-17 0 -2.92E-17 0 -2.92E-17 0 0.0242913 0 -5.84E-17 0 0.1614174 0 -0.259622 2.644E-16 0 -4.3E-16 0 -3.94E-16 0 0.7962501 0 -2.68E-16 0 1.5620689 0
14 138 DeptD Male Accept 161.69576 5.0857166 0.0524112 161.69576 145.91036 179.18892 -23.69576 -1.863466 -1.912007 -2.564589 -2.49948 -2.535876 0.444168 0.3840251 -0.014872 7.395E-17 5.546E-17 7.395E-17 7.395E-17 7.395E-17 7.395E-17 0.0242913 -7.4E-17 -3.7E-17 -7.4E-17 -0.102233 -7.4E-17 -0.259622 6.698E-16 2.677E-16 1.088E-15 9.883E-16 9.993E-16 9.871E-16 0.7962502 -5.83E-16 -1.7E-16 -7.26E-16 -0.989329 -6.39E-16
15 244 DeptD Female Reject 229.59014 5.4362957 0.0529774 229.59014 206.94685 254.71097 14.409858 0.9510056 0.941309 1.5784537 1.5947136 1.5889501 0.644368 0.3544502 -0.008952 9.714E-17 4.215E-17 7.729E-17 0.1080507 9.486E-17 1.562E-16 0.0146225 -1.63E-16 -8.21E-17 -1.21E-16 -0.108051 -1.39E-16 -0.156284 8.797E-16 2.034E-16 1.137E-15 1.4439895 1.282E-15 2.085E-15 0.479316 -1.29E-15 -3.77E-16 -1.19E-15 -1.045628 -1.2E-15
16 131 DeptD Female Accept 145.40986 4.9795564 0.0549138 145.40986 130.57234 161.93344 -14.40986 -1.194986 -1.215586 -1.622204 -1.594714 -1.610208 0.4384866 0.152763 -0.008953 -1.11E-17 1.113E-17 -1.11E-17 -0.068433 -2.23E-17 -5.56E-17 0.0146225 5.565E-17 1.113E-17 3.339E-17 0.0684334 4.452E-17 -0.156284 -1.01E-16 5.371E-17 -1.64E-16 -0.914544 -3.01E-16 -7.43E-16 0.4793161 4.388E-16 5.112E-17 3.278E-16 0.6622441 3.847E-16
17 138 DeptE Male Reject 116.93791 4.7616431 0.0733181 116.93791 101.28541 135.00933 21.062088 1.9477076 1.8932348 3.106604 3.1959883 3.1630862 0.6286041 1.3298634 -0.01253 -2.46E-17 0 -4.92E-17 0 -9.84E-17 -4.92E-17 0.0204658 4.919E-17 0 9.838E-17 4.919E-17 0.2969142 -0.218736 -2.23E-16 0 -7.24E-16 0 -1.33E-15 -6.57E-16 0.6708529 3.878E-16 0 9.657E-16 4.76E-16 2.5655562
18 53 DeptE Male Accept 74.062089 4.3049038 0.0747292 74.062089 63.971469 85.744365 -21.06209 -2.447392 -2.579791 -3.368885 -3.195988 -3.298475 0.4135965 0.5541756 -0.01253 9.346E-17 4.673E-17 7.788E-17 6.231E-17 9.346E-17 9.346E-17 0.0204658 -9.35E-17 -6.23E-17 -1.25E-16 -9.35E-17 -0.188049 -0.218736 8.464E-16 2.255E-16 1.146E-15 8.327E-16 1.263E-15 1.247E-15 0.6708529 -7.37E-16 -2.86E-16 -1.22E-15 -9.04E-16 -1.624883
19 299 DeptE Female Reject 240.61047 5.4831793 0.0518118 240.61047 217.37632 266.32799 58.389531 3.7642437 3.6255985 6.0928851 6.3258808 6.2443737 0.6459102 5.6150979 -0.036434 3.776E-17 -1.86E-16 2.849E-17 4.013E-17 0.4195937 1.351E-16 0.0595093 -1.63E-16 9.503E-17 -6.36E-17 -1.36E-16 -0.419594 -0.636029 3.419E-16 -8.98E-16 4.193E-16 5.362E-16 5.669627 1.803E-15 1.9506733 -1.29E-15 4.365E-16 -6.24E-16 -1.32E-15 -3.625598
20 94 DeptE Female Accept 152.38953 5.02644 0.0537902 152.38953 137.14149 169.33293 -58.38953 -4.72996 -5.093896 -6.812612 -6.325881 -6.602426 0.4409215 2.4276557 -0.036434 1.812E-16 2.718E-16 1.359E-16 1.359E-16 -0.265748 9.059E-17 0.0595093 -9.06E-17 -2.26E-16 -1.36E-16 -9.06E-17 0.2657478 -0.636029 1.641E-15 1.312E-15 1.999E-15 1.816E-15 -3.590832 1.209E-15 1.9506734 -7.14E-16 -1.04E-15 -1.33E-15 -8.77E-16 2.2962557
21 351 DeptF Male Reject 228.36589 5.4309491 0.0531121 228.36589 205.78898 253.41969 122.63411 8.1151336 7.5151464 12.598896 13.604754 13.255617 0.6441967 25.777847 -0.076153 1.196E-15 -4.48E-16 0 7.474E-16 2.99E-16 0.9240428 0.1243841 -0.924043 -0.924043 -0.924043 -0.924043 -0.924043 -1.329404 1.083E-14 -2.16E-15 0 9.988E-15 4.039E-15 12.333171 4.077222 -7.285899 -4.244304 -9.071119 -8.942148 -7.984408
22 22 DeptF Male Accept 144.63448 4.9742098 0.0550438 144.63448 129.84299 161.111 -122.6345 -10.1971 -12.744 -17.00283 -13.6048 -15.6051 0.4382161 11.106051 -0.076153 -2.84E-16 5.68E-16 2.84E-16 -9.47E-17 9.467E-17 -0.58524 0.1243845 0.58524 0.58524 0.58524 0.58524 0.58524 -1.329408 -2.57E-15 2.741E-15 4.179E-15 -1.27E-15 1.279E-15 -7.811181 4.0772344 4.6145047 2.6881185 5.745169 5.6634853 5.0569034
23 317 DeptF Female Reject 208.77408 5.3412527 0.05543 208.77408 187.28133 232.73339 108.22592 7.4901924 6.9525291 11.611039 12.508961 12.194621 0.6414552 21.533862 0.8184913 -0.885183 -0.885183 -0.885183 -0.885183 -0.885183 -0.885183 0.1089309 0.8851832 0.8851832 0.8851832 0.8851832 0.8851832 14.288419 -8.016767 -4.272084 -13.0246 -11.82959 -11.96076 -11.81451 3.5706788 6.9794986 4.0658143 8.6896432 8.5660954 7.6486325
24 24 DeptF Female Accept 132.22611 4.8845134 0.0572835 132.22611 118.18364 147.93708 -108.2261 -9.411816 -11.59923 -15.41621 -12.50898 -14.22795 0.4338873 9.225178 -0.62732 0.5606277 0.5606277 0.5606277 0.5606277 0.5606277 0.5606277 0.1089311 -0.560628 -0.560628 -0.560628 -0.560628 -0.560628 -10.95113 5.0773916 2.7057097 8.2490839 7.4922269 7.5753036 7.4826807 3.5706851 -4.420441 -2.575069 -5.503555 -5.425307 -4.844235

In R, one way to fit the (A, DS) model is with the following command:

berk.jnt = glm(Freq~Admit+Gender+Dept+Gender*Dept, family=poisson(), data=berk.data)

This model implies that association between D and S does NOT depend on level of the variable A. That is, the association between department and sex is independent of the rejection/acceptance decision.

\(\theta_{DS(A=reject)}=\theta_{DS(A=accept)}=\dfrac{\text{exp}(\lambda_{ij}^{DS})\text{exp}(\lambda_{i'j'}^{DS})}{\text{exp}(\lambda_{i'j}^{DS})\text{exp}(\lambda_{ij'}^{DS})}\)

Since we are assuming that the (DS) distribution is independent of A, then we are assuming that the conditional distribution of DS, given A, is the same as the unconditional distribution of DS. And equivalently, the conditional distribution of A, given DS, is the same as the unconditional distribution of A. In other words, if this model fits well, neither department nor sex tells us anything significant about admission status.

The first estimated coefficient 1.9436 for the DS associations implies that the estimated odds ratio between sex and department (specifically, A versus F) is \(\exp(1.9436) = 6.98\) with 95% CI

\((\exp(1.695), \exp(2.192))= (5.45, 8.96)\)

> summary(berk.jnt)

Call:
glm(formula = Freq ~ Admit + Gender + Dept + Gender * Dept, family = poisson(), 
    data = berk.data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-12.744   -3.208   -0.058    2.495    9.869  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)       4.88451    0.05728  85.269  < 2e-16 ***
AdmitRejected     0.45674    0.03051  14.972  < 2e-16 ***
GenderMale        0.08970    0.07492   1.197   0.2312    
DeptA            -1.14975    0.11042 -10.413  < 2e-16 ***
DeptB            -2.61301    0.20720 -12.611  < 2e-16 ***
DeptC             0.55331    0.06796   8.141 3.91e-16 ***
DeptD             0.09504    0.07483   1.270   0.2040    
DeptE             0.14193    0.07401   1.918   0.0551 .  
GenderMale:DeptA  1.94356    0.12683  15.325  < 2e-16 ***
GenderMale:DeptB  3.01937    0.21771  13.869  < 2e-16 ***
GenderMale:DeptC -0.69107    0.10187  -6.784 1.17e-11 ***
GenderMale:DeptD  0.01646    0.10334   0.159   0.8734    
GenderMale:DeptE -0.81123    0.11573  -7.010 2.39e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2650.10  on 23  degrees of freedom
Residual deviance:  877.06  on 11  degrees of freedom
AIC: 1062.1

But, we should really check the overall fit of the model first, to determine if these estimates are meaningful.

Model Fit

The goodness-of-fit statistic (in the output as "Residual deviance") 877.06 indicates that the model does not fit since the "Value/df" is much larger than 1.

How did we get these degrees of freedom? As usual, we're comparing this model to the saturated model, and the df is the difference in the numbers of parameters involved:

\(DF = (IJK-1) - [(I-1)+(J-1)+(K-1)+(J-1)(K-1)] = (I-1)(JK-1) \)

With \(I=2\), \(J=2\), and \(K=6\) corresponding to the levels of A, S, and D, respectively, this works out to be 11.

As before, we can look at the residuals for more information about why this model fits poorly. Recall that adjusted residuals have approximately a N(0, 1) distribution. In general, we have a lack of fit if (1) we have a large number of cells and adjusted residuals are greater than 3, or (2) we have a small number of cells and adjusted residuals are greater than 2. Here is only part of the output. Notice that the absolute value of the standardized residual for the first six cells are all large (e.g., in cell 1, the value is \(-15.18\)). Many other such residuals are rather large as well, indicating that this model fits poorly.

Evaluate the residuals

> fits = fitted(berk.jnt)
> resids = residuals(berk.jnt,type="pearson")
> h = lm.influence(berk.jnt)$hat
> adjresids = resids/sqrt(1-h)
> round(cbind(berk.data$Freq,fits,adjresids),2)
         fits adjresids
1  512 319.90     15.18
2  313 505.10    -15.18
3   89  41.88      9.42
4   19  66.12     -9.42
5  353 217.15     12.59
6  207 342.85    -12.59
...

Next, let us look at what is often the most interesting model of conditional independence.