nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Logistic Regression in ArcGIS is this a possibility?

AnthoniaOnyeahialam

Is binary logistic regression in ArcGIS a possibility ? I know one can do linear regression but binary response regression is what am not sure of. Because I have this sort of data which is also mapped so am wondering what I can do with it in ArcGIS.

Find more posts tagged with

Accepted answers

All comments

LaurenRosenshein

This is a great question that we get a lot! At this time logistic regression is not available in ArcGIS, but we do have a sample script available that helps you run logistic regression using the R statistical package right from inside ArcMap. You can find the sample script here.

GeorgeChaka

Hi,

I was wondering if its possible to post some tutorial on how to use the Logit Regression sample script in ArcGIS 10. Is it realistic to expect Logistic regression to be made available in ArcGIS anytime soon?

Thank you for all the answers you and your team provide to the numerous questions.

George

PaulLeonard

we do have a sample script available that helps you run logistic regression using the R statistical package right from inside ArcMap. here.

Two of the required packages for R are no longer on CRAN and are outdated, 'clustTool' and 'Design'. Besides that, there have been some user issues with this script; read the comments. I haven't tested with 10.0 but can confirm issues with 10.1

Cheers

JeffreyEvans

The "clustTool" library is in fact not available for versions of R > 2.12 however, the "Design" library has been folded into the "rms" (Regression Modeling Strategies) library without change in function names. The clustTool library is only required for clustering models. The implementation of logistic regression uses the rms library so, if you only need the logistic implementation you should be fine. I think that I will rewrite the code and post a new version with some new bells-and-whistles. Hopefully ESRI will not have a problem with this.

AllisonWhite

The "clustTool" library is in fact not available for versions of R > 2.12 however, the "Design" library has been folded into the "rms" (Regression Modeling Strategies) library without change in function names. The clustTool library is only required for clustering models. The implementation of logistic regression uses the rms library so, if you only need the logistic implementation you should be fine. I think that I will rewrite the code and post a new version with some new bells-and-whistles. Hopefully ESRI will not have a problem with this.

Hi Dr. Evans-

I tried to install your Arc-script creating an R Toolbox inside ArcMap to perform logistic regression but I am having difficulty getting it to work. I am using ArcGIS 10.1 and R version 3.0.1., but I see in the installation instructions (and in the description) that you were using ArcGIS 10 and R 2.11.1.. Do you know how I can get it to work using the newest versions of ArcGIS and R? I can see the R Toolbox in the ArcToolbox but then I get error messages when I actually try to run the logit models. Maybe I need different R packages?

I called ESRI today and they said that their development team is working on creating logistic regression functionality in ArcGIS but wouldn't tell me when it would be available. Any suggestions you may have are greatly appreciated! Thank you.

JeffreyEvans

This is not my toolbox, I was just planning on updating the clustering functionality.

In taking a quick look at the R code ESRI has built in a very specific (2.14) version dependency. Since this condition is false then it is trying to add the depreciated library "Design". You can fix this by opening the file"LogitWithR.r", located in the Scripts directory, in a text editor.

Delete these lines

versionBool = checkRVersion(2, 14)
if (versionBool){
    library(rms)
}else{
    require(Design)
}

And replace with

require(rms)

AllisonWhite

This is not my toolbox, I was just planning on updating the clustering functionality.

In taking a quick look at the R code ESRI has built in a very specific (2.14) version dependency. Since this condition is false then it is trying to add the depreciated library "Design". You can fix this by opening the file"LogitWithR.r", located in the Scripts directory, in a text editor.

Delete these lines
versionBool = checkRVersion(2, 14)
if (versionBool){
    library(rms)
}else{
    require(Design)
}
And replace with
require(rms)

Hi Dr. Evans-

I looked at LogitWithR.r in WordPad and I don't see the code that you referenced above. This is what I see:

#### Load Libraries ####
print("Loading Libraries....")
library(maptools)
require(Design)
library(sm)

#### Get Arguments ####
Args = commandArgs()
print(Args)

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(".dbf", "", Args[10], ignore.case = TRUE)
diagTable = sub(".dbf", "", Args[11], ignore.case = TRUE)

#### Get Ind Var Names ####
independentVars = strsplit(independentVarString, ";")
independentVars = c(unlist(independentVars))

### Make Formula ####
form = as.formula(paste(dependentVar, paste(independentVars, collapse='+'), sep='~'))

print("Begin Calculations....")
### Using Maptools ####
shp = readShapeSpatial(inputFC)

### Do Logit ####
print("Logit....")
fit = lrm(form, shp, x = TRUE, y = TRUE)

