glm2: Fitting Generalized Linear Models

Description

Fits generalized linear models using the same model specification as glm in the stats package, but with a modified default fitting method. The method provides greater stability for models that may fail to converge using glm.

Usage

glm2(formula, family = gaussian, data, weights, subset, na.action, 
     start = NULL, etastart, mustart, offset, control = list(...), 
     model = TRUE, method = "glm.fit2", x = FALSE, y = TRUE, 
     singular.ok = TRUE, contrasts = NULL, ...)

Arguments

formula

as for glm

family

as for glm

data

as for glm

weights

as for glm

subset

as for glm

na.action

as for glm

start

as for glm

etastart

as for glm

mustart

as for glm

offset

as for glm

control

as for glm except by default control is passed to glm.fit2 instead of glm.fit

model

as for glm

method

the method used in fitting the model. The default method "glm.fit2" uses iteratively reweighted least squares with modified step-halving that forces the deviance to decrease at each iteration; see help documentation for glm.fit2. As in glm, the alternative method "model.frame" returns the model frame and does no fitting.

as for glm

singular.ok

as for glm. NB this is ignored (and defaults to TRUE) for R versions < 3.5.0.

contrasts

as for glm

…

as for glm

Value

The value returned by glm2 has exactly the same documentation as the value returned by glm, except for:

method

the name of the fitter function used, which by default is "glm.fit2".

Details

glm2 is a modified version of glm in the stats package. It fits generalized linear models using the same model specification as glm. It is identical to glm except for minor modifications to change the default fitting method. The default method uses a stricter form of step-halving to force the deviance to decrease at each iteration and is implemented in glm.fit2. Like glm, user-supplied fitting functions can be used with glm2 by passing a function or a character string naming a function to the method argument. See Marschner (2011) for a discussion of the need for a modified fitting method.

References

Marschner, I.C. (2011) glm2: Fitting generalized linear models with convergence problems. The R Journal, Vol. 3/2, pp.12-15.

Examples

Run this code

# NOT RUN {
library(glm2)
data(crabs)
data(heart)

#==========================================================
# EXAMPLE 1: logistic regression null model
# (behaviour of glm and glm2 for different starting values)
#==========================================================

y <- c(1,1,1,0)
# intercept estimate = log(0.75/0.25) = 1.098612

#--- identical behaviour ---#
fit1 <- glm(y ~ 1, family=binomial(link="logit"),
	control=glm.control(trace=TRUE))
fit2 <- glm2(y ~ 1, family=binomial(link="logit"),
	control=glm.control(trace=TRUE))
print.noquote(c(fit1$coef,fit2$coef))

#--- convergence via different paths ---#
fit1 <- glm(y ~ 1, family=binomial(link="logit"),start=-1.75,
	control=glm.control(trace=TRUE))
fit2 <- glm2(y ~ 1, family=binomial(link="logit"),start=-1.75,
	control=glm.control(trace=TRUE))
print.noquote(c(fit1$coef,fit2$coef))

#--- divergence of glm to infinite estimate ---#
fit1 <- glm(y ~ 1, family=binomial(link="logit"),start=-1.81)
fit2 <- glm2(y ~ 1, family=binomial(link="logit"),start=-1.81)
print.noquote(c(fit1$coef,fit2$coef))


#=======================================================================
# EXAMPLE 2: identity link Poisson (successful boundary convergence
# using 4 identical approaches in glm and glm2 with the method argument) 
#=======================================================================

satellites <- crabs$Satellites
width.shifted <- crabs$Width - min(crabs$Width)
dark <- crabs$Dark
goodspine <- crabs$GoodSpine

fit1 <- glm(satellites ~ width.shifted + factor(dark) + factor(goodspine), 
 family = poisson(link="identity"), start = rep(1,4))

fit2 <- glm2(satellites ~ width.shifted + factor(dark) + factor(goodspine), 
 family = poisson(link="identity"), start = rep(1,4))

fit1.eq <- glm2(satellites ~ width.shifted + factor(dark) + factor(goodspine), 
 family = poisson(link="identity"), start = rep(1,4), method = "glm.fit")

fit2.eq <- glm(satellites ~ width.shifted + factor(dark) + factor(goodspine), 
 family = poisson(link="identity"), start = rep(1,4), method = "glm.fit2")

noquote(c("deviances: ",fit1$dev,fit2$dev,fit1.eq$dev,fit2.eq$dev))
noquote(c("converged: ",fit1$conv,fit2$conv,fit1.eq$conv,fit2.eq$conv))
noquote(c("boundary:  ",fit1$bound,fit2$bound,fit1.eq$bound,fit2.eq$bound))

#===================================================================
# EXAMPLE 3: identity link Poisson (periodic non-convergence in glm)
#===================================================================

R1 <- crabs$Rep1
satellites <- crabs$Satellites[R1]
width.shifted <- crabs$Width[R1] - min(crabs$Width)
dark <- crabs$Dark[R1]
goodspine <- crabs$GoodSpine[R1]

fit1 <- glm(satellites ~ width.shifted + factor(dark) + factor(goodspine), 
 family = poisson(link="identity"), start = rep(1,4), 
 control = glm.control(trace=TRUE))

fit2 <- glm2(satellites ~ width.shifted + factor(dark) + factor(goodspine), 
 family = poisson(link="identity"), start = rep(1,4), 
 control = glm.control(trace=TRUE))

noquote(c("deviances: ",fit1$dev,fit2$dev))
noquote(c("converged: ",fit1$conv,fit2$conv))

#===============================================================
# EXAMPLE 4: log link binomial (periodic non-convergence in glm)
#===============================================================

patients <- heart$Patients
deaths <- heart$Deaths
agegroup <- heart$AgeGroup
severity <-heart$Severity
delay <- heart$Delay
region <- heart$Region
start.p <- sum(deaths)/sum(patients)

fit1 <- glm(cbind(deaths,patients-deaths) ~ factor(agegroup) + factor(severity)
 + factor(delay) + factor(region), family = binomial(link="log"), 
 start = c(log(start.p), rep(0,8)), control = glm.control(trace=TRUE,maxit=100))

fit2 <- glm2(cbind(deaths,patients-deaths) ~ factor(agegroup) + factor(severity)
 + factor(delay) + factor(region), family = binomial(link="log"), 
 start = c(log(start.p), rep(0,8)), control = glm.control(trace=TRUE))

noquote(c("deviances: ",fit1$dev,fit2$dev))
noquote(c("converged: ",fit1$conv,fit2$conv))

#====================================================================
# EXAMPLE 5: identity link Poisson (aperiodic non-convergence in glm)
#====================================================================

R2 <- crabs$Rep2
satellites <- crabs$Satellites[R2]
width.shifted <- crabs$Width[R2] - min(crabs$Width)
dark <- crabs$Dark[R2]
goodspine <- crabs$GoodSpine[R2]

fit1 <- glm(satellites ~ width.shifted + factor(dark) + factor(goodspine), 
 family = poisson(link="identity"), start = rep(1,4), 
 control = glm.control(trace=TRUE))

fit2 <- glm2(satellites ~ width.shifted + factor(dark) + factor(goodspine), 
 family = poisson(link="identity"), start = rep(1,4), 
 control = glm.control(trace=TRUE))

noquote(c("deviances: ",fit1$dev,fit2$dev))
noquote(c("converged: ",fit1$conv,fit2$conv))

# }

Run the code above in your browser using DataCamp Workspace

Description

Usage

Arguments

Value

Details

References

See Also

Examples