R Language
Lineaire modellen (regressie)

Syntaxis

• lm (formule, gegevens, subset, gewichten, na.actie, methode = "qr", model = WAAR, x = ONWAAR, y = ONWAAR, qr = WAAR, enkelvoud.ok = WAAR, contrasten = NULL, offset, .. .)

parameters

Lineaire regressie op de dataset mtcars

Het ingebouwde mtcars -dataframe bevat informatie over 32 auto's, inclusief hun gewicht, brandstofefficiëntie (in mijl per gallon), snelheid, etc. (gebruik help(mtcars) meer informatie over de dataset).

Als we geïnteresseerd zijn in de relatie tussen brandstofefficiëntie ( mpg ) en gewicht ( wt ), kunnen we die variabelen gaan plotten met:

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

De plots tonen een (lineaire) relatie !. Als we dan lineaire regressie willen uitvoeren om de coëfficiënten van een lineair model te bepalen, zouden we de functie lm :

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

De ~ hier betekent "verklaard door", dus de formule mpg ~ wt betekent dat we mpg voorspellen zoals uitgelegd door wt. De handigste manier om de uitvoer te bekijken is met:

summary(fit)

Welke geeft de output:

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

Dit geeft informatie over:

• de geschatte helling van elke coëfficiënt ( wt en de y-intercept), wat suggereert dat de best passende voorspelling van mpg 37.2851 + (-5.3445) * wt
• De p-waarde van elke coëfficiënt, wat suggereert dat het onderscheppen en gewicht waarschijnlijk niet te wijten zijn aan toeval
• Algemene schattingen van fit zoals R ^ 2 en aangepaste R ^ 2, die laten zien hoeveel van de variatie in mpg wordt verklaard door het model

We kunnen een regel toevoegen aan onze eerste plot om de voorspelde mpg :

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

Het is ook mogelijk om de vergelijking aan die plot toe te voegen. Krijg eerst de coëfficiënten met coef . Vervolgens paste0 we met behulp van paste0 de coëfficiënten samen met de juiste variabelen en +/- om de vergelijking te bouwen. Ten slotte voegen we het toe aan de plot met behulp van mtext :

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

Het resultaat is:

De regressie plotten (basis)

Doorgaan met het voorbeeld van mtcars , hier is een eenvoudige manier om een plot van uw lineaire regressie te produceren die mogelijk geschikt is voor publicatie.

Monteer eerst het lineaire model en

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

Teken vervolgens de twee relevante variabelen en voeg de regressielijn toe binnen het definitiedomein:

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

Ik ben er bijna! De laatste stap is het toevoegen van de plot, de regressievergelijking, het rsquare en de correlatiecoëfficiënt. Dit wordt gedaan met behulp van de 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')

Merk op dat u elke andere parameter zoals de RMSE kunt toevoegen door de vectorfunctie aan te passen. Stel je voor dat je een legende met 10 elementen wilt. De vectordefinitie zou de volgende zijn:

rp = vector('expression',10)

en je moet r[1] definiëren .... tot r[10]

Hier is de output:

weging

Soms willen we dat het model sommige datapunten of voorbeelden meer gewicht geeft dan andere. Dit is mogelijk door het gewicht voor de invoergegevens op te geven terwijl u het model leert. Er zijn over het algemeen twee soorten scenario's waarbij we niet-uniforme gewichten kunnen gebruiken boven de voorbeelden:

• Analytische gewichten: reflecteer de verschillende nauwkeurigheidsniveaus van verschillende waarnemingen. Als bijvoorbeeld het analyseren van gegevens waarbij elke waarneming het gemiddelde is, het resultaat is van een geografisch gebied, is het analytische gewicht evenredig met het omgekeerde van de geschatte variantie. Handig bij het omgaan met gemiddelden in gegevens door een proportioneel gewicht te geven, gegeven het aantal waarnemingen. Bron

• Bemonsteringsgewichten (Inverse Probability Weights - IPW): een statistische techniek voor het berekenen van statistieken die zijn gestandaardiseerd voor een populatie die verschilt van die waarin de gegevens zijn verzameld. Studieontwerpen met een verschillende steekproefpopulatie en populatie van doelafleiding (doelpopulatie) komen veel voor in de toepassing. Handig bij het omgaan met gegevens die waarden missen. Bron

De functie lm() voert analytische weging uit. Voor steekproefgewichten wordt het survey gebruikt om een enquêteontwerpobject te bouwen en svyglm() . Standaard gebruikt het survey steekproefgewichten. (OPMERKING: lm() en svyglm() met family gaussian() produceren allemaal dezelfde puntschattingen, omdat ze allebei oplossen voor de coëfficiënten door de gewogen kleinste kwadraten te minimaliseren. Ze verschillen in hoe standaardfouten worden berekend.)

