Binomial Family Function

Family function for fitting generalized linear models to binomial responses

binomialff(link = "logitlink", multiple.responses = FALSE,
       parallel = FALSE, zero = NULL, bred = FALSE, earg.link = FALSE)

Arguments

Link function; see Links and CommonVGAMffArguments for more information.

multiple.responses

Multivariate response? If TRUE, then the response is interpreted as \(M\) independent binary responses, where \(M\) is the number of columns of the response matrix. In this case, the response matrix should have \(Q\) columns consisting of counts (successes), and the weights argument should have \(Q\) columns consisting of the number of trials (successes plus failures).

If FALSE and the response is a (2-column) matrix, then the number of successes is given in the first column, and the second column is the number of failures.

parallel

A logical or formula. Used only if multiple.responses is TRUE. This argument allows for the parallelism assumption whereby the regression coefficients for a variable is constrained to be equal over the \(M\) linear/additive predictors. If parallel = TRUE then the constraint is not applied to the intercepts.

zero

An integer-valued vector specifying which linear/additive predictors are modelled as intercepts only. The values must be from the set {1,2,...,\(M\)}, where \(M\) is the number of columns of the matrix response. See CommonVGAMffArguments for more information.

Details at CommonVGAMffArguments.

bred

Details at CommonVGAMffArguments. Setting bred = TRUE should work for multiple responses (multiple.responses = TRUE) and all VGAM link functions; it has been tested for logitlink only (and it gives similar results to brglm but not identical), and further testing is required. One result from fitting bias reduced binary regression is that finite regression coefficients occur when the data is separable (see example below). Currently hdeff.vglm does not work when bred = TRUE.

Details

This function is largely to mimic binomial, however there are some differences.

When used with cqo and cao, it may be preferable to use the clogloglink link.

Value

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

References

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

Altman, M. and Gill, J. and McDonald, M. P. (2004). Numerical Issues in Statistical Computing for the Social Scientist, Hoboken, NJ, USA: Wiley-Interscience.

Ridout, M. S. (1990). Non-convergence of Fisher's method of scoring—a simple example. GLIM Newsletter, 20(6).

Author

Thomas W. Yee

Note

If multiple.responses is FALSE (default) then the response can be of one of two formats: a factor (first level taken as failure), or a 2-column matrix (first column = successes) of counts. The argument weights in the modelling function can also be specified as any vector of positive values. In general, 1 means success and 0 means failure (to check, see the y slot of the fitted object). Note that a general vector of proportions of success is no longer accepted.

The notation \(M\) is used to denote the number of linear/additive predictors.

If multiple.responses is TRUE, then the matrix response can only be of one format: a matrix of 1's and 0's (1 = success).

Fisher scoring is used. This can sometimes fail to converge by oscillating between successive iterations (Ridout, 1990). See the example below.

Warning

See the above note regarding bred.

The maximum likelihood estimate will not exist if the data is completely separable or quasi-completely separable. See Chapter 10 of Altman et al. (2004) for more details, and safeBinaryRegression and hdeff.vglm. Yet to do: add a sepcheck = TRUE, say, argument to further detect this problem and give an appropriate warning.

Examples

shunua <- hunua[sort.list(with(hunua, altitude)), ]  # Sort by altitude
fit <- vglm(agaaus ~ poly(altitude, 2), binomialff(link = clogloglink),
            data = shunua)
if (FALSE) { # \dontrun{
plot(agaaus ~ jitter(altitude), shunua, ylab = "Pr(Agaaus = 1)",
     main = "Presence/absence of Agathis australis", col = 4, las = 1)
with(shunua, lines(altitude, fitted(fit), col = "orange", lwd = 2)) } # }

# Fit two species simultaneously
fit2 <- vgam(cbind(agaaus, kniexc) ~ s(altitude),
             binomialff(multiple.responses = TRUE), data = shunua)
if (FALSE) { # \dontrun{
with(shunua, matplot(altitude, fitted(fit2), type = "l",
     main = "Two species response curves", las = 1)) } # }

# Shows that Fisher scoring can sometime fail. See Ridout (1990).
ridout <- data.frame(v = c(1000, 100, 10), r = c(4, 3, 3), n = rep(5, 3))
(ridout <- transform(ridout, logv = log(v)))
#>      v r n     logv
#> 1 1000 4 5 6.907755
#> 2  100 3 5 4.605170
#> 3   10 3 5 2.302585
# The iterations oscillates between two local solutions:
glm.fail <- glm(r / n ~ offset(logv) + 1, weight = n,
               binomial(link = 'cloglog'), ridout, trace = TRUE)
