จากตัวอย่างปัญหาข้างต้น เป็นเรื่องการศึกษาความสัมพันธ์ระหว่างตัวแปรต้นกับตัวแปรตาม อาจเป็นการศึกษาไปเพื่อดู ความสัมพันธ์กันแบบเชิงเส้นเพื่อนำไปพยากรณ์เหตุการณ์ในอนาคตก็ได้

การวิเคราะห์การถดถอยเบื้องต้น(An introduction of regresdsion analysis)

การวิเคราะห์การถดถอยถอย คือการหาความสัมพันธ์ ระหว่างตัวอิสระกับตัวแปรตาม เช่นถ้าเรามีข้อมูลน้ำหนัก และส่วนสูงเฉลี่ยของเพศหญิงชาวอเมริกา

Height = ความสูง(นิ้ว)

Weight = น้ำหนัก(ปอนด์)

# ดูตัวแปรและข้อมูลเบื้องต้นของกรอบข้อมูล women
str(women)
## 'data.frame':    15 obs. of  2 variables:
##  $ height: num  58 59 60 61 62 63 64 65 66 67 ...
##  $ weight: num  115 117 120 123 126 129 132 135 139 142 ...
# วาดกราฟ
ggplot(data=women,   # ข้อมูล women
       aes(x=weight, # แกน x คือ weight
           y=height  # แกน y คือ height
           ))+
geom_point()         # วาดกราฟแบบจุด

จากกราฟ สามารถเขียนความสัมพันธ์ในของแบบฟังก์ชันได้คือ \[height =f(weight) \] หรือจากกราฟจะเห็นว่า เส้นที่จะใช้อธิบายความสัมพันธ์ได้ดีก็คือเส้นตรง

# วาดกราฟ
ggplot(data=women,   # ข้อมูล women
       aes(x=weight, # แกน x คือ weight
           y=height  # แกน y คือ height
           ))+
geom_point()+       # วาดกราฟแบบจุด  
geom_smooth(method = lm, # แสดงความสัมพันธ์ด้วยเส้นตรง  
            se = FALSE   # ไม่แสดงช่วงความเชื่อมั่น
            )

สามารถเขียนได้ในรูปแบบฟังก์ชันเชิงเส้นหรือสมการเส้นตรงดังนี้ \[height =a+b\times weight \] ค่า \(a\) จะเรียกว่าค่า intercept และค่า \(b\) เรียกว่าค่า coefficient หรือค่า slope ค่า \(a\) และ\(b\) เป็นที่เราต้องการหาจากข้อมูลที่มีนั้นเอง จากเส้นตรงที่ได้ จะเห็นไม่ได้ว่า เส้นตรงนี้เป็นเพียงเส้นที่ประมาณการขึ้น เพราะไม่ได้ลากผ่านคู่อันดับ \((weight,height)\) ทั้งหมด ดังนี้จะมีึค่าคาดเคลื่อนเกิดขึ้นดังนี้ คือ \[height_i =a+b\times weight_i+\varepsilon_i \] ค่า \(\varepsilon\) คืิอคาดเคลื่อนที่เกิดขึ้นจากการประมาณค่าหรือการพยากรณ์ส่วนสูงจากน้ำหนัก ซึ่งถ้าเราสมใจจะใช้สมการเส้นตรงในการศึกษาความสัมพันธ์ระหว่างน้ำหนักกับส่วนสูงนี้ ค่า \(a\) และ \(b\) ควรจะเป็นเท่าไหร่ และค่าคาดเคลื่อน(\(\varepsilon\)) ที่เกิดขึ้นควรจะมีคุณสมบัติอย่างไร การประมาณ \(a\) และ \(b\) ที่นิยมใช้กันในทางสถิติ โดยเบื้องต้นคือวิธีกำลังสองน้อยที่สุด(ordinary lease square) คือการหา \(a\) และ \(b\) ที่ทำให้ \(\sum_{i}^n \varepsilon_i^2\) มีึค่่าน้อยที่สุด ซึ่งสามารถเขียนได้ดังนี้

\[\begin{align*} \varepsilon_i &=height_i -a-b\times weight_i,\forall i=1,2,\cdots,15\\ \varepsilon_i^2&=(height_i -a-b\times weight_i)^2,\forall i=1,2,\cdots,15\\ \sum_{i=1}^{15}\varepsilon_i^2&=\sum_{i=1}^{15}(height_i -a-b\times weight_i)^2 \end{align*}\] จากสมการสุดท้ายสามารถหาค่า \(a\) และ \(b\) ได้โดยอนุพันธ์ของ \(\sum_{i=1}^{15}\varepsilon_i^2\) เทียบกับ \(a\) และ เทียบกับ \(b\) ตามลำดับได้ผลลัพธ์ดังนี้คือ

\[\begin{align*} \dfrac{\partial}{\partial a}\sum_{i=1}^{15}\varepsilon_i^2 &=-2\sum_{i=1}^{15}(height_i-a-b\times weight_i)\\ \dfrac{\partial}{\partial b}\sum_{i=1}^{15}\varepsilon_i^2 &=-2\sum_{i=1}^{15} weight_i\times(height_i-a-b\times weight_i)\\ \end{align*}\]

คุณสมบัติของเครื่องผลรวม(summation symbol)

\[\begin{align*} 1)& \sum_{i=1}^n c =n\times c\\ 2)& \sum_{i=1}^nx_i=x_1+x_2+\cdots+x_n \\ 3)& \sum_{i=1}^n c\times x_i=c\sum_{i=1}^nx_i\\ 3)& \sum_{i=1}^n (x_i+y_i)=\sum_{i=1}^nx_i+\sum_{i=1}^ny_i\\ 5)& \sum_{i=1}^n x_i y_i=x_1y_1+x_2y_2+\cdots x_ny_n\\ 6)& \sum_{i=1}^n x_i^k= x_1^k+x_2^k+\cdots+x_n^k,k\geq 0\\ 7)& \sum_{i=1}^n(x_i-y_i)^2=\sum_{i=1}^n x_i^2-2\sum_{i=1}^nx_iy_i+\sum_{i=1}^ny_i^2 \end{align*}\] และค่าเฉลี่ยคือ \[\bar{x}=\dfrac{\sum_{i=1}^n x_i}{n} \] ดังนั้นเมื่อกำหนดให้ \(\displaystyle\dfrac{\partial}{\partial a}\sum_{i=1}^{15}\varepsilon_i^2=0\) \(\displaystyle\dfrac{\partial}{\partial a}\sum_{i=1}^{15}\varepsilon_i^2=0\) จะได้

