R Language
रैखिक मॉडल (प्रतिगमन)

वाक्य - विन्यास

• lm (सूत्र, डेटा, सबसेट, वज़न, na.action, विधि = "qr", मॉडल = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, इसके विपरीत = NULL, ऑफसेट, ।। ।)

पैरामीटर

Mtcars डेटासेट पर रैखिक प्रतिगमन

अंतर्निहित mtcars डेटा फ़्रेम में 32 कारों के बारे में जानकारी होती है, जिसमें उनके वजन, ईंधन दक्षता (मील-प्रति-गैलन में), गति आदि शामिल हैं (डेटासेट के बारे में अधिक जानने के लिए, help(mtcars) उपयोग करें)।

यदि हम ईंधन दक्षता ( mpg ) और वजन ( wt ) के बीच संबंधों में रुचि रखते हैं, तो हम उन चरों के साथ साजिश शुरू कर सकते हैं:

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

भूखंड (रैखिक) संबंध दिखाता है! फिर अगर हम एक रैखिक मॉडल के गुणांक निर्धारित करने के लिए रैखिक प्रतिगमन करना चाहते हैं, तो हम lm फ़ंक्शन का उपयोग करेंगे:

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

~ का अर्थ "द्वारा समझाया गया" है, इसलिए सूत्र mpg ~ wt अर्थ है कि हम mpg की भविष्यवाणी कर रहे हैं जैसा कि wt द्वारा समझाया गया है। आउटपुट को देखने का सबसे सहायक तरीका निम्न है:

summary(fit)

जो आउटपुट देता है:

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

इसके बारे में जानकारी प्रदान करता है:

• प्रत्येक गुणांक ( wt और y- इंटरसेप्ट) की अनुमानित ढलान, जो mpg का सबसे अच्छा फिट होने का अनुमान 37.2851 + (-5.3445) * wt
• प्रत्येक गुणांक का पी-मूल्य, जो बताता है कि अवरोधन और भार संभवतः संयोग के कारण नहीं हैं
• फिट के समग्र अनुमान जैसे कि R ^ 2 और समायोजित R ^ 2, जो बताते हैं कि mpg में कितनी भिन्नता मॉडल द्वारा बताई गई है

हम पूर्वानुमानित mpg दिखाने के लिए अपने पहले प्लॉट में एक पंक्ति जोड़ सकते हैं:

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

उस भूखंड के समीकरण को जोड़ना भी संभव है। सबसे पहले, के साथ गुणांक पाने coef । तब paste0 का उपयोग paste0 हम समीकरण बनाने के लिए उपयुक्त चर और +/- साथ गुणांक को ढहते हैं। अंत में, हम इसे mtext का उपयोग करके प्लॉट में 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) 

परिणाम है:

प्रतिगमन (आधार) प्लॉटिंग

mtcars उदाहरण पर जारी रखते हुए, यहां अपने रैखिक प्रतिगमन के एक भूखंड का उत्पादन करने का एक सरल तरीका है जो प्रकाशन के लिए संभवतः उपयुक्त है।

पहले लीनियर मॉडल को फिट करें और

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

फिर ब्याज के दो चर को प्लॉट करें और परिभाषा डोमेन के भीतर प्रतिगमन लाइन जोड़ें:

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

लगभग वहाँ पहुँच गया! अंतिम चरण भूखंड, प्रतिगमन समीकरण, rsquare के साथ-साथ सहसंबंध गुणांक को जोड़ना है। यह vector फ़ंक्शन का उपयोग करके किया 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')

ध्यान दें कि आप वेक्टर फ़ंक्शन को स्वीकार करके RMSE जैसे किसी अन्य पैरामीटर को जोड़ सकते हैं। कल्पना कीजिए कि आप 10 तत्वों के साथ एक किंवदंती चाहते हैं। वेक्टर परिभाषा निम्नलिखित होगी:

rp = vector('expression',10)

और आपको परिभाषित r[1] .... to r[10]

यहाँ उत्पादन है:

भार

