R：来自经过训练的 GLM 模型的分类公式 [提供了可重现的示例]

Question

问题

(1) 以下示例代码中名为“model1”的拟合模型的 class化公式是什么？（是公式 A、B 还是两者都不是？）

(2) 'model1' 如何判断 class == 1 vs. 2?

公式A：

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

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

用例

将 R 用于 fit/train 二元 class化模型，然后解释模型以手动计算 class化 Excel，而不是 R。

模型系数

>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 代码示例

# 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

Answer 1

GLM 模型具有 link 函数和线性预测器。您没有在上面指定 link 函数。

令 Y = {0,1} 且 X 为 n x p 矩阵。（使用伪 LaTeX）这导致 \hat Y= \phi(X \hat B) = \eta

其中
- \eta 是线性预测器
- \phi() 是 link 函数

线性预测器只是 X %*% \hat B，分类返回 P(Y=1|X) = \phi^{-1}(\eta)——即反 link 函数。反 link 函数显然取决于 link 的选择。对于 logit，您有逆 logit P(Y=1|X) = exp(eta) / (1+ exp(eta))

Answer 2

您的模型定义为

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

family=binomial(link='logit')) 位表示响应 y 是一系列伯努利试验，即根据参数取值 1 或 0 的变量 p，并且 p = exp(m) / (1 + exp(m ))，其中 m 是数据的函数，称为线性预测器。

公式y~.表示m = a + b PetalLength + c 花瓣宽度，其中 a、b、c是模型系数。

因此y = 1的概率是

> 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

我们可以检查这与 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

最后，根据 P(Y = 1) 是高于还是低于某个阈值，可以将响应分为两类。例如，

> 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：来自经过训练的 GLM 模型的分类公式 [提供了可重现的示例]

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

r

classification

logistic-regression