Here is logistic regression R code that will work with the current toolbox (all the dependencies are correct). Copy the code into a text editor then save and replace "LogitWithR.r" in the "Scripts" directory of the toolbox.##############################################################################
# PROGRAM: LogitRegression (logistic regression called form ArcGIS)
# USE: LOGISTIC (BIONOMIAL) REGRESSION
# REQUIRES: R > 2.15.0, rms
#
# ARGUMENTS:
# ldata DATAFRAME OBJECT CONTANING VARIABLES
# y RESPONSE VARIABLE IN ldata
# x INDEPENDENT VARIABLES(S) IN ldata
# penelty APPLY REGRESSION PENELTY (TRUE/FALSE)
# ... ADDITIONAL ARGUMENTS PASSED TO lrm
#
# VALUE:
# A LIST OBJECT CONTANING OBJECTS:
# model lrm MODEL OBJECT
# diagTable DATAFRAME OF REGRESSION DIAGNOSTICS
# coefTable DATAFRAME OF REGRESSION COEFFICENTS
# Residuals DATAFRAME OF RESIDUALS AND STANDARDIZED RESIDUALS
#
# REFERENCES:
# Le Cessie S, Van Houwelingen JC: Ridge estimators in logistic regression.
# Applied Statistics 41:191�??201, 1992.
#
# Shao J: Linear model selection by cross-validation. JASA 88:486�??494, 1993.
#
# EXAMPLES:
# require(sp)
# require(rms)
#
# data(meuse)
# coordinates(meuse) <- ~x+y
# meuse@data <- data.frame(DepVar=rbinom(dim(meuse)[1], 1, 0.5), meuse@data)
# names(meuse@data)
#
# # RUN Bernoulli Trials TEST
# binom.test( c(length(meuse[meuse@data$DepVar == 1 ,]$DepVar),
# length(meuse[meuse@data$DepVar == 0 ,]$DepVar)), 0.65)
#
# # RUN LOGISTIC MODEL
# lmodel <- LogitRegression(meuse@data, y="DepVar", x=c("dist","cadmium","copper"))
# lmodel$model
# lmodel$diagTable
# lmodel$coefTable
#
# # ADD RESIDUALS, STANDARDIZED RESIDUALS AND PROBABILITIES
# meuse@data <- data.frame(meuse@data, Residual=lmodel$Residuals[,1],
# StdResid=lmodel$Residuals[,2],
# Probs=predict(lmodel$model,meuse@data[,rownames(lmodel$model$var)[-1]],
# type="fitted") )
#
# # PLOT PROBABILITIES
# library(lattice)
# trellis.par.set(sp.theme())
# spplot(meuse, c("Probs"))
#
# CONTACT:
# Jeffrey S. Evans
# Senior Landscape Ecologist
# The Nature Conservancy - Central Science
# Adjunct Assistant Professor
# Zoology and Physiology
# University of Wyoming
# Laramie, WY
# (970)672-6766
# jeffrey_evans@tnc.org
######################################################################################
LogitRegression <- function (ldata, y, x, penalty=TRUE, ...) {
if (!require (rms)) stop("sp PACKAGE MISSING")
if ( is.na(match(y,names(ldata))) )
stop("Dependent Variable not present in data")
xNames <- intersect(x,names(ldata))
if (length(xNames) < length(x))
stop("Mismatch in Independent Variable Names")
form=as.formula(paste(y, paste(x, collapse='+'), sep='~'))
fit <- lrm(form, data=ldata, x=TRUE, y=TRUE, ...)
bf <- pentrace(fit, seq(0.2,1,by=0.05))
if (penalty) {
pen = bf$penalty
} else {
pen = 0.0
}
allPens=bf$results.all[,1]
allAICs=bf$results.all[,3]
for (i in 1:length(allPens)){
penValue=allPens
if (penValue == pen){
aic = allAICs
}
}
if (penalty){ fit = update(fit, penalty=bf$penalty) }
res <- residuals.lrm(fit)
resSTD <- (res - mean(res)) / sqrt(var(res))
allIndVars <- c("Intercept")
allIndVars <- append(allIndVars, x)
k <- length(allIndVars)
d <- matrix(0, k, 4)
d[,1] <- fit$coefficients
d[,2] <- sqrt(diag(fit$var))
d[,3] <- d[,1] / d[,2]
d[,4] <- pnorm(abs(d[,3]), lower.tail = FALSE) * 2.0
coefList=list("Variable"=allIndVars, "Coef"=d[,1],
"StdError"=d[,2], "Wald"=d[,3], "Prob"=d[,4])
coefFrame=data.frame(coefList)
diagFrame=data.frame(Names=c(names(fit$stats),"PEN","AIC"),
Value=c(as.vector(fit$stats), pen, aic))
return(list(model=fit, diagTable=diagFrame, coefTable=coefFrame,
Residuals=data.frame(res=res,resSTD=resSTD)))
}
### Version check
Rv <- paste(round(as.numeric(R.version$major),digits=0),
as.numeric(R.version$minor),sep=".")
if (Rv < 2.14) stop("NEED TO UPGRADE VERSION OF R TO >= 2.14")
#### Load required packages
print("Loading Libraries....")
if (!require(rms)) stop("rms PACKAGE MISSING")
if (!require(sp)) stop("sp PACKAGE MISSING")
if (!require(rgdal)) stop("rgdal PACKAGE MISSING")
if (!require(foreign)) stop("foreign PACKAGE MISSING")
#### Get Arguments
Args = commandArgs()
inputFC = sub(".shp", "", Args[5], ignore.case=TRUE)
outputFC = sub(".shp", "", Args[6], ignore.case=TRUE)
dependentVar = Args[7]
independentVarString = Args[8]
usePenalty = as.integer(Args[9])
usePenalty <- usePenalty == 1
coefTable = sub(".csv", "", Args[10], ignore.case=TRUE)
diagTable = sub(".csv", "", Args[11], ignore.case=TRUE)
print(paste(commandArgs(), collapse=" "))
### Extract path and shapefile names
tmp <- unlist(strsplit(inputFC, "/"))[-length(unlist(strsplit(inputFC, "/")))]
path <- paste(tmp,collapse="/")
in.shp.name <- unlist(strsplit(inputFC, "/"))[length(unlist(strsplit(inputFC, "/")))]
out.shp.name <- unlist(strsplit(outputFC, "/"))[length(unlist(strsplit(outputFC, "/")))]
#### Get variable Names
independentVars = strsplit(independentVarString, ";")
independentVars = c(unlist(independentVars))
### Read shapefile using rgdal
shp <- readOGR(path, in.shp.name)
### Logistic regression model
print("Running model....")
lmodel <- LogitRegression(shp@data, y=dependentVar, x=independentVars,
penalty=usePenalty)
print( lmodel$model )
### Add columns to data and write shapefile
print("Writing Output....")
shp@data$Probs <- predict(lmodel$model,shp@data[,rownames(lmodel$model$var)[-1]],
type="fitted")
shp@data$Residual = lmodel$Residuals[,1]
shp@data$StdResid = lmodel$Residuals[,2]
writeOGR(shp, path, out.shp.name, driver="ESRI Shapefile", check_exists=TRUE,
overwrite_layer=TRUE)
### Write comma seperated ASCII flatfile of model coefficients and diagnositics
### model lrm MODEL OBJECT
### diagTable DATAFRAME OF REGRESSION DIAGNOSTICS
### coefTable DATAFRAME OF REGRESSION COEFFICENTS
### Residuals DATAFRAME OF RESIDUALS AND STANDARDIZED RESIDUALS
write.csv(lmodel$coefTable, coefTable)
write.csv(lmodel$diagTable, diagTable)
#if (AIF == TRUE) {
# p <- pentrace(lmodel$model, penalty=seq(0.1,1,by=0.05), which="aic.c")
# plot(p)
# }
print("Model Complete...")