\[\begin{align*} \sum_{i=1}^{15} height_i-15a-b\sum_{i=1}^{15}weight_i&=0\\ \sum_{i=1}^{15}height_i\times weight_i-a\sum_{i=1}^{15}height_i-b\sum_{i=1}^{15}height_i^2&=0 \end{align*}\] สามารถจัดอยู่ในรูปของเมตริกซ์ \(Ax=b\) ได้คือ \[\displaystyle\begin{bmatrix} 15&\sum_{i=1}^{15}weight_i\\ \sum_{i=1}^{15}weight_i&\sum_{i=1}^{15}weight_i^2 \end{bmatrix}\begin{bmatrix}a\\b\end{bmatrix}= \begin{bmatrix}\sum_{i=1}^{15}height_i\\ \sum_{i=1}^{15}height_i \times weight_i\end{bmatrix}\] หลังจากนี้สามารถใช้วิธีทางพีชคณิตเชิงเส้นเพื่อค่า \(a\) และ \(b\) ได้ดังนี้คือ \[\begin{align*} a&=\dfrac{(\sum_{i=1}^{15}height_i)(\sum_{i=1}^{15}weight_i^2)-\sum_{i=1}^{15}weight_i\times (\sum_{i=1}^{15}height_i\times weight_i)}{n\sum_{i=1}^{15}weight_i^2-(\sum_{i=1}^{15}weight_i)^2}, \\ b&=\dfrac{15(\sum_{i=1}^{15}height_i\times weight_i)-(\sum_{i=1}^{15}height_i)(\sum_{i=1}^{15}weight_i)}{n\sum_{i=1}^{15}weight_i^2-(\sum_{i=1}^{15}weight_i)^2} \end{align*}\] สามารถใช้ภาษาอาร์ในการคำนวณได้ดังนี้

# หาผลรวมของตัวแปร weight
sum_h<-sum(women$height)
sum_h
## [1] 975
# หาผลรวมของตัวแปร height
sum_w<-sum(women$weight)
sum_w
## [1] 2051
# หาผลรวมของตัวแปร weight คูณกับ height
sum_hw<-sum(women$height*women$weight)
sum_hw
## [1] 134281
# หาผลรวมของตัวแปร weight กำลังสอง
sum_w2<-sum(women$weight^2)
sum_w2
## [1] 283803
# ประมาณค่า a จากสูตร
a<-(sum_h*sum_w2-sum_w*sum_hw)/(15*sum_w2-(sum_w)^2)
a 
## [1] 25.72346
# ประมาณค่า b จากสูตร
b<-(15*sum_hw-sum_h*sum_w)/(15*sum_w2-(sum_w)^2)
b 
## [1] 0.2872492

สรุปสำหรับกรณีทั่วไปคือ ถ้าพิจารณาความสัมพันธ์ระหว่างตัวแปรอิสระ \(x_i\) และตัวแปรตาม\(y_i\) ด้วยสมการการถดภอยอย่างง่ายโดยมีข้อมูลจำนวน \(n\) ตัวอย่างจะได้ \[y_i=a+bx_i+\varepsilon_i,i=1,2,\cdots,n\] โดยที่ \[\begin{align*} a&=\frac{\left(\Sigma y_i\right)\left(\Sigma x_i^{2}\right)-\left(\Sigma x_i\right)(\Sigma x_i y_i)}{n\left(\Sigma x_i^{2}\right)-(\Sigma x_i)^{2}}\\ b&=\frac{n\left(\Sigma x_i y_i\right)-\left(\Sigma x_i\right)\left(\Sigma y_i\right)}{n\left(\Sigma x_i^{2}\right)-\left(\Sigma x_i\right)^{2}} \end{align*}\] ถ้าจัดรูปให้ดีจะได้ \[\begin{align*} a&=\bar{y}-b\bar{x}\\ b&=\dfrac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^n(x_i-\bar{x})^2}=\dfrac{cov(x,y)}{var(x)}=r_{xy}\sqrt{\dfrac{\sum_{i=1}^n(y_i-\bar{y})^2}{\sum_{i=1}^n(x_i-\bar{x})^2}}\\ &=r_{x,y}\dfrac{\sigma_y}{\sigma_x} \end{align*}\] โดยที่ \[\begin{align*} cov(x,y)&=\dfrac{1}{n}\sum_{i=1}^n (x_i-\bar{x})(y-\bar{y})\\ var(x)&=\sigma_x^2=\dfrac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2\\ var(y)&=\sigma_y^2=\dfrac{1}{n}\sum_{i=1}^n(y_i-\bar{y})^2\\ r_{x,y}&=cor(x,y)=\dfrac{cov(x,y)}{\sigma_x \sigma_y}=\dfrac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})} {\sqrt{\sum_{i=1}^n(x_i-\bar{x})^2}\sqrt{\sum_{i=1}^n(y_i-\bar{y})^2}} \end{align*}\] สรุปได้ว่า

# คำนวณค่า b โดยการคำนวณโดยตรงจากสูตร
b<-cov(women$height,women$weight)/var(women$weight)
b
## [1] 0.2872492
# คำนวณค่า b โดยการคำนวณโดยตรงจากอีกสูตร
b<-cor(women$height,women$weight)*sd(women$height)/sd(women$weight)
b
## [1] 0.2872492
# คำนวณค่า a จากสูตร
a<-mean(women$height)-b*mean(women$weight)
a
## [1] 25.72346

การพยากรณ์(predict)ค่าที่สนใจก็สามารถทำได้โดยใช้่ \[y=a+bx \] ได้เลยเช่น

y<-a+b*150
y
## [1] 68.81084