print("Adjustment....")
### AIC Alternative Measure of Model Performance ####
bf = pentrace(fit, seq(.2,1,by=.05))
if (usePenalty) {
    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 (usePenalty){
    fit = update(fit, penalty = bf$penalty)
    }

### Residuals ####
res = residuals.lrm(fit)
resOut = c(res)
resSTD = (resOut - mean(resOut)) / sqrt(var(resOut))

### Create Output Shape File ####
print("Writing Output....")
shp$Residual = resOut
shp$StdResid = resSTD
writeSpatialShape(shp, outputFC)

### Write Coefficient DBF Table ####
allIndVars = c("Intercept")
allIndVars = append(allIndVars, independentVars)
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)
write.dbf(coefFrame, coefTable)

### Write Diagnostic DBF Table ####
diagNames = c("Obs", "Max Deriv", "Model L.R.", "d.f.", "P", "C",
              "Dxy", "Gamma", "Tau-a", "R2", "Brier",
              "Penalty", "AIC")
allStats = c(as.vector(fit$stats), pen, aic)
diagValues = matrix(allStats, length(diagNames), 1)
diagList = list("Diag_Name" = diagNames, "Diag_Value" = diagValues)
diagFrame = data.frame(diagList)
write.dbf(diagFrame, diagTable)

print("Calculations Complete...")

JeffreyEvans

You are using an older version of the toolbox. Download the "more current" 10-10.1 compliant version here:
http://www.arcgis.com/home/item.html?id=a5736544d97a4544aa47d06baf910f6d

AllisonWhite

You are using an older version of the toolbox. Download the "more current" 10-10.1 compliant version here:
http://www.arcgis.com/home/item.html?id=a5736544d97a4544aa47d06baf910f6d

Hi (again) Dr. Evans-

I installed the most recent toolbox and am seeing the following error messages when I try to run a logit: (Do you know what this is about?). Thank you very much for all of your help thus far.

