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

Posted 2023-03-13

【中文标题】R：来自经过训练的 GLM 模型的分类公式 [提供可重现的示例]

【英文标题】：R: Classification formula from trained GLM model [reproducible example provided]

【发布时间】：2016-10-02 13:07:47

【问题描述】：

问题

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

(2) 'model1' 如何确定 class== 1 vs. 2？

公式 A：

类(物种1:2) = (-31.938998) + (-7.501714 * [花瓣长度]) + (63.670583 * [花瓣宽度])

公式 B：

类(物种1:2) = 1.346075e-14 + (5.521371e-04 * [PetalLength]) + (4.485211e+27 * [PetalWidth])

用例

使用 R 拟合/训练二元分类模型，然后解释模型以在 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

【参考方案1】：

您的模型定义为

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

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

公式y~.表示m = a + b PetalLength + c PetalWidth，其中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 

【参考方案2】：

GLM 模型具有链接函数和线性预测器。你没有在上面指定你的链接函数。

让 Y = 0,1 和 X 是一个 n x p 矩阵。 （使用pseudo-LaTeX）这导致\hat Y= \phi(X \hat B) = \eta

在哪里 - \eta 是线性预测器 - \phi()是链接函数

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

【讨论】：

感谢@Alex 的评论，但是OP 的示例代码中提供了链接功能。见第 95 行model1 = glm(y~., data=as.data.frame(x[, impF]), family = binomial(link='logit'))

@Webby 很公平......我略读了你的问题。我没看到。我的建议是您查看 GLM 模型，因为您应该了解您使用的模型。

我同意，但这并不能回答我的任何一个问题。

我不明白你的回答。是公式A、B还是C（给定OP中提供的模型系数）？