Lektion 7 – eigene Funktion schreiben

In der Lektion 6 haben wir angeschaut, wie man eine logistische Regression schätzen kann. In dieser Lektion schreibe wir eine eigene Funktion. Dabei schauen wir uns an, was man dazu alles beachten muss und wir man vorgeht. Das anhand eine Beispiels einer Simulation mit einer logistischen Regression.

Das folgende Video erklärt alles Schritt für Schritt. Wie immer folgt unterhalb des Videos eine Zusammenfassung der wichtigsten Punkte.

Eine Funktion ist eine Handlungsanweisung mit Input und Output. Eine eigene Funktion schreiben wir, wenn wir etwas mehr als zwei Mal machen. Damit verhindern wir häufiges Repetieren und vermeiden Fehler durch übersichtlicheren Code. Dabei gehen wir wie folgt vor:

  1. ein Beispiel schreiben
  2. relevante Variablen identifizieren und durch temporäre Variablen ersetzen
  3. den Funktionskopf hinzufügen mit diesen Variablen als Argumente
  4. die Funktion testen
  5. den Input kontrollieren mit Überprüfungen

Beispiel

Zuerst erstellen wir ein Beispiel in dem wir vorausgesagte Wahrscheinlichkeiten mit Konfidenzintervall berechnen (wenn du das Beispiel nicht verstehst, empfehle ich dir Lektion 6, wo das ganze erklärt wird):

df_selects = schlegel::selects2015
model_logit = glm(participation ~ age * gender, 
                  family = "binomial", data = df_selects)
betas = coef(model_logit)
vcov = vcov(model_logit)
betas_sim = MASS::mvrnorm(1000, betas, vcov)
x = c(1, 24, 0, 0)
yhat = betas_sim %*% x
pred = exp(yhat) / (1 + exp(yhat))
conf.int = quantile(pred, probs = c(0.025, 0.975))
c(conf.int[1], mean(pred), conf.int[2])

Als nächstes identifizieren wir die relevanten Variablen und ersetzen sie durch temporäre Variablen:

betas = coef(model) # model
vcov = vcov(model) # model
betas_sim = MASS::mvrnorm(n, betas, vcov) # n
x = x_values # x
yhat = betas_sim %*% x
pred = exp(yhat) / (1 + exp(yhat))
conf.int = quantile(pred, probs = c(0.025, 0.975))
c(conf.int[1], mean(pred), conf.int[2])

Nun müssen wir noch den Funktionskopf hinzufügen:

predProb = function(model, x_values, n){
  betas = coef(model) # model
  vcov = vcov(model) # model
  betas_sim = MASS::mvrnorm(n, betas, vcov) # n
  x = x_values # x
  yhat = betas_sim %*% x
  pred = exp(yhat) / (1 + exp(yhat))
  conf.int = quantile(pred, probs = c(0.025, 0.975))
  c(conf.int[1], mean(pred), conf.int[2])
}

Als nächstes testen wir die Funktion:

predProb(model_logit, c(1, 24, 0, 0), 1000)

Als nächstes überprüfen wir die Argumente (das ist für den eigenen gebraucht nicht zwingend notwendig, Hilfe aber Fehler zu erkennen):

predProb = function(model, x_values, n){
  if(!("glm" %in% class(model))){
    stop("Das Modell ist kein glm()-Modell!")
  }
  if(model$family$link != "logit"){
    stop("Das Modell ist kein Logit Modell!")
  }
  if(!is.numeric(x_values)){
    stop("Die x_values müssen numerisch sein!")
  }
  if(length(x_values) != length(coef(model))){
    stop("x_values hat nicht die gleiche Länge
          wie es Koeffizienten hat!")
  }
  if(!is.numeric(n) || round(n) != n){
    stop("n muss eine ganze Zahl sein!")
  }
  betas = coef(model) # model
  vcov = vcov(model) # model
  betas_sim = MASS::mvrnorm(n, betas, vcov) # n
  x = x_values # x
  yhat = betas_sim %*% x
  pred = exp(yhat) / (1 + exp(yhat))
  conf.int = quantile(pred, probs = c(0.025, 0.975))
  c(conf.int[1], mean(pred), conf.int[2])
}

Nun können wir noch einen Standardwert für n hinzufügen, damit die Angabe freiwillig ist.

predProb = function(model, x_values, n = 1000){
  if(!("glm" %in% class(model))){
    stop("Das Modell ist kein glm()-Modell!")
  }
  if(model$family$link != "logit"){
    stop("Das Modell ist kein Logit Modell!")
  }
  if(!is.numeric(x_values)){
    stop("Die x_values müssen numerisch sein!")
  }
  if(length(x_values) != length(coef(model))){
    stop("x_values hat nicht die gleiche Länge 
          wie es Koeffizienten hat!")
  }
  if(!is.numeric(n) || round(n) != n){
    stop("n muss eine ganze Zahl sein!")
  }
  betas = coef(model) # model
  vcov = vcov(model) # model
  betas_sim = MASS::mvrnorm(n, betas, vcov) # n
  x = x_values # x
  yhat = betas_sim %*% x
  pred = exp(yhat) / (1 + exp(yhat))
  conf.int = quantile(pred, probs = c(0.025, 0.975))
  c(conf.int[1], mean(pred), conf.int[2])
}

Mit der Standardeinstellung von 1000:

predProb(model_logit, c(1, 24, 0, 0))

Mit 10’000:

predProb(model_logit, c(1, 24, 0, 0), 10000)

Nun ist die Funktion fertig.