R Language
Lineare Modelle (Regression)

Syntax

• lm (Formel, Daten, Teilmenge, Gewichte, na.Aktion, Methode = "qr", Modell = TRUE, x = FALSCH, y = FALSE, qr = TRUE, Singular.ok = TRUE, Kontraste = NULL, Offset, .. .)

Parameter

Lineare Regression im mtcars-Dataset

Der eingebaute in mtcars Datenrahmen enthält Informationen über 32 Autos, einschließlich ihres Gewichts, ihrer Kraftstoffeffizienz (in Meilen pro Gallone), Geschwindigkeit usw. (Um mehr über den Datensatz, Verwendung finden heraus help(mtcars) ).

Wenn wir an der Beziehung zwischen Kraftstoffeffizienz ( mpg ) und Gewicht ( wt ) interessiert sind, beginnen wir diese Variablen mit:

plot(mpg ~ wt, data = mtcars, col=2)

Die Diagramme zeigen eine (lineare) Beziehung! Wenn wir dann die lineare Regression durchführen wollen, um die Koeffizienten eines linearen Modells zu bestimmen, würden wir die Funktion lm verwenden:

fit <- lm(mpg ~ wt, data = mtcars)

Das ~ bedeutet hier "erklärt durch", also bedeutet die Formel mpg ~ wt dass wir mpg voraussagen, wie es durch wt erklärt wird. Der hilfreichste Weg zum Anzeigen der Ausgabe ist mit:

summary(fit)

Was gibt die Ausgabe:

