9  Geradenmodelle 1

Schlüsselwörter

Statistik, Prognose, Modellierung, R, Datenanalyse, Regression

9.1 Lernsteuerung

Abb. Abbildung 1.3 zeigt den Standort dieses Kapitels im Lernpfad und gibt damit einen Überblick über das Thema dieses Kapitels im Kontext aller Kapitel.

9.1.1 Lernziele

• Sie können ein Punktmodell von einem Geradenmodell begrifflich unterscheiden.
• Sie können die Bestandteile eines Geradenmodells aufzählen und erläutern.
• Sie können die Güte eines Geradenmodells anhand von Kennzahlen bestimmen.
• Sie können Geradenmodelle sowie ihre Modellgüte in R berechnen.

9.1.2 Benötigte R-Pakete

library(tidyverse)
library(easystats)

\[ \definecolor{ycol}{RGB}{230,159,0} \definecolor{modelcol}{RGB}{86,180,233} \definecolor{errorcol}{RGB}{0,158,115} \definecolor{beta0col}{RGB}{213,94,0} \definecolor{beta1col}{RGB}{0,114,178} \definecolor{xcol}{RGB}{204,121,167} \]

9.1.3 Benötigte Daten

mariokart <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/openintro/mariokart.csv")

9.2 Vorhersagen

Vorhersagen sind eine nützliche Sache, unter (mindestens) folgenden Voraussetzungen:

1. Sie sind präzise
2. Wir kennen die Präzision
3. Jemand interessiert sich für die Vorhersage

Die Methode des Vorhersagens, die wir hier betrachten, nennt man auch lineare Regression.

9.2.1 Vorhersagen ohne Prädiktor

Beispiel 9.1 Nach intensiver Beschäftigung mit Statistik sind Sie allgemein als Checker bekannt. Viele jüngere Studis fragen Sie um Rat. eines Tages kommt eine Studentin, Toni, und fragt: “Welche Statistiknote kann ich in der Klausur erwarten?” Sie entgegnen: “Wie viel hast du denn gelernt?”. Die Antwort: “Sag ich nicht.”

Nach kurzem Überlegen geben sie den Notenschnitt der letzten Klausur als Prognose für diese Person. Dazu rechnen Sie schnell den Notenschnitt (Mittelwert) aus.

Zuerst importieren Sie die Daten der letzten Klausur. Die Syntax in Listing 9.1 wird bei Ihnen nur funktionieren, wenn auf Ihrem Computer dieser Ordner mit dieser Datei existiert. Andernfalls müssen Sie die Daten erst herunterladen¹:

Listing 9.1: Der Datensatz ‘noten2’ liegt im Unterordner ‘Noten.’

noten2 <- read.csv("daten/noten2.csv")

Dann rechnen Sie den Mittelwert aus:

noten2 %>% 
  summarise(mw = mean(y))  # y ist der Punktwert in der Klausur

Ihre Antwort lautet also: “Im Schnitt haben die Studis bei der letzten Klausur gut 70% der Punkte erzielt. Diesen Wert kannst du erwarten. Solange ich keine genaueren Infos habe, z.B. wieviel du gelernt hast, kann ich dir keine genauere Vorhersage machen, sorry!”\(\square\)

Hinweis

Ohne Kenntnis eines Prädiktors (UV) (wie z.B. Lernzeit) ist der Mittelwert ein geeigneter Vorhersagewert für jede Beobachtung, s. Abbildung 9.1. Wir nutzen den Mittelwert als Punktmodell für den Klausurerfolg.\(\square\)

Abbildung 9.1: Mittelwert als Vorhersagewert, bzw. Mittelwert als Punktmodell

9.2.2 Nullmodell (Punktmodell)

Modelle ohne Prädiktor, Punktmodelle also, kann man so bezeichnen: y ~ 1. Da das Modell null Prädiktoren hat, nennt man es auch manchmal “Nullmodell”.

Auf Errisch kann man dieses Nullmodell so spezifizieren:

lm0 <- lm(y ~ 1, data = noten2)
lm0
## 
## Call:
## lm(formula = y ~ 1, data = noten2)
## 
## Coefficients:
## (Intercept)  
##        71.1

lm steht für “lineares Modell”, die 1 sagt, dass es keine Prädiktoren gibt. In dem Fall wird der Mittelwert als Gerade verwendet. Der zurückgemeldete Koeffizient (Intercept) ist hier der Modell des Punktmodells. Da es ein Punktmodell ist, sagt es für alle Beobachtungen (hier Studentis) den gleichen Wert vorher.

Die Regressionsgleichung lautet demnach: y_pred = 71.08. In Worten: “Wir sagen für jede Beobachtung einen Wert von ca. 71 vorher”.

9.2.3 Vorhersagen mit Prädiktor

Beispiel 9.2 (Toni verrät die Lernzeit) Toni entschließt sich dann doch noch, die Lernzeit zu verraten: “Okay, also ich hab insgesamt 42 Stunden gelernt, insgesamt.” Jetzt müssen Sie erstmal nachdenken: “Wie viele Klausurpunkte sag ich vorher, wenn Toni 42 Stunden gelernt hat?”

Sie visualisieren sich zur Hilfe die vorliegenden Daten, s. Abbildung 9.2, a).²

library(DataExplorer)
noten2 %>% 
  plot_scatterplot(by = "y")  # Y-Variable muss angegeben werden

Auf dieser Basis antworten Sie Toni: “Bei 42 Stunden Lernzeit solltest du so 46 Punkte bekommen. Könnte mit dem Bestehen eng werden.” Toni ist nicht begeistert von Ihrer Prognose und zieht von dannen.\(\square\)

Der Trend (im Sinne eines linearen Zusammenhangs) von Lernzeit und Klausurpunkte ist deutlich zu erkennen. Mit einem Lineal könnte man eine entsprechende Gerade in das Streudiagramm einzeichnen, s. Abbildung 9.2, b).

(a) Gemeinsame Verteilung (Zusammenhang) von Lernzeit (X) und Noten (Y)

(b) Eine ‘Trendgerade’ (blau) im Datensatz noten2. Mittelwerte sind mit gestrichelten Linien eingezeichnet. Tonis Vorhersage ist mit einem Punkt gekennzeichnet.

Abbildung 9.2: Noten und Lernzeit: Rohdaten und Modell

Eine Gerade eignet sich, um einen linearen Trend zusammenzufassen.

9.3 Geradenmodelle

9.3.1 Achsenabschnitt und Steigung definieren eine Gerade

Wir verwenden eine Gerade als Modell für die Daten, s. Abbildung 9.2, b. Anders gesagt: Wir modellieren die Daten (bzw. deren Zusammenhang) mit einer Geraden.

Ein Geradenmodell ist eine Verallgemeinerung des Punktmodells: Ein Punktmodell sagt für alle Beobachtungen den gleichen Wert vorher. Abbildung 9.1 und Abbildung 9.2 stellen ein Punktmodell einem Geradenmodell gegenüber.

In einem Geradenmodell wird nicht mehr (notwendig) für jede Beobachtung die gleiche Vorhersage \(\hat{y}\) gemacht (wie das bei einem Punktmodell der Fall ist).

Definition 9.1 (Gerade) Eine Gerade ist das, was man bekommt, wenn man eine lineare Funktion in ein Koordinatensystem einzeichnet. Man kann sie durch durch zwei Koeffizienten festlegen: Achsenabschnitt (engl. intercept), und Steigung (engl. slope). Häufig wird (z.B. im Schulunterricht) der Achsenabschnitt mit \(t\) und die Steigung mit \(m\) bezeichnet:

\(f(\color{xcol}{x})=\color{ycol}{y}={m} \color{xcol}{x} + \color{beta0col}{t}\).

In der Statistik wird folgende Nomenklatur bevorzugt: \(f(\color{xcol}{x})=\color{ycol}{\hat{y}}=\color{beta0col}{\beta_0} + {\beta_1} \color{xcol}{x}\) oder \(f(\color{xcol}{x})=\color{ycol}{\hat{y}}= \color{beta0col}{b_0} + {b_1} \color{xcol}{x}\) .

Die Nomenklatur mit \(b_0, b_1\) hat den Vorteil, dass man das Modell einfach erweitern kann: \(b_2, b_3, ...\). Anstelle von \(b\) liest man auch oft \(\beta\). Griechische Buchstaben werden meist verwendet, um zu zeigen, dass man an einer Aussage über eine Population, nicht nur über eine Stichprobe, machen möchte.