ค่าคาดเคลื่อน \(\varepsilon\)มีก่ีแจงแจงแบบปกติ ที่ค่าคาดหวังเท่ากับ 0 และความแปรปรวนเท่ากับ \(\sigma^2\) หรือ \(\varepsilon\sim N(0,\sigma^2)\) จากคุณสมบัตินี้ทำให้ได้ว่าหวังใช้ประมาณค่า \(a\) และ \(b\) ด้วยวิธีกำลังสองน้อยที่สุดแล้วจะได้ \(\sum{i=1}^n \varepsilon_i=0\) หรือ \(\bar{\varepsilon}=0\) ดังนั้น \(\sigma^2=var(\varepsilon)=\dfrac{\sum_{i=1}^n(\varepsilon_i-\bar{\varepsilon})^2}{n}=\dfrac{\sum_{i=1}^n\varepsilon_i^2}{n}\) แต่ถ้าต้องการคุณสมบัติความไม่เอนเอียง การคำนวณค่าแปรปรวนของค่าความคาดเคลื่อนจะใช้สูตร \[\sigma_{\varepsilon}^2=var(\varepsilon)=\dfrac{\sum_{i=1}^n\varepsilon_i^2}{n-2} \] เราเรียกค่าความแปรปรวน \(\sigma_{\varepsilon}^2\) นี้ว่า mean square error (MSE)

# คำนวณค่าคาดเคลื่อน
# error = y-a-bx
error <-women$height-a-b*women$weight
# จำนวนข้อมูลทั้งหมด
n<-length(women[,1])
n
## [1] 15
#คำนวณค่า mean square error จากสูตร
MSE<-sum(error^2)/(n-2)
MSE
## [1] 0.1936344

ระวังสับสนกับสูตรการความแปรปรวนจากข้อมูลที่เก็บได้โดยทั่วไป ในกรณีนี้เป็นสูตรประมาณค่าความแปรปรวนของค่าคาดเคลื่อนที่เกิดขึ้นจากการประมาณค่าด้วยวิธีกำลังสองน้อยที่สุด (ถ้ามีจำนวน \(n\) ใหญ่มากพอ การหารด้วย \(n\) หรือ \(n-1\) หรือ \(n-2\) จะไม่แตกต่างกันอย่างมีนัยยะสำคัญ) ซึ่งจะถูกใช้ในการวิเคราะห์ความแปรปรวนเพื่อพิจารณาความเหมาะสมตัวแปรการถดถอยนี้

ความแปรปรวนของพารามิเตอร์ \(a\) และ \(b\)

สามารถคำนวณได้จากสูตรนี้ ความแปรปรวนของ \(a\) สามารถคำนวณได้จากสูตร \[\sigma_{a}^2=\sigma_{\varepsilon}^2\left( \dfrac{1}{n}+\dfrac{\bar{x}^2}{\sum_{i=1}^n(x_i-\bar{x})^2} \right) \]

# ความแปรปรวณของตัวแปร a โดยสูตร
variance_a <- MSE*(1/n+mean(women$weight)^2/(sum((women$weight-mean(women$weight))^2)))
variance_a
## [1] 1.089406

ความแปรปรวนของ \(b\) คำนวณได้จาก \[\sigma_{b}^2=\dfrac{\sigma_{\varepsilon}^2}{\sum_{i=1}^n(x_i-\bar{x})^2} \]

# ความแปรปรวณของตัวแปร b โดยสูตร
variance_b<-MSE/sum((women$weight-mean(women$weight))^2)
variance_b
## [1] 5.757901e-05

การพืสูจน์ที่มาของสูตรดูได้จาก… ส่วนเบี่ยงเบนมาตราฐาน\(\sigma_{a}\) และ \(\sigma_b\) คำนวณได้จาก \[\begin{align*} \sigma_a=\sqrt{\sigma_a^2}\\ \sigma_b=\sqrt{\sigma_b^2} \end{align*}\]

# ส่วนเบี่ยงเบนมาตรฐานของตัวแปร a
SD_a <-sqrt(variance_a)
SD_a
## [1] 1.043746
# ส่วนเบี่ยงเบนมาตรฐานของตัวแปร b
SD_b <-sqrt(variance_b)
SD_b
## [1] 0.007588083

เมื่อคำนวณความแปรปรวนหรือส่วนเบี้ยงเบียนได้แล้วกัน ก็จะสามารถคำนวณช่วงความเชื่อมั่นได้ต่อไป ถ้ากำหนดช่วงความเชื่อมั่นไว้ เท่ากับ \(1-\alpha\) ก็จะสามาถคำนวณได้จาก

ข่วงความช่วงของ \(a\) คำนวณได้จากสูตร

\[a\pm t^{-1}_{n-2}(1-\frac{\alpha}{2})\times\sigma_{a}=\left[a- t^{-1}_{n-2}(1-\frac{\alpha}{2})\times\sigma_{a},a+ t^{-1}_{n-2}(1-\frac{\alpha}{2})\times\sigma_{a}\right] \] ถ้าต้องการช่วงความเชื่อมั่นไว้ $95% $ หรือ \(\alpha=5\%\) จะได้

# ช่วงความเชื่อมั่นต่ำสุดของ a
lower_a<- a-qt(1-.05/2,n-2)*SD_a
lower_a
## [1] 23.46858
# ช่วงความเชื่อมั่นสูงสุดของ a
upper_a<- a+qt(1-.05/2,n-2)*SD_a
upper_a
## [1] 27.97833

ในทำนองเดียวกัน สามารถคำนวณช่วงของเชื่อของ \(b\) ได้เท่ากับ \[b\pm t^{-1}_{n-2}(1-\frac{\alpha}{2})\times\sigma_{b}=\left[b- t^{-1}_{n-2}(1-\frac{\alpha}{2})\times\sigma_{b},b+ t^{-1}_{n-2}(1-\frac{\alpha}{2})\times\sigma_{b}\right] \] ถ้าต้องการช่วงความเชื่อมั่นไว้ $95% $ หรือ \(\alpha=5\%\) จะได้

# ช่วงความเชื่อมั่นต่ำสุดของ b
lower_b<- b-qt(1-.05/2,n-2)*SD_b
lower_b
## [1] 0.2708562
# ช่วงความเชื่อมั่นสูงสุดของ b
upper_b<- b+qt(1-.05/2,n-2)*SD_b
upper_b
## [1] 0.3036423

ช่วงความเชื่อมั่นของค่าพยากรณ์ \(y^\ast\) เมื่อทราบค่า \(x^\ast\) สามารถคำนวณได้จากสูตร \[a+bx^\ast \pm t_{df=n-2}^{-1}(1-\frac{\alpha}{2})\sigma_\varepsilon\sqrt{\dfrac{1}{n}+\dfrac{(x^\ast-\bar{x})^2}{\sum_{i=1}^n(x_i-\bar{x})^2}} \] การกระจายให้ผู้อ่านลองทำด้วยตนเอง สมมุติให้ \(x^\ast= 145\) จะได้ช่วงความเชื่อมั่นของ \(y^\ast\) คือ