#> Deviance = 22.72791 Iterations - 1
#> Deviance = 28.29819 Iterations - 2
#> Deviance = 19.01799 Iterations - 3
#> Deviance = 21.07849 Iterations - 4
#> Deviance = 17.97632 Iterations - 5
#> Deviance = 18.59905 Iterations - 6
#> Deviance = 17.90326 Iterations - 7
#> Deviance = 18.42591 Iterations - 8
#> Deviance = 17.88036 Iterations - 9
#> Deviance = 18.37207 Iterations - 10
#> Deviance = 17.87175 Iterations - 11
#> Deviance = 18.3519 Iterations - 12
#> Deviance = 17.86832 Iterations - 13
#> Deviance = 18.34386 Iterations - 14
#> Deviance = 17.86692 Iterations - 15
#> Deviance = 18.34059 Iterations - 16
#> Deviance = 17.86634 Iterations - 17
#> Deviance = 18.33924 Iterations - 18
#> Deviance = 17.8661 Iterations - 19
#> Deviance = 18.33868 Iterations - 20
#> Deviance = 17.866 Iterations - 21
#> Deviance = 18.33845 Iterations - 22
#> Deviance = 17.86596 Iterations - 23
#> Deviance = 18.33835 Iterations - 24
#> Deviance = 17.86595 Iterations - 25
#> Deviance = 18.33831 Iterations - 1
#> Deviance = 17.86594 Iterations - 2
#> Deviance = 18.3383 Iterations - 3
#> Deviance = 17.86594 Iterations - 4
#> Deviance = 18.33829 Iterations - 5
#> Deviance = 17.86594 Iterations - 6
#> Deviance = 18.33829 Iterations - 7
#> Deviance = 17.86593 Iterations - 8
#> Deviance = 18.33829 Iterations - 9
#> Deviance = 17.86593 Iterations - 10
#> Deviance = 18.33829 Iterations - 11
#> Deviance = 17.86593 Iterations - 12
#> Deviance = 18.33829 Iterations - 13
#> Deviance = 17.86593 Iterations - 14
#> Deviance = 18.33829 Iterations - 15
#> Deviance = 17.86593 Iterations - 16
#> Deviance = 18.33829 Iterations - 17
#> Deviance = 17.86593 Iterations - 18
#> Deviance = 18.33829 Iterations - 19
#> Deviance = 17.86593 Iterations - 20
#> Deviance = 18.33829 Iterations - 21
#> Deviance = 17.86593 Iterations - 22
#> Deviance = 18.33829 Iterations - 23
#> Deviance = 17.86593 Iterations - 24
#> Deviance = 18.33829 Iterations - 25
coef(glm.fail)
#> (Intercept) 
#>   -5.157362 
# vglm()'s half-stepping ensures the MLE of -5.4007 is obtained:
vglm.ok <- vglm(cbind(r, n-r) ~ offset(logv) + 1,
               binomialff(link = clogloglink), ridout, trace = TRUE)
#> Iteration 1: deviance = 22.72791
#> Iteration 2: deviance = 28.29819
#> Taking a modified step.
#> Iteration  2 :  deviance = 18.08526
#> Iteration 3: deviance = 17.80972
#> Iteration 4: deviance = 18.20771
#> Taking a modified step.
#> Iteration  4 :  deviance = 17.46081
#> Iteration 5: deviance = 17.46624
#> Taking a modified step.
#> Iteration  5 :  deviance = 17.43666
#> Iteration 6: deviance = 17.43668
#> Taking a modified step.
#> Iteration  6 :  deviance = 17.43661
#> Iteration 7: deviance = 17.43661
#> Taking a modified step.
#> Iteration  7 :  deviance = 17.43661
coef(vglm.ok)
#> (Intercept) 
#>   -5.400638 
wsdm(vglm.ok)
#> (Intercept) 
#>   0.7105305 
#> attr(,"seems.okay")
#> [1] TRUE

# Separable data
set.seed(123)
threshold <- 0
bdata <- data.frame(x2 = sort(rnorm(nn <- 100)))
bdata <- transform(bdata, y1 = ifelse(x2 < threshold, 0, 1))
fit <- vglm(y1 ~ x2, binomialff(bred = TRUE),
            data = bdata, criter = "coef", trace = TRUE)
#> Iteration 1: coefficients = -0.10015601,  2.10784287
#> Iteration 2: coefficients = -0.11063311,  3.70413974
#> Iteration 3: coefficients = -0.10881416,  5.58548660
#> Iteration 4: coefficients = -0.11736263,  7.34812525
#> Iteration 5: coefficients = -0.13796225,  8.32167765
#> Iteration 6: coefficients = -0.14943194,  8.49745083
#> Iteration 7: coefficients = -0.15084587,  8.49588971
#> Iteration 8: coefficients = -0.15083084,  8.49572183
#> Iteration 9: coefficients = -0.15082982,  8.49572911
#> Iteration 10: coefficients = -0.15082987,  8.49572905
#> Iteration 11: coefficients = -0.15082987,  8.49572904
coef(fit, matrix = TRUE)  # Finite!!
#>             logitlink(prob)
#> (Intercept)      -0.1508299
#> x2                8.4957290
summary(fit, wsdm = TRUE, hdiff = 0.0001)
#> 
#> Call:
#> vglm(formula = y1 ~ x2, family = binomialff(bred = TRUE), data = bdata, 
#>     criter = "coef", trace = TRUE)
#> 
#> Coefficients: 
#>             Estimate Std. Error z value Pr(>|z|) WSDM    
#> (Intercept)  -0.1508     0.4626  -0.326    0.744 0.06    
#> x2            8.4957     2.0883   4.068 4.74e-05 0.78 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Name of linear predictor: logitlink(prob) 
#> 
#> Residual deviance: 15.0176 on 98 degrees of freedom
#> 
#> Log-likelihood: -7.5088 on 98 degrees of freedom
#> 
#> Number of Fisher scoring iterations: 11 
#> 
#> No Hauck-Donner effect found in any of the estimates
#> 
#> Max-WSDM (excluding intercepts): 0.780
#> 

if (FALSE)  plot(depvar(fit) ~ x2, data = bdata, col = "blue", las = 1)
lines(fitted(fit) ~ x2, data = bdata, col = "orange")
#> Error in plot.xy(xy.coords(x, y), type = type, ...): plot.new has not been called yet
abline(v = threshold, col = "gray", lty = "dashed")  # \dontrun{}
#> Error in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): plot.new has not been called yet