Das “Dach” über y, \(\color{modelcol}{\hat{y}}\), drückt aus, dass es sich den den geschätzten, bzw. vom Modell vorhergesagten (“modellierten”) Wert für \(\color{ycol}{y}\) handelt, nicht das tatsächliche (empirische, beobachtete) \(\color{ycol}{y}\). \(\square\)

Abbildung 9.3 skizziert die Elemente einer Regression. (Bildquelle: Basierend auf TikZ-Quellcode von Henri Menke.)

Abbildung 9.3: Achsenabschnitt und Steigung einer Regressionsgeraden (Menk, 2014)

Das einfache lineare Modell

Das einfache lineare Modell nimmt den Wert einer abhängigen metrischen Variablen, y als lineare Funktion von unabhängigen Variablen, x an, plus einem Fehlerterm, e. \(\square\)

\[\begin{aligned} \color{ycol}{y} &= f(\color{xcol}{x}) + \color{errorcol}{\epsilon} \\ \color{ycol}{y_i} &= \color{beta0col}{\beta_0} + \color{beta1col}{\beta_1} \cdot \color{modelcol}{x_i} + \color{errorcol}{\epsilon_i} \square \end{aligned}\]

Mit:

• \(\color{beta0col}{\beta_0}\): geschätzter y-Achsenabschnitt laut Modell
• \(\color{beta1col}{\beta_1}\): geschätzte Steigung laut Modell
• \(\color{errorcol}{\epsilon}\): Fehler des Modells

Je nach Datenlage können sich Regressionsgerade in Steigung oder Achsenabschnitt unterscheiden, s. Abbildung 9.4.

(a) Datensatz 1

(b) Datensatz 2

Abbildung 9.4: Regressionsanalysen mit verschiedenen Koeffizienten, aber gleicher Modellgüte

Abbildung 9.5 zeigt ein interaktives Beispiel einer linearen Funktion. Sie können Punkte per Klick/Touch hinzufügen.

resetButton = {
  const resetButton = Inputs.button("Reset");

  d3.select(resetButton).on("input", () => {

    regressionPlot.reset();
  });

  return resetButton;
}

viewOptions = {
  const viewOptions = Inputs.radio(
    ["None", "Absolute Error", "Squared Error"],
    { label: "View", value: "Absolute Error" }
  );

  d3.select(viewOptions).on("input", () => {
    regressionPlot.updateView(viewOptions.value);
  });

  return viewOptions;
}

rSquaredPlot = RSquaredPlot({ width: width })

regressionPlot = {
  const regressionPlot = RegressionPlot(data.slice(0, 1), {
    width: width,
    xDomain: [0, d3.max(data, ([x]) => x) + 5],
    yDomain: [0, d3.max(data, ([_, y]) => y) + 5],
    r: 6,
    showGrid: true
  });
  // Attach listener
  d3.select(regressionPlot).on("input", function () {
    rSquaredPlot.update(this.value, this.transition);
  });

  rSquaredPlot.update(regressionPlot.value);

  return regressionPlot;
}

width = 800
data = {
  const numPoints = 15;
  const xScale = 10;
  const xShift = 5;
  const yScale = 6;
  const yShift = 0;

  const data = d3.range(numPoints).map((i) => {
    const xCoord = xShift + xScale * d3.randomUniform()();
    const yCoord = yShift + xCoord + yScale * d3.randomUniform()();

    return [xCoord, yCoord];
  });

  return data;
}

animation = {
  // Perform opening animation
  const totalTime = 6_000;
  const ease = d3.easeQuadOut;

  d3.range(1, data.length).forEach((i) => {
    setTimeout(() => {
      const [xCoord, yCoord] = data[i];

      regressionPlot.updateDatapoint(xCoord, yCoord);
    }, totalTime * ease(i / data.length));
  });
}

function RSquaredPlot({
  data,
  marginTop = 0, // top margin, in pixels
  marginRight = 30, // right margin, in pixels
  marginBottom = 0, // bottom margin, in pixels
  marginLeft = 53, // left margin, in pixels
  width = 640, // outer width, in pixels
  height = 20, // outer height, in pixels
  barHeight = 25, // height of bar, in pixels
  titleSize = 20, // fontsize of title text
  labelSize = 12 // fontsize of label text
} = {}) {
  if (data === undefined) data = 0;

  const xScale = d3
    .scaleLinear()
    .domain([0, 1])
    .range([marginLeft, width - marginRight]);

  const svg = d3
    .create("svg")
    .attr("width", width)
    .attr("height", height)
    .attr("viewBox", [0, 0, width, height])
    .attr("style", "max-width: 100%; height: auto; height: intrinsic;");

  // Draw title.
  const titleGroup = svg.append("g");

  // titleGroup
  //   .append("line")
  //   .attr("stroke", "black")
  //   .attr("x1", xScale(0))
  //   .attr("x2", xScale(0))
  //   .attr("y1", 0)
  //   .attr("y2", height);

  titleGroup
    .append("text")
    .attr("fill", "black")
    .attr("text-anchor", "end")
    .attr("dominant-baseline", "middle")
    .attr("font-family", "serif")
    .attr("font-size", titleSize)
    .attr("x", xScale(0))
    .attr("dx", -10)
    .attr("y", height / 2)
    .text("R²");

  // Draw whole bar.
  svg
    .append("g")
    .append("rect")
    .attr("fill", "gainsboro")
    .attr("x", xScale(0))
    .attr("y", marginTop)
    .attr("width", width - marginLeft - marginRight)
    .attr("height", height - marginTop - marginBottom);

  // Draw bar.
  const bar = updateBar(svg.append("g").selectAll("rect"), 0);

  // Draw label
  const labelFormat = d3.format(".2f");
  const label = updateLabel(svg.append("g").selectAll("text"), 0);

  // Update according to data.
  update(data);

  function updateBar(rect, data, transition = true) {
    return rect.data([data]).join(
      (enter) =>
        enter
          .append("rect")
          .attr("fill", "black")
          .attr("x", xScale(0))
          .attr("y", marginTop)
          .attr("width", 0)
          .attr("height", height - marginTop - marginBottom),
      (update) =>
        update.call((update) => {
          if (transition) update = update.transition();

          update.attr("width", (d) => xScale(isNaN(d) ? 0 : d) - xScale(0));
        })
    );
  }

  function updateLabel(text, data, transition = true) {
    return text.data([data]).join(
      (enter) =>
        enter
          .append("text")
          .attr("dominant-baseline", "middle")
          .attr("font-family", "sans-serif")
          .attr("font-size", labelSize)
          .attr("font-weight", "bold")
          .attr("fill", "black")
          .attr("text-anchor", "start")
          .attr("dx", 10)
          .attr("x", xScale(0))
          .attr("y", height / 2)
          .text(labelFormat("0")),
      (update) =>
        update.call((update) => {
          // Check if bar is too short
          const check = (d) => d < 0.1;

          if (transition) update = update.transition();

          update
            .attr("x", (d) => xScale(d))
            .text((d) => labelFormat(d))
            .attr("fill", (d) => (check(d) ? "black" : "white"))
            .attr("text-anchor", (d) => (check(d) ? "start" : "end"))
            .attr("dx", (d) => (check(d) ? 10 : -10));
        })
    );
  }

  // Main function for data updates.
  function update(data, transition = true) {
    updateBar(bar, data, transition);
    updateLabel(label, data, transition);
  }

  return Object.assign(svg.node(), { update });
}

// to draw lines at the origin (to show intercept and stuff)
// create voronoi overlay? for dragging points
// viz error / squared error
// make the error lines look like a weight
// show other metrics like R^2 etc.
function RegressionPlot(
  data,
  {
    x = ([x]) => x, // accessor function for x-coordinate
    y = ([, y]) => y, // accessor function for y-coordinate
    r = 6, // radius of dots, in pixels
    marginTop = 20, // top margin, in pixels
    marginRight = 30, // right margin, in pixels
    marginBottom = 30, // bottom margin, in pixels
    marginLeft = 40, // left margin, in pixels
    inset = r * 2, // inset the default range, in pixels
    insetTop = inset, // inset the default y-range
    insetRight = inset, // inset the default x-range
    insetBottom = inset, // inset the default y-range
    insetLeft = inset, // inset the default x-range
    width = 640, // outer width, in pixels
    height = 500, // outer height, in pixels
    xType = d3.scaleLinear, // type of x-scale
    xDomain, // [xmin, xmax]
    xRange = [marginLeft + insetLeft, width - marginRight - insetRight], // [left, right]
    yType = d3.scaleLinear, // type of y-scale
    yDomain, // [ymin, ymax]
    yRange = [height - marginBottom - insetBottom, marginTop + insetTop], // [bottom, top]
    showGrid = true // whether to show grid lines
  } = {}
) {
  const errorColor = d3.schemeSet3[3];

  let errorLinesOpacity = 1;
  let errorSquaresOpacity = 0;

  // Init drag object.
  const drag = d3
    .drag()
    // .on("start", dragstarted)
    .on("drag", dragged);
  // .on("end", dragended);

  // To give each data point a unique id
  let idCounter = 0;

  // Compute data values.
  data = data.map((d) => ({
    xCoord: x(d),
    yCoord: y(d),
    id: idCounter++
  }));

  // Store copy of initial data object
  const initData = data.map((d) => ({ ...d }));

  // Compute default domains.
  // if (xDomain === undefined) xDomain = d3.extent(data, (d) => d.xCoord);
  if (xDomain === undefined) xDomain = [0, d3.max(data, (d) => d.xCoord)];
  if (yDomain === undefined) yDomain = [0, d3.max(data, (d) => d.yCoord)];

  // Init linear regressor.
  const linearRegression = d3
    .regressionLinear()
    .x((d) => d.xCoord)
    .y((d) => d.yCoord)
    .domain(xDomain);

  // Construct scales and axes.
  const xScale = xType(xDomain, xRange);
  const yScale = yType(yDomain, yRange);
  const xAxis = d3.axisBottom(xScale).ticks(width / 80);
  const yAxis = d3.axisLeft(yScale).ticks(height / 80);

  // Draw svg.
  const svg = d3
    .create("svg")
    .attr("width", width)
    .attr("height", height)
    .attr("viewBox", [0, 0, width, height])
    .attr("style", "max-width: 100%; height: auto; height: intrinsic;");

  // Draw xAxis.
  const axisOpacity = 1;

  const xGroup = svg
    .append("g")
    .attr("transform", `translate(0,${height - marginBottom})`)
    .attr("opacity", axisOpacity)
    .call(xAxis)
    .call((g) => g.select(".domain").remove());
  // Draw grid lines
  xGroup
    .selectAll(".tick line")
    .clone()
    .attr("y2", marginTop + marginBottom - height)
    // Draw line at origin
    .attr("stroke-opacity", (d) => {
      if (d == 0) {
        return axisOpacity;
      } else {
        return showGrid ? 0.1 : 0;
      }
    });

  // Draw yAxis.
  const yGroup = svg
    .append("g")
    .attr("transform", `translate(${marginLeft},0)`)
    .attr("opacity", axisOpacity)
    .call(yAxis)
    .call((g) => g.select(".domain").remove());
  yGroup
    .selectAll(".tick line")
    .clone()
    .attr("x2", width - marginLeft - marginRight)
    // Draw line at origin
    .attr("stroke-opacity", (d) => {
      if (d == 0) {
        return axisOpacity;
      } else {
        return showGrid ? 0.1 : 0;
      }
    });

  // Draw error squares.
  let errorSquares = updateErrorSquares(
    svg.append("g").selectAll("rect"),
    data
  );

  // Draw error lines.
  let errorLines = updateErrorLines(svg.append("g").selectAll("line"), data);

  // Draw regression line.
  const regressionLineGroup = svg
    .append("g")
    .attr("stroke", "black")
    .attr("stroke-width", 2);
  // .attr("stroke-dasharray", "20,20");

  const regressionLine = updateRegressionLine(
    regressionLineGroup.selectAll("line"),
    data
  );

  // Draw space for plot interactions.
  const plotRect = svg
    .append("rect")
    .attr("width", width)
    .attr("height", height)
    .attr("opacity", 0)
    .on("click", addDatapoint);

  // Draw data points.
  const circlesGroup = svg.append("g").attr("fill", "black");
  let circles = updateCircles(circlesGroup.selectAll("circle"), data);

  // Drag interactions for circles.
  function dragstarted(event, d) {
    d3.select(this).raise().attr("stroke", "red");
  }
  function dragged(event, d) {
    d3.select(this)
      // Update data point, as well as its position on the plot
      .attr("cx", () => {
        d.xCoord = xScale.invert(event.x);
        return event.x;
      })
      .attr("cy", () => {
        d.yCoord = yScale.invert(event.y);
        return event.y;
      });

    updateRegressionLine(regressionLine, data, false);
    errorSquares = updateErrorSquares(errorSquares, data, false);
    errorLines = updateErrorLines(errorLines, data, false);
  }
  function dragended(event, i) {
    d3.select(this).attr("stroke", null);
  }

  // Click interaction for circles.
  function removeDatapoint(event, dCurr) {
    if (event.defaultPrevented) return; // dragged

    // Remove data point; faster way to do this?
    data = data.filter((d) => d.id !== dCurr.id);

    update(data);
  }

  // Click interaction for plot.
  function addDatapoint(event) {
    const [xm, ym] = d3.pointer(event);

    updateDatapoint(xScale.invert(xm), yScale.invert(ym));
  }

  // Helper function for regression line update.
  function updateRegressionLine(line, data, transition = true) {
    const regressionData = linearRegression(data);

    // Dispatch event and R^2 value
    svg.node().value = regressionData.rSquared;
    svg.node().transition = transition;
    svg.dispatch("input", { bubbles: true });

    return line
      .data([regressionData])
      .join("line")
      .call((line) => {
        if (transition) line = line.transition();
        line
          .attr("x1", (d) => xScale(d[0][0]))
          .attr("x2", (d) => xScale(d[1][0]))
          .attr("y1", (d) => yScale(d[0][1]))
          .attr("y2", (d) => yScale(d[1][1]));
      });
  }

  // Helper function to update circles based on new data.
  function updateCircles(circles, data) {
    return circles
      .data(data, (d) => d.id)
      .join(
        (enter) =>
          enter
            .append("circle")
            .attr("cx", (d) => xScale(d.xCoord))
            .attr("cy", (d) => yScale(d.yCoord))
            // To transition from 0 radius
            .attr("r", 0)
            // Attach interactions
            .call(drag)
            .on("click", removeDatapoint)
            // Add transition
            .call((enter) =>
              enter
                .transition()
                .ease(d3.easeBackOut.overshoot(1.7))
                .attr("r", r)
            ),
        (update) =>
          update
            .transition()
            .attr("cx", (d) => xScale(d.xCoord))
            .attr("cy", (d) => yScale(d.yCoord)),
        (exit) =>
          exit
            .transition()
            .ease(d3.easeBackIn.overshoot(1.7))
            .attr("r", 0)
            .remove()
      );
  }

  // Helper function to update error lines based on new data.
  function updateErrorLines(lines, data, transition = true) {
    const regressionData = linearRegression(data);

    return lines
      .data(data, (d) => d.id)
      .join(
        (enter) =>
          enter
            .append("line")
            .attr("stroke", errorColor)
            .attr("stroke-width", 2)
            .attr("opacity", errorLinesOpacity)
            // Start at data point
            .attr("x1", (d) => xScale(d.xCoord))
            .attr("y1", (d) => yScale(d.yCoord))
            .attr("x2", (d) => xScale(d.xCoord))
            .attr("y2", (d) => yScale(d.yCoord))
            // Add transition
            .call((enter) =>
              enter
                .transition()
                .attr("y2", (d) => yScale(regressionData.predict(d.xCoord)))
            ),
        (update) =>
          update.call((update) => {
            if (transition) update = update.transition();
            update
              .attr("x1", (d) => xScale(d.xCoord))
              .attr("y1", (d) => yScale(d.yCoord))
              .attr("x2", (d) => xScale(d.xCoord))
              .attr("y2", (d) => yScale(regressionData.predict(d.xCoord)));
          }),
        (exit) =>
          exit
            .transition()
            .attr("y2", (d) => yScale(d.yCoord))
            .remove()
      );
  }

  // Helper function to update error squares based on new data.
  function updateErrorSquares(rects, data, transition = true) {
    const regressionData = linearRegression(data);

    const computeHeight = (d) => {
      const yCoordPred = regressionData.predict(d.xCoord);

      return Math.abs(yScale(yCoordPred) - yScale(d.yCoord));
    };
    // Compute which direction the box should face
    const computeTransform = (d) => {
      const yCoordPred = regressionData.predict(d.xCoord);
      const check = yScale(yCoordPred) - yScale(d.yCoord) < 0;
      const deg = check ? -90 : 90;

      return `rotate(${deg}, ${xScale(d.xCoord)}, ${yScale(d.yCoord)})`;
    };

    return rects
      .data(data, (d) => d.id)
      .join(
        (enter) =>
          enter
            .append("rect")
            .attr("stroke", "none")
            .attr("fill", errorColor)
            .attr("opacity", errorSquaresOpacity)
            .attr("transform", computeTransform)
            // Start at data point
            .attr("x", (d) => xScale(d.xCoord))
            .attr("y", (d) => yScale(d.yCoord))
            .attr("width", 0)
            .attr("height", 0)
            // Add transition
            .call((enter) => {
              enter
                .transition()
                .attr("width", computeHeight)
                .attr("height", computeHeight);
            }),
        (update) =>
          update.call((update) => {
            update.attr("transform", computeTransform);

            if (transition) update = update.transition();

            update
              .attr("x", (d) => xScale(d.xCoord))
              .attr("y", (d) => yScale(d.yCoord))
              .attr("width", computeHeight)
              .attr("height", computeHeight);
          }),
        (exit) => exit.transition().attr("width", 0).attr("height", 0).remove()
      );
  }

  // Resets the plot to the initial data
  function reset() {
    update(initData, true);
  }

  // Updates which error type to show
  function updateView(option) {
    if (option === "None") {
      errorLinesOpacity = 0;
      errorSquaresOpacity = 0;

      errorLines.transition().attr("opacity", 0);
      errorSquares.transition().attr("opacity", 0);
    } else if (option === "Absolute Error") {
      errorLinesOpacity = 1;
      errorSquaresOpacity = 0;

      errorLines.transition().attr("opacity", 1);
      errorSquares.transition().attr("opacity", 0);
    } else if (option === "Squared Error") {
      errorLinesOpacity = 0;
      errorSquaresOpacity = 0.5;

      errorLines.transition().attr("opacity", 0);
      errorSquares.transition().attr("opacity", 0.5);
    }
  }

  // Adds a new datapoint and updates the plot
  function updateDatapoint(xCoord, yCoord) {
    // Add datapoint
    data = [...data, { xCoord, yCoord, id: idCounter++ }];

    update(data);
  }

  // Main function that updates the plot based on new data
  function update(newData, transition = true) {
    // Upate local data object
    data = newData.map((d) => ({ ...d }));

    updateRegressionLine(regressionLine, data, transition);
    circles = updateCircles(circles, data);
    errorSquares = updateErrorSquares(errorSquares, data, transition);
    errorLines = updateErrorLines(errorLines, data, transition);
  }

  return Object.assign(svg.node(), {
    update,
    updateDatapoint,
    updateView,
    reset
  });
}

