Skip to content

Instantly share code, notes, and snippets.

@ricardocarvalhods
Forked from MrFlick/makeglm.R
Last active August 29, 2015 14:25
Show Gist options
  • Save ricardocarvalhods/bc339ef349bb2ea61768 to your computer and use it in GitHub Desktop.
Save ricardocarvalhods/bc339ef349bb2ea61768 to your computer and use it in GitHub Desktop.
makeglm <- function(formula, ..., family, data=NULL) {
dots <- list(...)
out<-list()
tt <- terms(formula, data=data)
if(!is.null(data)) {
mf <- model.frame(tt, data)
vn <- sapply(attr(tt, "variables")[-1], deparse)
if((yvar <- attr(tt, "response"))>0)
vn <- vn[-yvar]
xlvl <- lapply(data[vn], function(x) if (is.factor(x))
levels(x)
else if (is.character(x))
levels(as.factor(x))
else
NULL)
attr(out, "xlevels") <- xlvl[!vapply(xlvl,is.null,NA)]
attr(tt, "dataClasses") <- sapply(data[vn], stats:::.MFclass)
}
out$terms <- tt
coef <- numeric(0)
stopifnot(length(dots)>1 & !is.null(names(dots)))
for(i in seq_along(dots)) {
if((n<-names(dots)[i]) != "") {
v <- dots[[i]]
if(!is.null(names(v))) {
coef[paste0(n, names(v))] <- v
} else {
stopifnot(length(v)==1)
coef[n] <- v
}
} else {
coef["(Intercept)"] <- dots[[i]]
}
}
out$coefficients <- coef
out$rank <- length(coef)
if (!missing(family)) {
out$family <- if (class(family) == "family") {
family
} else if (class(family) == "function") {
family()
} else if (class(family) == "character") {
get(family)()
} else {
stop(paste("invalid family class:", class(family)))
}
out$qr <- list(pivot=seq_len(out$rank))
out$deviance <- 1
out$null.deviance <- 1
out$aic <- 1
class(out) <- c("glm","lm")
} else {
class(out) <- "lm"
out$fitted.values <- predict(out, newdata=dd)
out$residuals <- out$mf[attr(tt, "response")] - out$fitted.values
out$df.residual <- nrow(data) - out$rank
out$model <- data
#QR doesn't work
}
out
}
set.seed(15)
dd <- data.frame(
X1=runif(50),
X2=factor(sample(letters[1:4], 50, replace=T)),
X3=rpois(50, 5),
Outcome = sample(0:1, 50, replace=T)
)
mymodel<-glm(Outcome~X1+X2+X3, data=dd, family=binomial)
predict(mymodel, type="response")
newmodel <- makeglm(Outcome~X1+X2+X3, family=binomial, data=dd,
-.5, X1=1, X2=c(b=1.5, c=1, d=1.5), X3=-.15)
predict(newmodel, newdata=dd, type="response")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment