loglinb2.RdFits a loglinear model to two binary responses.
loglinb2(exchangeable = FALSE, zero = "u12")Logical.
If TRUE, the two marginal probabilities are constrained
to be equal. Should be set TRUE for ears, eyes,
etc. data.
Which linear/additive predictors are modelled as
intercept-only?
A NULL means none of them.
See CommonVGAMffArguments for more information.
The model is $$P(Y_1=y_1,Y_2=y_2) = \exp(u_0+ u_1 y_1+u_2 y_2+u_{12} y_1 y_2)$$ where \(y_1\) and \(y_2\) are 0 or 1, and the parameters are \(u_1\), \(u_2\), \(u_{12}\). The normalizing parameter \(u_0\) can be expressed as a function of the other parameters, viz., $$u_0 = -\log[1 + \exp(u_1) + \exp(u_2) + \exp(u_1 + u_2 + u_{12})].$$ The linear/additive predictors are \((\eta_1,\eta_2,\eta_3)^T = (u_1,u_2,u_{12})^T\).
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 four joint probabilities, labelled as
\((Y_1,Y_2)\) = (0,0), (0,1), (1,0), (1,1),
respectively.
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.
The response must be a two-column matrix of ones and zeros only.
This is more restrictive than binom2.or,
which can handle more types of input formats.
Note that each of the 4 combinations of the multivariate response
need to appear in the data set.
After estimation, the response attached to the object is also a
two-column matrix; possibly in the future it might change into a
four-column matrix.
coalminers <- transform(coalminers, Age = (age - 42) / 5)
# Get the n x 4 matrix of counts
fit0 <- vglm(cbind(nBnW,nBW,BnW,BW) ~ Age, binom2.or, coalminers)
counts <- round(c(weights(fit0, type = "prior")) * depvar(fit0))
# Create a n x 2 matrix response for loglinb2()
# bwmat <- matrix(c(0,0, 0,1, 1,0, 1,1), 4, 2, byrow = TRUE)
bwmat <- cbind(bln = c(0,0,1,1), wheeze = c(0,1,0,1))
matof1 <- matrix(1, nrow(counts), 1)
newminers <-
data.frame(bln = kronecker(matof1, bwmat[, 1]),
wheeze = kronecker(matof1, bwmat[, 2]),
wt = c(t(counts)),
Age = with(coalminers, rep(age, rep(4, length(age)))))
newminers <- newminers[with(newminers, wt) > 0,]
fit <- vglm(cbind(bln,wheeze) ~ Age, loglinb2(zero = NULL),
weight = wt, data = newminers)
coef(fit, matrix = TRUE) # Same! (at least for the log odds-ratio)
#> u1 u2 u12
#> (Intercept) -7.8071500 -3.69405946 4.4551596
#> Age 0.1030801 0.04012099 -0.0332305
summary(fit)
#>
#> Call:
#> vglm(formula = cbind(bln, wheeze) ~ Age, family = loglinb2(zero = NULL),
#> data = newminers, weights = wt)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept):1 -7.807150 0.229187 -34.065 < 2e-16 ***
#> (Intercept):2 -3.694059 0.099871 -36.988 < 2e-16 ***
#> (Intercept):3 4.455160 0.278700 15.986 < 2e-16 ***
#> Age:1 0.103080 0.004476 23.030 < 2e-16 ***
#> Age:2 0.040121 0.002232 17.974 < 2e-16 ***
#> Age:3 -0.033230 0.005492 -6.051 1.44e-09 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Names of linear predictors: u1, u2, u12
#>
#> Log-likelihood: -12863.55 on 102 degrees of freedom
#>
#> Number of Fisher scoring iterations: 5
#>
#> Warning: Hauck-Donner effect detected in the following estimate(s):
#> '(Intercept):1', '(Intercept):2'
#>
# Try reconcile this with McCullagh and Nelder (1989), p.234
(0.166-0.131) / 0.027458 # 1.275 is approximately 1.25
#> [1] 1.274674