d3 = require("d3-regression", "d3")

Abbildung 9.5: Interaktives Beispiel für eines lineares Modell. Fügen Sie Punkte per Klick/Touch hinzu.

Beispiel 9.3 (Toni will es genau wissen) Da Toni Sie als Statistik-Profi abgespeichert hat, w erden Sie wieder konsultiert. “Okay, ich hab noch zwei Fragen. Erstens: Wie viele Punkte bekomme ich, wenn ich gar nicht lerne? Zweitens, wie viele Punkte bekomme ich pro gelernte Stunde? Ist immerhin meine Lebenszeit, krieg ich nicht zurück!”

Das sind gute Fragen. Den \(\color{ycol}{Y}\)-Wert (Klausurpunkte) bei \(\color{xcol}{X}=0\) gibt der Achsenabschnitt zurück.

Schnell skizzieren Sie dazu ein Diagramm, s. Abbildung 9.6. Puh, die Antwort wird Toni nicht gefallen …\(\square\)

Abbildung 9.6: Der Achsenabschnitt: Wie viele Punkt kann Toni erwarten bei 0 Lernstunden? (roter Punkt bei x=0)

Anstelle auf Abbildung 9.6 zu schauen, können Sie sich auch von R Tonis Klausurerfolg vorhersagen (to predict) lassen:

🧑‍🏫 Hey R, predicte mir mal auf Basis vom Modell “lm1” den Lernerfolg für Toni, wenn der x=0 Stunden lernt.

🤖 Okay, ich predicte mit Modell “lm1” und nehme als neue Datentabelle Tonis Lernzeit (x=0)!

tonis_lernzeit <- tibble(x = 0)  # `tibble` erstellt eine Tabelle

predict(lm1, newdata = tonis_lernzeit)
##   1 
## 8.6

9.3.2 Spezifikation eines Geradenmodells

Ein Geradenmodell kann man im einfachsten Fall so spezifizieren, s. Gleichung 9.2 :

\[\color{ycol}{\hat{y}} \sim \color{xcol}{\text{x}} \tag{9.1}\]

Lies: “Laut meinem Modell ist mein (geschätztes) \(\color{ycol}{\hat{y}}\) irgendeine Funktion von \(\color{xcol}{\text{x}}\)”.

Wir erinnern uns, dass \(\color{ycol}{Y}\) die \(\color{ycol}{AV}\) und \(\color{xcol}{X}\) die \(\color{xcol}{UV}\) ist:

\[\color{ycol}{AV} \sim \color{xcol}{UV} \tag{9.2}\]

Wir werden als Funktion (erstmal) nur Geraden verwenden. Die genauen Werte der Gerade lassen wir uns (erstmal) vom Computer ausrechnen.

Gleichung 9.2 können Sie so ins Errische übersetzen:

lm(y ~ x, data = meine_daten)

lm steht für “lineares Modell”, also eine Gerade als Modell. Die Gerade nennt man auch Regressionsgerade (an anderer Stelle in diesem Buch unscharf als “Trendgerade” bezeichnet).

Beispiel 9.4 (Zahlen für Toni) Toni ist nicht zufrieden mit Ihren Vorhersagen: “Jetzt hör mal auf mit deinem Lineal hier herum zu malen. Ich will es genau wissen, sag mir präzise Zahlen!”.

lm1 <- lm(y ~ x, data = noten2)
lm1
## 
## Call:
## lm(formula = y ~ x, data = noten2)
## 
## Coefficients:
## (Intercept)            x  
##       8.603        0.879

R gibt Ihnen die beiden Koeffizienten für die Gerade aus. Den Namen des Objekts können Sie frei aussuchen, z.B. mein_erstes_lm.

Die Regressionsgleichung lautet demnach: y_pred = 8.6 + 0.88*x

8.6 ist der Achsenabschnitt, d.h. der Wert von \(\color{ycol}{Y}\) wenn \(\color{xcol}{x}=0\). 0.88 ist das Regressionsgewicht, d.h. die Steigung der Regressionsgeraden: Für jede Stunde Lernzeit steigt der vorhergesagte Klausurerfolg um 0.88 Punkte.

Mit Kenntnis der beiden Koeffizienten kann man beliebige \(\color{ycol}{Y}\)-Werte ausrechnen gegeben bestimmte \(\color{xcol}{X}\)-Werte. Hat jemand zum Beispiel 10 Stunden gelernt, würden wir folgendes Klausurergebnis vorhersagen:

lernzeit <- 10
y_pred <- 8.6 + 0.88*lernzeit
y_pred
## [1] 17

Beispiel 9.5 (Vorhersage für Klausurerfolg, nächster Versuch) Sie versuchen, noch etwas Gutes für Toni zu tun. R hilft Ihnen dabei und rechnet die erwartete Punktzahl aus, wenn Toni 73 Stunden lernt. Sie dürfen es aber auch selber rechnen, wenn Ihnen das lieber ist.

tonis_lernzeit2 <- tibble(x = 73)  # Der Befehl `tibble` erstellt eine Tabelle in R.

tonis_lernzeit2 ist eine Tabelle mit einer Zeile und einer Spalte:

tonis_lernzeit2

predict(lm1, newdata = tonis_lernzeit2)
##  1 
## 73

Die Syntax von predict lautet:

predict(name_des_objekts, newdata = tabelle_mit_prädiktorwerten)

Hinweis

Mit predict bekommt man eine Vorhersage; im Standard eine “Punkt-Vorhersage”, eine einzelne Zahl.\(\square\)

9.3.3 Vorhersagefehler

Die Differenz zwischen vorhergesagten Wert für eine (neue) Beobachtung, \(\color{modelcol}{\hat{y_0}}\) und ihrem tatsächlichen Wert nennt man Vorhersagefehler (error, \(e_i\)) oder Residuum: \(\color{errorcol}{e_i} = \color{ycol}{y_i} - \color{modelcol}{\hat{y}_i}\).

(a) Residuen beim Geradenmodell (lm1)

(b) Residuen beim Punktmodell (lm0)

Abbildung 9.7: Vorhersagefehler als Abweichungsbalken

Wie ist es mit den Vorhersagefehlern von beiden Modellen bestellt?

Lassen wir uns von R die Streuung (Residuen) in Form der mittleren Absolutabweichung (MAE) ausgeben (aus dem Paket easystats):

mae(lm0)
## [1] 11
mae(lm1)
## [1] 8