# ช่วงความเชื่อมั่นต่ำสุดของ y เมื่อให้ค่า x
lower_y<-a+b*145-qt(1-.05/2,n-2)*
         sqrt(MSE)*sqrt(1/n+(145-mean(women$weight))^2/((n-1)*var(women$weight)))
lower_y
## [1] 67.09421
# ช่วงความเชื่อมั่นสูงสุดของ y เมื่อให้ค่า x
upper_y<-a+b*145+qt(1-.05/2,n-2)*
         sqrt(MSE)*sqrt(1/n+(145-mean(women$weight))^2/((n-1)*var(women$weight)))   
upper_y
## [1] 67.65497

ในหนังสิอบางเล่มกล่าวว่า ถ้า \(n-2>30\) สามารถเปลี่ยนจากการใช้การแจกแจง \(t\) ไปเป็นการแจกแจงปกติ \(z\) แทนได้ ผู้อ่านลองพิจารณา ค่าจากการคำนวณข้างล่างแล้วพิจารณาดูว่า แทนได้จริงหรือไม่ โดยการพิจาณา ความแตกต่างระหว่าง ค่าจากการแจกแจงปกติ กับการแจกแจง t ที่เพิ่มค่า degree of freedom ขึ้นไปเรื่อยๆ

qnorm(.9) - qt(.9,30)
## [1] -0.02886346
qnorm(.9) - qt(.9,100)
## [1] -0.008523196
qnorm(.9) - qt(.9,500)
## [1] -0.001695455
qnorm(.9) - qt(.9,1000)
## [1] -0.0008471559
qnorm(.9) - qt(.9,10000)
## [1] -8.466419e-05

จะเห็นได้ว่ายิ่งค่า degree of freedom เพิ่มขึ้นความคาดเคลื่อนการจากประมาณค่าสถิติ \(t\) ด้วยค่า \(z\) จะมีค่าคาดเคลื่อนลดลงเรื่อยๆ เมื่อ ค่า degree of freedom เพิ่มขึ้นหรือเมื่อมีจำนวนข้อมูลมากขึ้น เพราะจากทฤษฎีเราทราบว่า \[\lim_{df\rightarrow \infty}t(x,df)=z(x) \] ### สมมุติฐานและคุณสมบัติที่ได้จากวิธีกำลังสองน้อยที่สุด สมมุติฐานสำหรับการวิเคราะห์การถดถอยแบบเชิงเส้นอย่างง่าย 1) ตัวแปรต้นและตัวแปรเมื่อทำการ scatter point ควรมีรูปร่างคล้ายเส้นตรง เช่นข้อมูลชุด cars เป็นข้อมูลระหว่างความเร็วที่ใช้สูงสุด และระยะทางที่ใช้เบรค

str(cars)
## 'data.frame':    50 obs. of  2 variables:
##  $ speed: num  4 4 7 7 8 9 10 10 10 11 ...
##  $ dist : num  2 10 4 22 16 10 18 26 34 17 ...

เมื่อทำการ scatter plot ระหว่าง \((speed,dist)\)

ggplot(data=cars,aes(x=speed,y=dist))+
geom_point()

หรือข้อมูลขุด mtcars

str(mtcars)
## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

เมื่อทำการ scatterplot ระหว่างที่ได้ต่อน้ำมัน 1 แกลลอน(miles/(US) gallon) mpg กับขนาดของลูกสูบ(disp)

# วาดกราฟ
ggplot(data=mtcars, # ข้อมูล mtcars
       aes(x=mpg,y=disp))+ # แกน x คือ mpg แกน y คือ disp
geom_point()+     #วาด scatter plot
geom_smooth(method = "auto",se=FALSE)+ # วาดกราฟเส้นแบบอัตโนมัติ ไม่แสดงช่วงความเชื่อมั่น
geom_smooth(method = "lm",se=FALSE,color="red")  # วาดกราฟเส้นแบบเส้นตรง ไม่แสดงช่วงความเชื่อมั่น
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

จะเห็นว่าความสัมพันธ์อาจจะไม่ใช่เส้นตรง

การพิจารณาว่าค่าคาดเคลื่อนมีการแจกแจงแบบปกติหรือไม่สามารถพิจารณาได้จากการกราฟ histogram และกราฟ QQ plot

# คำนวณหาค่าความคาดเคลื่อน
error<-women$height-a-b*women$weight
# เปลี่ยน error เป็นกรอบข้อมูล
error<-data.frame(error=error)
# วาดกราฟ
ggplot(error, aes(x=error)) + 
geom_histogram(aes(y=..density..)) # วาดกราฟ histogram
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

จะเป็นว่า ไม่สามารถใช้ได้การเพราะ histogram จะใช้การได้ดีถ้าข้อมูลมีจำนวนมาก จึงควรเลื่อกใช้ QQ plot โดยคำสั่งเฉพาะดังนั้

ggplot(error, # ข้อมูลคือ error
       aes(sample=error))+ # sample คือ error
stat_qq()+ # วาด qqplot โดยใช้จุด
stat_qq_line() # วาด qq plot โดยเพิ่มเส้น

ถ้าแต่ละจุดมีอยู่ห่างจากเส้นตรงไม่มาก แสดงว่าค่าคาดเคลื่อนอาจจะมีการแจกแจงแบบปกติ ถ้าต้องการยืนยันว่าค่าคลานเคลื่อนมีการจากแจกแจงแบบปกติหรือไม่ให้วิธีการทดสอบสมมมุติฐาน โดยอาจจะใช้วิธี Jarque-Bera Test จาก package tseries โดยJarque-Bera TestZ(JB test)นี้มีสมมุติฐานหลักคือ ข้อมูลมีการแจกแจงแบบปกติ

library(tseries)
## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
# Jarque-Bera Test 
 jarque.bera.test(error$error) 
## 
##  Jarque Bera Test
## 
## data:  error$error
## X-squared = 1.4855, df = 2, p-value = 0.4758

จะเห็นถ้ากำหนดค่า p-value = 5% เราจะยอมรับค่าคาดเคลื่อนมีการแจกแจงแบบปกติ

การถดสอบสมมุติฐาน