Executing: LogitRegressionR "F:\Russia Shapefiles\Export_Output.shp" C:\Users\Allison\Documents\ArcGIS\Export_Output_LogitRegressio.shp 2007Tata_8 2007Tata_3;2007Tata_4;2007Tata_5 USE_PENALTY C:\Users\Allison\Documents\ArcGIS\Export_Output_LogitRegressio.dbf C:\Users\Allison\Documents\ArcGIS\Export_Output_LogitRegressio.dbf
Start Time: Mon Jun 10 18:36:25 2013
Running script LogitRegressionR...
Loading required package: foreign
Loading required package: sp
Loading required package: grid
Loading required package: lattice
Checking rgeos availability: FALSE
Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
which has a restricted licence. It is disabled by default;
to enable gpclib, type gpclibPermit()
Package `sm', version 2.2-5: type help(sm) for summary information
Loading required package: rms
Loading required package: Hmisc
Loading required package: survival
Loading required package: splines
Hmisc library by Frank E Harrell Jr

Type library(help='Hmisc'), ?Overview, or ?Hmisc.Overview')
to see overall documentation.

NOTE:Hmisc no longer redefines [.factor to drop unused levels when
subsetting. To get the old behavior of Hmisc type dropUnusedLevels().

Attaching package: 'Hmisc'

The following object is masked from 'package:survival':

    untangle.specials

The following object is masked from 'package:maptools':

    label

The following object is masked from 'package:base':

    format.pval, round.POSIXt, trunc.POSIXt, units

Attaching package: 'rms'

The following object is masked from 'package:survival':

    Surv

Error in parse(text = x) : <text>:1:5: unexpected symbol
1: 2007Tata_8
       ^
Calls: as.formula ... formula -> formula.character -> formula -> eval -> parse
Execution halted

Completed script LogitRegressionR...
Failed to execute (LogitRegressionR).
Failed at Mon Jun 10 18:36:27 2013 (Elapsed Time: 2.00 seconds)

JeffreyEvans

I would recommend cutting out the middle-man and running this directly in R. I rewrote the R code so you can run it directly. Set arguments in the "USER ARGUMENTS" section and the rest of the code is automated. The OutShp will be written to the working directory. This will require the sp, rgdal and rms libraries.

####################################################
# USER ARGUMENTS 
####################################################
# SET WORKING DIRECTORY
setwd("D:/TMP")

# SET VARIABLE NAMES
dependentVar = "SppDen"
independentVars = c("Var1","Var2","Var3")

usePenalty=TRUE         # USE PENALTY
InShp="MyShape"         # IN SHAPEFILE (NO .shp NEEDED)
OutShp="LogitModel"     # OUT SHAPEFILE (NO .shp NEEDED) 
diagTable = "Diag.csv"  # DIAGNOSTIC TABLE
coefTable = "Coef.csv"  # COEFFICENT TABLE

####################################################
# START MODEL
####################################################
require(sp)
require(rgdal)
require(rms)

# READ SHAPEFILE 
shp <- readOGR(getwd(), InShp)

# CHECK VARIABLE NAMES 
if ( is.na(match(dependentVar,names(shp@data))) )
  stop("Dependent Variable not present in data")  

xNames <- intersect(independentVars,names(shp@data)) 
  if (length(xNames) < length(independentVars)) 
    stop("Mismatch in Independent Variable Names")

# CREATE FORMULA
form=as.formula(paste(dependentVar, paste(independentVars, collapse='+'), sep='~'))

# LOGISTIC REGRESSION WITH AIC
fit <- lrm(form, data=shp@data, x = TRUE, y = TRUE) 
  bf <- pentrace(fit, seq(.2,1,by=.05))
  if (usePenalty) {
      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 (usePenalty){ fit = update(fit, penalty=bf$penalty) }

# RESIDUAL ERROR
res <- residuals.lrm(fit)
  resSTD <- (res - mean(res)) / sqrt(var(res))

# ADD RESIDUALS, STANDARDIZED RESIDUALS AND PROBABILITIES 
shp@data <- data.frame(shp@data,Residual=res, StdResid=resSTD, 
                Probs=predict(fit,shp@data[,rownames(fit$var)[-1]],
    type="fitted") )  

# WRITE SHAPEFILE    
writeOGR(shp, dsn=getwd(), OutShp, driver="ESRI Shapefile", 
         check_exists=TRUE, overwrite_layer=TRUE)
   
# WRITE COEFFICENT AND DIAGNOSTIC TABLES 
allIndVars = c("Intercept")
allIndVars = append(allIndVars, independentVars)
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)
  write.csv(coefFrame, coefTable, row.names=FALSE)  

daigFrame <- data.frame(Names=c(names(fit$stats),"PEN","AIC"),
                      Value=c(as.vector(fit$stats), pen, aic))
  write.csv(daigFrame, diagTable, row.names=FALSE)

paolapasino

https://github.com/Esri/R-toolbox-py

Try this updated tool box it works ( not fully) for me

Paola

JeffreyEvans

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...")

MarilynMontgomery

Hi Dr. Evans,

thank you for posting the script for R for logistic regression. I started using R only a few weeks ago and I got it to work without the "middle man" (arcgis).

However I wanted to make sure, this is not an autologistic regression, correct? I want the spatial covariate for a autologistic/spatial logistic model; and it doesn't look like this script takes the geography of the observations into account. Am I right?

Thanks,

Marilyn

JeffreyEvans

Marilyn,
You are correct this is not an autologistic model. However, all you have to do to specify an autologistic regression using this code is to spatially lag your dependent variable and include it on the x side of the model as an independent variable (autocovariate). You can easily spatially lag a variable using the spdep package.

The caveat with this type of model is that there is an assumption that the autocovariate does not correlate with the other independent variables. If this is the case then the model becomes quite biased and the coefficients often misleading. It is highly recommend that you explore mixed effect models and treat the autocorrelation as a random effect. This type of model can be specified using the GLMM function in the nlme package. If you would like to correctly specify an autologistic model you could put it in an MCMC context and specify the autocorrelation function as a design matrix.

Just in case you have not read them I attached the Dormann (2007) paper on autologistic regression and well as the He et al., (2003) simulation paper.

Best,
Jeff

He2003.pdf

Dormann2007.pdf

MarilynMontgomery

Hi Jeff,

Thank you very much for your helpful reply. Does the lag.listw command in R do the same thing as autocov_dist, if I use a distance-based weight matrix with lag.listw and specify the same distances for lag and autocovariate? It seems like I could use either of these commands to calculate an autocovariate for autologistic regression.

And could you give me any more guidance on how to "correctly specify an autologistic model ... put it in an MCMC context and specify the autocorrelation function as a design matrix"? When you say autocorrelation function as a design matrix, are you referring to the autocovariate?

Thanks again!

Marilyn

JeffreyEvans

The "autocov_dist" and "lag.listw" functions do not result in the same thing. For a autologistic model you want to use the "autocov_dist". The "lag.listw" results in a sparse matrix based on the spatial lag of the defined neighbor relationship. In concert with knearneigh (where k = (n-1)) you could trick this function to return lagged weights but the "autocov_dist" function is specifically designed to produce a autocovariate.

I have to say that an MCMC specification of a model is a very complicated procedure and requires either writing your own MCMC function, which is a huge undertaking if you are not familiar with the mathematical/statistical theory, or using Jags/WinBUGS through an R interface.

The design matrix is defining a specific autocovariance function in the parameter portion of the Bayesian model rather than depending on an autocovariate, which can bias the model. This type of model is really beyond the scope of this forum, as is much of the R advice thus far. If you have further questions I would recommend contacting me off list. jeffrey_evans<at>tnc.org

Join the Conversation

Signed in members can post, follow updates, and more. New here? Register a free account.

Quick Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.