Vergleichen wir MAE im Nullmodell mit MAE in lm1:

verhaeltnis_fehler_mae <- mae(lm1) / mae(lm0)
verhaeltnis_fehler_mae
## [1] 0.71

Ah! Das Geradenmodell ist viel besser: Von lm0 zu lm1 haben die mittlere (Absolut-)Länge des Fehlerbalkens auf 71 Prozent verbessert. Nicht schlecht!

Definition 9.2 (Fehlerstreuung) Als Fehlerstreuung bezeichnen wir die Gesamtheit der Abweichungen der beobachteten Werte (\(y_i\)) vom vorhergesagten Wert (\(\hat{y}_i\)).\(\square\)

Zur Berechnung der Fehlerstreuung gibt es mehrere Kenngrößen wie MAE oder MSE.

Hinweis

Ein Geradenmodell ist immer besser als ein Punktmodell (im Hinblick auf die Verringerung der Fehlerstreuung), solange X mit Y korreliert ist.\(\square\)

Natürlich können wir - in Analogie zur Varianz - auch den mittleren Quadratfehlerbalken (Mean Squared Error, MSE) berechnen. Wer mag, kann den MSE auch von Hand berechnen: mean((noten2$y-mean(noten2$y))^2).

mse(lm0)
## [1] 193
mse(lm1)
## [1] 106

verhaeltnis_fehler_mse <- mse(lm1)/mse(lm0)
verhaeltnis_fehler_mse
## [1] 0.55

9.3.4 Berechnung der Modellkoeffizienten

Aber wie legt man die Regressionsgerade in das Streudiagramm, bildlich gesprochen?

Die Regressionskoeffizienten (hier synonym: Modellparameter) b0 und b1 wählt man so, dass die Residuen minimal sind,

Abbildung 9.8 veranschaulicht die Minimierung der Residuen (Vorhersagefehler).

Minimierung der Residuen

Berechnung der Modellkoeffizienten durch Minimierung der Residuen

Minimierung der quadrierten Residuen

Minimierung der quadrierten Residuen

Abbildung 9.8: Bildquelle: Karsten Lübke, FOM Hochschule

Genauer gesagt wird die Summe der quadrierten Residuen minimiert, s. Gleichung 9.3.

\[\text{min}\sum_i \color{errorcol}{e_i}^2 \tag{9.3}\]

Es gibt verschiedene Möglichkeiten, um die Koeffizienten zu berechnen (die sind aber nicht in diesem Buch zu finden). Eine schöne Darstellung dazu findet sich bei Kaplan (2009).

“Von Hand” können Sie die Optimierung von b0 und b1 in dieser App der FOM-Hochschule³ ausprobieren.

9.4 R-Quadrat als Maß der Modellgüte

Anders gesagt, wir haben uns um \(1 - 0.55\) verbessert:

1 - verhaeltnis_fehler_mse
## [1] 0.45

Definition 9.3 (R-Quadrat) Die Verringerung (als Anteil) der Fehlerstreuung der Zielvariablen von lm0 zum gerade untersuchten Modell nennt man R-Quadrat (\(R^2\)). R-Quadrat (\(R^2\)) e ines Modells \(m\) ist definiert als die Verringerung der Streuung, wenn man das Modell \(m\) mit dem Nullmodell \(m_0\) vergleicht: \(R^2 =1- \frac{\text{MSE}_{m}}{\text{MSE}_{m0}}\). R-Quadrat ist ein Maß der Modellgüte: Je größer \(R^2\), desto besser die Vorhersage. Da es ein Anteilsmaß ist, liegt der Wertebereich zwischen 0 uns 1. Im Nullmodell liegt R-Quadrat per Definition bei 0. Im Fall von Modellen des Typs \(y\sim x\) gilt: \(R^2 = r_{xy}^2\). \(\square\)

Einfach gesagt: \(R^2\) gibt an, wie gut (zu welchem Anteil) ein Modell die Zielvariable erklärt.

Wir können R-Quadrat (\(R^2\)) uns von R z.B. so ausgeben lassen:

r2(lm1)
## # R2 for Linear Regression
##        R2: 0.448
##   adj. R2: 0.442

Bei einer perfekten Korrelation ist \(r=1\), daher ist dann auch \(R^2 = 1\). Das gilt bei Modellen mit einem Prädiktor; gibt es mehrere Prädiktoren gilt die Beziehung nur, wenn die Prädiktoren alle paarweise unabhängig sind, vgl. Abbildung 9.9.

(a) Keine Korrelation, r ≅ 0 und R2 ≅ 0. Prognose durch Mittelwert; die Regressionsgerade ist (ungefähr) parallel zur X-Achse

(b) Perfekte Korrelation, r = 1 und R2 = 1. Prognose gleich beobachtetem Wert

Abbildung 9.9: Extremfälle von R-Quadrat: 0 und 1

Bei einer perfekten Korrelation \(R^2=1\) liegen die Punkte auf der Geraden. Im gegenteiligen Extremfall von \(R^2=0\) ist die Vorhersage genauso gut, wie wenn man für jedes \(y\) den Mittelwert, \(\color{ycol}{\bar{y}}\) , vorhersagen würde.

Hinweis

Je größer R-Quadrat, desto besser passt das Modell zu den Daten; desto besser “erklärt” das Modell die Daten (desto besser der “Fit”, sagt man).

Diese App der FOM-Hochschule erlaubt es Ihnen mit der Größe der Residuen eines linearen Modells zu spielen.

9.5 Interpretation eines Regressionsmodells

9.5.1 Modellgüte

Die Residuen (Vorhersagefehler) bestimmen die Modellgüte: Sind die Residuen im Schnitt groß, so ist die Modellgüte gering (schlecht), und umgekehrt. Verschiedenen Koeffizienten stehen zur Verfügung: R-Quadrat, r (als Korrelation von tatsächlichem \(y\) und vorhergesagten \(\hat{y}\)), MSE, RMSE, MAE, …

9.5.2 Koeffizienten

Die Modellkoeffizienten, also Achsenabschnitt (\(\beta_0\);lies: “beta Null”) und Steigung (\(\beta_1\)) sind nur eingeschränkt zu interpretieren, wenn man die zugrundeliegenden kausalen Abhängigkeiten nicht kennt. Nur aufgrund eines statistischen Zusammenhangs darf man keine kausalen Abhängigkeiten annehmen. Ohne eine guten Grund für eine Kausalbehauptung kann man kann nur deskriptiv argumentieren. Oder sich mit der Modellgüte und den Vorhersagen begnügen. Was auch was wert ist.

9.5.2.1 Achsenabschnitt (b0)

Im Modell lm1 liegt der Achsenabschnitt bei \(\textcolor{ycol}{y}=8.6\). Beobachtungen mit \(\color{xcol}{x}=0\) können also diesen \(\textcolor{ycol}{Y}\)-Wert erwarten.

Leider ist es häufig so, dass Prädiktorwerte von 0 in der Praxis nicht realistisch sind, so dass der Achsenabschnitt dann wenig nützt.

Beispiel 9.6 (Regression Größe und Gewicht) Nutzt man Körpergröße und das Gewicht von Menschen vorherzusagen, ist der Achsenabschnitt von Körpergröße wenig nützlich, da es keine Menschen gibt der Größe 0.\(\square\)

9.5.2.2 Geradensteigung (b1)

“Im Modell lm1 beträgt der Regressionskoeffizient b1 \(0.88\). Zwei Studenti, deren Lernzeit sich um eine Stunde unterscheidet, unterscheiden sich laut Modell um den Wert von b1.”

Vorsicht

Häufig liest man, der “Effekt des Prädiktors” auf die AV betrage z.B. \(0.88\). “Effekt” ist aber ein Wort, dass man kausal verstehen kann. Ohne weitere Absicherung kann man aber Regressionskoeffizienten nicht kausal verstehen. Daher sollte man das Wort “Effekt” mit Vorsicht genießen. Manche sprechen daher auch von einem “statistischen Effekt”.\(\square\).

9.6 Wie man mit Statistik lügt

Der Unterschied in Modellgüte zwischen, sagen wir, \(r=.1\) und \(r=.2\) ist viel kleiner als zwischen \(r=.7\) und \(r=.8\). \(R^2\) ist ein (lineares) Maß der Modellgüte und da \(r = \sqrt{R^2}\), darf \(r\) nicht wie \(R^2\) als Maß der Modellgüte interpretiert werden. Abbildung 9.10 zeigt den Zusammenhang von \(r\) und \(R^2\).

Abbildung 9.10: Der Zusammenhang von r und R-Quadrat ist nicht linear.

Vorsicht