Call:
lm(formula = mpg ~ wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5432 -2.3647 -0.1252  1.4096  6.8727 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
wt           -5.3445     0.5591  -9.559 1.29e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared:  0.7528,    Adjusted R-squared:  0.7446 
F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

Hier finden Sie Informationen zu:

• Die geschätzte Steigung jedes Koeffizienten ( wt und y- 37.2851 + (-5.3445) * wt ), die auf die best-fit-Vorhersage von mpg schließen lässt, beträgt 37.2851 + (-5.3445) * wt
• Der p-Wert jedes Koeffizienten legt nahe, dass der Achsenabschnitt und das Gewicht wahrscheinlich nicht zufällig sind
• Gesamteinschätzungen der Anpassung wie R ^ 2 und bereinigtes R ^ 2, die zeigen, wie viel von der Variation von mpg durch das Modell erklärt wird

Wir könnten unserem ersten Plot eine Linie hinzufügen, um das vorhergesagte mpg :

abline(fit,col=3,lwd=2)

Es ist auch möglich, die Gleichung zu diesem Diagramm hinzuzufügen. coef zuerst die Koeffizienten mit coef . Mit paste0 wir die Koeffizienten mit geeigneten Variablen und +/- , um die Gleichung paste0 . Zum Schluss fügen wir es mit mtext zum Plot:

bs <- round(coef(fit), 3) 
lmlab <- paste0("mpg = ", bs[1],
             ifelse(sign(bs[2])==1, " + ", " - "), abs(bs[2]), " wt ")
mtext(lmlab, 3, line=-2) 

Das Ergebnis ist:

Die Regression plotten (Basis)

mtcars dem mtcars Beispiel fort. Hier können Sie auf einfache Weise eine mtcars Darstellung Ihrer linearen Regression erstellen, die möglicherweise für die Veröffentlichung geeignet ist.

Passen Sie zuerst das lineare Modell und an

fit <- lm(mpg ~ wt, data = mtcars)

Plotten Sie dann die beiden interessierenden Variablen und fügen Sie die Regressionslinie innerhalb der Definitionsdomäne hinzu:

plot(mtcars$wt,mtcars$mpg,pch=18, xlab = 'wt',ylab = 'mpg')
lines(c(min(mtcars$wt),max(mtcars$wt)),
as.numeric(predict(fit, data.frame(wt=c(min(mtcars$wt),max(mtcars$wt))))))

Fast dort! Der letzte Schritt ist das Hinzufügen der Regressionsgleichung, des Rechtecks ​​sowie des Korrelationskoeffizienten. Dies geschieht mit der vector :

rp = vector('expression',3)
rp[1] = substitute(expression(italic(y) == MYOTHERVALUE3 + MYOTHERVALUE4 %*% x), 
          list(MYOTHERVALUE3 = format(fit$coefficients[1], digits = 2),
                        MYOTHERVALUE4 = format(fit$coefficients[2], digits = 2)))[2]
rp[2] = substitute(expression(italic(R)^2 == MYVALUE), 
             list(MYVALUE = format(summary(fit)$adj.r.squared,dig=3)))[2]
rp[3] = substitute(expression(Pearson-R == MYOTHERVALUE2), 
             list(MYOTHERVALUE2 = format(cor(mtcars$wt,mtcars$mpg), digits = 2)))[2]

legend("topright", legend = rp, bty = 'n')

Beachten Sie, dass Sie durch Anpassen der Vektorfunktion beliebige andere Parameter wie RMSE hinzufügen können. Stellen Sie sich vor, Sie möchten eine Legende mit 10 Elementen. Die Vektordefinition wäre die folgende:

rp = vector('expression',10)

und Sie müssen r[1] .... bis r[10]

Hier ist die Ausgabe:

Gewichtung

Manchmal möchten wir, dass das Modell einigen Datenpunkten oder Beispielen mehr Gewicht verleiht als anderen. Dies ist möglich, indem Sie das Gewicht für die Eingabedaten angeben, während Sie das Modell lernen. Es gibt im Allgemeinen zwei Arten von Szenarien, in denen wir in den Beispielen möglicherweise nicht einheitliche Gewichte verwenden:

• Analytische Gewichte: Reflektieren Sie die unterschiedlichen Genauigkeitsgrade verschiedener Beobachtungen. Wenn zum Beispiel Daten analysiert werden, bei denen jede Beobachtung die Durchschnittsergebnisse aus einem geografischen Gebiet ist, ist das Analysegewicht proportional zum Inversen der geschätzten Varianz. Nützlich beim Umgang mit Durchschnittswerten in Daten durch Angabe eines proportionalen Gewichts aufgrund der Anzahl der Beobachtungen. Quelle

• Stichprobengewichte (inverse Wahrscheinlichkeitsgewichte - IPW): Eine statistische Methode zur Berechnung von Statistiken, die auf eine Grundgesamtheit normiert sind, die sich von derjenigen unterscheidet, in der die Daten erfasst wurden. Studiendesigns mit einer unterschiedlichen Stichprobenpopulation und einer Zielpopulation (Zielpopulation) sind häufig in der Anwendung. Nützlich bei Daten mit fehlenden Werten. Quelle

Die Funktion lm() führt eine analytische Gewichtung durch. Für Stichprobengewichte wird das survey verwendet, um ein Vermessungsobjekt zu erstellen und svyglm() . Standardmäßig verwendet das survey Stichprobengewichte. (ANMERKUNG: lm() und svyglm() mit der Familie gaussian() . Werden alle die gleichen Punktschätzungen zu produzieren, weil sie sowohl für die Koeffizienten zu lösen , indem die gewichteten kleinsten Quadrate zu minimieren Sie unterscheiden sich darin , wie Standardfehler berechnet.)

Testdaten

data <- structure(list(lexptot = c(9.1595012302023, 9.86330744180814, 
8.92372556833205, 8.58202430280175, 10.1133857229336), progvillm = c(1L, 
1L, 1L, 1L, 0L), sexhead = c(1L, 1L, 0L, 1L, 1L), agehead = c(79L, 
43L, 52L, 48L, 35L), weight = c(1.04273509979248, 1.01139605045319, 
1.01139605045319, 1.01139605045319, 0.76305216550827)), .Names = c("lexptot", 
"progvillm", "sexhead", "agehead", "weight"), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -5L))

Analytische Gewichte

lm.analytic <- lm(lexptot ~ progvillm + sexhead + agehead, 
                            data = data, weight = weight)
summary(lm.analytic)

Ausgabe

Call:
lm(formula = lexptot ~ progvillm + sexhead + agehead, data = data, 
    weights = weight)

Weighted Residuals:
         1          2          3          4          5 
 9.249e-02  5.823e-01  0.000e+00 -6.762e-01 -1.527e-16

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.016054   1.744293   5.742    0.110
progvillm   -0.781204   1.344974  -0.581    0.665
sexhead      0.306742   1.040625   0.295    0.818
agehead     -0.005983   0.032024  -0.187    0.882