ในการวิเคราะห์การถดถอยเบื้องต้น จะสนใจทดสอบสมสมมติฐานหลักๆ คือ

  1. ค่าพารามิตอร์ \(a\) หรือ \(b\) เท่ากับ 0 หรือไม่

  2. ค่าคาดเคลื่อนที่ประมาณได้มีการแจงแจกแบบปกติหรือไม่(แสดงไปแล้วในหัวข้อก่อนหน้านี้)

ในการทดสอบสมมุติฐานว่า \(a\) และ\(b\) เท่ากับ 0 หรือไม่นั้นจะมีการตั้งสมมุติฐานหลักและสมมุติฐานรองดังนี้ สำหรับ \(a\) \[\begin{align*} H_0:& a=0\\ H_a:& a\neq 0\\ \end{align*}\] สำหรับ \(b\) \[\begin{align*} H_0:& b=0\\ H_a:& b\neq 0\\ \end{align*}\] จากหัวข้อก่อนหน้า ทำให้เราคาดได้ว่าค่า \(a\) และ \(b\) มีการแจกแจงแบบ \(t\) ดังนี้คือ \[\begin{align*} \dfrac{a}{\sigma_a}\sim t_{df=n-2}()\\ \dfrac{b}{\sigma_b}\sim t_{df=n-2}()\\ \end{align*}\] ที่ระดับนัยสำคัญ \(\alpha\) จะยอมรับสมมุติฐานถ้า สำหรับ \(a\) \[t_{df=n-2}^{-1}(\dfrac{\alpha}{2})\leq \dfrac{a}{\sigma_a}\leq t_{df=n-2}^{-1}(1-\dfrac{\alpha}{2})\] จากข้อมูลชุดเดิม จะทำเช็คค่า สถิติที่คำนวณนี้อยู่ระหว่างต่ำสุดกับค่าสูงสุดหรือไม่

# ค่าที่คำนวณได้อยู่ระหว่างค่าต่ำสุดกับสูงสุดหรือไม่
(a/SD_a > qt(.05/2,n-2))&(a/SD_a < qt(1-.05/2,n-2))
## [1] FALSE

แสดงว่ายอมรับสมมุติฐานรอง ที่ค่า \(\alpha= 5\%\) เพราะว่าค่าสถิติ ที่คำนวณได้มีค่ามากกว่าคู่สูงสุด นั่นคือค่า \(a\neq 0\) สำหรับ \(b\) \[t_{df=n-2}^{-1}(\dfrac{\alpha}{2})\leq \dfrac{b}{\sigma_b}\leq t_{df=n-2}^{-1}(1-\dfrac{\alpha}{2})\] จากข้อมูลชุดเดิม จะทำเช็คค่า สถิติที่คำนวณนี้อยู่ระหว่างต่ำสุดกับค่าสูงสุดหรือไม่

# ค่าที่คำนวณได้อยู่ระหว่างค่าต่ำสุดกับสูงสุดหรือไม่
(b/SD_b > qt(.05/2,n-2))&(b/SD_b < qt(1-.05/2,n-2))
## [1] FALSE

แสดงว่ายอมรับสมมุติฐานรอง ที่ค่า \(\alpha=5\%\) เพราะค่าสถิติที่คำนวณได้มีค่ามากกว่าค่าสูงสุด นั่นคือยอมรับว่าค่า \(b\neq 0\)

การพิจารณาว่าสมการถดถอยอย่างง่ายที่ใช้ เหมาะสมดีไม่

นอกจากทดสอบสมมุติฐาน ค่า \(a\) และ \(b\) แล้วอีกวิธีหนึ่งที่สามารถใช้ได้ คือการพิจารณาค่าสัมประสิทธิ์ตัวกำหนด \(R^2\) จากค่าสัมประสิทธิ์สหสัมพันธ์ \(r\) ที่คำนวณได้โดยมีความสัมพันธ์กันดังนี้ \[R^2=r^2=cor(x,y)^2=\dfrac{cov(x,y)^2}{var(x)\times var(y)} \] จากตัวอย่างจะได้

# คำนวณค่าส้มประสิทธิ์ตัวกำหนดโดยสูตร
R2<-cor(women$weight,women$height)^2
R2
## [1] 0.9910098

ค่า \(R^2=.9910\) หมายความว่า สมการการถดถอยนี้สามารถอธิบาย ส่วนสูงจากน้ำหนักของผู้หญิงอเมริกาได้ถึง \(99.10\%\) สำหรับค่า adjust \(R^2\) จะพิจารณาภายหลัง เมื่อเราทำศึกษาวการวิเคราะห์ถดถอยที่ตัวแปรต้นหลายตัว

คำสั่งใน R

โดยปกติแล้ว เราสามารถใช้โปรแกรมอาร์หาสมการการถดถอยออกมาให้ได้เลย โดยใช้คำสั่งดังนี้

# ประมาณค่าพารามิเตอร์ในตัวแบบการถดถอย
fit.lm <-lm(height~weight, # height= a+b*weight
            data=women)    # ข้อมูลคือ women
fit.lm
## 
## Call:
## lm(formula = height ~ weight, data = women)
## 
## Coefficients:
## (Intercept)       weight  
##     25.7235       0.2872

และใช้คำสั่ง summary() เพื่อดูค่าทั้งหมดจะได้

# แสดงค่าทั้งหมดของสมการการถดถอย
summary(fit.lm)
## 
## Call:
## lm(formula = height ~ weight, data = women)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83233 -0.26249  0.08314  0.34353  0.49790 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 25.723456   1.043746   24.64 2.68e-12 ***
## weight       0.287249   0.007588   37.85 1.09e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.44 on 13 degrees of freedom
## Multiple R-squared:  0.991,  Adjusted R-squared:  0.9903 
## F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14

จากตารางข้่างบนสามารถบอกได้ว่า - ค่า \(a\) เท่ากับ 25.723456 โดยมีค่า \(\sigma_a=1.043746\) จากทดสอบสมมุติฐานที่ว่า \(H_0: a=0\) ที่ระดับความเชื่อมั่น \(5\%\) จะปฎิเสธสมมุติฐานหลัก เพราะมีค่า P-value \(=2.68e-12<.05\)

ส่วนการพยากรณ์ด้วยค่าที่สนใจสามารถใช้คำสั่ง ได้เช่น ถ้าสนใจค่า \(weight =150\) ต้องสร้างกรอบข้อมูลขึ้นมาใหม่โดยมีตัวแปรชื่อว่า weight ไม่เช่นนั้นจะไม่สามารถใช้คำสั่งได้