Unterschiede zwischen Korrelationsdifferenzen dürfen nicht linear interpretiert werden. \(\square\)

9.7 Fallbeispiel Mariokart

9.7.1 Der Datenwahrsager legt los

Als mittlerweile anerkannter Extrem-Datenanalyst in dem Online-Auktionshaus, in dem Sie arbeiten, haben Sie sich neue Ziele gesetzt. Sie möchten eine genaue Vorhersage von Verkaufspreisen erzielen. Als Sie von diesem Plan berichteten, leuchteten die Augen Ihrer Chefin. Genaue Vorhersagen, das ist etwas von hoher betriebswirtschaftlicher Relevanz. Auf geht’s!

Daten laden (und die üblichen Pakete starten, nicht vergessen):

mariokart <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/openintro/mariokart.csv")

lm2 <- lm(total_pr ~ start_pr, data = mariokart)
r2(lm2)
## # R2 for Linear Regression
##        R2: 0.005
##   adj. R2: -0.002

Oh nein! Unterirdisch schlecht. Anstelle von bloßem Rumprobieren überlegen Sie und schauen dann in Abbildung 8.10 nach, welche Variable am stärksten korreliert mit total_pr; es resultiert lm3:

lm3 <- lm(total_pr ~ ship_pr, data = mariokart)
parameters(lm3)

Tabelle 9.1: Modellparameter von lm3

Parameter	Coefficient	SE	95% CI	t(141)	p
(Intercept)	36.25	2.54	(31.23, 41.26)	14.28	< .001
ship pr	4.34	0.57	(3.22, 5.46)	7.67	< .001

Der Achsenabschnitt liegt bei ca. 36 Euro, wie man in Tabelle 9.1 sieht: Ein Spiel, das mit Null Euro Preis startet, kann laut lm3 etwa 36 Euro finaler Verkaufspreis erwarten. Pro Euro an Versandkosten (ship_pr) steigt der zu erwartende finale Verkaufspreis um ca. 4 Euro. (Die Spalte 95 CI gibt einen Schätzbereich für den jeweiligen Modellkoeffizienten an, denn es handelt sich bei den Koeffizienten um Schätzwerte; der wahre Wert in der Population ist unbekannt. Wir kennen schließlich nur eine Stichprobe der Größe \(n=143\).)

Die Regressionsgleichung von lm3 lautet demnach:

total_pr_pred = 36.25 + 4.34*ship_pr.

In Worten:

Der vorhergesagte Gesamptreis eines Spiels liegt bei 36.25€ “Sockelbetrag” plus 4.34 mal die Versandkosten.

9.7.2 Vertiefung

Man kann sich die erwarteten Werte (“expectations”) des Verkaufspreises in Abhängigkeit vom Wert der UV (ship_pr) auch schätzen (“to estimate”) lassen, und zwar so:

estimate_expectation(lm3) %>%  # aus dem Paket 'easystats'
  head()  # nur die ersten paar vorhergesagten Werte

Model-based Expectation

ship_pr	Predicted	SE	95% CI	Residuals
4.00	53.59	1.87	(49.89, 57.30)	-2.04
3.99	53.55	1.87	(49.85, 57.25)	-16.51
3.50	51.43	1.82	(47.82, 55.03)	-5.93
0.00	36.25	2.54	(31.23, 41.26)	7.75
0.00	36.25	2.54	(31.23, 41.26)	34.75
4.00	53.59	1.87	(49.89, 57.30)	-8.59

Variable predicted: total_pr

Ah, bei 4 Euro Versandkosten ist laut dem Modell knapp 54 Euro Verkaufspreis zu erwarten, fassen Sie sich die Ausgabe zusammen.

🤖 Das sieht man in der Spalte Predicted, dort steht der vorhersagte Wert für total_pr für einen bestimmten Wert von ship_pr.

🧑‍🎓 Kann ich auch predict benutzen? Ich würde gerne den Verkaufspreis wissen, wenn die Versandkosten bei 1 und bei 4 Euro liegen.

🤖 Ja, klar!

neue_daten <- tibble(
  ship_pr = c(1, 4)) # zwei Werte zum Vorhersagen

predict(lm3, newdata = neue_daten)
##  1  2 
## 41 54

Aber nützlich wäre noch, das Modell (bzw. die Schätzung der erwarteten Werte) als Diagramm zu bekommen. Das erreicht man z.B. so, s. Abbildung 10.10.

estimate_expectation(lm3) %>% plot()

Abbildung 9.11: Verbildlichung der erwarteteten Werte laut lm3

estimate_expectation heißt sinngemäß “schätze den zu erwartenden Wert”. Kurz gesagt: Wir wollen eine Vorhersage von R.

Am wichtigsten ist Ihnen aber im Moment die Frage, wie “gut” das Modell ist, spricht wie lang oder kurz die (absoluten) Vorhersagefehler-Balken sind:

mae(lm3)
## [1] 13

Das Modell erklärt einen Anteil von ca. 0.29 der Gesamtstreuung.

r2(lm3)
## # R2 for Linear Regression
##        R2: 0.294
##   adj. R2: 0.289

mae(lm3)
## [1] 13

Im nächsten Meeting erzählen Sie Ihrem Chef “Ich kann den Verkaufspreis von Mariokart-Spielen im Schnitt auf 13 Dollar genau vorhersagen!”. Hört sich gut an. Allerdings hätte ihr Chef es gerne genauer. Kann man da noch was machen?

9.8 Fallstudie Immobilienpreise

Vorsicht

Diese Fallstudie stellt die Prüfungsleistung “Prognosewettbewerb” einführend dar. Es empfiehlt sich für Sie, diese Fallstudie sorgsam zu bearbeiten.\(\square\)

9.8.1 Hintergrund

In dieser Fallstudie geht es darum, die Preise von Immobilien vorherzusagen. Kurz gesagt: Sagen Sie die Hauspreise vorher, und reichen Sie Ihre Vorhersagen als CSV bei Kaggle ein.

Kaggle ist eine Webseite, die Prognosewettbewerbe veranstaltet.

In dieser Fallstudie nehmen Sie teil an der Kaggle-Competition Ames House Prices.⁴

• Beschreibung⁵
• Ziel/Aufgabe⁶
• Spielregeln⁷

9.8.2 Benötigte R-Pakete

library(tidyverse)
library(easystats)

9.8.3 Daten

Wenn Sie sich nicht bei Kaggle einloggen möchten, können Sie die Daten von Kaggle herunterladen und zwar hier.

Im Einzelnen müssen Sie folgende Dateien herunterladen:

• Data_description.txt: Code book, d.h. Beschreibung der Variablen im Datensatz
• train.csv: Daten von Häusern, die Sie nutzen, um Modelle zu erstellen
• test.csv: Daten von Häusern, von denen Sie den Kaufpreis vorhersagen sollen
• sample_submission.csv: Beispielhafte Prognosedatei, die Datei also, mit der Sie Ihre Vorhersagen einreichen

Sie können auch so auf die Daten zugreifen:

d_train_path_online <- paste0(
    "https://raw.githubusercontent.com/sebastiansauer/",
    "Lehre/main/data/ames-kaggle/train.csv")

d_test_path_online <- paste0(
  "https://raw.githubusercontent.com/sebastiansauer/",
  "Lehre/main/data/ames-kaggle/test.csv")

d_train <- read.csv(d_train_path_online)
d_test <- read.csv(d_test_path_online)

Laden Sie diese Daten am besten herunter und speichern Sie sie in einem passenden Unterverzeichnis (Ihres Projektordners in RStudio) ab.

Das Code Book können Sie hier einsehen und herunterladen.⁸

9.8.4 Prognosedatei

Die Prognosedatei ist die Datei, die Ihre Vorhersagen (Prognosen) enthält. Sie soll prinzipiell so aussehen wie in Tabelle 9.2 dargestellt.

Tabelle 9.2: Beispiel den Aufbau der Prognose-Datei

id	SalePrice
1461	169277
1462	187758
1463	183584

Die Prognosedatei besteht also aus zwei Spalten: der Spalte id und der Spalte Saleprice. Die Spalte id gibt an, welches Haus in einer bestimmten Zeile Ihrer Prognosedatei gemeint ist - für welches Haus Sie also gerade einen Kaufpreis vorhersagen. die Spalte SalePrice ist Ihre Vorhersage für den Kaufpreis das Hauses mit der Id, die in der betreffenden Zeile steht. Insgesamt soll die Prognosedatei genau so viele Zeilen haben wie der Test-Datensatz, also die Tabelle, die die vorherzusagenden Werte angibt.

Alles klar?

Los geht’s!

9.8.5 Daten importieren von der Festplatte