Residual standard error: 0.8971 on 1 degrees of freedom
Multiple R-squared:  0.467, Adjusted R-squared:  -1.132 
F-statistic: 0.2921 on 3 and 1 DF,  p-value: 0.8386

Probenahmegewichte (IPW)

library(survey)
data$X <- 1:nrow(data)             # Create unique id
        
# Build survey design object with unique id, ipw, and data.frame
des1 <- svydesign(id = ~X,  weights = ~weight, data = data)
        
# Run glm with survey design object
prog.lm <- svyglm(lexptot ~ progvillm + sexhead + agehead, design=des1)

Ausgabe

Call:
svyglm(formula = lexptot ~ progvillm + sexhead + agehead, design = des1)

Survey design:
svydesign(id = ~X, weights = ~weight, data = data)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 10.016054   0.183942  54.452   0.0117 *
progvillm   -0.781204   0.640372  -1.220   0.4371  
sexhead      0.306742   0.397089   0.772   0.5813  
agehead     -0.005983   0.014747  -0.406   0.7546  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.2078647)

Number of Fisher Scoring iterations: 2

Prüfung auf Nichtlinearität mit polynomialer Regression

Wenn Sie mit linearer Regression arbeiten, müssen wir manchmal die Nichtlinearität der Daten überprüfen. Eine Möglichkeit, dies zu tun, besteht darin, ein Polynommodell anzupassen und zu prüfen, ob es besser zu den Daten passt als ein lineares Modell. Es gibt andere Gründe, z. B. theoretische Gründe, die auf ein quadratisches oder ein Modell höherer Ordnung schließen lassen, da davon ausgegangen wird, dass die Variablenbeziehung von Natur aus polynomisch ist.

Passen Sie ein quadratisches Modell für das mtcars Dataset an. Für ein lineares Modell siehe Lineare Regression im mtcars-Dataset .

Zuerst disp wir ein Streudiagramm der Variablen mpg (Meilen / Gallone), disp (Displacement (cu.in.)) und wt (Weight (1000 lbs)). Die Beziehung zwischen mpg und disp erscheint nicht linear.

plot(mtcars[,c("mpg","disp","wt")])

Eine lineare Anpassung zeigt, dass disp nicht signifikant ist.

fit0 = lm(mpg ~ wt+disp, mtcars)
summary(fit0)

# Coefficients:
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 34.96055    2.16454  16.151 4.91e-16 ***
#wt          -3.35082    1.16413  -2.878  0.00743 ** 
#disp        -0.01773    0.00919  -1.929  0.06362 .  
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 2.917 on 29 degrees of freedom
#Multiple R-squared:  0.7809,    Adjusted R-squared:  0.7658

Um das Ergebnis eines quadratischen Modells zu erhalten, haben wir I(disp^2) hinzugefügt. Das neue Modell erscheint bei Betrachtung von R^2 besser und alle Variablen sind signifikant.

fit1 = lm(mpg ~ wt+disp+I(disp^2), mtcars)
summary(fit1)

# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 41.4019837  2.4266906  17.061  2.5e-16 ***
#wt          -3.4179165  0.9545642  -3.581 0.001278 ** 
#disp        -0.0823950  0.0182460  -4.516 0.000104 ***
#I(disp^2)    0.0001277  0.0000328   3.892 0.000561 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 2.391 on 28 degrees of freedom
#Multiple R-squared:  0.8578,    Adjusted R-squared:  0.8426 

Da wir drei Variablen haben, ist das angepasste Modell eine Fläche, die dargestellt wird durch:

mpg = 41.4020-3.4179*wt-0.0824*disp+0.0001277*disp^2

Eine andere Möglichkeit, die Regression von Polynomen anzugeben, ist die Verwendung von poly mit dem Parameter raw=TRUE . Andernfalls werden orthogonale Polynome berücksichtigt (weitere Informationen finden Sie in der help(ploy) ). Wir erhalten das gleiche Ergebnis mit:

summary(lm(mpg ~ wt+poly(disp, 2, raw=TRUE),mtcars))

Was ist schließlich, wenn wir eine Darstellung der geschätzten Oberfläche anzeigen müssen? Es gibt viele Möglichkeiten, 3D Darstellungen in R erstellen. Hier verwenden wir Fit3d aus dem p3d Paket.