कभी-कभी हम चाहते हैं कि मॉडल दूसरों की तुलना में कुछ डेटा बिंदुओं या उदाहरणों को अधिक वजन दे। मॉडल सीखते समय इनपुट डेटा के लिए वजन निर्दिष्ट करके यह संभव है। आम तौर पर दो प्रकार के परिदृश्य होते हैं जहाँ हम उदाहरणों पर गैर-समान भार का उपयोग कर सकते हैं:

• विश्लेषणात्मक भार: विभिन्न अवलोकनों की सटीकता के विभिन्न स्तरों को दर्शाते हैं। उदाहरण के लिए, यदि डेटा का विश्लेषण जहां प्रत्येक अवलोकन भौगोलिक क्षेत्र से औसत परिणाम है, तो विश्लेषणात्मक भार अनुमानित विचरण के व्युत्क्रमानुपाती होता है। प्रेक्षणों की संख्या को देखते हुए आनुपातिक भार प्रदान करके डेटा में औसत से निपटने के लिए उपयोगी। स्रोत

• सैम्पलिंग वेट (उलटा प्रायिकता वज़न - IPW): आँकड़ों की गणना के लिए एक सांख्यिकीय तकनीक जो उस डेटा से एकत्र की गई जनसंख्या से अलग है। एक विषम नमूना जनसंख्या और लक्ष्य निर्धारण की आबादी (लक्ष्य जनसंख्या) के साथ अध्ययन डिजाइन आवेदन में आम हैं। लापता मान वाले डेटा के साथ काम करते समय उपयोगी। स्रोत

lm() फ़ंक्शन एनालिटिकल वेटिंग करता है। सैंपलिंग वेट के लिए survey पैकेज का उपयोग एक सर्वेक्षण डिज़ाइन ऑब्जेक्ट बनाने और svyglm() चलाने के लिए किया जाता है। डिफ़ॉल्ट रूप से, survey पैकेज नमूना भार का उपयोग करता है। (नोट: lm() , और svyglm() साथ gaussian() सभी एक ही बिंदु अनुमान का उत्पादन करेंगे, क्योंकि वे दोनों गुणांक के लिए कम से कम भार वर्ग को हल करते हैं। वे मानक त्रुटियों की गणना कैसे की जाती हैं, इसमें भिन्नता है।)

परीक्षण डेटा

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

विश्लेषणात्मक वजन

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

उत्पादन

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

सैंपलिंग वेट (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)

उत्पादन

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

बहुपद प्रतिगमन के साथ nonlinearity के लिए जाँच

कभी-कभी रैखिक प्रतिगमन के साथ काम करते समय हमें डेटा में गैर-रैखिकता की जांच करने की आवश्यकता होती है। ऐसा करने का एक तरीका एक बहुपद मॉडल को फिट करना है और यह जांचना है कि क्या यह एक रैखिक मॉडल की तुलना में बेहतर डेटा फिट बैठता है। ऐसे अन्य कारण हैं, जैसे कि सैद्धांतिक, जो एक द्विघात या उच्चतर क्रम मॉडल को इंगित करने के लिए संकेत देता है क्योंकि यह माना जाता है कि चर संबंध स्वाभाविक रूप से बहुपद है।

mtcars डेटासेट के लिए एक द्विघात मॉडल फिट करते हैं। रैखिक मॉडल के लिए mtcars डेटासेट पर रैखिक प्रतिगमन देखें।

पहले हम चर का एक बिखराव साजिश बनाने के mpg (मील / गैलन), disp (विस्थापन (cu.in.)), और wt (वजन (1000 पाउंड))। mpg और disp के बीच संबंध गैर-रैखिक दिखाई देता है।

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

एक रेखीय फिट दिखाएगा कि disp महत्वपूर्ण नहीं है।

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

फिर, एक द्विघात मॉडल का परिणाम प्राप्त करने के लिए, हमने I(disp^2) जोड़ा। R^2 समय नया मॉडल बेहतर दिखाई देता है और सभी चर महत्वपूर्ण होते हैं।

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 