Wir können die Daten auch von der Festplatte importieren; oft müssen wir das auch - wenn die Daten nämlich nicht öffentlich zugreifbar auf einem Server liegen.

d_train_path <- "daten/ames-kaggle/train.csv"
d_test_path <- "daten/ames-kaggle/test.csv"
d_train <- read.csv(d_train_path)
d_test <- read.csv(d_test_path)

In diesem Beispiel gehen wir davon aus, dass die Dateien train.csv und test.csv in einem Unterordner namens daten/ames-kaggle liegen. Sie müssen sie dort abspeichern. Dieser Ordner muss ein Unterordner Ihres aktuellen R-Projekts sein.\(\square\)

Vorsicht

Wenn das Importieren von der Festplatte nicht klappt … Es ist hilfreich, wenn man Daten von der eigenen Festplatte importieren kann. Aber fürs Erste können Sie die Daten auch von oben angegeben Online-Pfad importieren.\(\square\)

9.8.6 Ein erster Blick in die Daten

Schauen wir uns einmal die Verteilung der metrischen Variablen an, Tabelle 9.3.

describe_distribution(d_train)

Tabelle 9.3: Verteilung der metrischen Variablen im ames-Datensatz

Variable	Mean	SD	IQR	Range	Skewness	Kurtosis	n	n_Missing
Id	730.50	421.61	730.50	(1.00, 1460.00)	0.00	-1.20	1460	0
MSSubClass	56.90	42.30	50.00	(20.00, 190.00)	1.41	1.58	1460	0
LotFrontage	70.05	24.28	21.00	(21.00, 313.00)	2.16	17.45	1201	259
LotArea	10516.83	9981.26	4060.00	(1300.00, 2.15e+05)	12.21	203.24	1460	0
OverallQual	6.10	1.38	2.00	(1.00, 10.00)	0.22	0.10	1460	0
OverallCond	5.58	1.11	1.00	(1.00, 9.00)	0.69	1.11	1460	0
YearBuilt	1971.27	30.20	46.00	(1872.00, 2010.00)	-0.61	-0.44	1460	0
YearRemodAdd	1984.87	20.65	37.00	(1950.00, 2010.00)	-0.50	-1.27	1460	0
MasVnrArea	103.69	181.07	166.00	(0.00, 1600.00)	2.67	10.08	1452	8
BsmtFinSF1	443.64	456.10	712.75	(0.00, 5644.00)	1.69	11.12	1460	0
BsmtFinSF2	46.55	161.32	0.00	(0.00, 1474.00)	4.26	20.11	1460	0
BsmtUnfSF	567.24	441.87	585.00	(0.00, 2336.00)	0.92	0.47	1460	0
TotalBsmtSF	1057.43	438.71	503.50	(0.00, 6110.00)	1.52	13.25	1460	0
X1stFlrSF	1162.63	386.59	509.75	(334.00, 4692.00)	1.38	5.75	1460	0
X2ndFlrSF	346.99	436.53	728.00	(0.00, 2065.00)	0.81	-0.55	1460	0
LowQualFinSF	5.84	48.62	0.00	(0.00, 572.00)	9.01	83.23	1460	0
GrLivArea	1515.46	525.48	649.75	(334.00, 5642.00)	1.37	4.90	1460	0
BsmtFullBath	0.43	0.52	1.00	(0.00, 3.00)	0.60	-0.84	1460	0
BsmtHalfBath	0.06	0.24	0.00	(0.00, 2.00)	4.10	16.40	1460	0
FullBath	1.57	0.55	1.00	(0.00, 3.00)	0.04	-0.86	1460	0
HalfBath	0.38	0.50	1.00	(0.00, 2.00)	0.68	-1.08	1460	0
BedroomAbvGr	2.87	0.82	1.00	(0.00, 8.00)	0.21	2.23	1460	0
KitchenAbvGr	1.05	0.22	0.00	(0.00, 3.00)	4.49	21.53	1460	0
TotRmsAbvGrd	6.52	1.63	2.00	(2.00, 14.00)	0.68	0.88	1460	0
Fireplaces	0.61	0.64	1.00	(0.00, 3.00)	0.65	-0.22	1460	0
GarageYrBlt	1978.51	24.69	41.00	(1900.00, 2010.00)	-0.65	-0.42	1379	81
GarageCars	1.77	0.75	1.00	(0.00, 4.00)	-0.34	0.22	1460	0
GarageArea	472.98	213.80	244.50	(0.00, 1418.00)	0.18	0.92	1460	0
WoodDeckSF	94.24	125.34	168.00	(0.00, 857.00)	1.54	2.99	1460	0
OpenPorchSF	46.66	66.26	68.00	(0.00, 547.00)	2.36	8.49	1460	0
EnclosedPorch	21.95	61.12	0.00	(0.00, 552.00)	3.09	10.43	1460	0
X3SsnPorch	3.41	29.32	0.00	(0.00, 508.00)	10.30	123.66	1460	0
ScreenPorch	15.06	55.76	0.00	(0.00, 480.00)	4.12	18.44	1460	0
PoolArea	2.76	40.18	0.00	(0.00, 738.00)	14.83	223.27	1460	0
MiscVal	43.49	496.12	0.00	(0.00, 15500.00)	24.48	701.00	1460	0
MoSold	6.32	2.70	3.00	(1.00, 12.00)	0.21	-0.40	1460	0
YrSold	2007.82	1.33	2.00	(2006.00, 2010.00)	0.10	-1.19	1460	0
SalePrice	1.81e+05	79442.50	84075.00	(34900.00, 7.55e+05)	1.88	6.54	1460	0

9.8.7 Ein erstes Vorhersagemodell

9.8.7.1 Welche Variablen eignen sich zur Vorhersage?

Eine einfache Antwort auf die Frage, welche Variablen sich zur Vorhersage eignen, ist, die Korrelation aller Prädiktoren mit der abhängigen Variablen⁹ zu berechnen, s. Listing 9.2.

Listing 9.2: Welche Variablen korrelieren stärker als .3?

d_train %>% 
  select(-Id) %>% 
  correlation() %>%  # berechne Korrelationen
  filter(Parameter2 == "SalePrice") %>%   # aber nur, wo die zweite Variable "SalesPrice" ist
  arrange(-abs(r)) %>%   # sortiere absteigend nach der Höhe des Korrelationskoeffizienten r
  filter(abs(r) > .3)  # nur |r| > .3

Tabelle 9.4: Korrelation der Prädiktoren (UV) mit der AV

Aha! Ein Menge Information.¹⁰

Diese Variablen sind einigermaßen stark mit unserer Zielvariablen SalePrice korreliert. Nutzen wir also diese Variablen (oder einige von ihnen) zur Vorhersage.

9.8.7.2 Modell 1

Im ersten Modell gehen wir davon aus, dass der Verkaufspreis im Großen und Ganzen durch den Zustand der Immobilie (OverallQual) vorhergesagt werden kann. Diese Variable ist am stärksten mit der Zielvariable korreliert und ist daher ein guter Kandidat für die Vorhersage.

m1 <- lm(SalePrice ~ OverallQual, data = d_train)
parameters(m1)  # aus easystats

Parameter	Coefficient	SE	95% CI	t(1458)	p
(Intercept)	-96206.08	5756.41	(-1.07e+05, -84914.35)	-16.71	< .001
OverallQual	45435.80	920.43	(43630.29, 47241.31)	49.36	< .001

Wie gut ist das Modell?

rmse(m1)  # aus easystats
## [1] 48589

Im Schnitt liegen wir 4.54^{4} Dollar daneben. Ob das viel oder weniger ist, wird sich im Vergleich mit anderen Modellen zeigen.

R-Quadrat liefert einen anderen Blick auf die Modellgüte:

r2(m1)  # aus easystats
## # R2 for Linear Regression
##        R2: 0.626
##   adj. R2: 0.625

9.8.7.3 Model 2

Berechnen wir als nächstes ein Modell mit mehreren UV, m2.

Hinweis

Mann kann mehrere UV (Prädiktorvariablen) in ein Regressionsmodell aufnehmen. Dazu trennt man sie mit einem Pluszeichen in lm():

mein_modell <- lm(av ~ uv1 + uv2 + ... + uv_n, data = meine_daten)

Dabei ist das Pluszeichen kein arithmetischer Operator, sondern sagt nur “als UV nimm UV1 und UV2 und …”. \(\square\)

m2 <- lm(SalePrice ~ OverallQual + GrLivArea + GarageCars, data = d_train)
parameters(m2)

Tabelle 9.5 zeigt die Koeffizienten von m2.

Tabelle 9.5: Modellparameter von m1