# สร้างกรอบข้อมูล newdata โดยมีตัวแปรภายในเป็น weight=150
newdata<-data.frame(weight=150)
# พยากรณ์ height ที่ weight=150
predict(fit.lm,newdata)
##        1 
## 68.81084

ถ้าต้องการช่วงความเชื่อมั่น \(95\%\) ด้วยก็สามารถทำได้โดยเพิ่มคำสั่ง

# ช่วงความเชื่อ 95% ของ height ที่ weight= 150
predict(fit.lm,newdata,interval = "prediction")
##        fit      lwr      upr
## 1 68.81084 67.80522 69.81646

หรือ

# ช่วงความเชื่อ 95% ของ height ที่ weight= 150
predict(fit.lm,newdata,interval = "confidence")
##        fit     lwr      upr
## 1 68.81084 68.4829 69.13878

ความแตกต่าง ระหว่าง การเลือกใช้ “prediction” กับ “interval” ก็คือ ช่วงความเชื่อมั่นที่เกิดจาก prediction จะเกิดจากส่วนเบี่ยงเบนมาตราฐาน \(\sigma_\varepsilon\) ของ \(\varepsilon\) ส่วน ช่วงความเชื่อมั่นด้วยคำสั่ง interval จะเกิดจากสูตร \[a+bx^\ast \pm t_{df=n-2}^{-1}(1-\frac{\alpha}{2})\sigma_\varepsilon\sqrt{\dfrac{1}{n}+\dfrac{(x^\ast-\bar{x})^2}{\sum_{i=1}^n(x_i-\bar{x})^2}} \] ถ้าค่า \(x^\ast\) ห่างจาก \(\bar{x}\) มาก ช่วงความเชื่อมั่นก็จะกว้างขึ้น

# กรอบข้อมูล height.prediction จะมีค่า พยากรณ์ ช่วงความเชื่อมั่นต่ำสุด และช่วงความเชื่อมั่นสูงสุด
height.prediction <- predict(fit.lm,interval="prediction")
## Warning in predict.lm(fit.lm, interval = "prediction"): predictions on current data refer to _future_ responses
# รวมกรอบข้อมูล height.prediction เข้ากับ women
women<-cbind(women,height.prediction)
str(women)
## 'data.frame':    15 obs. of  5 variables:
##  $ height: num  58 59 60 61 62 63 64 65 66 67 ...
##  $ weight: num  115 117 120 123 126 129 132 135 139 142 ...
##  $ fit   : num  58.8 59.3 60.2 61.1 61.9 ...
##  $ lwr   : num  57.7 58.3 59.2 60 60.9 ...
##  $ upr   : num  59.8 60.4 61.2 62.1 62.9 ...

สร้างกราฟ เพื่อให้เห็นภาพชัดเจนขึ้น

ggplot(data=women,aes(x=weight,y=height))+
geom_point()+
geom_smooth(method="lm",se =TRUE)+ #วาดสมการถดถอย
geom_line(aes(y = lwr), color = "red", linetype = "dashed")+ #วาดเส้นช่วงความเชื่อมั่นต่ำสุด
geom_line(aes(y = upr), color = "red", linetype = "dashed")  #วาดเส้นช่วงความเชื่อมั่นสูงสุด

คำสั่ง geom_smooth(method=“lm”,se =TRUE) เป็นวาดเส้นสมการการถดถอย และช่วงความเชื่อมั่นที่เกิดจากการใช้ค่า\(weight\) เอง เส้นสีแดงขีด คือ ช่วงความเชื่อมั่นที่เกิดส่วนเบี่ยงเบนมาตราฐานของ \(\varepsilon\) หรือค่าคาดเคลื่อนนั่นเอง

การวิเคราะห์การถดถอยหลายตัวแปร

ก็คือการการถดถอยที่มีตัวแปรต้นหรือตัวแปรอิสระตั้งแต่ 2 ตัวเป็นไปนั้นเอง หรือสามารถเขียนได้ในรูป \[Y=\beta_0+\beta_1 X_1+\beta_2 X_2+\cdots+\beta_n X_n+\varepsilon \] เพื่อให้เข้าใจได้ง่ายขึ้น จะขอยก ตัวอย่างต่อไปเป็นยอดขาย Album ที่ขึ้นอยู่กับโฆษณา และระยะเวลาการออกอากาศ สามารถ โหลดไฟล์ได้จาก สามารถโหลดไฟล์ที่นี่

Album<- read.csv(file="/Users/somsak_mac2/Documents/thaiXelatex/Statistics with R/Album Sales.csv")
str(Album)
## 'data.frame':    200 obs. of  3 variables:
##  $ adverts: num  10.3 985.7 1445.6 1188.2 574.5 ...
##  $ sales  : int  330 120 360 270 220 170 70 210 200 300 ...
##  $ airplay: int  43 28 35 33 44 19 20 22 21 40 ...

ดังน้ั้นถ้าเราสนใจความสัมพันธ์ดังนี้ \[ \text{sales }=a+b_1*\text{adverts}+b_2*\text{airplay}+\varepsilon\] พิจารณาค่า correlation โดย

cor(Album)
##           adverts     sales   airplay
## adverts 1.0000000 0.5784877 0.1018828
## sales   0.5784877 1.0000000 0.5989188
## airplay 0.1018828 0.5989188 1.0000000

สามารถสร้างกราฟ scatter plot ของตัวแปรแต่คู่ทั้งหมดได้โดยใช้คำสั่ง

pairs(Album)

ถ้าต้องการให้เรียงรูปแบบของกราฟ โดยกำหนดตัวแปรเอง ทำได้ดังนี้

pairs(~sales+adverts+airplay,data=Album)

คำถาม ถ้่าต้องการวาดกราฟนี้ โดยไม่ใช้ตัวแปร airplay จะทำอย่างไร?

ถ้าต้องการใช้สี หรือเปลี่ยนรูปทรงของจุดในกราฟ ก็สามารถทำได้เหมือนกับคำสั่ง plot

pairs(~sales+adverts+airplay,data=Album,
      col ="blue",
      pch= 8)

ถ้าต้องการให้กราฟมีความสวยงามขึ้นโดยใช้ ggplot ขอแนะนำให้ใช้ชุดคำสั่ง GGally และใช้คำสั่งต่อไปนี้

library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
## 
## Attaching package: 'GGally'
## The following object is masked from 'package:dplyr':
## 
##     nasa
ggpairs(Album)