जैसा कि हमारे पास तीन चर हैं, सज्जित मॉडल एक सतह है जिसका प्रतिनिधित्व करता है:

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

बहुपद प्रतिगमन को निर्दिष्ट करने का एक अन्य तरीका पैरामीटर raw=TRUE साथ poly का उपयोग कर रहा है, अन्यथा ऑर्थोगोनल बहुपद help(ploy) अधिक जानकारी के लिए help(ploy) देखें) पर विचार किया जाएगा। हमें एक ही परिणाम प्राप्त होता है:

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

अंत में, क्या होगा यदि हमें अनुमानित सतह के एक भूखंड को दिखाने की आवश्यकता है? वैसे R में 3D प्लॉट बनाने के लिए कई विकल्प हैं। यहाँ हम का उपयोग Fit3d से p3d पैकेज।

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

गुणवत्ता मूल्यांकन

एक प्रतिगमन मॉडल के निर्माण के बाद परिणाम की जांच करना और यह तय करना महत्वपूर्ण है कि क्या मॉडल उपयुक्त है और हाथ में डेटा के साथ अच्छी तरह से काम करता है। यह अवशेष प्लॉट और साथ ही अन्य नैदानिक भूखंडों की जांच करके किया जा सकता है।

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

ये भूखंड दो मान्यताओं की जाँच करते हैं जो मॉडल बनाते समय बनाई गई थीं:

1. कि अनुमानित चर (इस मामले में mpg ) में अनुमानित मूल्य भविष्यवाणियों के एक रैखिक संयोजन (इस मामले में wt ) द्वारा दिया गया है। हम उम्मीद करते हैं कि यह अनुमान निष्पक्ष होगा। तो अवशिष्ट भविष्यवक्ताओं के सभी मूल्यों के लिए अर्थ के आसपास केंद्रित होना चाहिए। इस मामले में हम देखते हैं कि अवशिष्ट अंत में सकारात्मक होते हैं और बीच में नकारात्मक, चर के बीच एक गैर-रैखिक संबंध का सुझाव देते हैं।
2. कि वास्तविक अनुमानित चर सामान्य रूप से अपने अनुमान के आसपास वितरित किया जाता है। इस प्रकार, अवशिष्टों को सामान्य रूप से वितरित किया जाना चाहिए। सामान्य रूप से वितरित डेटा के लिए, एक सामान्य क्यूक्यू भूखंड में अंक विकर्ण पर या उसके करीब होना चाहिए। यहाँ सिरों पर तिरछी कुछ मात्रा है।

'प्रेडिक्ट' फ़ंक्शन का उपयोग करना

एक बार एक मॉडल बनाया गया है predict नए डेटा के साथ परीक्षण करने के लिए मुख्य कार्य है। हमारे उदाहरण का उपयोग करेगा mtcars निर्मित विस्थापन के खिलाफ प्रति गैलन वापसी मील डाटासेट:

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

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

Coefficients:
(Intercept)         disp  
   29.59985     -0.04122

यदि मेरे पास विस्थापन के साथ एक नया डेटा स्रोत था तो मैं प्रति गैलन अनुमानित मील देख सकता था।

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

प्रक्रिया का सबसे महत्वपूर्ण हिस्सा मूल डेटा के समान कॉलम नाम के साथ एक नया डेटा फ़्रेम बनाना है। इस मामले में, मूल डेटा में एक स्तंभ लेबल disp था, मुझे उसी नाम के नए डेटा को कॉल करना सुनिश्चित था।

सावधान

आइए कुछ सामान्य नुकसान देखें:

1. नए ऑब्जेक्ट में data.frame का उपयोग नहीं कर रहा है:

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

2. नए डेटा फ़्रेम में समान नामों का उपयोग नहीं करना:

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

शुद्धता

भविष्यवाणी की सटीकता की जांच करने के लिए आपको नए डेटा के वास्तविक y मूल्यों की आवश्यकता होगी। इस उदाहरण में, newdf को 'mpg' और '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