Parameter	Coefficient	SE	95% CI	t(1456)	p
(Intercept)	-98832.49	4842.90	(-1.08e+05, -89332.69)	-20.41	< .001
OverallQual	27104.83	1072.18	(25001.64, 29208.01)	25.28	< .001
GrLivArea	50.67	2.55	(45.67, 55.68)	19.86	< .001
GarageCars	21298.96	1807.06	(17754.23, 24843.69)	11.79	< .001

Wie gut sind die Vorhersagen des Modells m2 für die Daten von d_train?

rmse(m2)
## [1] 40566

Im Schnitt liegen unsere Vorhersagen 2.71^{4} Dollar daneben. Ist das gut?

Betrachten wir noch \(R^2\):

r2(m2)
## # R2 for Linear Regression
##        R2: 0.739
##   adj. R2: 0.739

Hinweis

Ob die Modellgüte (R-Quadrat, RMSE, etc.) “gut” oder “hoch” ist, beantwortet man am besten relativ, also im Vergleich zu anderen Modellen. \(\square\)

9.8.7.4 Nullmodell

Zum Vergleich berechnen wir das maximal einfache Modell: ohne Prädiktoren. Man nennt es das “Nullmodell”. In diesem Modell sagen wir für jedes Haus einfach den mittleren Preis aller Häuser vorher.

m0 <- lm(SalePrice ~ 1, data = d_train)

Wie gut ist die Vorhersage des Nullnomdells?

rmse(m0)
## [1] 79415

Beim Nullmodell liegen wir ca. 80 Tausend Dollar daneben.

Das R-Quadrat der Nullmodells ist per Definition Null:

r2(m0)
## # R2 for Linear Regression
##        R2: 0.000
##   adj. R2: 0.000

9.8.8 Vorhersagen im Test-Datensatz mit m2

Wir haben jetzt unseren Champion, m2. Alle Hoffnung ruht auf diesem Modell. Ob die Vorhersagen im Test-Sample präzise sein werden? Oder himmelweit daneben? Enttäusche uns nicht!

Hier sind die Vorhersagen:

m2_pred <- predict(m2, newdata = d_test)
head(m2_pred)
##      1      2      3      4      5      6 
## 103395 152441 161838 187676 225467 190260

1

predicte anhand der Regressionsgerade von m1 und zwar anhand der Daten aus d_test

2

zeige den “Kopf” der Vorhersagen (m1_pred), d.h. die ersten paar Vorhersagen

Die Vohersagen fügen wir jetzt dem Test-Sample hinzu:

d_test <- 
  d_test %>% 
  mutate(SalePrice = m2_pred)

9.8.9 Einreichen!

9.8.9.1 Wir brauchen zwei Spalten: Id und SalePrice

So, wir haben unsere Vorhersagen! Jetzt reichen wir diese Vorhersagen ein.

Für die Prognosedatei (submission file) zum Einreichen brauchen wir nur die Spalten id und SalePrice:

m2_subm <-
  d_test %>% 
  select(Id, SalePrice)

Kaggle möchte keine fehlenden Werten in den Vorhersagen, also prüfen wir das mal:

m2_subm %>% 
  drop_na() %>%
  nrow()
## [1] 1458

1

Lass alle Zeilen mit NAs (fehlenden Werten in irgendeiner Spalte) fallen, filtere diese Zeilen also raus

2

zähle die Anzahl der Zeilen (die noch verbleiben)

Die Anzahl der Zeilen, die wir hier erhalten, ist gleich zu den Anzahl der Zeilen von d_test. Es gibt also keine fehlenden Werte.

nrow(d_test)
## [1] 1459

9.8.9.2 Hochladen

Diesen Tibble speichern wir als CSV-Datei an geeigneter Stelle ab.¹¹

write_csv(m2_subm, "daten/ames-kaggle/m1-subm.csv")

Und dann laden Sie diese Datei, m1_subm.csv bei Kaggle hoch und hoffen auf einen Hauptgewinn.

Das Modell erzielte einen Score von 0.55521.

9.8.10 Fazit

Diese Fallstudie hat ein einfaches Prognosemodell vorgestellt. Sicherlich gibt es viele Ansätze, dieses Modell zu verbessern.

Hier sind einige Fragen, die Sie sich dazu stellen können:

• Welche Prädiktoren sollte ich in das Modell aufnehmen?
• Wie gehe ich mit fehlenden Werten um?
• Wenn ein Prädiktor schief ist, sollte ich ihn dann log-transformieren?
• Vielleicht sollte man manche Prädiktoren quadrieren?
• Wie gehe ich mit nominalskalierten Variablen um, wenn diese viele Stufen haben?
• …

Viel Spielraum für Ihre Kreativität!

9.9 Aufgaben

Eine Aufgabe, die eine Einführung zum Kaggle-Wettbewerb Ames House Prices bietet¹², finden Sie im Datenwerk.¹³

Die Webseite datenwerk.netlify.app stellt eine Reihe von einschlägigen Übungsaufgaben bereit. Sie können die Suchfunktion der Webseite nutzen, um die Aufgaben mit den folgenden Namen zu suchen:

• Aussagen-einfache-Regr
• interpret-koeff-lm
• korr-als-regr
• Linearitaet1a
• lm1
• mtcars-regr01
• nichtlineare-regr1
• penguins-regr02
• regression1
• regression1b
• Regression3
• Regression4
• Regression5
• Regression6

Schauen Sie sich die Aufgaben beim Datenwerk an, vor allem die Tags regression und lm.

Nicht alle Aufgaben aus dieser Sammlung passen zum Stoff; vielleicht können Sie einige Aufgaben nicht lösen. Ignorieren Sie einfach diese Aufgaben.

Beachten Sie die Hinweise zu den Aufgaben.¹⁴

9.10 Literaturhinweise

Gelman et al. (2021) liefert eine deutlich umfassendere Einführung in die Regressionsanalyse als dieses Kapitel es tut. Eine moderne, R-orientierte Einführung in Statistik inklusive der Regressionsanalyse findet sich bei Cetinkaya-Rundel & Hardin (2021). Ein Klassiker mit viel Aha-Potenzial ist Cohen et al. (2003)

Cetinkaya-Rundel, M., & Hardin, J. (2021). Introduction to Modern Statistics. https://openintro-ims.netlify.app/

Cohen, J., Cohen, P., West, S. G., & Aiken, L. S. (2003). Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences, 3rd Ed (S. xxviii, 703). Lawrence Erlbaum Associates Publishers.

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and Other Stories. Cambridge University Press.

Kaplan, D. T. (2009). Statistical Modeling: A Fresh Approach. CreateSpace. https://dtkaplan.github.io/SM2-bookdown/

Menk. (2014, Juli 29). Linear Regression. https://texample.net/tikz/examples/linear-regression/

---

1. https://raw.githubusercontent.com/sebastiansauer/statistik1/main/daten/noten.csv

2. Die Daten stehen hier zum Download bereit: https://raw.githubusercontent.com/sebastiansauer/statistik1/main/daten/noten2.csv

3. https://fomshinyapps.shinyapps.io/KleinsteQuadrate/

4. https://www.kaggle.com/competitions/house-prices-advanced-regression-techniques/overview

5. https://www.kaggle.com/competitions/house-prices-advanced-regression-techniques/overview/description

6. https://www.kaggle.com/competitions/house-prices-advanced-regression-techniques/overview/evaluation

7. https://www.kaggle.com/competitions/house-prices-advanced-regression-techniques/rules

8. <ttps://github.com/sebastiansauer/Lehre/blob/main/data/ames-kaggle/data_description.txt>

9. die vorherzusagende Variable, auch Ziel- oder Outcome-Variable genannt

10. Wenn Sie Teile der Ausgabe der Tabelle nicht verstehen: Im Zweifel einfach ignorieren. Wenn Sie die R-Syntax nicht verstehen: Führen Sie die Syntax schrittweise aus. Zuerst d_train ausführen und das Ergebnis betrachten. Dann d_train %>% select(-Id) ausführen, wieder die Ausgabe betrachten, usw.

11. Es bietet sich an write_csv zu verwenden, da write.csv automatisch (ungefragt) noch eine Id-Spalte ohne Namen einfügt (mit den Zeilennummern), das mag aber Kaggle nicht. Kaggle erwartet exakt zwei Spalten und zwar mit den Namen Id und SalePrice.

12. https://www.kaggle.com/competitions/house-prices-advanced-regression-techniques/overview

13. https://sebastiansauer.github.io/Datenwerk/posts/ames-kaggle1/ames-kaggle1.html

14. https://sebastiansauer.github.io/Datenwerk/hinweise