R: Classification Formel aus trainierten GLM-Modell [reproduzierbares Beispiel vorgesehen]

Question

FRAGT

(1), was die Einstufung Formel aus dem fit-Modell in Beispiel Code unter dem Namen 'model1'? (Ist es Formel A, B oder Weder?)

(2) Wie kann 'model1' bestimmen, ob Klasse == 1 vs. 2?

• Formel A:

  Klasse (Species {1: 2}) = (-31,938998) + (-7,501714 * [PetalLength]) + (63,670583 * [PetalWidth])

• Formel B:

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

USE CASE

Verwenden R ein binäres Klassifikationsmodell zu passen/Zug, dann das Modell zum Zwecke der manuellen Berechnung Einstufungen in Excel interpretieren, nicht R.

Modellkoeffizienten

>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-BEISPIEL

# 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

ist Ihr Modell als

definiert

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

Das family=binomial(link='logit')) Bit sagt, dass die Antwort y eine Reihe von Bernoulli-Versuchen ist, das heißteine Variable, die Werte 1 oder 0 in Abhängigkeit von einem Parameter p ist eine Funktion der Daten nimmt und daß p = exp (m)/(1 + exp (m)), wobei m , genannt der lineare Prädiktor.

Die Formel y~. bedeutet, dass m = a + b PetalLength + c PetalWidth, wo ein , b, c die Modellkoeffizienten sind.

Daher ist die Wahrscheinlichkeit des y = 1 ist

> 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 

wir überprüfen können, dass diese das gleiche wie die Ausgabe von 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 

Schließlich ist, kann die Reaktion in zwei klassifiziert werden, Kategorien abhängig davon, ob P (Y = 1) über oder unter einem bestimmten Schwellenwert liegt. Zum Beispiel

> 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 

Answer 2

Ein GLM-Modell hat eine Link-Funktion und einen linearen Prädiktor. Sie haben Ihre Link-Funktion oben nicht angegeben.

Sei Y = {0,1} und X sei eine n x p Matrix. (Unter Verwendung von pseudo-latex) Dies führt zu \hat Y= \phi(X \hat B) = \eta

wo
- \eta der linearen Prädiktor
- \phi() ist die Link-Funktion

Der lineare Prädiktor nur X %*% \hat B ist und die Klassifizierung zurück zu P(Y=1|X) = \phi^{-1}(\eta) - dh die inverse Verbindungsfunktion. Die inverse Link-Funktion hängt offensichtlich von der Wahl der Verbindung ab. Für eine logit haben Sie die inverse logit P(Y=1|X) = exp(eta)/(1+ exp(eta))