R Language
Linjära modeller (regression)

Syntax

• lm (formel, data, delmängd, vikter, na.action, metod = "qr", modell = SANT, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, kontraster = NULL, offset, .. .)

parametrar

Linjär regression på mtcars-datasättet

Den inbyggda mtcars- dataramen innehåller information om 32 bilar, inklusive deras vikt, bränsleeffektivitet (i mil per gallon), hastighet osv. (För att veta mer om datasatsen, använd help(mtcars) ).

Om vi är intresserade av förhållandet mellan bränsleeffektivitet ( mpg ) och vikt ( wt ) kan vi börja plotta dessa variabler med:

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

Tomterna visar ett (linjärt) förhållande !. Om vi sedan vill utföra linjär regression för att bestämma koefficienterna för en linjär modell, skulle vi använda lm funktionen:

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

~ betyder "förklarat av", så formeln mpg ~ wt betyder att vi förutsäger mpg som förklaras av wt. Det mest användbara sättet att visa utgången är med:

summary(fit)

Vilket ger utgången:

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

Detta ger information om:

• den uppskattade lutningen för varje koefficient ( wt och y-intercept), vilket antyder att den bästa passningen av mpg är 37.2851 + (-5.3445) * wt
• P-värdet för varje koefficient, vilket antyder att fånget och vikten förmodligen inte beror på slump
• Övergripande uppskattningar av passform som R ^ 2 och justerad R ^ 2, som visar hur mycket av variationen i mpg förklaras av modellen

Vi kan lägga till en linje till vår första plot som visar den förutsagda mpg :

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

Det är också möjligt att lägga till ekvationen till denna tomt. Först får du koefficienterna med coef . Sedan använder vi paste0 kollapsar vi koefficienterna med lämpliga variabler och +/- , för att bygga ekvationen. Slutligen lägger vi till den till mtext med 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) 

Resultatet är:

Plotta regressionen (bas)

Fortsätt med mtcars exemplet, här är ett enkelt sätt att producera ett diagram över din linjära regression som är potentiellt lämplig för publicering.

Montera först den linjära modellen och

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

Plotta sedan de två variablerna av intresse och lägg till regressionslinjen inom definitionsdomänen:

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

Nästan där! Det sista steget är att lägga till plottet, regressionsekvationen, rsquare såväl som korrelationskoefficienten. Detta görs med 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')

Observera att du kan lägga till alla andra parametrar som RMSE genom att anpassa vektorfunktionen. Föreställ dig att du vill ha en legend med 10 element. Vektordefinitionen skulle vara följande:

rp = vector('expression',10)

och du måste definiera r[1] .... till r[10]

Här är utgången:

viktning

Ibland vill vi att modellen lägger större vikt vid vissa datapunkter eller exempel än andra. Detta är möjligt genom att ange vikten för inmatningsdata medan du lär dig modellen. Det finns i allmänhet två slags scenarier där vi kan använda ojämna vikter jämfört med exemplen:

• Analytiska vikter: återspegla de olika precisionsnivåerna för olika observationer. Om man till exempel analyserar data där varje observation är det genomsnittliga resultatet från ett geografiskt område, är den analytiska vikten proportionell mot den inversa av den uppskattade variationen. Användbart vid hantering av medelvärden i data genom att tillhandahålla en proportionell vikt med antalet observationer. Källa

• Provtagningsvikter (invers sannolikhetsvikter - IPW): en statistisk teknik för att beräkna statistik standardiserad till en annan population än den där data samlades in. Studieutformningar med en olika provtagningspopulation och population av målinferens (målpopulation) är vanliga vid tillämpning. Användbart när man hanterar data som har saknade värden. Källa

Funktionen lm() gör analytisk viktning. För samplingsvikter används survey för att bygga ett undersökningsdesignobjekt och köra svyglm() . Som standard använder survey provvikt. (OBS: lm() och svyglm() med familjegussisk gaussian() ger alla samma poänguppskattningar, eftersom de båda löser för koefficienterna genom att minimera de vägda minsta rutorna. De skiljer sig i hur standardfel beräknas.)

