Loglinear Model for Three Binary Responses

Fits a loglinear model to three binary responses.

Usage

loglinb3(exchangeable = FALSE, zero = c("u12", "u13", "u23",
         if (u123.arg) "u123" else NULL),
         u123.arg = FALSE)

Arguments

exchangeable

Logical. If TRUE, the three marginal probabilities are constrained to be equal, as well as their second-order interaction terms.

zero

Which linear/additive predictors are modelled as intercept-only? A NULL means none. See CommonVGAMffArguments for further information.

u123.arg

Logical. Include the 3rd-order interaction u123? Note that u123.arg = TRUE might become the default some time in the future.

Details

The full model is $P(Y_1=y_1,Y_2=y_2,Y_3=y_3) =$ $$\exp(u_0+u_1 y_1+u_2 y_2+u_3 y_3+u_{12} y_1 y_2+ u_{13} y_1 y_3+u_{23} y_2 y_3 + u_{123} y_1 y_2 y_3 )$$ where $y_1$, $y_2$ and $y_3$ are 0 or 1, and the parameters are $u_1$, $u_2$, $u_3$, $u_{12}$, $u_{13}$, $u_{23}$, and if u123.arg then $u_{123}$ too. The normalizing parameter $u_0$ can be expressed as a function of the other parameters. The the parameters are estimated by identitylink. Note that a third-order association parameter, $u_{123}$ for the product $y_1 y_2 y_3$, is assumed to be zero for this family function by default; it is estimated if u123.arg is TRUE. Note the default for this argument might change in the future.

The linear/additive predictors are, for the full model, $(\eta_1,\eta_2,\ldots,\eta_6,\eta_7)^T = (u_1,u_2,u_3,u_{12},u_{13},u_{23},u_{123})^T$. By default, the last element is not there since u123.arg = FALSE.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, rrvglm and vgam.

When fitted, the fitted.values slot of the object contains the eight joint probabilities, labelled as $(Y_1,Y_2,Y_3)$ = (0,0,0), (0,0,1), (0,1,0), (0,1,1), (1,0,0), (1,0,1), (1,1,0), (1,1,1), respectively.

References

Yee, T. W. and Wild, C. J. (2001). Discussion to: “Smoothing spline ANOVA for multivariate Bernoulli observations, with application to ophthalmology data (with discussion)” by Gao, F., Wahba, G., Klein, R., Klein, B. Journal of the American Statistical Association, 96, 127–160.

McCullagh, P. and Nelder, J. A. (1989). Generalized Linear Models, 2nd ed. London: Chapman & Hall.

Author

Thomas W. Yee and Yunhao (Harry) Han who added u123.

Note

The response must be a 3-column matrix of ones and zeros only. Note that each of the 8 combinations of the multivariate response need to appear in the data set, therefore data sets will need to be large in order for this family function to work. After estimation, the response attached to the object is also a 3-column matrix; possibly in the future it might change into a 8-column matrix.

Examples

lfit1 <- vglm(cbind(cyadea, beitaw, kniexc) ~ altitude,
             loglinb3, data = hunua, trace = TRUE)
#> Iteration 1: loglikelihood = -748.25937
#> Iteration 2: loglikelihood = -746.63602
#> Iteration 3: loglikelihood = -746.63166
#> Iteration 4: loglikelihood = -746.63166
coef(lfit1, matrix = TRUE)
#>                       u1           u2           u3       u12       u13      u23
#> (Intercept) -0.977471200 -1.890117628 -0.377188302 0.6079725 0.1550935 1.117167
#> altitude    -0.000570146  0.003850299  0.001611092 0.0000000 0.0000000 0.000000
lfit2 <- vglm(cbind(cyadea, beitaw, kniexc) ~ altitude,
             loglinb3(u123 = TRUE), hunua, trace = TRUE)
#> Iteration 1: loglikelihood = -745.42318
#> Iteration 2: loglikelihood = -745.30255
#> Iteration 3: loglikelihood = -745.3025
#> Iteration 4: loglikelihood = -745.3025
coef(lfit2, matrix = TRUE)
#>                        u1           u2           u3        u12        u13
#> (Intercept) -0.8326528535 -1.723757823 -0.298068984 0.09095367 -0.1418004
#> altitude    -0.0006060001  0.003856669  0.001624325 0.00000000  0.0000000
#>                   u23      u123
#> (Intercept) 0.8564277 0.7907418
#> altitude    0.0000000 0.0000000
head(fitted(lfit2))
#>         000       001        010       011       100        101        110
#> 1 0.2567062 0.2205342 0.06479800 0.1310820 0.1057142 0.07881153 0.02922532
#> 2 0.2618034 0.2212893 0.06358448 0.1265546 0.1084686 0.07956207 0.02885232
#> 3 0.2769515 0.2229594 0.05991447 0.1135784 0.1168498 0.08163323 0.02768578
#> 4 0.2819412 0.2233193 0.05868636 0.1094578 0.1196781 0.08226198 0.02728312
#> 5 0.2819412 0.2233193 0.05868636 0.1094578 0.1196781 0.08226198 0.02728312
#> 6 0.2769515 0.2229594 0.05991447 0.1135784 0.1168498 0.08163323 0.02768578
#>          111
#> 1 0.11312857
#> 2 0.10988522
#> 3 0.10042742
#> 4 0.09737225
#> 5 0.09737225
#> 6 0.10042742
summary(lfit2)
#> 
#> Call:
#> vglm(formula = cbind(cyadea, beitaw, kniexc) ~ altitude, family = loglinb3(u123 = TRUE), 
#>     data = hunua, trace = TRUE)
#> 
#> Coefficients: 
#>                 Estimate Std. Error z value Pr(>|z|)    
#> (Intercept):1 -0.8326529  0.2320303  -3.589 0.000333 ***
#> (Intercept):2 -1.7237578  0.2703873  -6.375 1.83e-10 ***
#> (Intercept):3 -0.2980690  0.2053965  -1.451 0.146728    
#> (Intercept):4  0.0909537  0.4012030   0.227 0.820655    
#> (Intercept):5 -0.1418004  0.2972991  -0.477 0.633389    
#> (Intercept):6  0.8564277  0.2773943   3.087 0.002019 ** 
#> (Intercept):7  0.7907418  0.4906812   1.612 0.107067    
#> altitude:1    -0.0006060  0.0009225  -0.657 0.511250    
#> altitude:2     0.0038567  0.0009619   4.010 6.08e-05 ***
#> altitude:3     0.0016243  0.0009668   1.680 0.092937 .  
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Number of linear predictors:  7 
#> 
#> Names of linear predictors: u1, u2, u3, u12, u13, u23, u123
#> 
#> Log-likelihood: -745.3025 on 2734 degrees of freedom
#> 
#> Number of Fisher scoring iterations: 4 
#>