library(p3d)
Init3d(family="serif", cex = 1)
Plot3d(mpg ~ disp+wt, mtcars)
Axes3d()
Fit3d(fit1)

Qualitätsprüfung

Nach dem Erstellen eines Regressionsmodells ist es wichtig, das Ergebnis zu überprüfen und zu entscheiden, ob das Modell geeignet ist und mit den vorliegenden Daten gut funktioniert. Dies kann durch Untersuchen des Residuendiagramms sowie anderer Diagnoseplots erfolgen.

# fit the model
fit <- lm(mpg ~ wt, data = mtcars)
# 
par(mfrow=c(2,1))
# plot model object
plot(fit, which =1:2)

Diese Diagramme überprüfen zwei Annahmen, die beim Erstellen des Modells getroffen wurden:

1. Dass der erwartete Wert der vorhergesagten Variablen (in diesem Fall mpg ) durch eine Linearkombination der Prädiktoren (in diesem Fall wt ) gegeben ist. Wir erwarten, dass diese Einschätzung unvoreingenommen ist. Daher sollten die Residuen um den Mittelwert für alle Werte der Prädiktoren zentriert sein. In diesem Fall sehen wir, dass die Residuen an den Enden tendenziell positiv und in der Mitte negativ sind, was auf eine nichtlineare Beziehung zwischen den Variablen schließen lässt.
2. Die tatsächliche vorhergesagte Variable ist normalerweise um ihre Schätzung verteilt. Daher sollten die Residuen normal verteilt sein. Bei normalverteilten Daten sollten die Punkte in einem normalen QQ-Diagramm auf oder nahe der Diagonale liegen. An den Enden ist etwas Schräglauf vorhanden.

Verwenden der Funktion "Vorhersagen"

Sobald ein Modell erstellt ist predict ist predict die Hauptfunktion zum Testen mit neuen Daten. In unserem Beispiel wird der eingebaute Datensatz von mtcars , um Meilen pro Gallone gegen Verschiebung zu mtcars :

my_mdl <- lm(mpg ~ disp, data=mtcars)
my_mdl

Call:
lm(formula = mpg ~ disp, data = mtcars)

Coefficients:
(Intercept)         disp  
   29.59985     -0.04122

Wenn ich eine neue Datenquelle mit Verschiebung hätte, könnte ich die geschätzten Meilen pro Gallone sehen.

set.seed(1234)
newdata <- sample(mtcars$disp, 5)
newdata
[1] 258.0  71.1  75.7 145.0 400.0

newdf <- data.frame(disp=newdata)
predict(my_mdl, newdf)
       1        2        3        4        5 
18.96635 26.66946 26.47987 23.62366 13.11381

Der wichtigste Teil des Prozesses ist das Erstellen eines neuen Datenrahmens mit den gleichen Spaltennamen wie die Originaldaten. In diesem Fall hatten die ursprünglichen Daten eine Spalte mit der Bezeichnung disp . Ich war sicher, die neuen Daten denselben Namen zu nennen.

Vorsicht

Schauen wir uns ein paar häufige Fallstricke an:

1. kein data.frame im neuen Objekt verwenden:

   predict(my_mdl, newdata)
   Error in eval(predvars, data, env) : 
      numeric 'envir' arg not of length one
   

2. In neuen Datenrahmen nicht dieselben Namen verwenden:

   newdf2 <- data.frame(newdata)
   predict(my_mdl, newdf2)
   Error in eval(expr, envir, enclos) : object 'disp' not found
   

Richtigkeit

Um die Genauigkeit der Vorhersage zu überprüfen, benötigen Sie die tatsächlichen y-Werte der neuen Daten. In diesem Beispiel benötigt newdf eine Spalte für 'mpg' und 'disp'.

newdf <- data.frame(mpg=mtcars$mpg[1:10], disp=mtcars$disp[1:10])
#     mpg  disp
# 1  21.0 160.0
# 2  21.0 160.0
# 3  22.8 108.0
# 4  21.4 258.0
# 5  18.7 360.0
# 6  18.1 225.0
# 7  14.3 360.0
# 8  24.4 146.7
# 9  22.8 140.8
# 10 19.2 167.6

p <- predict(my_mdl, newdf)

#root mean square error
sqrt(mean((p - newdf$mpg)^2, na.rm=TRUE))
[1] 2.325148