Testdata

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

Analytiska vikter

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

Produktion

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

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

Produktion

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

Kontroll av olinjäritet med polynomregression

Ibland när vi arbetar med linjär regression måste vi kontrollera om icke-linearitet i data. Ett sätt att göra detta är att anpassa en polynomisk modell och kontrollera om den passar uppgifterna bättre än en linjär modell. Det finns andra skäl, såsom teoretiska, som indikerar att passa in i en kvadratisk eller högre ordning, eftersom det tros att variablernas förhållande är i sig polynomiskt.

Låt oss anpassa en kvadratisk modell för mtcars . För en linjär modell se Linjär regression i mtcars-datasättet .

Först skapar vi en spridningsdiagram över variablerna mpg (Miles / gallon), disp (Displacement (cu.in.)) och wt (Weight (1000 lbs)). Förhållandet mellan mpg och disp verkar olinjär.

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

En linjär passning visar att disp inte är signifikant.

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

Sedan, för att få resultatet av en kvadratisk modell, la vi till I(disp^2) . Den nya modellen verkar bättre när man tittar på R^2 och alla variabler är betydande.

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 

Eftersom vi har tre variabler är den monterade modellen en yta som representeras av:

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

Ett annat sätt att specificera polynomregression är att använda poly med parameter raw=TRUE , annars kommer ortogonala polynom att övervägas (se help(ploy) för mer information). Vi får samma resultat med:

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

Slutligen, vad händer om vi måste visa ett diagram över den uppskattade ytan? Det finns väl många alternativ att göra 3D tomter i R Här använder vi Fit3d från p3d paketet.

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

Kvalitetsbedömning

Efter att ha byggt en regressionsmodell är det viktigt att kontrollera resultatet och bestämma om modellen är lämplig och fungerar bra med uppgifterna. Detta kan göras genom att undersöka restplottet såväl som andra diagnostiska diagram.

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

Dessa tomter kontrollerar för två antaganden som gjordes medan modellen byggdes:

1. Att det förväntade värdet på den förutsagda variabeln (i detta fall mpg ) ges av en linjär kombination av prediktorerna (i detta fall wt ). Vi räknar med att denna uppskattning är opartisk. Så resterna bör centreras kring medelvärdet för alla värden på prediktorerna. I detta fall ser vi att resterna tenderar att vara positiva i ändarna och negativa i mitten, vilket antyder ett icke-linjärt samband mellan variablerna.
2. Att den faktiska förutsagda variabeln normalt fördelas runt dess uppskattning. Således bör resterna normalt fördelas. För normalt distribuerad data bör punkterna i ett normalt QQ-diagram ligga på eller nära diagonalen. Det finns en del skev i ändarna här.

Använd "förutsäga" -funktionen

När en modell har byggts är predict den viktigaste funktionen att testa med nya data. Vårt exempel kommer att använda de inbyggda mtcars att regressera miles per gallon mot förskjutning:

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

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

Coefficients:
(Intercept)         disp  
   29.59985     -0.04122

Om jag hade en ny datakälla med förskjutning kunde jag se de uppskattade milerna per gallon.

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

Den viktigaste delen av processen är att skapa en ny dataram med samma kolumnnamn som originaldata. I det här fallet hade de ursprungliga uppgifterna en kolumn märkt disp , jag var säker på att kalla den nya informationen samma namn.

Varning

Låt oss titta på några vanliga fallgropar:

1. använder inte ett data.frame i det nya objektet:

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

2. använder inte samma namn i ny dataram:

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

Noggrannhet

För att kontrollera förutsägelsens noggrannhet behöver du de verkliga y-värdena för den nya datan. I det här exemplet newdf en kolumn för 'mpg' och '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