Testgegevens

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 gewichten

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

uitgang

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

Bemonsteringsgewichten (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)

uitgang

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

Controleren op niet-lineariteit met polynoomregressie

Soms moeten we bij het werken met lineaire regressie controleren op niet-lineariteit in de gegevens. Een manier om dit te doen is om een polynoommodel te passen en te controleren of het beter bij de gegevens past dan een lineair model. Er zijn andere redenen, zoals theoretisch, die aangeven dat ze passen in een kwadratisch of hoger orde model omdat wordt aangenomen dat de variabelenrelatie inherent polynoom van aard is.

Laten we een kwadratisch model voor de mtcars gegevensset passen. Zie Lineaire regressie op de mtcars-gegevensset voor een lineair model.

Eerst maken we een spreidingsdiagram van de variabelen mpg (Miles / gallon), disp (Verplaatsing (cu.in.)) en wt (Gewicht (1000 lbs)). De relatie tussen mpg en disp lijkt niet-lineair.

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

Een lineaire aanpassing zal aantonen dat disp niet significant is.

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

Vervolgens hebben we I(disp^2) toegevoegd om het resultaat van een kwadratisch model te krijgen. Het nieuwe model lijkt beter wanneer we kijken naar R^2 en alle variabelen zijn significant.

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 

Omdat we drie variabelen hebben, is het gepaste model een oppervlak dat wordt voorgesteld door:

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

Een andere manier om polynoomregressie op te geven, is poly met de parameter raw=TRUE , anders worden orthogonale polynomen in overweging genomen (zie de help(ploy) voor meer informatie). We krijgen hetzelfde resultaat met:

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

Tot slot, wat als we een plot van het geschatte oppervlak moeten laten zien? Welnu, er zijn veel opties om 3D plots in R . Hier gebruiken we Fit3d van p3d pakket.

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

Kwaliteitsbeoordeling

Na het bouwen van een regressiemodel is het belangrijk om het resultaat te controleren en te beslissen of het model geschikt is en goed werkt met de beschikbare gegevens. Dit kan worden gedaan door de residualplot en andere diagnostische plots te onderzoeken.

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

Deze plots controleren op twee veronderstellingen die zijn gemaakt tijdens het bouwen van het model:

1. Dat de verwachte waarde van de voorspelde variabele (in dit geval mpg ) wordt gegeven door een lineaire combinatie van de voorspellers (in dit geval wt ). We verwachten dat deze schatting onbevooroordeeld is. De residuen moeten dus gecentreerd zijn rond het gemiddelde voor alle waarden van de voorspellers. In dit geval zien we dat de residuen de neiging hebben om positief te zijn aan het uiteinde en negatief in het midden, wat een niet-lineair verband tussen de variabelen suggereert.
2. Dat de werkelijke voorspelde variabele normaal rond zijn schatting wordt verdeeld. De residuen moeten dus normaal worden verdeeld. Voor normaal verdeelde gegevens moeten de punten in een normale QQ-plot op of dichtbij de diagonaal liggen. Er zit een beetje scheef aan de uiteinden hier.

Gebruik van de functie 'voorspellen'

Zodra een model is gebouwd, is predict de belangrijkste functie om met nieuwe gegevens te testen. Ons voorbeeld zal de ingebouwde dataset van mtcars gebruiken om mijlen per gallon tegen verplaatsing te mtcars :

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

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

Coefficients:
(Intercept)         disp  
   29.59985     -0.04122

Als ik een nieuwe gegevensbron met verplaatsing had, kon ik de geschatte mijlen per gallon zien.

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

Het belangrijkste onderdeel van het proces is het maken van een nieuw gegevensframe met dezelfde kolomnamen als de oorspronkelijke gegevens. In dit geval hadden de oorspronkelijke gegevens een kolom met het label disp , ik was er zeker van dat ik de nieuwe gegevens dezelfde naam zou geven.

Voorzichtigheid

Laten we een paar veel voorkomende valkuilen bekijken:

1. geen data.frame gebruiken in het nieuwe object:

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

2. geen gebruik maken van dezelfde namen in nieuw gegevensframe:

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

Nauwkeurigheid

Om de nauwkeurigheid van de voorspelling te controleren, hebt u de werkelijke y-waarden van de nieuwe gegevens nodig. In dit voorbeeld heeft newdf een kolom nodig voor 'mpg' en '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