4 Kategoriale Regression
Mit der Klasse der generalisierten linearen Regressionmodelle steht uns bereits eine große Bandbreite an Modellen für verschiedene Verteilungsformen der Zielvariable zur Verfügung. Was uns bis jetzt noch nicht möglich ist, ist die Modellierung einer Zielvariablen mit kategorialen Ausprägungen. Hierbei unterscheiden wir zwischen Kategorien ohne sinnvolle Interpretation der Reihenfolge (nominal) und mit sinnvoller Interpretation der Reihenfolge, jedoch ohne sinnvoller Interpretation der Abstände (ordinal). In beiden Fällen hat die Zielvariable die Form:
\[ Y_i = k, ~~~k = 1, ..., K~~~k, K \in \mathbb{N} \] wobei \(k_i\) für die Kategorie der \(i-\)ten Beobachtung steht und \(c\) verschiedene Kategorien einnehmen kann.
4.1 Nominal skalierte Zielvariable
Mit der nominal skalierten Zielvariable beginnend, können wir feststellen, dass wir bereits einen Spezialfall dieses Modells kennen gelernt haben, das logistische Regressionsmodell. Hierbei handelt es sich um eine Zielvariable mit zwei Kategorien, bei dem wir die Wahrscheinlichkeit modelliert haben, dass eine Kategorie der Zielvariablen auftritt. Hierbei steht \(\pi_i\) für die (bedingte) Wahrscheinlichkeit, dass \(Y_i\) den Wert \(1\) annimmt und entsprechend liegt die Wahrscheinlichkeit, dass wir \(Y_i = 0\) beobachten bei \(1-\pi_i\). Der lineare Prediktor des logistischen Regressionsmodells entspricht dem Quotienten der logarithmierten Wahrscheinlichkeiten der Kategorie \(1\) und \(0\):
\[ \ln \frac{\pi_i}{1-\pi_i} = \beta_0 + \beta_1 x_{i1} + ... + \beta_p x_{ip} = \boldsymbol{\beta}^T \boldsymbol{x_i} \]
Bei dieser Betrachtung vergleichen wir somit die Wahrscheinlichkeit, dass \(Y_i\) der Kategorie \(1\) entspricht mit der Wahrscheinlichkeit, dass \(Y_i\) der Kategorie \(0\) entspricht. Wir setzen somit Kategorie \(1\) in Referenz zu Kategorie \(0\), weshalb wir letztere auch als Referenzkategorie bezeichnen.
4.1.1 Modellierung der Wahrscheinlichkeit je Kategorie
Haben wir es nun mit einer Variablen mit mehreren Kategorien zu tun, können wir den grundsätzlichen Gedanken des Modells mit zwei Kategorien analog fortsetzen. Angenommen \(Y_i\) hat drei Kategorien, z.B.,
\[ 1: \text{Pop} \] \[ 2: \text{Rock} \]
\[ 3: \text{Klassik} \]
Als erklärende Variablen \(\boldsymbol{x_i}\) stellen wir uns Variablen wie die Existenz verzerrter Instrumente, die Anzahl der gesungenen Wörter und ähnliches vor. Wir möchten wie bei der logistischen Regression nun bedingte Wahrscheinlichkeiten schätzen, dass es sich bei einem Song \(Y_i\) und einen Pop, Rock oder Klassik Song handelt. Als erste Idee könnten wir zunächst versuchen für jede Kategorie \(k = \lbrace1, 2, 3\rbrace\) und jeden Song \(i\) eine Wahrscheinlichkeit zu schätzen, also \(\pi_{ik}\). Da es sich um disjunkte Kategorien und Wahrscheinlichkeiten handelt, muss \(\pi_{ik} >= 0\) und für alle \(i\):
\[ \sum_{k = 1}^{K} \pi_{ik} = 1 \] gelten. Durch kurze Überlegung stellen wir jedoch fest, dass es aufgrund dieser Anforderungen nicht notwendig ist jede einzelne Wahrscheinlichkeit zu modellieren, da wir eine Kategorie auch als Gegenwahrscheinlichkeit zur Summe der Wahrscheinlichkeiten der übrigen Kategorien bestimmen können. Sei z.B. \(\pi_{i1} = 0.20\) und \(\pi_{ik} = 0.50\), so muss
\[\pi_{i3} = 1 - \pi_{i1} - \pi_{i2} = 1 - 0.20 - 0.50 = 0.30\]
gelten. Gegeben wir entscheiden uns daher Kategorie \(3\) als Referenzkategorie zu verwenden, so können wir erneut einen Vergleich der (logarithmierten) Kategoriewahrscheinlichkeiten betrachten. In unserem Beispiel:
\[ \ln \frac{\pi_{i1}}{\pi_{i3}} = \beta_{0, 1} + \beta_{1, 1} x_{i1} + ... + \beta_{p, 1} x_{ip} = \boldsymbol{\beta}_1^T \boldsymbol{x_i} \]
\[ \ln \frac{\pi_{i2}}{\pi_{i3}} = \beta_{0, 2} + \beta_{1, 2} x_{i1} + ... + \beta_{p, 2} x_{ip} = \boldsymbol{\beta}_2^T \boldsymbol{x_i} \]
Was dabei (hoffentlich) ins Auge springt, ist dass jeder Vergleich bzw. jede Gleichung eigene Betakoeffizienten \(\boldsymbol{\beta}_k, ~k = 1, 2\) hat. Dies ist notwendig, da ansonsten die Wahrscheinlichkeiten für die ersten beiden Kategorien identisch sein müssten. Allgemein ist es also notwendig eine Referenzkategorie \(K\) mit Wahrscheinlichkeit \(\pi_{iK} = 1 - \sum_{k = 1}^{K-1} \pi_{ik}\) zu bestimmen und dann je Kategorie \(r = 1, ..., K-1\) die zugehörige Wahrscheinlichkeit \(\pi_{ik}\) mittels eines linearen Prediktors je Kategorie und einer geeigneten Transformation zu modellieren. Wie könnte die geeignete Transformation hier aussehen. Als erste Bedingung haben wir \(\pi_{ik} >= 0\), weshalb es sich anbietet eine Transformation auf den nicht-negativen bzw. positiven Wertebereich mittels der \(e-\)Funktion vorzunehmen, also:
\[ e^{\boldsymbol{\beta}_k^T \boldsymbol{x_i}} \]
Zusätzlich muss die Summe aller Einzelwahrscheinlichkeiten auf \(1\) normiert werden, d.h. wir dividieren alle transformierten linearen Prediktoren durch einen Wert, so dass die Summe der Einzelwahrscheinlichkeiten wieder \(1\) ergibt. Dies gelingt durch die Normierung:
\[ 1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}} \]
und zwar weil:
\[ \begin{split} \pi_{iK} & = 1 - \sum_{k = 1}^{K-1} \pi_{ik} = \frac{1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}}}{1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}}} - \frac{e^{\boldsymbol{\beta}_1^T \boldsymbol{x_i}}}{1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}}} - ... - \frac{e^{\boldsymbol{\beta}_{K-1}^T \boldsymbol{x_i}}}{1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}}} = \\ & = \frac{1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}} - \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}}}{1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}}} = \frac{1}{1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}}} \end{split} \] und
\[ \pi_{ik} = \frac{e^{\boldsymbol{\beta}_k^T \boldsymbol{x_i}}}{1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}}}\]
womit
\[ \pi_{iK} + \sum_{k = 1}^{K-1} \pi_{ik} = \frac{1}{1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}}} + \frac{e^{\boldsymbol{\beta}_1^T \boldsymbol{x_i}}}{1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}}} + ... + \frac{e^{\boldsymbol{\beta}_{K-1}^T \boldsymbol{x_i}}}{1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}}} = \frac{1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}}}{1 + \sum_{s = 1}^{K-1} e^{\boldsymbol{\beta}_s^T \boldsymbol{x_i}}} = 1\]
4.1.2 Schätzen des Modells
Wir wissen nun also, wie unser Modell aufgebaut ist. Zum Schätzen des Modells greifen wir auf die Maximum-Likelihood-Methode zurück. Die abhängige Variable ist kategorial verteilt, dies bedeutet, je Datenpunkt entspricht die Wahrscheinlichkeit:
\[ f(\boldsymbol{y}_i | \boldsymbol{\pi}) = \pi_{i1}^{y_{i1}} \cdot \pi_{i2}^{y_{i2}} \cdot ... \cdot \pi_{iK-1}^{y_{iK-1}} \cdot \underbrace{\left( 1 - \pi_{i1}^{y_{i1}} \cdot \pi_{i2}^{y_{i2}} \cdot ... \cdot \pi_{iK-1}^{y_{iK-1}} \right)}_{\pi_{iK}} {}^{1 - y_{i1} - ... - y_{iK-1}} \] d.h. der Wert der Likelihoodfunktion eines Datenpunktes \(y_i\) entspricht der Wahrscheinlichkeit der Kategorie des Datenpunktes. Wir können dies auch durch:
\[ f ( y_i | \boldsymbol { \pi }_i ) = \pi _ { i1 } ^ { y _ { i1 } } \cdot \ldots \cdot \pi _ { K } ^ { y _ { iK } } = \prod_{j = 1}^{K} \pi_{ij}^{e_{ij}} \] mit
\[ e_{ij} = \begin{cases} 1 & y_i = j \\ 0 & \text{sonst} \end{cases} \]
darstellen, womit die Log-Likelihood-Funktion für \(i = 1, ..., n\) Datenpunkte sich folgendermaßen darstellt:
\[ \ell \left( \boldsymbol{\beta}_1, ..., \boldsymbol{\beta}_{K-1} \right) = \sum_{i =1}^{n} \sum_{j = 1}^{K} e_{ij} \cdot \ln \left( \pi_{ij} \right), ~~~ i = 1, ..., n, ~~~j = 1, ..., K \]
Der Likelihood-Schätzer ist:
\[ \left( \hat{\boldsymbol{\beta}}_1, ..., \hat{\boldsymbol{\beta}}_{K-1} \right) = \underset{\boldsymbol{\beta}_1, ..., \boldsymbol{\beta}_{K-1}}{\mathrm{argmax}} \ln \mathcal{L} \left( \boldsymbol{\beta}_1, ..., \boldsymbol{\beta}_{K-1} \right) \]
4.1.3 Simulation eines kategorialen Modells
Dem Motto was man simulieren kann, hat man verstanden zunächst ein simuliertes Beispiel mit drei Kategorien und einer erkärenden Variable.
#Beta-Koeffizienten der ersten Kategorie
beta1 <- matrix(c(0.2, 1), nrow = 1)
#Beta-Koeffizienten der zweiten Kategorie
beta2 <- matrix(c(-0.3, -0.7), nrow = 1)
#Sitchprobenumfang
n <- 1000
#erklärende Variable
set.seed(42)
x <- rnorm(n)
#Daten mit Konstante
X <- cbind(1, x)
#lineare Prediktoren für Kategorie 1 und 2
lin_pred1 <- X %*% t(beta1)
lin_pred2 <- X %*% t(beta2)
#Umwandlung der linearen Prediktoren in Wahrscheinlichkeiten mittels eine Funktion
prob_transform <- function(lp1, lp2){
denominator <- 1 + exp(lp1) + exp(lp2)
pi1 <- exp(lp1) / denominator
pi2 <- exp(lp2) / denominator
pi3 <- 1 - pi1 - pi2
return(c(pi1, pi2, pi3))
}
#Matrix mit den Wahrscheinlichkeiten für Kategorie 1, 2 und 3 je Zeile
prob <- matrix(nrow = nrow(X), ncol = 3)
#Y Variable in Kodierung (1, 0, 0) für Kategorie 1, (0, 1, 0) für Kategorie 2 und (0, 0, 1) für Kategorie 3
y_dummy <- matrix(nrow = nrow(X), ncol = 3)
#Y Variable mit Kategorie 1, 2 oder 3
y <- c()
for(i in 1:nrow(X)){
prob[i, ] <- prob_transform(lin_pred1[i,1], lin_pred2[i, 1])
y_dummy[i, ] <- rmultinom(1, size = 1, prob = prob[i,])
y <- c(y, match(1, y_dummy[i,]))
}
head(y)## [1] 3 3 2 2 1 3#die kategoriale Regression wird mit dem nnet Package umgesetzt
library(nnet)
#vereinfachter Datensatz mit erklärender und abhängiger Variable
dataset <- data.frame(x = X[,2], y = as.factor(y))
#für eine leichtere Handhabung wird Kategorie 3 als Referenzkategorie festgelegt
dataset$y <- relevel(dataset$y, ref = 3)
#Die Syntax ist wie bei den anderen Regressionsmodellen
model <- multinom(y ~ x, data = dataset)## # weights: 9 (4 variable)
## initial value 1098.612289
## iter 10 value 953.158936
## final value 953.158674
## converged## Call:
## multinom(formula = y ~ x, data = dataset)
##
## Coefficients:
## (Intercept) x
## 1 0.31150699 0.8820988
## 2 -0.04055658 -0.5429803
##
## Std. Errors:
## (Intercept) x
## 1 0.08443403 0.09933589
## 2 0.09297055 0.09854583
##
## Residual Deviance: 1906.317
## AIC: 1914.3174.1.4 Reale Beispiele
Nach der Simulation betrachten wir ein paar reale Beispiel. Zunächst den bekannten Iris Datensatz, bei dem anhand der Blütenblätter und der Blüten Pflanzen in drei verschiedene Pflanzenkategorien eingeteilt werden. Zunächste schätzen wir das Modell für gesamten Datensatz.
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa#Die Syntax ist wie bei den anderen Regressionsmodellen
model <- multinom(Species ~ Sepal.Length + Sepal.Width + Petal.Length, data = iris)## # weights: 15 (8 variable)
## initial value 164.791843
## iter 10 value 15.825451
## iter 20 value 12.188082
## iter 30 value 11.888231
## final value 11.885927
## converged## Call:
## multinom(formula = Species ~ Sepal.Length + Sepal.Width + Petal.Length,
## data = iris)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## setosa -0.1013477 2.295778 9.3685323 -15.40205
## virginica -38.2328882 -3.853459 -0.6372467 13.15161
##
## Std. Errors:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## setosa 39.93598 76.648841 106.51757 274.384874
## virginica 14.06089 1.709709 2.29092 3.908092
##
## Residual Deviance: 23.77185
## AIC: 39.77185Mit der hier gewählten Referenzkategorie sind die Einflüsse der erklärenden Variablen auf die anderen Kategorien konträr. Gegeben alle übrigen Variablen, erhöht sich die Wahrscheinlichkeit für Kategorie “setose” bei größerer Petal.Length im Vergleich zur Referenzkategorie (“versicolor”), während diese sich im Vergleich zur Referenzkategorie für Kategorie “virginica” reduziert. Dieses Paket gibt keine p-Values aus, jedoch können wir diese leicht selbst bestimmen, da die Standardfehler mit in der Ausgabe des Modells enthalten sind. Die standardisierten realisierten Teststatistiken für die Nullhypothesen \(H_0: \beta_j = 0\) erhalten wir durch:
Da die Schätzer asymptotisch normal verteilt sind, können die p-Values approximativ durch die Normalverteilung bestimmt werden.
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## setosa 0.997975168 0.97610542 0.9299141 0.9552358467
## virginica 0.006546084 0.02420444 0.7808881 0.0007648122Bezieht man die Analyse ausschließlich darauf, ob die \(\beta-\)Koeffizienten statistisch signifikant von Null verschieden sind, würde man hier bei einem Signifikanzniveau von \(\alpha = 0.05\) zu dem Ergebnis kommen, dass lediglich die Sepal.Length und Petal.Length einen von Null verschiedenen Einfluss auf die Differenzierung zwischen der Referenzkategorie und der Kategorie “virginica” hat. Wie wir gleich sehen werden bedeutet dies jedoch nicht, dass es sich insgesamt um ein schlechtes Modell zur Unterscheidung der verschiedenen Pflanzenkategorien handelt. Hierzu unterteilen wir den Datensatz in zufällig ausgewählte Trainings- und Testdaten. Wir schätzen das Modell auf Basis der Trainingsdaten und vergleichen die Fähigkeit des Modells für beide Datensätze die Pflanzenkategorien anhand der erlärenden Variablen zu unterscheiden. Wir betrachten hierzu wieder eine Confusion Matrix, in der wir sehen können, wo das Modell richtig und wie das Modell falsch zuordnet.
#Paket zur Bestimmung der Confusion Matrix
library(caret)
#Könnten wir
idx <- 1:dim(iris)[1]
#teile die Daten zufällig auf, 100 Trainingsdaten, 50 Testdaten
set.seed(43)
train_idx <- sample(idx, 100)
test_idx <- idx[-train_idx]
train_data <- iris[train_idx, ]
test_data <- iris[test_idx, ]
#Regresssionmodell auf Basis der Trainingsdaten
train_model <- multinom(Species ~ Sepal.Length + Sepal.Width + Petal.Length, data = train_data)## # weights: 15 (8 variable)
## initial value 109.861229
## iter 10 value 13.226434
## iter 20 value 9.249979
## iter 30 value 9.143287
## iter 40 value 9.109571
## iter 50 value 9.108833
## iter 60 value 9.108654
## iter 70 value 9.108608
## iter 80 value 9.108588
## iter 90 value 9.108542
## iter 100 value 9.108513
## final value 9.108513
## stopped after 100 iterations#In sample Prognose
train_predict <- predict(train_model, train_data)
#Out of sample Prognose
test_predict <- predict(train_model, test_data)
#Confusion Matrix in sample
confusionMatrix(train_predict, train_data$Species)## Confusion Matrix and Statistics
##
## Reference
## Prediction versicolor setosa virginica
## versicolor 34 0 2
## setosa 0 32 0
## virginica 1 0 31
##
## Overall Statistics
##
## Accuracy : 0.97
## 95% CI : (0.9148, 0.9938)
## No Information Rate : 0.35
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.955
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: versicolor Class: setosa Class: virginica
## Sensitivity 0.9714 1.00 0.9394
## Specificity 0.9692 1.00 0.9851
## Pos Pred Value 0.9444 1.00 0.9688
## Neg Pred Value 0.9844 1.00 0.9706
## Prevalence 0.3500 0.32 0.3300
## Detection Rate 0.3400 0.32 0.3100
## Detection Prevalence 0.3600 0.32 0.3200
## Balanced Accuracy 0.9703 1.00 0.9622## Confusion Matrix and Statistics
##
## Reference
## Prediction versicolor setosa virginica
## versicolor 14 0 1
## setosa 0 18 0
## virginica 1 0 16
##
## Overall Statistics
##
## Accuracy : 0.96
## 95% CI : (0.8629, 0.9951)
## No Information Rate : 0.36
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9398
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: versicolor Class: setosa Class: virginica
## Sensitivity 0.9333 1.00 0.9412
## Specificity 0.9714 1.00 0.9697
## Pos Pred Value 0.9333 1.00 0.9412
## Neg Pred Value 0.9714 1.00 0.9697
## Prevalence 0.3000 0.36 0.3400
## Detection Rate 0.2800 0.36 0.3200
## Detection Prevalence 0.3000 0.36 0.3400
## Balanced Accuracy 0.9524 1.00 0.9554Wir können hier beobachten, dass es sowohl in wie auch out of sample zu einer sehr hohen Genauigkeit kommt und kaum Fehler bei der Einteilung der Kategorien gemacht werden.
Zuletzt noch ein kleines Spaßprojekt. Folgenden Datensatz habe ich mit Hilfe der Spotify API erzeugt. Spotify speichert zu allen Songs Attribute wie beispielsweise die Loudness oder Speechiness bzw. Dinge wie die Tonart. Bei der Erzeugung der Daten habe ich mir Songs aus Playlists ausgeben lassen, die einer von mehreren von mir vorgegebenen Kategorien entsprechen. Zu jedem Song habe ich die Kategorie und die Atrribute der Songs im Datensatz abgespeichert. Wir gehen wie zuletzt vor und betrachten, wie gut Kategorien mit dem Modell prognostiziert werden können.
spotify <- read.csv('Data/spotify_playlists.csv')
spotify$genre <- relevel(spotify$genre, ref = 'Pop')
#Könnten wir
idx <- 1:dim(spotify)[1]
#teile die Daten zufällig auf, 75% der Daten als Trainingsdaten
set.seed(43)
train_idx <- sample(idx, floor(0.75 * dim(spotify)[1]))
test_idx <- idx[-train_idx]
train_data <- spotify[train_idx, ]
test_data <- spotify[test_idx, ]
#Regresssionmodell auf Basis der Trainingsdaten
train_model <- multinom(genre ~ popularity + danceability + energy + key + loudness + mode + speechiness + acousticness + instrumentalness + liveness + valence + tempo + duration + time_signature, data = train_data)## # weights: 112 (90 variable)
## initial value 7053.924290
## iter 10 value 6705.699623
## iter 20 value 6463.071346
## iter 30 value 5719.993192
## iter 40 value 4886.182988
## iter 50 value 4516.660998
## iter 60 value 4420.100829
## iter 70 value 4408.668583
## iter 80 value 4407.879703
## iter 90 value 4407.827663
## final value 4407.826898
## converged## Call:
## multinom(formula = genre ~ popularity + danceability + energy +
## key + loudness + mode + speechiness + acousticness + instrumentalness +
## liveness + valence + tempo + duration + time_signature, data = train_data)
##
## Coefficients:
## (Intercept) popularity danceability energy key
## Alternative 0.6806983 -0.04699888 -6.514242 2.1075697 -0.008168461
## Electro -4.7325959 -0.05626319 6.959047 5.2686530 -0.015564976
## Hip-Hop -7.4434533 -0.02210377 7.893357 -0.4436977 0.031406255
## Metal -1.3113586 -0.05763150 -8.497366 10.1581048 -0.031272681
## Rap -9.8279961 -0.04677611 9.498826 1.1396654 -0.003812815
## Rock 1.2222894 -0.03402294 -7.131107 1.5785348 -0.004937185
## loudness mode speechiness acousticness instrumentalness
## Alternative -0.328777031 0.56778496 -4.219367 -2.2371862 5.152105
## Electro -0.006561962 0.02774649 -2.354864 -2.1421418 7.313404
## Hip-Hop -0.136311612 0.05429582 11.398454 -1.2471146 6.947054
## Metal -0.361847775 -0.10548310 -3.952947 -4.4632504 4.510148
## Rap -0.163915795 -0.04236722 11.390545 0.5124731 2.663859
## Rock -0.281081895 0.51746655 -7.581207 -1.2765763 2.006243
## liveness valence tempo duration time_signature
## Alternative 0.7096854 1.678355 -0.005486187 5.919771e-06 0.2736189
## Electro 0.8244114 -2.611466 0.002707091 -2.002072e-06 0.1786506
## Hip-Hop 1.4528968 -1.929671 0.001104803 1.651699e-06 0.3204392
## Metal 0.6092870 -2.383509 -0.010624768 1.354584e-05 -0.3135209
## Rap 0.7859081 -1.759950 0.006069261 7.063760e-06 0.1901984
## Rock 0.3698131 1.258438 -0.007452544 1.051553e-05 0.1485007
##
## Std. Errors:
## (Intercept) popularity danceability energy key
## Alternative 3.605697e-06 2.234834e-04 1.670070e-06 3.005845e-06 1.746761e-05
## Electro 4.135225e-07 3.562046e-05 3.246526e-07 2.283845e-07 2.519638e-06
## Hip-Hop 3.605613e-06 2.975914e-04 2.505460e-06 2.068595e-06 1.869486e-05
## Metal 3.608574e-06 2.101942e-04 1.420990e-06 3.513874e-06 1.730734e-05
## Rap 2.843289e-06 2.370238e-04 1.971767e-06 1.657647e-06 1.332859e-05
## Rock 3.660434e-06 2.486486e-04 1.794923e-06 2.954323e-06 1.758753e-05
## loudness mode speechiness acousticness
## Alternative 1.844856e-05 2.685899e-06 3.655764e-07 4.681903e-07
## Electro 2.516155e-06 2.829583e-07 3.876081e-08 1.479618e-07
## Hip-Hop 2.349910e-05 1.995799e-06 8.315430e-07 7.580111e-07
## Metal 1.197346e-05 2.158409e-06 4.159913e-07 8.730307e-09
## Rap 1.833332e-05 1.621623e-06 7.211951e-07 5.476871e-07
## Rock 1.834965e-05 2.731633e-06 3.021527e-07 5.277203e-07
## instrumentalness liveness valence tempo
## Alternative 1.361751e-07 6.698191e-07 2.167923e-06 0.0008055734
## Electro 1.185741e-07 7.273359e-08 3.155960e-07 0.0000710470
## Hip-Hop 1.936622e-07 4.360330e-07 1.472482e-06 0.0008361002
## Metal 1.034728e-07 7.441829e-07 1.651807e-06 0.0007834543
## Rap 1.021610e-08 2.864002e-07 1.082806e-06 0.0007023945
## Rock 1.913517e-07 6.341795e-07 2.235605e-06 0.0008007902
## duration time_signature
## Alternative 5.506016e-07 1.419443e-05
## Electro 3.657873e-07 1.663013e-06
## Hip-Hop 6.465757e-07 1.427949e-05
## Metal 5.031936e-07 1.399302e-05
## Rap 5.579644e-07 1.123101e-05
## Rock 5.370671e-07 1.429618e-05
##
## Residual Deviance: 8815.654
## AIC: 8995.654#p-values
z <- summary(train_model)$coefficients / summary(train_model)$standard.errors
2 * (1 - pnorm(abs(z)))## (Intercept) popularity danceability energy key loudness mode
## Alternative 0 0 0 0 0 0 0
## Electro 0 0 0 0 0 0 0
## Hip-Hop 0 0 0 0 0 0 0
## Metal 0 0 0 0 0 0 0
## Rap 0 0 0 0 0 0 0
## Rock 0 0 0 0 0 0 0
## speechiness acousticness instrumentalness liveness valence
## Alternative 0 0 0 0 0
## Electro 0 0 0 0 0
## Hip-Hop 0 0 0 0 0
## Metal 0 0 0 0 0
## Rap 0 0 0 0 0
## Rock 0 0 0 0 0
## tempo duration time_signature
## Alternative 9.740431e-12 0.000000e+00 0
## Electro 0.000000e+00 4.416761e-08 0
## Hip-Hop 1.863758e-01 1.063304e-02 0
## Metal 0.000000e+00 0.000000e+00 0
## Rap 0.000000e+00 0.000000e+00 0
## Rock 0.000000e+00 0.000000e+00 0#In sample Prognose
train_predict <- predict(train_model, train_data)
#Out of sample Prognose
test_predict <- predict(train_model, test_data)
#Confusion Matrix in sample
confusionMatrix(train_predict, train_data$genre)## Confusion Matrix and Statistics
##
## Reference
## Prediction Pop Alternative Electro Hip-Hop Metal Rap Rock
## Pop 255 43 40 75 3 47 67
## Alternative 22 243 45 26 56 20 137
## Electro 51 35 205 28 18 22 15
## Hip-Hop 39 31 21 247 0 128 10
## Metal 6 111 32 2 437 2 96
## Rap 20 26 27 123 1 192 10
## Rock 90 158 17 7 47 9 283
##
## Overall Statistics
##
## Accuracy : 0.5137
## 95% CI : (0.4972, 0.53)
## No Information Rate : 0.1785
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.43
##
## Mcnemar's Test P-Value : 3.488e-09
##
## Statistics by Class:
##
## Class: Pop Class: Alternative Class: Electro
## Sensitivity 0.52795 0.37558 0.52972
## Specificity 0.91248 0.89725 0.94781
## Pos Pred Value 0.48113 0.44262 0.54813
## Neg Pred Value 0.92633 0.86866 0.94402
## Prevalence 0.13324 0.17848 0.10676
## Detection Rate 0.07034 0.06703 0.05655
## Detection Prevalence 0.14621 0.15145 0.10317
## Balanced Accuracy 0.72021 0.63641 0.73876
## Class: Hip-Hop Class: Metal Class: Rap Class: Rock
## Sensitivity 0.48622 0.7776 0.45714 0.45793
## Specificity 0.92653 0.9187 0.93541 0.89092
## Pos Pred Value 0.51891 0.6370 0.48120 0.46318
## Neg Pred Value 0.91712 0.9575 0.92932 0.88885
## Prevalence 0.14014 0.1550 0.11586 0.17048
## Detection Rate 0.06814 0.1206 0.05297 0.07807
## Detection Prevalence 0.13131 0.1892 0.11007 0.16855
## Balanced Accuracy 0.70638 0.8481 0.69628 0.67442## Confusion Matrix and Statistics
##
## Reference
## Prediction Pop Alternative Electro Hip-Hop Metal Rap Rock
## Pop 82 18 17 31 2 11 28
## Alternative 5 77 14 7 18 6 60
## Electro 14 12 74 6 6 9 5
## Hip-Hop 18 9 12 89 0 38 3
## Metal 3 42 10 0 138 3 21
## Rap 3 8 10 30 0 57 1
## Rock 38 46 7 3 22 6 90
##
## Overall Statistics
##
## Accuracy : 0.5021
## 95% CI : (0.4735, 0.5306)
## No Information Rate : 0.1754
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.4161
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: Pop Class: Alternative Class: Electro
## Sensitivity 0.50307 0.36321 0.51389
## Specificity 0.89771 0.88967 0.95117
## Pos Pred Value 0.43386 0.41176 0.58730
## Neg Pred Value 0.92059 0.86791 0.93536
## Prevalence 0.13482 0.17535 0.11911
## Detection Rate 0.06782 0.06369 0.06121
## Detection Prevalence 0.15633 0.15467 0.10422
## Balanced Accuracy 0.70039 0.62644 0.73253
## Class: Hip-Hop Class: Metal Class: Rap Class: Rock
## Sensitivity 0.53614 0.7419 0.43846 0.43269
## Specificity 0.92330 0.9228 0.95181 0.87812
## Pos Pred Value 0.52663 0.6359 0.52294 0.42453
## Neg Pred Value 0.92596 0.9516 0.93364 0.88164
## Prevalence 0.13730 0.1538 0.10753 0.17204
## Detection Rate 0.07361 0.1141 0.04715 0.07444
## Detection Prevalence 0.13978 0.1795 0.09016 0.17535
## Balanced Accuracy 0.72972 0.8324 0.69513 0.65541Ein isolierter Blick auf die Accuracy des Modells würde auf keine gute Performance hindeuten, jedoch können wir mit Hilfe der Confusion Matrix gut erkennen, dass die schlechte Performance eher durch sehr ähnliche Genres beeinflusst wird. Beispielsweise sind die Grenzen zwischen den Kategorien Metal und Rock sicherlich fließend, was so auch vom Modell wahrgenommen wird. Zur Übung könnte man ein wenig experimentieren, was mit der Modellperformance passiert, wenn man ähnliche Kategorien, wie Metal und Rock zusammen in einer Kategorie abbildet.
4.2 Ordinal skalierte Zielvariable
Beim Modell ordinalen Regressionmodell gehen wir weiterhin von einer abhängigen Variablen mit meherern Kategorien aus, wobei nun eine Reihenfolge der Kategorien sinnvoll gebildet werden kann. Dennoch ist die Zielvariable weiterhin:
\[ Y_i = k, ~~~k = 1, ..., K~~~k, K \in \mathbb{N} \]
Entsprechend ändert sich im Vergleich zum kategorialen Regressionsmodell nichts an der Wahrscheinlichkeitsverteilung von \(Y_i\), was sich jedoch grundlegend ändert, ist die Frage, wie wir die Kategorienwahrscheinlichkeiten modellieren. Um eine Reihung zu berücksichtigen, brauchen wir einen Zusammenhang zwischen Höhe einer Variable und Kategorie der Zielvariable. Sinngemäß, je größer (oder kleiner) etwas ist, umso größer (oder kleiner) wird die Wahrscheinlichkeit einer in eine höhere Kategorie zu fallen.
Dies wird im ordinalen Regressionsmodell durch eine sogenannte latente (also nicht beobachtbare) Variable erreicht. Wir können uns hierunter z.B. vorstellen, dass gewisse Variablen den nicht beobachtbaren Gefallen einer Person approximieren, welcher Einfluss auf die Bewertung des Produktes hat (beispielsweise in Form eines fünf Sterne Systems). Gegeben erklärende Variablen \(\boldsymbol{x}_i^T\) und zugehörigen \(\beta-\)Koeffizienten resultiert der nicht beobachtbare latente Wert (\(U_i\))
\[ U_i = -\boldsymbol{x}_i^T \boldsymbol{\beta} + \epsilon_i, \]
wobei neben den deterministischen Zusammenhang \(\boldsymbol{x}_i^T \boldsymbol{\beta}\) der nicht beobachtbare Gefallen auch einer gewissen Zufälligkeit unterlegen ist, die durch \(\epsilon_i\) ausgedrückt wird. Nehmen wir nun an, \(U_i\) bzw. \(\epsilon_i\) ist standard normal verteilt. Fällt der Gefallen unter eine Grenze \(\theta_1\), so gibt der Kunde lediglich einen Stern. Die Wahrscheinlichkeit wäre:
\[ P(Y_i = 1) = P(U_i \leq \theta_1) = P(-\boldsymbol{x}_i^T \boldsymbol{\beta} + \epsilon_i \leq \theta_1) = \Phi(\theta_1 + \boldsymbol{x}_i^T \boldsymbol{\beta}) \]
Je kleiner \(\theta_1 + \boldsymbol{x}_i^T \boldsymbol{\beta}\), umso geringer die Wahrscheinlichkeit, dass \(Y_i = 1\). Weiter gilt für die nächste Kategorie, dass \(Y_i = 2\), falls \(U_i\) zwar über der Grenze \(\theta_1\), jedoch unter der nächsten Grenze \(\theta_2\) liegt, entsprechend:
\[ P(Y_i = 2) = P(U_i \leq \theta_2) - P(U_i \leq \theta_1) = \Phi(\theta_2 + \boldsymbol{x}_i^T \boldsymbol{\beta}) - \Phi(\theta_1 + \boldsymbol{x}_i^T \boldsymbol{\beta}) \] Und so geht es immer weiter, wobei die letzte Kategorie wieder als Gegenwahrscheinlichkeit dargestellt werden kann. Somit bei \(K\) Kategorien:
\[ P(Y_i = K) = 1 - P(U_i \leq \theta_{K-1}) \]
Im Gegensatz zum kategorialen Regressionsmodell haben wir es wieder “nur” mit einem Koeffizientenvektor \(\boldsymbol{\beta}\) zu tun. Hierbei setzen wir \(\beta_0 = 0\). Wenn man so möchte sind die Schwellwerte \(\theta_r\) kategorienspezifische Konstanten der Regressionsgerade.
Wird \(\epsilon_i \sim N(0,1)\) angenommen spricht man ordinalen probit Modell. Der Standard ist jedoch, dass man an Stelle der Normalverteilung die logistische Verteilung verwendet, was dem ordinalen logit Modell entspricht.
4.2.1 Simulation eines kategorialen Modells
Dem Motto was man simulieren kann, hat man verstanden zunächst ein simuliertes Beispiel mit drei Kategorien und einer erkärenden Variable.
library(ordinal)
#Simulation
#Datenumfang
n <- 200
#erster Cut Point
theta1 <- -1
#zweiter Cut Point
theta2 <- 1
#beta Koeffizient
beta <- 1.5
#erklärende Variable
set.seed(123)
x <- rnorm(n)
#zufällige Abweichungen
eps <- rnorm(n, mean = 0, sd = 1)
#latente Variable
U <- beta * x + eps
y <- c()
for(i in 1:n){
if(U[i] <= theta1){
y <- c(y, 1)
}else{
if(U[i] > theta1 & U[i] <= theta2){
y <- c(y, 2)
}else{
y <- c(y, 3)
}
}
}
table(y)## y
## 1 2 3
## 53 98 49## formula: as.factor(y) ~ x
##
## link threshold nobs logLik AIC niter max.grad cond.H
## probit flexible 200 -139.48 284.95 5(0) 1.46e-10 3.5e+00
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## x 1.3609 0.1392 9.779 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -1.0469 0.1286 -8.143
## 2|3 1.0654 0.1314 8.1114.2.2 Realbeispiel
Um einen realen Datensatz zu analysieren, verwenden wir einen Datensatz mit Imdb-Ratings, welche wir durch Runden in ordinale Kategorien bringen.
library(data.table)
library(plyr)
library(plotly)
library(ordinal)
library(DescTools)
dt <- data.table(read.csv('Data/imdb_movie_data.csv'))
#welche Datentypen haben wir
str(dt)## Classes 'data.table' and 'data.frame': 16744 obs. of 17 variables:
## $ X : int 0 1 2 3 4 5 6 7 8 9 ...
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Title : Factor w/ 16744 levels "¡Asu Mare! 2",..: 6728 14145 1486 1563 13558 11902 14430 4129 10469 6767 ...
## $ Year : int 2010 1999 2018 1985 1966 2018 2002 2012 1981 2009 ...
## $ Age : Factor w/ 6 levels "","13+","16+",..: 2 4 2 5 4 5 4 4 5 4 ...
## $ IMDb : num 8.8 8.7 8.5 8.5 8.8 8.4 8.5 8.4 8.4 8.3 ...
## $ Rotten.Tomatoes: Factor w/ 100 levels "","10%","100%",..: 87 87 84 97 98 98 96 87 96 89 ...
## $ Netflix : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Hulu : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Prime.Video : int 0 0 0 0 1 0 1 0 0 0 ...
## $ Disney. : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Type : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Directors : Factor w/ 11339 levels "","A. Bhimsingh",..: 1981 6287 803 9128 9753 1277 9219 8551 10210 8551 ...
## $ Genres : Factor w/ 1910 levels "","Action","Action,Adventure",..: 178 407 176 509 1906 663 941 1678 3 572 ...
## $ Country : Factor w/ 1304 levels "","Afghanistan,France",..: 1255 1061 1061 1061 583 1061 952 1061 1061 442 ...
## $ Language : Factor w/ 1103 levels "","Aboriginal,English",..: 352 103 103 103 785 493 279 261 268 261 ...
## $ Runtime : int 148 136 149 116 161 117 150 165 115 153 ...
## - attr(*, ".internal.selfref")=<externalptr>## X ID Title Year Age
## 0 0 0 0 0
## IMDb Rotten.Tomatoes Netflix Hulu Prime.Video
## 571 0 0 0 0
## Disney. Type Directors Genres Country
## 0 0 0 0 0
## Language Runtime
## 0 592#Wir löschen alle Zeilen, sobald eine NA auftaucht
dt <- dt[complete.cases(dt), ]
#Erzeuge eine ordinale Variable für die Alterskategorien
from_age_map <- c(sort(levels(dt$Age)))
to_age_map <- 1:length(from_age_map)
dt[, Age_ordinal := as.numeric(mapvalues(Age, from = from_age_map, to = to_age_map))]
#dt[, Age_ordinal := mapvalues(Age, from = from_age_map, to = to_age_map)]
head(dt)## X ID Title Year Age IMDb Rotten.Tomatoes Netflix
## 1: 0 1 Inception 2010 13+ 8.8 87% 1
## 2: 1 2 The Matrix 1999 18+ 8.7 87% 1
## 3: 2 3 Avengers: Infinity War 2018 13+ 8.5 84% 1
## 4: 3 4 Back to the Future 1985 7+ 8.5 96% 1
## 5: 4 5 The Good, the Bad and the Ugly 1966 18+ 8.8 97% 1
## 6: 5 6 Spider-Man: Into the Spider-Verse 2018 7+ 8.4 97% 1
## Hulu Prime.Video Disney. Type Directors
## 1: 0 0 0 0 Christopher Nolan
## 2: 0 0 0 0 Lana Wachowski,Lilly Wachowski
## 3: 0 0 0 0 Anthony Russo,Joe Russo
## 4: 0 0 0 0 Robert Zemeckis
## 5: 0 1 0 0 Sergio Leone
## 6: 0 0 0 0 Bob Persichetti,Peter Ramsey,Rodney Rothman
## Genres Country
## 1: Action,Adventure,Sci-Fi,Thriller United States,United Kingdom
## 2: Action,Sci-Fi United States
## 3: Action,Adventure,Sci-Fi United States
## 4: Adventure,Comedy,Sci-Fi United States
## 5: Western Italy,Spain,West Germany
## 6: Animation,Action,Adventure,Family,Sci-Fi United States
## Language Runtime Age_ordinal
## 1: English,Japanese,French 148 2
## 2: English 136 4
## 3: English 149 2
## 4: English 116 5
## 5: Italian 161 4
## 6: English,Spanish 117 5#Betrachte die Verteilung der Alterklassen
plot_ly()%>%
add_bars(x = sort(unique(dt$Age_ordinal)), y = table(dt$Age_ordinal)) %>%
layout(title = 'Verteilung Altersklassen')#Wir schätzen das Modell mit den erklärenden Variablen: Year, Age, Runtime
#Betrachten wir kurz die erklärenden Variablen
plot_ly()%>%
add_histogram(x = dt$Runtime)%>%
layout(title = 'Laufzeit')## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 83.00 92.00 93.77 104.00 1256.00## X ID Title Year Age IMDb Rotten.Tomatoes Netflix
## 1: 13179 13180 Colorado 1940 all 5.9 0
## 2: 4405 4406 Law of the Lawless 1964 6.1 0
## 3: 15295 15296 The Vatican Museums 2007 5.0 0
## 4: 12475 12476 Scarlett 2016 16+ 4.4 0
## 5: 14113 14114 The Inner Circle 2009 6.6 0
## ---
## 15815: 7613 7614 For the Birds 2016 7.9 0
## 15816: 14298 14299 Thanksgiving 2014 4.7 0
## 15817: 15142 15143 #LoveSwag 2015 13+ 5.4 0
## 15818: 16456 16457 Luxo Jr. 1986 all 7.3 0
## 15819: 2626 2627 Liefling The Movie 2010 6.3 1
## Hulu Prime.Video Disney. Type Directors
## 1: 0 1 0 0
## 2: 1 1 0 0
## 3: 0 1 0 0 Luca De Mata
## 4: 0 1 0 0
## 5: 0 1 0 0
## ---
## 15815: 0 1 0 0 Ralph Eggleston
## 15816: 0 1 0 0 Eli Roth
## 15817: 0 1 0 0 Austin Davoren
## 15818: 0 0 1 0 John Lasseter
## 15819: 0 0 0 0 Kendrick Hew
## Genres
## 1: Action,Adventure,Drama,Romance,Western
## 2: Action,Crime,Drama
## 3:
## 4: Drama,Romance
## 5: Drama
## ---
## 15815: Animation,Short,Comedy,Family
## 15816: Short,Comedy,Horror
## 15817: Short,Comedy,Drama,Romance
## 15818: Animation,Short,Family
## 15819: Short,Thriller
## Country Language
## 1: United States English
## 2: Russia Russian
## 3: Italy English
## 4: France,United States,Germany,Italy,United Kingdom,Spain,Austria English
## 5: Sweden Swedish
## ---
## 15815: United States None
## 15816: United States English
## 15817: United States English
## 15818: United States None
## 15819: Canada English
## Runtime Age_ordinal
## 1: 1256 6
## 2: 750 1
## 3: 410 1
## 4: 360 3
## 5: 360 1
## ---
## 15815: 3 1
## 15816: 2 1
## 15817: 2 2
## 15818: 2 6
## 15819: 1 1#ein kurzer Abgleich mit Imbd lässt vermuten, dass es sich bei den extremen Werten ggf. um Ausreisser handelt
#wir schätzen später auch ein Modell bei dem wir Filme mit extremen Laufzeiten heraus filtern
#Jahreszahl
plot_ly()%>%
add_histogram(x = dt$Year)%>%
layout(title = 'Jahreszahl')#offensichtlich sind in moderneren Zeiten mehr Filme gedreht worden
#als abhängige Variable haben wir das Imdb-Rating
plot_ly()%>%
add_histogram(x = dt$IMDb)%>%
layout(title = 'Imdb Rating')#für unser Modell erzeugen wir ordinale Kategorien
dt[, rating_class := as.factor(ceiling(dt$IMDb))]
table(dt$rating_class)##
## 0 1 2 3 4 5 6 7 8 9 10
## 4 2 58 466 1084 2225 3862 4831 2824 454 9plot_ly()%>%
add_bars(x = sort(unique(dt$rating_class)), y = table(dt$rating_class)) %>%
layout(title = 'Verteilung Imdb Ratingklassen')#nun schätzen wir mit der clm-Funktion
ord_model <- clm(rating_class ~ Age_ordinal + Year + Runtime, data = dt)
summary(ord_model)## formula: rating_class ~ Age_ordinal + Year + Runtime
## data: dt
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 15819 -26933.95 53893.91 11(0) 7.90e-09 3.4e+11
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## Age_ordinal -0.0389905 0.0082634 -4.718 2.38e-06 ***
## Year -0.0009758 0.0006500 -1.501 0.133
## Runtime 0.0064748 0.0005679 11.401 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1 -9.73633 1.38824 -7.013
## 1|2 -9.33071 1.35790 -6.871
## 2|3 -6.95978 1.30108 -5.349
## 3|4 -4.81500 1.29568 -3.716
## 4|5 -3.62581 1.29498 -2.800
## 5|6 -2.58482 1.29451 -1.997
## 6|7 -1.49213 1.29432 -1.153
## 7|8 -0.09152 1.29456 -0.071
## 8|9 2.07791 1.29548 1.604
## 9|10 6.04929 1.33691 4.525#aufgrund der Warnmeldung reskalieren wir die erkärenden Variablen
dt[, Age_ord_rescaled := (Age_ordinal - mean(Age_ordinal)) / sd(Age_ordinal)]
dt[, year_rescaled := (Year - mean(Year)) / sd(Year)]
dt[, rt_ord_rescaled := (Runtime - mean(Runtime)) / sd(Runtime)]
ord_model <- clm(rating_class ~ Age_ord_rescaled + year_rescaled + rt_ord_rescaled, data = dt)
summary(ord_model)## formula: rating_class ~ Age_ord_rescaled + year_rescaled + rt_ord_rescaled
## data: dt
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 15819 -26933.95 53893.91 11(0) 4.56e-12 3.9e+03
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## Age_ord_rescaled -0.06671 0.01414 -4.718 2.38e-06 ***
## year_rescaled -0.02046 0.01363 -1.501 0.133
## rt_ord_rescaled 0.18099 0.01588 11.401 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1 -8.29616 0.50007 -16.590
## 1|2 -7.89054 0.40833 -19.324
## 2|3 -5.51961 0.12527 -44.063
## 3|4 -3.37482 0.04422 -76.326
## 4|5 -2.18564 0.02632 -83.055
## 5|6 -1.14464 0.01860 -61.524
## 6|7 -0.05195 0.01597 -3.253
## 7|8 1.34866 0.01968 68.522
## 8|9 3.51808 0.04723 74.485
## 9|10 7.48946 0.33344 22.461#der Output spricht dafür, dass ältere Filme niederer Alterklassen mit langer Laufzeit ein besseres Rating haben
#die Aussagekraft dieser Ergebnisse steht hier allerdings nicht so sehr im Vordergrund
#Testen wir noch, ob sich etwas ändert, wenn wir extreme Daten der Runtime ausschließen
dt[, Runtime_win := Winsorize(Runtime)]
dt[, rt_ord_rescaled_win := (Runtime_win - mean(Runtime_win)) / sd(Runtime_win)]
ord_model <- clm(rating_class ~ Age_ord_rescaled + year_rescaled + rt_ord_rescaled_win, data = dt)
summary(ord_model)## formula: rating_class ~ Age_ord_rescaled + year_rescaled + rt_ord_rescaled_win
## data: dt
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 15819 -26901.31 53828.63 10(0) 7.34e-12 3.9e+03
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## Age_ord_rescaled -0.07110 0.01415 -5.024 5.05e-07 ***
## year_rescaled -0.03315 0.01368 -2.423 0.0154 *
## rt_ord_rescaled_win 0.20750 0.01467 14.145 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1 -8.30095 0.50007 -16.600
## 1|2 -7.89531 0.40833 -19.336
## 2|3 -5.52425 0.12527 -44.099
## 3|4 -3.37893 0.04423 -76.402
## 4|5 -2.18858 0.02633 -83.124
## 5|6 -1.14550 0.01862 -61.520
## 6|7 -0.04953 0.01599 -3.097
## 7|8 1.35497 0.01973 68.681
## 8|9 3.52633 0.04726 74.612
## 9|10 7.49755 0.33344 22.485