glim <- function(x, y, n, error = "gaussian", link = "identity", wt = 1, resid = "none", init, intercept = T, scale, offset = 0, sequence, eps = 0.0001, iter.max = 10) { error.int <- charmatch(error, c("gaussian", "poisson", "binomial", "gamma", "inverse gaussian")) if(is.na(error.int)) stop("Invalid error type") else if(error.int == 0) stop("Ambiguous error type") resid.int <- charmatch(resid, c("none", "pearson", "deviance")) if(is.na(resid.int)) stop("Invalid residual type") if(!missing(scale) && !is.numeric(scale)) { temp <- charmatch(scale, c("pearson", "deviance")) if(is.na(temp)) stop("Invalid scale option") } if(!(is.numeric(y) && is.numeric(x))) stop("Invalid y or x values") x <- as.matrix(x) nn <- nrow(x) nvar <- ncol(x) if(length(y) != nn) stop("Dimensions of x and y don't match") if(!missing(wt)) { if(length(wt) != nn) stop("Dimensions of x and wt don't match") if(any(wt < 0, na.rm = T)) stop("Weights must be >=0") } if(!missing(offset) && (length(offset) != nn)) stop( "Dimensions of x and offset don't match") # n is needed for the binomial, but is a dummy value for all others if(error.int != 3) n <- rep(1, nn) else { if(missing(n)) stop("Binomial fits require the n vector") else if(length(n) != nn) stop("Length of n vector is incorrect") n <- as.vector(n) } nomiss <- !(is.na(y) | is.na(wt) | is.na(n) | is.na(x %*% rep(1, nvar)) ) if(sum(nomiss) < nn) { warning(paste(nn - sum(nomiss), "deleted due to missing values" )) x <- x[nomiss, , drop = F] y <- y[nomiss] n <- n[nomiss] if(!missing(wt)) wt <- wt[nomiss] if(!missing(offset)) offset <- offset[nomiss] nn <- sum(nomiss) } if(missing(sequence)) sequence <- c(0, nvar) else { if(max(sequence) > nvar) stop("Invalid sequence argument") if(min(sequence) < 0) stop("Invalid sequence argument") if(any(diff(sequence) <= 0)) stop("Invalid sequence argument") } xn <- dimnames(x)[[2]] if(is.null(xn)) xn <- paste("X", 1:nvar, sep = "") if(intercept) { x <- cbind(rep(1, nn), x) xn <- c("Intercept", xn) nvar <- nvar + 1 sequence <- sequence + 1 } dimnames(x) <- list(dimnames(x)[[1]], xn) # Check for legal y values if(error.int != 1 && any(y < 0)) stop("Illegal y values") if(error.int >= 4 && any(y == 0)) stop("Illegal y values") if(error.int == 3) { y <- y/n if(any(y > 1)) stop("Illegal y values") } var <- switch(error.int, function(x, n) rep(1, length(x)), function(x, n) x, function(x, n) (x - x^2)/n, function(x, n) x^2, function(x, n) x^3) dev <- switch(error.int, function(y, mu, n) (y - mu)^2, function(y, mu, n) 2 * (y * log(ifelse(y == 0, 1, y/mu)) - (y - mu)), function(y, mu, n) 2 * n * (y * log(ifelse(y == 0, 1, y/mu)) + (1 - y) * log( ifelse(y == 1, 1, (1 - y)/(1 - mu)))), function(y, mu, n) -2 * (log(y/mu) - (y - mu)/mu), function(y, mu, n) (y - mu)^2/(y * mu^2)) # # Now create the link function-- this list could be easily expanded # f will be the inverse of the link, g the link, and deriv the # derivative of f. G is only needed for initial values, so # defer defining it until I need it. # link.int <- charmatch(link, c("identity", "log", "logit", "sqrt", "inverse", "probit", "loglog")) if(is.na(link.int)) stop("Invalid link type") else if(link.int == 0) stop("Ambiguous link type") f <- switch(link.int, function(x) x, function(x) exp(x), function(x) { z <- exp(ifelse(x > 80, 80, ifelse(x < -80, -80, x))) z/(z + 1) } , function(x) x^2, function(x) 1/x, pnorm, function(x) { z <- exp(ifelse(x > 80, 80, ifelse(x < -80, -80, x))) 1 - exp( - z) } ) deriv <- switch(link.int, function(x) rep(1, length(x)), function(x) exp(x), function(x) exp(x)/((1 + exp(x))^2), function(x) 2 * x, function(x) -1/(x^2), function(x) 0.3989422 * exp(-0.5 * x^2), function(x) exp(x) * exp( - exp(x))) temp.y <- switch(error.int, y, y + 0.01, (n * y + 0.5)/(n + 1), y, y) # here I need g, but only once, so never really define it eta <- switch(link.int, temp.y, log(temp.y), log(temp.y/(1 - temp.y)), sqrt(temp.y), 1/temp.y, qnorm(temp.y), log( - log(1 - temp.y))) if(sum(is.na(eta)) > 0) stop("Illegal value(s) of y for this link function") if(missing(init)) { tempd <- deriv(eta) z <- eta + (y - temp.y)/tempd - offset w <- sqrt((wt * tempd^2)/var(temp.y, n)) fit <- qr(x * w) if(fit$rank < nvar) stop("X matrix is not full rank") init.qr <- fit$qr[1:nvar, ] init.qty <- qr.qty(fit, z * w)[1:nvar] } else if(length(init) != nvar) stop("Argument 'init' is the wrong length") if(!missing(wt)) nn <- sum(wt > 0) #adjust df for obs with 0 weight # # Time to actually do the regression(s) # # Models with zero covars or only a mean are a special case # (but intercept + offset is too hard-- iterate it) deviance <- df <- NULL if(sequence[1] == 1 && intercept && missing(offset)) { if(missing(wt)) yhat <- sum(y * n)/sum(n) else yhat <- sum((y * n * wt)[wt > 0])/sum((n * wt)[wt > 0]) if((yhat > 1) && (link.int == 3 || link.int == 6 || link.int == 7)) yhat <- 1 if((yhat < 0) && link.int > 1 && link.int != 5) yhat <- 0 deviance <- sum(dev(y, yhat, n) * wt) df <- nn - 1 sequence <- sequence[-1] } else if(sequence[1] == 0) { deviance <- sum(dev(y, f(offset), n) * wt) df <- nn sequence <- sequence[-1] } for(modnum in sequence) { model <- 1:modnum if(missing(init)) coef <- backsolve(init.qr, init.qty, k = modnum)[model] else coef <- init[model] tempx <- as.matrix(x[, model, drop = F]) nvar <- length(model) eta <- c(tempx %*% coef) #make it a vector yhat <- f(eta + offset) olddev <- sum(dev(y, yhat, n) * wt) fini <- F for(i in seq(length = iter.max)) { tempd <- deriv(eta + offset) if(any(is.na(yhat) | is.na(tempd))) { warning( "Coef vector is diverging, iteration truncated" ) break } leave.in <- ((1 + tempd) != 1) z <- eta + (y - yhat)/tempd w <- sqrt((wt * tempd^2)/var(yhat, n)) fit <- qr(tempx[leave.in, ] * w[leave.in]) if(fit$rank < nvar) { warning(paste( "Weighted X matrix no longer of full rank at iteration", i)) break } coef <- qr.coef(fit, (z * w)[leave.in]) eta <- c(tempx %*% coef) yhat <- f(eta + offset) newdev <- sum(dev(y, yhat, n) * wt) if(abs((olddev - newdev)/(newdev + 1)) < eps) { fini <- T break } else olddev <- newdev } if(fini == F) warning(paste("Model with", length(model), "variables did not converge")) deviance <- c(deviance, newdev) df <- c(df, nn - nvar) } # Put together the results vector coef <- as.vector(coef) names(coef) <- dimnames(tempx)[[2]] if(missing(scale)) scale <- switch(error.int, newdev/(nn - nvar), 1, 1, newdev/(nn - nvar), newdev/(nn - nvar)) else if(!is.numeric(scale)) scale <- switch(charmatch(scale, c("pearson", "deviance")), sum((y - yhat)^2/(var(yhat) * (nn - nvar))), newdev/(nn - nvar)) varmat <- backsolve(fit$qr[1:nvar, ], diag(nvar)) varmat <- varmat %*% t(varmat) temp <- list(coef = coef, var = varmat * scale, deviance = deviance, df = df, scale = scale, intercept = intercept) if(resid.int == 1) temp else if(resid.int == 3) { resid <- rep(NA, length(nomiss)) resid.good <- sqrt(dev(y, yhat, n)) resid[nomiss] <- ifelse(y > yhat, resid.good, - resid.good) c(temp, list(resid = resid)) } else { resid <- rep(NA, length(nomiss)) resid[nomiss] <- (y - yhat)/ifelse(y == yhat, 1, sqrt(var(yhat, n))) c(temp, list(resid = resid)) } }