fast generalized linear model fitting

fast generalized linear model fitting

fastglmPure(x, y, family = gaussian(), weights = rep(1, NROW(y)),
  offset = rep(0, NROW(y)), start = NULL, etastart = NULL,
  mustart = NULL, method = 0L, tol = 1e-07, maxit = 100L)

Arguments

x	input model matrix. Must be a matrix object
y	numeric response vector of length nobs.
family	a description of the error distribution and link function to be used in the model. For fastglmPure this can only be the result of a call to a family function. (See family for details of family functions.)
weights	an optional vector of 'prior weights' to be used in the fitting process. Should be a numeric vector.
offset	this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be a numeric vector of length equal to the number of cases
start	starting values for the parameters in the linear predictor.
etastart	starting values for the linear predictor.
mustart	values for the vector of means.
method	an integer scalar with value 0 for the column-pivoted QR decomposition, 1 for the unpivoted QR decomposition, 2 for the LLT Cholesky, 3 for the LDLT Cholesky, 4 for the full pivoted QR decomposition, 5 for the Bidiagonal Divide and Conquer SVD
tol	threshold tolerance for convergence. Should be a positive real number
maxit	maximum number of IRLS iterations. Should be an integer

Value

A list with the elements

coefficients

a vector of coefficients

se

a vector of the standard errors of the coefficient estimates

rank

a scalar denoting the computed rank of the model matrix

df.residual

a scalar denoting the degrees of freedom in the model

residuals

the vector of residuals

s

a numeric scalar - the root mean square for residuals

fitted.values

the vector of fitted values

Examples

set.seed(1)
x <- matrix(rnorm(1000 * 25), ncol = 25)
eta <- 0.1 + 0.25 * x[,1] - 0.25 * x[,3] + 0.75 * x[,5] -0.35 * x[,6] #0.25 * x[,1] - 0.25 * x[,3]
y <- 1 * (eta > rnorm(1000))

yp <- rpois(1000, eta ^ 2)
yg <- rgamma(1000, exp(eta) * 1.75, 1.75)

# binomial
system.time(gl1 <- glm.fit(x, y, family = binomial()))
#>    user  system elapsed 
#>   0.007   0.000   0.008 

system.time(gf1 <- fastglmPure(x, y, family = binomial(), tol = 1e-8))
#>    user  system elapsed 
#>   0.003   0.000   0.004 

system.time(gf2 <- fastglmPure(x, y, family = binomial(), method = 1, tol = 1e-8))
#>    user  system elapsed 
#>   0.004   0.001   0.003 

system.time(gf3 <- fastglmPure(x, y, family = binomial(), method = 2, tol = 1e-8))
#>    user  system elapsed 
#>   0.002   0.000   0.003 

system.time(gf4 <- fastglmPure(x, y, family = binomial(), method = 3, tol = 1e-8))
#>    user  system elapsed 
#>   0.002   0.000   0.002 

max(abs(coef(gl1) - gf1$coef))
#> [1] 1.998401e-15
max(abs(coef(gl1) - gf2$coef))
#> [1] 8.881784e-16
max(abs(coef(gl1) - gf3$coef))
#> [1] 1.776357e-15
max(abs(coef(gl1) - gf4$coef))
#> [1] 1.776357e-15

# poisson
system.time(gl1 <- glm.fit(x, yp, family = poisson(link = "log")))
#>    user  system elapsed 
#>   0.007   0.000   0.008 

system.time(gf1 <- fastglmPure(x, yp, family = poisson(link = "log"), tol = 1e-8))
#>    user  system elapsed 
#>   0.003   0.000   0.003 

system.time(gf2 <- fastglmPure(x, yp, family = poisson(link = "log"), method = 1, tol = 1e-8))
#>    user  system elapsed 
#>   0.003   0.000   0.003 

system.time(gf3 <- fastglmPure(x, yp, family = poisson(link = "log"), method = 2, tol = 1e-8))
#>    user  system elapsed 
#>   0.002   0.000   0.002 

system.time(gf4 <- fastglmPure(x, yp, family = poisson(link = "log"), method = 3, tol = 1e-8))
#>    user  system elapsed 
#>   0.002   0.000   0.002 

max(abs(coef(gl1) - gf1$coef))
#> [1] 3.538836e-16
max(abs(coef(gl1) - gf2$coef))
#> [1] 3.816392e-16
max(abs(coef(gl1) - gf3$coef))
#> [1] 2.63678e-16
max(abs(coef(gl1) - gf4$coef))
#> [1] 2.428613e-16

# gamma
system.time(gl1 <- glm.fit(x, yg, family = Gamma(link = "log")))
#>    user  system elapsed 
#>   0.020   0.001   0.022 

system.time(gf1 <- fastglmPure(x, yg, family = Gamma(link = "log"), tol = 1e-8))
#>    user  system elapsed 
#>   0.008   0.000   0.010 

system.time(gf2 <- fastglmPure(x, yg, family = Gamma(link = "log"), method = 1, tol = 1e-8))
#>    user  system elapsed 
#>   0.008   0.000   0.008 

system.time(gf3 <- fastglmPure(x, yg, family = Gamma(link = "log"), method = 2, tol = 1e-8))
#>    user  system elapsed 
#>   0.007   0.000   0.006 

system.time(gf4 <- fastglmPure(x, yg, family = Gamma(link = "log"), method = 3, tol = 1e-8))
#>    user  system elapsed 
#>   0.006   0.000   0.005 

max(abs(coef(gl1) - gf1$coef))
#> [1] 6.661338e-16
max(abs(coef(gl1) - gf2$coef))
#> [1] 5.551115e-16
max(abs(coef(gl1) - gf3$coef))
#> [1] 8.881784e-16
max(abs(coef(gl1) - gf4$coef))
#> [1] 8.881784e-16