R: Classification formula from trained GLM model [reproducible example provided]

Question

QUESTIONS

(1) what is the classification formula from the fit model in example code below named 'model1'? (is it formula A, B or Neither?)

(2) how does 'model1' determine if class == 1 vs. 2?

Formula A:

class(Species{1:2}) = (-31.938998) + (-7.501714 * [PetalLength]) + (63.670583 * [PetalWidth])
Formula B:

class(Species{1:2}) = 1.346075e-14 + (5.521371e-04 * [PetalLength]) + (4.485211e+27 * [PetalWidth])

USE CASE

Use R to fit/train a binary classification model, then interpret the model for the purpose of manual calculating classifications in Excel, not R.

MODEL COEFFICIENTS

>coef(model1)
#(Intercept) PetalLength  PetalWidth 
#-31.938998   -7.501714   63.670583 

>exp(coef(model1))
#(Intercept)  PetalLength   PetalWidth 
#1.346075e-14 5.521371e-04 4.485211e+27

R CODE EXAMPLE

# Load data (using iris dataset from Google Drive because uci.edu link wasn't working for me today)
#iris <- read.csv(url("http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data"), header = FALSE)
iris <- read.csv(url("https://docs.google.com/spreadsheets/d/1ovz31Y6PrV5OwpqFI_wvNHlMTf9IiPfVy1c3fiQJMcg/pub?gid=811038462&single=true&output=csv"), header = FALSE)
dataSet <- iris

#assign column names
names(dataSet) <- c("SepalLength", "SepalWidth", "PetalLength", "PetalWidth", "Species")

#col names
dsColNames <- as.character(names(dataSet))

#num of columns and rows
dsColCount <- as.integer(ncol(dataSet))
dsRowCount <- as.integer(nrow(dataSet))

#class ordinality and name
classColumn <- 5 
classColumnName <- dsColNames[classColumn]
y_col_pos <- classColumn

#features ordinality
x_col_start_pos <- 1
x_col_end_pos <- 4

# % of [dataset] reserved for training/test and validation  
set.seed(10)
sampleAmt <- 0.25
mainSplit <- sample(2, dsRowCount, replace=TRUE, prob=c(sampleAmt, 1-sampleAmt))

#split [dataSet] into two sets
dsTrainingTest <- dataSet[mainSplit==1, 1:5] 
dsValidation <- dataSet[mainSplit==2, 1:5]
nrow(dsTrainingTest);nrow(dsValidation);

# % of [dsTrainingTest] reserved for training
sampleAmt <- 0.5
secondarySplit <- sample(2, nrow(dsTrainingTest), replace=TRUE, prob=c(sampleAmt, 1-sampleAmt))

#split [dsTrainingTest] into two sets 
dsTraining <- dsTrainingTest[secondarySplit==1, 1:5]
dsTest <- dsTrainingTest[secondarySplit==2, 1:5]
nrow(dsTraining);nrow(dsTest);

nrow(dataSet) == nrow(dsTrainingTest)+nrow(dsValidation)
nrow(dsTrainingTest) == nrow(dsTraining)+nrow(dsTest)

library(randomGLM)

dataSetEnum <- dsTraining[,1:5]
dataSetEnum[,5] <- as.character(dataSetEnum[,5])
dataSetEnum[,5][dataSetEnum[,5]=="Iris-setosa"] <- 1 
dataSetEnum[,5][dataSetEnum[,5]=="Iris-versicolor"] <- 2 
dataSetEnum[,5][dataSetEnum[,5]=="Iris-virginica"] <- 2 
dataSetEnum[,5] <- as.integer(dataSetEnum[,5])

x <- as.matrix(dataSetEnum[,1:4])
y <- as.factor(dataSetEnum[,5:5])

# number of features
N <- ncol(x)

# define function misclassification.rate
if (exists("misclassification.rate") ) rm(misclassification.rate);
misclassification.rate<-function(tab){
  num1<-sum(diag(tab))
  denom1<-sum(tab)
  signif(1-num1/denom1,3)
}

#Fit randomGLM model - Ensemble predictor comprised of individual generalized linear model predictors
RGLM <- randomGLM(x, y, classify=TRUE, keepModels=TRUE,randomSeed=1002)

RGLM$thresholdClassProb

tab1 <- table(y, RGLM$predictedOOB)
tab1
# y  1  2
# 1  2  0
# 2  0 12

# accuracy
1-misclassification.rate(tab1)

# variable importance measure
varImp = RGLM$timesSelectedByForwardRegression
sum(varImp>=0)

table(varImp)

# select most important features
impF = colnames(x)[varImp>=5]
impF

# build single GLM model with most important features
model1 = glm(y~., data=as.data.frame(x[, impF]), family = binomial(link='logit'))

coef(model1)
#(Intercept) PetalLength  PetalWidth 
#-31.938998   -7.501714   63.670583 

exp(coef(model1))
#(Intercept)  PetalLength   PetalWidth 
#1.346075e-14 5.521371e-04 4.485211e+27 

confint.default(model1)
#                2.5 %   97.5 %
#(Intercept) -363922.5 363858.6
#PetalLength -360479.0 360464.0
#PetalWidth  -916432.0 916559.4

Ernest A · Accepted Answer

Your model is defined as

model1 <- glm(y~., data=as.data.frame(x[, impF]), family=binomial(link='logit'))

The family=binomial(link='logit')) bit is saying that the response y is a series of Bernoulli trials, i.e. a variable that takes values 1 or 0 depending on a parameter p, and that p = exp(m) / (1 + exp(m)), where m is a function of the data, called the linear predictor.

The formula y~. means that m = a + b PetalLength + c PetalWidth, where a, b, c are the model coefficients.

Therefore the probability of y = 1 is

> m <- model.matrix(model1) %*% coef(model1)
> exp(m) / (1+exp(m))
            [,1]
20  3.448852e-11
50  1.253983e-13
65  1.000000e+00
66  1.000000e+00
87  1.000000e+00
105 1.000000e+00
106 1.000000e+00
107 1.000000e+00
111 1.000000e+00
112 1.000000e+00
116 1.000000e+00
118 1.000000e+00
129 1.000000e+00
130 1.000000e+00

We can check that this is the same as the output of fitted.values

> fitted.values(model1)
          20           50           65           66           87          105 
3.448852e-11 1.253983e-13 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
         106          107          111          112          116          118 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
         129          130 
1.000000e+00 1.000000e+00

Finally, the response can be classified in two categories depending on whether P(Y = 1) is above or below a certain threshold. For example,

> ifelse(fitted.values(model1) > 0.5, 1, 0)
 20  50  65  66  87 105 106 107 111 112 116 118 129 130 
  0   0   1   1   1   1   1   1   1   1   1   1   1   1

R: Classification formula from trained GLM model [reproducible example provided]

Answers (2)

Related Questions