Why homogeneous association?
These models can be described by an undirected graph.
Mutual independence of \(Y\) and \((X,Z)\)
Models for contingency tables
Estimation via pseudo-likelihood.
\(N_{ij} \sim \text{Poisson}(\mu_{ij})\), \(1 \leq i \leq I, 1 \leq j \leq J\)
Independence model \[ \log(\mu_{ij}) = \lambda + \lambda_i^X + \lambda_j^Y \]
Saturated model \[ \log(\mu_{ij}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_{ij}^{XY} \]
\[ P(X=i,Y=j) = \frac{\mu_{ij}}{\sum_{ab} \mu_{ab}}. \]
\[ \begin{aligned} \lambda^X_I &= 0 \\ \lambda^Y_J &= 0 \\ \lambda^{XY}_{iJ} &= \lambda^{XY}_{Ij} = 0 \end{aligned} \]
\[ \begin{aligned} \log(\mu_{ijk}) &= \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k \\ & \qquad + \lambda^{XY}_{ij} +\lambda^{YZ}_{jk} +\lambda^{XZ}_{ik} \\ & \qquad + \lambda^{XYZ}_{ijk} \end{aligned} \]
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k \]
\[ \mu_{ijk} \propto e^{\lambda^X_i} \cdot e^{\lambda^Y_j} \cdot e^{\lambda^Z_k}. \]
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} \]
\[ \mu_{ijk} \propto e^{\lambda^Y_j} \cdot e^{\lambda^X_i + \lambda^{Z}_k + \lambda^{XZ}_{ik}}. \]
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \]
\[ \mu_{ijk} \propto e^{\lambda^Y_j + \lambda^{YZ}_{jk}} \cdot e^{\lambda^X_i + \lambda^{Z}_k + \lambda^{XZ}_{ik}}. \]
\[ H_0: \lambda^{XY}_{ij} \equiv 0, \lambda^{XYZ}_{ijk} \equiv 0 \]
\[ H_a: \lambda^{XYZ}_{ijk} \equiv 0 \]
## Loading required package: vcd## Loading required package: grid## Loading required package: gnmdata(DaytonSurvey)
Dayton.aggregated = aggregate(Freq ~ cigarette + alcohol + marijuana,
data=DaytonSurvey, FUN=sum)
Dayton.aggregated## cigarette alcohol marijuana Freq
## 1 Yes Yes Yes 911
## 2 No Yes Yes 44
## 3 Yes No Yes 3
## 4 No No Yes 2
## 5 Yes Yes No 538
## 6 No Yes No 456
## 7 Yes No No 43
## 8 No No No 279full.model = glm(Freq ~ .^2, data=Dayton.aggregated, family=poisson) # homogeneous association: only up to 2way interactions
summary(full.model)##
## Call:
## glm(formula = Freq ~ .^2, family = poisson, data = Dayton.aggregated)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## 0.02044 -0.09256 -0.33428 0.49134 -0.02658 0.02890 0.09452 -0.03690
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.81387 0.03313 205.699 < 2e-16 ***
## cigaretteNo -3.01575 0.15162 -19.891 < 2e-16 ***
## alcoholNo -5.52827 0.45221 -12.225 < 2e-16 ***
## marijuanaNo -0.52486 0.05428 -9.669 < 2e-16 ***
## cigaretteNo:alcoholNo 2.05453 0.17406 11.803 < 2e-16 ***
## cigaretteNo:marijuanaNo 2.84789 0.16384 17.382 < 2e-16 ***
## alcoholNo:marijuanaNo 2.98601 0.46468 6.426 1.31e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2851.46098 on 7 degrees of freedom
## Residual deviance: 0.37399 on 1 degrees of freedom
## AIC: 63.417
##
## Number of Fisher Scoring iterations: 4reduced.model = glm(Freq ~ alcohol + marijuana + cigarette +
alcohol:marijuana +
cigarette:marijuana,
data=Dayton.aggregated, family=poisson)
summary(reduced.model)##
## Call:
## glm(formula = Freq ~ alcohol + marijuana + cigarette + alcohol:marijuana +
## cigarette:marijuana, family = poisson, data = Dayton.aggregated)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## 0.0584 -0.2619 -0.8663 2.2287 4.5702 -4.3441 -9.7716 6.8353
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.81261 0.03316 205.450 <2e-16 ***
## alcoholNo -5.25227 0.44837 -11.714 <2e-16 ***
## marijuanaNo -0.72847 0.05538 -13.154 <2e-16 ***
## cigaretteNo -2.98919 0.15111 -19.782 <2e-16 ***
## alcoholNo:marijuanaNo 4.12509 0.45294 9.107 <2e-16 ***
## marijuanaNo:cigaretteNo 3.22431 0.16098 20.029 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2851.46 on 7 degrees of freedom
## Residual deviance: 187.75 on 2 degrees of freedom
## AIC: 248.8
##
## Number of Fisher Scoring iterations: 5anova(reduced.model, full.model) # null: conditional independence of alcohol & cigarette | marijuana## Analysis of Deviance Table
##
## Model 1: Freq ~ alcohol + marijuana + cigarette + alcohol:marijuana +
## cigarette:marijuana
## Model 2: Freq ~ (cigarette + alcohol + marijuana)^2
## Resid. Df Resid. Dev Df Deviance
## 1 2 187.754
## 2 1 0.374 1 187.38\[ \theta_{ij(k)} = \frac{\mu_{ijk} \mu_{i+1,j+1,k}}{\mu_{i+1,j,k} \mu_{i,j+1,k}}. \]
\[ \log \theta_{ij(k)} = \lambda^{XY}_{ij} + \lambda^{XY}_{i+1,j+1} - \lambda^{XY}_{i+1,j} - \lambda^{XY}_{i,j+1} \]
These models can be described by an undirected graph.
\[ \begin{aligned} \log(\mu_{ijkl}) &= \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^W_l \\ & \qquad + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{XW}_{il} \\ & \qquad + \lambda^{YZ}_{jk} + \lambda^{YW}_{jl} + \lambda^{ZW}_{kl} \end{aligned} \]
Homogeneous model with this graph has same conditional independencies as \((XW, WYZ)\).
\(\implies\) conditional independencies are not going to tell us whether there are any 3-way interactions or not.
\(W\) is conditionally independent of \(Y\) given \((X,Z)\).
\(X\) is conditionally independent of \(Z\) given \((W,Y)\).
Model: \(\lambda^{XZ}_{ik} \equiv 0, \lambda^{YW}_{jl} \equiv 0\).
Generally, for vertex \(v\) set \(N(v)\) to be its neighbors in the graph \(G=(E,V)\).
The conditional independencies: \(v\) is independent of \((\{v\} \cup N(v))^c\) given \(N(v)\).
For simplicity let’s assume each of \((X,Y,Z,W)\) is binary.
In complete graph
\[ \begin{aligned} \log \left(\frac{P(Y=1|X=i,Z=k,W=l)}{P(Y=0|X=i,Z=k,W=l)}\right) &= \log\left(\frac{\mu_{i1kl}}{\mu_{i0kl}}\right) \\ &= \lambda^Y_1 - \lambda^Y_0 + \\ & \qquad \lambda^{XY}_{i1} - \lambda^{XY}_{i0} + \\ & \qquad \lambda^{YZ}_{1k} - \lambda^{YZ}_{0k} + \\ & \qquad \lambda^{YW}_{1l} - \lambda^{YW}_{0l} \end{aligned} \]
A logistic regression model with Binomial outcome \(Y_{ikl}\) and 3-column design matrix.
Number of rows is total count in the \(2x2x2x2\) table.
\[ \begin{aligned} \log \left(\frac{P(Y=1|X=i,Z=k,W=l)}{P(Y=0|X=i,Z=k,W=l)}\right) &= \log\left(\frac{\mu_{i1kl}}{\mu_{i0kl}}\right) \\ &= \lambda^Y_1 - \lambda^Y_0 + \\ & \qquad \lambda^{XY}_{i1} - \lambda^{XY}_{i0} + \\ & \qquad \lambda^{YZ}_{1k} - \lambda^{YZ}_{0k} \end{aligned} \]
Features for the logistic regression model for \(Y\) determined by \(N(Y)\).
Design matrix has 2 columns plus intercept.
\(\lambda^Y_1 - \lambda^Y_0\)
\(\lambda^{XY}_{11} - \lambda^{XY}_{10} - \lambda^{XY}_{01} + \lambda^{XY}_{00}\)
\(\lambda^{YW}_{11} - \lambda^{YW}_{01} - \lambda^{YW}_{10} + \lambda^{YW}_{00}\)
\(\lambda^X_1 - \lambda^X_0\)
\(\lambda^{XY}_{11} - \lambda^{XY}_{01} - \lambda^{XY}_{10} + \lambda^{XY}_{00}\)
\(\lambda^{XW}_{11} - \lambda^{XW}_{01} - \lambda^{XW}_{10} + \lambda^{XW}_{00}\)
Categorical models become baseline multinomial logits.
We estimate matrices of parameters – more complicated notation, conceptually the same.
Usual default constraints would set anything with a 0 in subscript to 0…
With these constraints, we get two estimates for \(\lambda^{XY}_{11} - \lambda^{XY}_{10} - \lambda^{XY}_{01} + \lambda^{XY}_{00}\) (equal to \(\lambda^{XY}_{11}\) under defaults), one from each logistic regression.
Actually three – one comes from the loglinear model.
Which interactions to include? (To be revisited later.)
Each of \(X,Y,Z,W\) has a corresponding logistic regression model formed by conditioning on its neighbours.
The pseudo-likelihood for parameters \(\Lambda = (\lambda, \lambda^X_i, \lambda^Y_j, \dots)\) is
\[ \begin{aligned} \log L_{\text{pseudo}}(\Lambda|X,Y,Z,W) &= \log L_{\text{logistic}}^Y(\Lambda|X,Z) + \\ & \qquad \log L_{\text{logistic}}^X(\Lambda|W,Y) + \\ & \qquad \log L_{\text{logistic}}^Z(\Lambda|Y,W) + \\ & \qquad \log L_{\text{logistic}}^W(\Lambda|X,Z) \end{aligned} \]
References:
Several more in between…
Suppose now we have vertices \(V = \{X_1, \dots, X_k\}\) where each \(X_i \in \{0,1\}\) and graph \(G=(E,V)\)
For binary data and usual constraints on parameters each vertex in \(V\) gets one parameter and each edge gets one parameter. One parameter \(\gamma\) for overall intensity.
Edge / vertex parameters can be represented by symmetric matrix \(\Theta\) with \(\Theta_{ij} \neq 0 \iff (i,j) \in E\).
Likelihood (with a little work, assuming we observe full table)
\[ \log L(\gamma,\Theta|X) = \gamma \cdot 2^k + \text{Tr}(\Theta S) - \text{sum}(N) \cdot \log \left( \sum_{x \in \{0,1\}^k} \exp(\gamma + x^T \Theta x) \right) \]
where
\[ S_{ij} = \begin{cases} \text{(1,1) entry in $X_i \times X_j$ marginal table} & i \neq j \\ \text{the number of 1's in $X_i$ marginal table.} & i=j \end{cases} \]
\[ S = B'B \]
Pseudo-likelihood is a sum of \(k\) different logistic regression log-likelihoods.
The parameters in the log-likelihoods are functions of \((\gamma, \Theta)\) – pseudolikelihood ties them all together.
Normalizing constant is simple.
Fit is very easy.
Since this is not a likelihood, cannot appeal to LR tests, Fisher information.
Bootstrap likely applicable.
If a variable \(X_i\) is not binary, then its pseudo-likelihood is baseline multinomial logit log-likelihood.
In this case, \(\Theta\) is not a matrix, every edge is represented by potentially a matrix of parameters.
\[ \log L_{\text{pseudo}}(\Theta) + \lambda \|\Theta_{\text{offdiag}}\|_1 \]
Besag’s original motivation was for inference in models similar to an Ising model and other models of spatial statistics.
Given graph \(G=(V,E)\) and a single observation \(X \in \{0,1\}^V\), we want to estimate parameters in a model
\[ P_{(\alpha,\beta)}(X) \propto \exp\left(\alpha \cdot \sum_{i \in V} X_i + \beta \cdot \sum_{(i,j) \in E} X_i X_j\right) \]
This is an exponential family:
Sufficient statistics: \(\left(\sum_{i \in V} X_i, \sum_{(i,j) \in E}X_i X_j \right)\)
Natural parameters: \((\alpha, \beta)\)
Reference measure: counting measure on \(\{0,1\}^V\).
\[ \begin{aligned} P_{(\alpha,\beta)}(X_i=1|X_{N(i)}) &\propto \exp\left(\alpha + \beta \sum_{j \in N(i)} X_j \right) \\ P_{(\alpha,\beta)}(X_i=0|X_{N(i)}) &\propto 1 \\ \end{aligned} \]
Differentiate w.r.t. \(\alpha\) \[ \sum_{i \in V} \left[X_i - \frac{\exp\left(\alpha + \beta \sum_{j \in N(i)} X_j \right)}{1 + \exp\left(\alpha + \beta \sum_{j \in N(i)} X_j \right)} \right] \]
Computing expectations at true \((\alpha^*,\beta^*)\):
\[ E_{(\alpha^*,\beta^*)} \left[\sum_{i \in V} \left[X_i - \frac{\exp\left(\alpha + \beta \sum_{j \in N(i)} X_j \right)}{1 + \exp\left(\alpha + \beta \sum_{j \in N(i)} X_j \right)} \right] \right] \]
At pseudo-MLE we set this equal to \[ 0= E_{(\alpha^*,\beta^*)} \left[\sum_{i \in V} \left[E_{(\alpha^*,\beta^*)}[X_i|X_{N(i)}] - \frac{\exp\left(\alpha^* + \beta^* \sum_{j \in N(i)} X_j \right)}{1 + \exp\left(\alpha^* + \beta^* \sum_{j \in N(i)} X_j \right)} \right] \right] \]
Pseudo-likelihood is an estimating equation…
Suppose we want to sample from \(P_{(\alpha,\beta)}\).
Distribution \(X_i | X_{-i}\) is the same as \(X_i | X_{N(i)}\).
\[ \begin{aligned} P_{(\alpha,\beta)}(X_i=1|X_{N(i)}) &\propto \exp\left(\alpha + \beta \sum_{j \in N(i)} X_j \right) \\ P_{(\alpha,\beta)}(X_i=0|X_{N(i)}) &\propto 1 \\ \end{aligned} \]
The “pieces” of pseudolikelihood are conditional likelihoods.
These conditional likelihoods are sometimes useful themselves following general principle of eliminating nuisance parameters.
In the homogeneous association model conditioning on all neighbours eliminate all nuisance parameters for all variables outside of \(\{v\} \cup N(v)\).
\[ \begin{aligned} -\log L(\beta, \sigma^2|X) &= \frac{1}{2 \sigma^2} \|Y-X\beta\|^2_2 + \frac{n}{2} \log(2\pi \sigma^2) \\ &= \frac{1}{2 \sigma^2} \|Y\|^2_2 - \frac{1}{\sigma^2}\beta^T(X^TY) + C(\sigma^2,\beta) \end{aligned} \]
An exponential family:
Sufficient statistic \((X^TY, \|Y\|^2_2)\).
Natural parameters \(\left(\beta / \sigma^2, -\frac{1}{2\sigma^2}\right)\)
Reference measure: \(dY\).
Consider the law \[P_{(\beta, \sigma^2)}\left(\|Y\|^2_2 \bigl \vert X^TY\right).\]
Not hard to see that this law is \(\|(I-H)Y\|^2_2 + \|HY\|^2_2\) where \(HY\) is measurable w.r.t. \(X^TY\).
Conditional likelihood for \(\sigma^2\) is based on \(P_{(\beta,\sigma^2)}(\|(I-H)Y\|^2_2)\) – a \(\chi^2_{n-p}\) likelihood.
Maximum (conditional) likelihood estimator is the usual unbiased estimate of \(\sigma^2\) not the biased MLE. This conditional likelihood estimator is better than the MLE!
This phenomenon has a life of its own in mixed effects model: REML (Residual Maximum Likelihood) (though elimination of nuisance parameter is not always justified by conditioning).
\[ \log L(\gamma,\Theta|X) = \gamma \cdot 2^k + \text{Tr}(\Theta S) - \sum_{x \in \{0,1\}^k} \exp(\gamma + x^T \Theta x) \]
Suppose data are continuous but marginally homogeneous.
A few choices yields Gaussian: \(X_s \sim N(0,\Sigma), 1 \leq s \leq n\)
\[ \log L(\Sigma) = - \frac{n}{2} \log |\Sigma| - \frac{1}{2} \text{Tr}(\Sigma^{-1}S) \]
\[ S = X'X = \sum_{s=1}^n X_sX_s' \]
\[ \log L(\Theta) = \frac{n}{2} \log |\Theta| - \frac{1}{2} \text{Tr}(\Theta S) \]
\[ \hat{\Theta} = \text{argmin}_{\Theta} \left[-\frac{n}{2} \log |\Theta| + \frac{1}{2} \text{Tr}(\Theta S) + \lambda \|\Theta_{\text{offdiag}}\|_1 \right] \]
\[ N_i \sim \text{Poisson}(\lambda_i) = \# \left\{s: X^D_s=i\right\} \]
Saturated: \(\lambda_i\) unconstrained
Ordinal with scores: \(\lambda_i = \beta_0 + \beta_1 \cdot i\)
\[ X^{\Delta} = i \iff X \in ((i-1)\Delta, i\Delta] \\ \]
\[ \exp(\lambda_i) \approx \Delta \cdot \exp \left( \gamma_0 + \gamma_1 (i \Delta) \right), \qquad \gamma_1 = \beta_1 / \Delta, \gamma_0 = \beta_0 - \log\Delta \]
\[ f_{\gamma}(x) \propto \exp(\gamma_0 + \gamma_1 \cdot x) \]
\[ \exp(-\gamma_0) = \int_{[0,1]} e^{\gamma_1 x} \; dx \]
No reason to restrict scores to just be \(i\)…
Can extend from \([0,1]\) to all of \(\mathbb{R}\) by adding “bins” of the form \([C,\infty)\) and \((-\infty,C]\) (with unrestricted parameters)
Thinking of index as \(i \cdot \Delta\) consider a model
\[ \exp(\lambda_i) = \Delta \cdot \exp \left(\gamma_0 + \sum_j \gamma_j h_j(i \cdot \Delta) \right) \]
\[ f(x) \propto \exp\left(\sum_j \gamma_j h_j(x) \right), \qquad \exp(-\gamma_0) = \int_{[-C,C]} \exp\left(\sum_j \gamma_j h_j(x) \right) \; dx \]
Bins need not be of same size – correspond to an offset in a GLM.
“Representative point” of a bin need not be the right hand endpoint, could be midpoint, left end point….
\[ \exp(\lambda_i) = \Delta_i \cdot \exp \left(\gamma_0 + \sum_j \gamma_j h_j(t_i) \right) \]
X = rnorm(20000)
bins = unique(c(seq(-4,-2,by=0.2),
seq(-2.1,2,by=0.1),
seq(2.2, 4, by=0.2)))
Y = table(cut(X, breaks=bins))
midpts = (bins[1:(length(bins)-1)] + bins[2:length(bins)])/2
widths = bins[2:length(bins)] - bins[1:(length(bins)-1)]
barplot(Y, widths)
df = cbind(Y, midpts, midpts^2)
colnames(df) = c('N', 'linear', 'quadratic')
df = data.frame(df)
M = glm(N ~ linear + quadratic,
family=poisson,
offset=log(widths),
data=df)## Warning in log(widths): NaNs produced##
## Call:
## glm(formula = N ~ linear + quadratic, family = poisson, data = df,
## offset = log(widths))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.2651 -0.5705 0.2531 1.0650 2.3851
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.994306 0.008699 1034.000 < 2e-16 ***
## linear 0.024856 0.007170 3.467 0.000526 ***
## quadratic -0.518815 0.005161 -100.524 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 27240.69 on 59 degrees of freedom
## Residual deviance: 145.03 on 57 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 558.3
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = N ~ linear + quadratic, family = poisson, data = df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2390 -0.9766 0.1222 1.3073 6.6948
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.649725 0.008816 754.252 <2e-16 ***
## linear 0.017415 0.007102 2.452 0.0142 *
## quadratic -0.446184 0.005311 -84.017 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 16794.95 on 60 degrees of freedom
## Residual deviance: 199.21 on 58 degrees of freedom
## AIC: 618.79
##
## Number of Fisher Scoring iterations: 4Suppose \((X,Y)\) are two discrete ordinal variables and \(N_{ij} \sim \text{Poisson}(\lambda_{ij})\) is an \(I\times J\) contingency table representing a sample from \((X,Y)\).
Ordinal (independence): \(\lambda_{ij} = \beta_0 + \beta_X i + \beta_Y j\).
Ordinal with dependence?: \(\lambda_{ij} = \beta_0 + \beta_X i + \beta_Y j + \beta^{XY}_{ij}\).
Suppose \((X^{\Delta}, Y^{\Delta})\) are such discretizations
Basis expansion (independence): \(\lambda_{ij} = \gamma_0 + \sum_j \gamma^X_j h^X_j(i\Delta) + \sum_l \gamma^Y_l h_l^Y(j\Delta)\)
\[ \sum_{j'l'} \gamma^{XY}_{j'l'} h^{XY}_{j'l'}(i\Delta,j\Delta). \]
Convenient to consider tensor products \(h^{XY}_{j'l'} = g^X_{j'l'}(x) g^Y_{j'l'}(y)\).
Example: if we set \(h^X_1(s)=h^Y_1(s)=s, h^X_2(s)=h^Y_2(s)=s^2\) and \(h^{XY}_{11}(x,y)=xy\) \(\implies\) Gaussian!
Use of tensor products can make conditionals nice (i.e. “simpler” pseudo-likelihood).
Suppose \(X\) is continuous and \(Y\) is discrete with \(K\) levels
For “main effects” of \(X\) we take \(h^X_1(x)=x, h^X_2(x)=x^2\).
For “main effects” of \(Y\) we take one-hot encoding \(h^Y_k(y) = \delta_j(y), 1 \leq k \leq K\).
For “interactions” we take \(h^{XY}_{1k}(x,y) = x \cdot \delta_k(y)\) – tensor product of (some) terms in the “main effect”. (Could also have the variance change as well.)
\(\implies\) conditionals \(X|Y\) and \(Y|X\) will be in the same family as implied by “main effect”.
\(X|Y\) is Gaussian (equivalent to lm(X ~ factor(Y)))
\(Y|X\) is multinomial (equivalent to a baseline multinomial logit nnet::multinom(Y ~ X)).
Suppose \(Z\) is discrete with \(L\) levels and \(Y\) is discrete with \(K\) levels
For “main effects” of \(Z\) we take one-hot encoding \(h^Z_l(l)=\delta_l(z) , 1 \leq l \leq L\).
For “main effects” of \(Y\) we take one-hot encoding \(h^Y_k(y) = \delta_j(y), 1 \leq k \leq K\).
For “interactions” we take \(h^{ZY}_{lk}(z,y) = \delta_l(z) \delta_k(y)\) – tensor product of terms in the “main effect”.
For all continuous variables use the Gaussian encoding with \(x \mapsto x, x \mapsto x^2\) as main effects.
For all discrete variables use the one-hot encoding \(x \mapsto \delta_l(x)\)
Use the tensor product formulation, but for any continuous interactions use only the \(x \mapsto x\) main effect.
Continuous / continuous: a single parameter \(\mathbb{R}\)
Continuous / discrete with \(L\) levels: a vector \(\mathbb{R}^L\)
Discrete with \(L\) levels / discrete with \(K\) levels: a matrix in \(\mathbb{R}^{L \times K}\).
For every continuous variable: Gaussian linear model
For every discrete variable: baseline multinomial logit
Want a penalty that can zero-out an entire vector (continuous / discrete) or matrix (discrete / discrete).
Consider
\[ \begin{aligned} \hat{\beta} &= \text{argmin}_{\beta} \frac{1}{2} \|Z-\beta\|^2_2 + \lambda \|\beta\|_2 \\ &= \begin{cases} 0 & \|Z\|_2 \leq \lambda \\ \frac{Z}{\|Z\|_2} \cdot (\|Z\|_2-\lambda, 0) & \text{otherwise.} \end{cases} \\ &= \max(\|Z\|_2-\lambda, 0) \cdot \frac{Z}{\|Z\|_2} \end{aligned} \]
\[ \|\beta\|_F = \text{Tr}(\beta'\beta)^{1/2} \]