Loglinear Model for Four Binary Responses
loglinb4.RdFits a loglinear model to four binary responses.
Arguments
- exchangeable
Logical. If
TRUE, the four marginal probabilities are constrained to be equal, as well as their higher-order interaction terms.- zero
Which linear/additive predictors are modelled as intercept-only? A
NULLmeans none. SeeCommonVGAMffArgumentsfor further information.- order4
Logical or either 2 or 3 or 4. If logical and
TRUEthen 4, else 3. It ends up an integer that is either 2 or 3 or 4. Any higher-order term is 0, e.g.,u1234 = 0iforder4 = 3.
Details
The full model is
\(P(Y_1=y_1,Y_2=y_2,Y_3=y_3,Y_4=y_4) =\)
$$\exp(u_0+u_1 y_1+u_2 y_2+u_3 y_3+u_4 y_4+
u_{12} y_1 y_2+
u_{13} y_1 y_3+
u_{14} y_1 y_4 +
u_{23} y_2 y_3 +
u_{24} y_2 y_4 +
u_{34} y_3 y_4 +
u_{123} y_1 y_2 y_3 +
u_{124} y_1 y_2 y_4 +
u_{134} y_1 y_3 y_4 +
u_{234} y_2 y_3 y_4 +
u_{1234} y_1 y_2 y_3 y_4)$$
where \(y_1\), \(y_2\) and
\(y_3\) and
\(y_4\) are 0
or 1, and the parameters are \(u_1\),
\(u_2\),
\(u_3\),
\(u_4\),
\(u_{12}\),
\(u_{13}\),
\(u_{14}\),
\(u_{23}\),
\(u_{24}\),
\(u_{34}\), and if
order4 >= 3 then
\(u_{123}\),
\(u_{124}\),
\(u_{134}\),
\(u_{234}\), too, and if
order4 == 4 then
\(u_{1234}\) too.
The normalizing parameter \(u_0\) can
be expressed as a function of the others.
The the parameters are estimated by
identitylink.
Unlike loglinb3,
a fourth-order (full) association parameter,
\(u_{1234}\) for the product
\(y_1 y_2 y_3 y_4\),
is not assumed to be zero for this
family function by default.
Note the default for this argument might change
in the future.
If the data cannot support such a high
order interaction term then reduce order4.
The linear/additive predictors are,
for the full model,
\((\eta_1,\eta_2,\ldots,\eta_{15})^T =
(u_1,u_2,u_3,u_4,u_{12},u_{13},\ldots,u_{1234})^T\).
The ordering agrees with combn
piecemeal.
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 joint probabilities,
labelled as
\((Y_1,Y_2,Y_3,Y_4)\) =
(0,0,0,0), (0,0,0,1), (0,0,1,0),
(0,0,1,1), (0,1,0,0), (0,1,0,1), (0,1,1,0),
(0,1,1,1),
(1,0,0,0), (1,0,0,1), (1,0,1,0),
(1,0,1,1), (1,1,0,0), (1,1,0,1), (1,1,1,0),
(1,1,1,1),
respectively.
Note
The response must be a 4-column matrix
of ones and zeros only.
By default, each of the 16 combinations of
the multivariate
response need to appear in the data set,
therefore data sets will usually
need to be large in order for this
family function to work most simply.
As stated above, reduce order4
if there are problems.
After estimation, the response attached
to the object is also
a 4-column matrix; possibly in the future
it might change into
a 16-column matrix.
Examples
lfit4 <- vglm(cbind(cyadea, beitaw, kniexc, vitluc) ~ altitude,
loglinb4, hunua, trace = TRUE)
#> Iteration 1: loglikelihood = -906.60564
#> Iteration 2: loglikelihood = -899.48683
#> Iteration 3: loglikelihood = -898.85592
#> Iteration 4: loglikelihood = -898.8472
#> Iteration 5: loglikelihood = -898.8472
coef(lfit4, matrix = TRUE)
#> u1 u2 u3 u4 u12
#> (Intercept) -0.721809701 -2.280219633 -0.321403594 -0.4123939 0.1916423
#> altitude -0.001068434 0.005287642 0.001657184 -0.0130585 0.0000000
#> u13 u14 u23 u24 u34 u123
#> (Intercept) -0.211515 -0.3812867 0.9101037 1.480132 0.1170015 0.8644671
#> altitude 0.000000 0.0000000 0.0000000 0.000000 0.0000000 0.0000000
#> u124 u134 u234 u1234
#> (Intercept) -0.09089224 0.5203574 -0.2296513 -0.5969292
#> altitude 0.00000000 0.0000000 0.0000000 0.0000000
head(fitted(lfit4))
#> 0000 0001 0010 0011 0100 0101 0110
#> 1 0.08127079 0.01661204 0.06841073 0.01571903 0.01337590 0.01201227 0.02797473
#> 2 0.08162848 0.01901265 0.06758252 0.01769491 0.01274284 0.01304009 0.02621272
#> 3 0.08094363 0.02789463 0.06376526 0.02470217 0.01078240 0.01632551 0.02110426
#> 4 0.08008875 0.03145008 0.06205488 0.02739298 0.01011907 0.01745838 0.01948040
#> 5 0.08008875 0.03145008 0.06205488 0.02739298 0.01011907 0.01745838 0.01948040
#> 6 0.08094363 0.02789463 0.06376526 0.02470217 0.01078240 0.01632551 0.02110426
#> 0111 1000 1001 1010 1011 1100
#> 1 0.02244629 0.03586701 0.005007177 0.02443572 0.006452445 0.007150114
#> 2 0.02396641 0.03641183 0.005792322 0.02439919 0.007341537 0.006884880
#> 3 0.02854947 0.03728241 0.008775081 0.02377090 0.010582646 0.006015421
#> 4 0.03002882 0.03728490 0.009999826 0.02338178 0.011861470 0.005705991
#> 5 0.03002882 0.03728490 0.009999826 0.02338178 0.011861470 0.005705991
#> 6 0.02854947 0.03728241 0.008775081 0.02377090 0.010582646 0.006015421
#> 1101 1110 1111
#> 1 0.004004516 0.02872957 0.01331645
#> 2 0.004393856 0.02720919 0.01437101
#> 3 0.005680054 0.02262007 0.01767676
#> 4 0.006139456 0.02110387 0.01879243
#> 5 0.006139456 0.02110387 0.01879243
#> 6 0.005680054 0.02262007 0.01767676
head(predict(lfit4))
#> u1 u2 u3 u4 u12 u13 u14
#> [1,] -0.8179688 -1.804332 -0.1722570 -1.5876588 0.1916423 -0.211515 -0.3812867
#> [2,] -0.8072845 -1.857208 -0.1888289 -1.4570738 0.1916423 -0.211515 -0.3812867
#> [3,] -0.7752314 -2.015838 -0.2385444 -1.0653189 0.1916423 -0.211515 -0.3812867
#> [4,] -0.7645471 -2.068714 -0.2551162 -0.9347339 0.1916423 -0.211515 -0.3812867
#> [5,] -0.7645471 -2.068714 -0.2551162 -0.9347339 0.1916423 -0.211515 -0.3812867
#> [6,] -0.7752314 -2.015838 -0.2385444 -1.0653189 0.1916423 -0.211515 -0.3812867
#> u23 u24 u34 u123 u124 u134 u234
#> [1,] 0.9101037 1.480132 0.1170015 0.8644671 -0.09089224 0.5203574 -0.2296513
#> [2,] 0.9101037 1.480132 0.1170015 0.8644671 -0.09089224 0.5203574 -0.2296513
#> [3,] 0.9101037 1.480132 0.1170015 0.8644671 -0.09089224 0.5203574 -0.2296513
#> [4,] 0.9101037 1.480132 0.1170015 0.8644671 -0.09089224 0.5203574 -0.2296513
#> [5,] 0.9101037 1.480132 0.1170015 0.8644671 -0.09089224 0.5203574 -0.2296513
#> [6,] 0.9101037 1.480132 0.1170015 0.8644671 -0.09089224 0.5203574 -0.2296513
#> u1234
#> [1,] -0.5969292
#> [2,] -0.5969292
#> [3,] -0.5969292
#> [4,] -0.5969292
#> [5,] -0.5969292
#> [6,] -0.5969292
summary(lfit4, HDEtest = FALSE)
#>
#> Call:
#> vglm(formula = cbind(cyadea, beitaw, kniexc, vitluc) ~ altitude,
#> family = loglinb4, data = hunua, trace = TRUE)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept):1 -0.7218097 0.2531669 -2.851 0.00436 **
#> (Intercept):2 -2.2802196 0.3332867 -6.842 7.83e-12 ***
#> (Intercept):3 -0.3214036 0.2303060 -1.396 0.16285
#> (Intercept):4 -0.4123939 0.3518526 -1.172 0.24117
#> (Intercept):5 0.1916423 0.4624818 0.414 0.67860
#> (Intercept):6 -0.2115150 0.3226579 -0.656 0.51212
#> (Intercept):7 -0.3812867 0.6141001 -0.621 0.53467
#> (Intercept):8 0.9101037 0.3259687 2.792 0.00524 **
#> (Intercept):9 1.4801323 0.5203235 2.845 0.00445 **
#> (Intercept):10 0.1170015 0.4413630 0.265 0.79094
#> (Intercept):11 0.8644671 0.5572768 1.551 0.12085
#> (Intercept):12 -0.0908922 0.9914524 -0.092 0.92696
#> (Intercept):13 0.5203574 0.8450575 0.616 0.53805
#> (Intercept):14 -0.2296513 0.6585910 -0.349 0.72731
#> (Intercept):15 -0.5969292 1.2344548 -0.484 0.62870
#> altitude:1 -0.0010684 0.0009909 -1.078 0.28091
#> altitude:2 0.0052876 0.0010777 4.907 9.27e-07 ***
#> altitude:3 0.0016572 0.0010384 1.596 0.11051
#> altitude:4 -0.0130585 0.0020796 -6.279 3.40e-10 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Number of linear predictors: 15
#>
#> Log-likelihood: -898.8472 on 5861 degrees of freedom
#>
#> Number of Fisher scoring iterations: 5
#>