10.1.3 - Saturated Model
10.1.3 - Saturated ModelWith the saturated model, the \(N = IJ\) counts in the cells are still assumed to be independent observations of a Poisson random variable, but no independence is assumed between the variables \(A\) and \(B\). The model is expressed as
\(\log(\mu_{ij})=\lambda+\lambda_i^A+\lambda_j^B+\lambda_{ij}^{AB}\)
Note the additional \(\lambda_{ij}^{AB}\) interaction term, compared with the independence model. We still have the same ANOVA-like constraints to avoid over-parameterization, and again, this changes with the software used. Ultimately, the inference, in the end, will not depend on the choice of constraints.
Parameter estimates and interpretation
The odds ratio is directly related to the interaction terms. For example, for a \(2 \times 2\) table,
\begin{align}
\log(\theta) &= \log\left(\dfrac{\mu_{11}\mu_{22}}{\mu_{12}\mu_{21}}\right)\\
&= \\
&= \lambda_{11}^{AB}+\lambda_{22}^{AB}-\lambda_{12}^{AB}-\lambda_{21}^{AB}
\end{align}
\(\begin{align*}\log(\theta)& = \log\left(\dfrac{\mu_{11}\mu_{22}}{\mu_{12}\mu_{21}}\right)\\ &= \lambda+\lambda_1^A+\lambda_1^B+\lambda_{11}^{AB}+\lambda+\lambda_2^A+\lambda_2^B+\lambda_{22}^{AB}-\lambda-\lambda_1^A-\lambda_2^B-\lambda_{12}^{AB}-\lambda-\lambda_2^A-\lambda_1^B-\lambda_{21}^{AB}\\ &= \lambda_{11}^{AB}+\lambda_{22}^{AB}-\lambda_{12}^{AB}-\lambda_{21}^{AB}\\ \end{align*}\)
How many odds ratios are there? There should be \((I-1) (J-1)\) which should be equal to the unique number of \(\lambda_{ij}\)s in the model.
How many unique parameters are there in this model?
| Term | # of terms | #of constraints | # of unique parameters |
|---|---|---|---|
| \(\lambda\) | \(1\) | \(0\) | \(1\) |
| \(\lambda_i^A\) | \(I\) | \(1\) | \(I - 1\) |
| \(\lambda_i^B\) | \(J\) | \(1\) | \(J - 1\) |
| \(\lambda_{ij}^{AB}\) | \(I J\) | \(I + J - 1\) | \((I - 1) (J - 1)\) |
With \(I J = N\) terms, the model is a perfect fit! Moreover, the saturated model includes the independence model as a special case, where the interaction terms are 0.
Example: Vitamin C
Let's evaluate the CATMOD and GENMOD output for VitaminCLoglin.sas and R code for VitaminCLoglin.R.
To fit the saturated models, we need to specify or add an interaction term in our code.
In CATMOD:
proc catmod data=ski order=data;
weight count;
model treatment*response=_response_;
loglin treatment|response;
run;
The LOGLIN statement (line) now specifies the saturated model. treatment|response specification adds an interaction term to the main effects model.
In GENMOD:
proc genmod data=ski order=data;
class treatment response;
model count = response*treatment /link=log dist=poisson lrci type3 obstats;
run;
The "*" sign in the response*treatment, indicates that we want the interaction term and ALL lower order terms; that is, the saturated model.
Notice, that we can also specify the saturated model by explicitly entering the main effects in the model statement:
model count = response treatment response*treatment/link=log dist=poisson lrci type3 obstats;
In these two model specifications under GENMOD, we are fitting different numbers of parameters, and the output sections on LR Statistics For Type 3 Analysis will be different, but the estimates of odds and odds-ratios will be the same. For the output, we have the same parts as we did for the independence model, except now we have more parameter estimates. We will consider the latter GENMOD parametrization and compare it to that for CATMOD.
So, what are the odds of getting cold given that a person took a placebo pill? Let's first start with log-odds:
\(\log(\mu_{11}/\mu_{12})=\log(\mu_{11})-\log(\mu_{12})=[\lambda+\lambda_1^A+\lambda_1^B+\lambda_{11}^{AB}]-[\lambda+\lambda_1^A+\lambda_2^B+\lambda_{12}^{AB}]\)
Based on GENMOD parameterization: \(\lambda_1^B + \lambda_{11}^{AB} - \lambda_2^B -\lambda_{12}^{AB} = -1.9708 + 0.7134 - 0 - 0\), and the odds are \(\exp(-1.2574) = 0.2844\).
Based on CATMOD parameterization: \([\lambda+\lambda_1^A + \lambda_1^B+\lambda_{11}^{AB} ]- [\lambda+\lambda_1^A- \lambda_1^B - \lambda_{11}^{AB}] = 2\lambda_1^B+2\lambda_{11}^{AB} = 2(-0.807 + 0.1784)\), and the odds are \(\exp(\log \mbox{odds}) = 0.2844\).
In R, the syntax for loglin() would be
loglin(ski, list(c(1, 2)), fit=TRUE, param=TRUE)
And the syntax for glm() would be
glm(Freq~Treatment*Cold, family=poisson(), data=ski.data)
Q1: What are the odds ratio based on CATMOD or loglin()? What about GENMOD or glm()?
For example, from the second way of fitting the model with GENMOD, the odds ratio is \(\exp(0.7134)=2.041\), with 95% CI \((\exp(0.079), \exp(1.3770)\). This corresponds to \(1/2.041=0.49\), which we saw before.
treatment*response placebo cold 1 0.7134 0.3293 0.0790 1.3770
From R, we get the same:
> summary(ski.sat)
...
TreatmentVitaminC:ColdNoCold 0.7134 0.3293 2.166 0.0303 *
Q2: Is the interaction of treatment and response significant (hint: look at LR Statistics for Type 3 Analysis). Where does df = 3 come from?
Q3: Model fits perfectly. Let's do some model selection. Compare the independence and saturated model, by calculating the LR statistics (e.g. difference in deviances?
\(2\times (973.0657-970.629)\approx 4.8717\)
Q4: How many df does this statistic have?
Hierarchical Models
These models include all lower-order terms that comprise higher-order terms in the model.
(A, B) is a simpler model than (AB)
Interpretation does not depend on how the variables are coded.
\(\log(\mu_{ij})=\lambda+\lambda_i^A+ \lambda_{ij}^{AB}\)
No, it is not a hierarchical model. It is missing \(\lambda^{B}_{j}\).
We want hierarchical models so that interaction terms represent associations. If there is a significant interaction term, we do NOT look at the lower-order terms, but only interpret the higher-order terms because the value of lower-order terms is coding dependent and can be misleading.
Next, we'll see how these models extend to three-way tables.