โดยเพิ่มการประมาณค่าการแจกแจงของแต่ละตัวแปร ด้วยวิธี kernel density ในแนวทแยงมุม และแสดงค่าสัมประสิทธิ์สหสัมพันธ์ระหว่างตัวแปรด้วย เช่น ค่าสัมประสิทธิ์สหสัมพันธ์ระหว่าง Sale กับ adverts มีค่าเท่ากับ .578

การวิเคราะห์การถดถอยโดยอาร์ คือ

 lm(sales~adverts+airplay, #สมการ (sales =a+b1*adverts+b2*airplay)
    data=Album)            # ข้อมูล Album
## 
## Call:
## lm(formula = sales ~ adverts + airplay, data = Album)
## 
## Coefficients:
## (Intercept)      adverts      airplay  
##    41.12381      0.08689      3.58879

ถ้าต้องการข้อมูลทางสถิติด้วย

 summary( lm(sales~adverts+airplay,data=Album))
## 
## Call:
## lm(formula = sales ~ adverts + airplay, data = Album)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -112.121  -30.027    3.952   32.072  155.498 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 41.123811   9.330952   4.407 1.72e-05 ***
## adverts      0.086887   0.007246  11.991  < 2e-16 ***
## airplay      3.588789   0.286807  12.513  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 49.38 on 197 degrees of freedom
## Multiple R-squared:  0.6293, Adjusted R-squared:  0.6255 
## F-statistic: 167.2 on 2 and 197 DF,  p-value: < 2.2e-16

การพยากรณ์ด้วยการถดถอยเดิมถ้ามีข้อมูลเพิ่มเข้ามา

ต้องสร้างกรอบข้อมูล โดยที่ ชื่อตัวแปรภายในกรอบข้อมูลต้องมีเขียนเดียวกับสมการการถดถอย และควรต้องเก็บค่าทั้งหมดที่ได้จากคำสั่ง lm ไว้ในตัวแปรใหม่

New.data <-data.frame(adverts= c(180,230, 250),
                      airplay=c(30,40, 45))

ก็จะพยากรณ์ได้โดยใช้คำสั่งนี้

Sale.fit<- lm(sales~adverts+airplay,data=Album)
predict(Sale.fit,         # ตัวแปรที่เก็บข้อมูลทั้งหมดจากการคำนวณหาสมการการถดถอย
        newdata= New.data)# ข้อมูลชุดใหม่ที่มี เมื่อต้องการพยากรณ์
##        1        2        3 
## 164.4272 204.6594 224.3411

ถ้าต้องการค่าคาดเคลื่อนมาพิจารณาว่ามีการแจกแจงแบบปกติหรือโดยใช้ QQ plot ก็ควรต้องสร้างตัวแปร สมมุติให้ชื่อ lm.Album

lm.Album<- lm(sales~adverts+airplay, #สมการ (sales =a+b1*adverts+b2*airplay)
    data=Album)            # ข้อมูล Album

การพิจาณา QQ plot ของค่า error ทำได้โดยใช้คำสั่ง

par(family = "TH Sarabun New") # เลือกใช้ font เพื่อแสดงภาษาไทยได้
qqnorm(lm.Album$residuals, pch = 1)
qqline(lm.Album$residuals, col = "steelblue", lwd = 2,
               main="QQ plot ของค่าคาดเคลื่อน")  # ชื่อกราฟ

ในกรณีที่ตัวแปรอิสระมีจำนวนมาก มีความจำเป็นต้องคัดเลือกตัวแปรอิสระเข้ามาในสมการการถดถอย โดยโปรแกรมอาร์แล้ว เราสามารถสร้างเงื่อนไขให้โปรแกรมคัดเลือกตัวแปรเข้าไปในสมการการถดถอยของเราได้ เช่น ข้อมูลชุด mtscars

str(mtcars)
## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

ถ้าเราสนใจ ตัวแปร (mpg) คือระยะทางที่รถเคลื่อนได้ต่อน้ำมัน 1 แกลลอน จะขึ้นอยู่กับปัจจัยใดบ้าง จะเห็นว่า มีตัวแปรอิสระทั้งหมด อีก 10 ตัว ชุดคำสั่งที่ต้องใช้เพิ่มเติมคือ MASS โดยเลือกสมการการถดถอยที่ดีที่สุดด้วยค่า AIC ที่น้่อยที่สุด

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
# ขั้นที่ 1 รันสมการด้วยคัวแปรอิสระทั้งหมด
full.model <- lm(mpg ~., data = mtcars)
# และรันสมการถดถอยเริ่มต้นสำหรับวิธี forward
initial.model <-lm(mpg ~wt, data = mtcars)

ในการคัดเลือกตัวแปรอิสระเข้าสมการการถดถอยมีทั้งหมด 3 วิธีคือ forward backward และ stepwise รายละเอียดเพิ่มเติมศึกษาได้จาก link ที่ให้ สำหรับคำสั่งในอาร์ จะเป็นวิธีึคัดเลือกแบบ stepwise ทั้งหมด

# Forward regression model
# กำหนดสมการเริ่มต้น และกำหนดสมการสุดท้าย
forward.model <- stepAIC(initial.model, #สมการเริ่มต้น
                         scope=list(lower=initial.model,upper=full.model),#สมการเริ่มต้น และสมการสุดท้าย
                      direction = "forward", 
                      trace = FALSE)
summary(forward.model)
## 
## Call:
## lm(formula = mpg ~ wt + cyl + hp, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9290 -1.5598 -0.5311  1.1850  5.8986 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 38.75179    1.78686  21.687  < 2e-16 ***
## wt          -3.16697    0.74058  -4.276 0.000199 ***
## cyl         -0.94162    0.55092  -1.709 0.098480 .  
## hp          -0.01804    0.01188  -1.519 0.140015    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.512 on 28 degrees of freedom
## Multiple R-squared:  0.8431, Adjusted R-squared:  0.8263 
## F-statistic: 50.17 on 3 and 28 DF,  p-value: 2.184e-11

ปัญหาของวิธีนี้ คือจะกำหนดให้สมการใดเป็นสมการเริ่มต้น ซึ่งมีทางเลือกมากกว่า

คำถาม ผู้อ่านลองเปลี่ยนสมการเริ่มต้น เป็นตัวแปรอื่นๆ ดู จะได้ผลเหมือนกันหรือไม่?

# Backaward regression model
# ใช้สมการที่มีตัวแปรอิสะมากที่สุกเป็นตัวเริ่มต้นได้เลย
backward.model <- stepAIC(full.model,
                      direction = "backward", 
                      trace = FALSE)
summary(backward.model)
## 
## Call:
## lm(formula = mpg ~ wt + qsec + am, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4811 -1.5555 -0.7257  1.4110  4.6610 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.6178     6.9596   1.382 0.177915    
## wt           -3.9165     0.7112  -5.507 6.95e-06 ***
## qsec          1.2259     0.2887   4.247 0.000216 ***
## am            2.9358     1.4109   2.081 0.046716 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.459 on 28 degrees of freedom
## Multiple R-squared:  0.8497, Adjusted R-squared:  0.8336 
## F-statistic: 52.75 on 3 and 28 DF,  p-value: 1.21e-11

จะเห็นวิธีการนี้จะง่ายกว่ามาก แต่ในทางเศรษฐศาสตร์ บ้างครั้งตัวแปรบ้างตัวก็มีความสำตัญมากในทางทฤษฎี จนไม่สามารถตัดออกจาสมการได้ เพราะฉะนั้นการคัดเลือกตัวแปรเข้ามาในสมการการถดถอย จำเป็นต้องเหตุผลอื่นๆ ประกอบด้วยนอกพิจารณาค่า AIC ที่ต่ำที่สุดเพียงอย่างเดียว

ถ้าต้องการพิจารณาเลือกสมการการถดถอยที่ดีที่สุดด้วยค่าอื่นๆ เช่น \(R^2\) หรือ adj. \(R^2\) ก็สามารถทำได้โดยใช้ชุดคำสั่ง caret และ leaps ในการทำงาน

library(MASS)
library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
## 
##     lift
library(leaps)

การใช้งานก็ไม่ต่างจากชุด คำสั่ง MASS มาก

Select.model <- train(mpg ~., # สมการการถดถอยที่ใส่อตัวแปรอิสระทั้งหมด 
                      data = mtcars, # ข้อมูล mtcars
                    method = "leapBackward",  # วิธีการตัดเลือกสมการถดถอยด้วยวิธี backward selection
                    tuneGrid = data.frame(nvmax = 1:10) #จำตัวแปรอิสระที่พิจารณาสูงสุด
                    )

ทำหรับคำสั่ง method นั้นสามารถเลือกใช้ได้อีก 2 แบบ คือ

เมื่อทำการคำนวณแล้วก็สามารถดูผลลัพธ์ของการคำนวณได้โดย

Select.model$results

จะเห็นได้ว่าจำนวนตัวแปรอิสระ 3 ตัวมีค่า \(R^2\) สูงที่สุด หรือพิจารณาจากคำสั่ง

Select.model$bestTune

คำสั่งต่อไป จะเป็นการบอกว่าตัว 3 ตัวใดที่ถูกคัดเลืิอก

summary(Select.model$finalModel)
## Subset selection object
## 10 Variables  (and intercept)
##      Forced in Forced out
## cyl      FALSE      FALSE
## disp     FALSE      FALSE
## hp       FALSE      FALSE
## drat     FALSE      FALSE
## wt       FALSE      FALSE
## qsec     FALSE      FALSE
## vs       FALSE      FALSE
## am       FALSE      FALSE
## gear     FALSE      FALSE
## carb     FALSE      FALSE
## 1 subsets of each size up to 3
## Selection Algorithm: backward
##          cyl disp hp  drat wt  qsec vs  am  gear carb
## 1  ( 1 ) " " " "  " " " "  "*" " "  " " " " " "  " " 
## 2  ( 1 ) " " " "  " " " "  "*" "*"  " " " " " "  " " 
## 3  ( 1 ) " " " "  " " " "  "*" "*"  " " "*" " "  " "

ก็คือ wt qsec และ am หลังจากนั้นก็ใช้คำสั่ง lm อีกครั้งก็จะได้สมการการถดถอยที่ดีที่สุดด้วยเกณฑ์ \(R^2\) ตามต้องการ จากวิธีการคัดเลือกสมการจากทั้งสอง package จะให้ผลไม่เหมือนกันเพราะเกณฑ์ในการคัดสิน มีความแตกต่างกันคือ ค่า AIC และค่า \(R^2\)

การแก้ปัญหาเรื่องข้อมูลสูญหายเบื้องต้น

ในบ้างครั้งข้อมูลที่เก็บมานั้นมีการสูญหาย เช่น

Data <-data.frame(x= c(1,2,3,NA,5),
                  y= c(2,4.35,6.1,8.1,9.8))
Data

จากตัวอย่างจะเห็น กรอบข้อมูล Data ตัวแปร x สูญหายไม่ 1 ค่า ซึ่งอาจเกิดขึ้นจากหลายสาเหตุ เช่น ลืมกรอก หรือไม่ทราบ หรือ ไม่สามารถสังเกตุ จึงไม่ได้กรอกไว้ จากตัวอย่าง จะมีเพียงค่าเดียวที่สูญหาย ค่าสมมุติว่าข้อมูลมีจำนวนหลักร้อยหรือหลักพันขึ้นไปการแก้ไขด้่วย มืออาจทำให้เสียเวลามาก ดังนั้นโปรแกรมดึงคำสั่งเพื่อจัดการปัญหานี้ คือ

na.omit(Data)

ถ้าในแถวใด มีที่ข้อมูลของบางตัวแปรสูญหาย คำสั่งนี้จะก็ตัดข้อมูลทิ้งไปแถว อย่าลืมสร้างตัวแปรใหม่ขึ้นมา หรือตัวใช้คำสั่งนี้ไปกำหนดค่าให้ตัวแปรเดิมก็ได้

หรือถ้าต้องการใช้วิธีแทนค่าที่สุญหายไปตัวเลขที่กำหนดขึ้นก็สามารถทำได้ ในกรณีที่มีจำนวนตัวแปรไม่มาก ก็สามารถทำได้โดย

Data[is.na(Data)] <- 0
Data

ในกรณีที่ข้อมูลจำตัวแปรมาก ตั้งแต่ 10 ตัวขึ้นไป ในการวิธีเขียนคำสั่งให้โปรแกรมอาร์แก้ไขโดยอัตโนมัติจะสะดวกแล้ว ซึ่งจะไม่พูดถึงในบทนี้

## บทอื่นๆ