การใช้ R language ในการคำนวณค่าทางสถิติ

ปราโมทย์ คูวิจิตรจารุ
เพิ่มเติมล่าสุด 18 เมษายน 2548


เป็น programming language สำหรับการคำนวณทางสถิติและแสดงกราฟฟิกโดยเฉพาะ ได้รับการพัฒนามาเพื่อทดแทน S language ที่พัฒนามาเป็นซอฟท์แวร์ทางการค้า S-PLUS (บริษัท Insightful) S-Plus เป็นโปรแกรมคำนวณทางสถิติที่ได้รับความนิยมเป็นอย่างมากอันหนึ่ง โดยภาษา S ได้รับการพัฒนาในช่วงปลายทศวรรษที่ 80 โดย AT&T lab ส่วนภาษา R นั้นเริ่มได้รับการพัฒนาโดย Robert Gentleman และ Ross Ihaka จากภาควิชาสถิติ มหาวิทยาลัยโอคแลนด์ ประเทศนิวซีแลนด์ในปี 1995 ในปัจจุบันอยู่ภายใต้การดูแลของทีมพัฒนาซึ่งเป็นอาสาสมัครจากทั่วโลก R พัฒนาขึ้นมาแบบ open-source ทำให้สามารถดาวโหลดโปรแกรมมาใช้ได้"ฟรี"ซึ่งมีอยู่บนหลายแพล็ตฟอร์มเช่น Unix, Linux, MS Window ซึ่งล่าสุดได้รับการพัฒนาถึงเวอร์ชั่น 2.0.1 แล้ว (stable version)โดย R language มาพร้อมกับเอกสารแนะนำเบื้องต้น และเอกสาร Help นอกจากนี้ยังสามารถอ่านเอกสารเฉพาะทางซึ่งเขียนโดยผู้เชี่ยวชาญในด้านนั้นๆได้จากเว็บไซต์ด้านบนเช่นกัน นอกจากนี้แล้วแหล่งข้อมูลที่สำคัญอีกที่คือ R-help mailing list ซึ่งสามารถค้นหาคำถามและคำตอบหรือจะสมัครเพื่อรับและส่งอีเมล์กับกลุ่มดังกล่าวก็ได้เช่นกัน และโดยทั่วไป R ถูกพัฒนาขึ้นมาทดแทน S ดังนั้นเอกสารของ S-PLUS จึงสามารถใช้ได้กับ R ได้เช่นกัน

ในปัจจุบันมีหนังสือที่เขียนแนะนำการใช้ R ในการคำนวณสถิติมากขึ้นเรื่อยๆ ดังเช่น

ตัวเลขและตัวอย่างที่ใช้ในการแนะนำในเอกสารนี้นำมาจาก

สารบัญ

R บน MS Windows

ในเอกสารนี้จะใช้ R เวอร์ชั่นบนแพล็ตฟอร์ม Microsoft Windows โดย R จะมาพร้อมกับ RGui ซึ่งเมื่อคลิ้กที่โปรแกรมก็จะปรากฏหน้าต่าง R console ขึ้นมาอย่างในภาพด้านล่าง และต่อไปในส่วนที่เป็นโค้ดหรือผลที่ได้จากการรันใน R จะแสดงอยู่บนพื้นหลังสีดำเพื่อให้ต่างกับเนื้อหาอื่น

จากเมนู Help > Html help จะทำการเปิดเอกสารที่มากับโปรแกรมด้วยเว็บบราวเซอร์ เช่นเอกสารแนะนำการใช้เบื้องต้น หรือค้นหาคำสั่งและรูปแบบการใช้งานและคำอธิบายพร้อมตัวอย่าง

จุดที่เห็นเครื่องหมาย > นั้นคือหมายถึงว่าเป็นจุดสำหรับการป้อนคำสั่งต่างๆ เริ่มต้นเราจะดูการใช้ built-in function สำหรับแสดงวิธีการเขียนอ้างอิงหากนำ R มาใช้ในการคำนวณ โดยใช้คำสั่ง citation()

> citation()
To cite R in publications, use

  R Development Core Team (2003). R: A language and environment for
  statistical computing. R Foundation for Statistical Computing,
  Vienna, Austria. ISBN 3-900051-00-3, URL http://www.R-project.org.

คำนวณค่าเฉลี่ยและค่าเบี่ยงเบนมาตรฐาน

คำนวนค่าเฉลี่ยและค่าเบี่ยงเบนมาตรฐานของผลการไตเตรตดังนี้ 10.08, 10.11, 10.09, 10.10, 10.12 mL

> titrant <-c(10.08,10.11,10.09,10.10,10.12)
> mean(titrant)
[1] 10.1
> sd(titrant)
[1] 0.01581139
> titrant2 =c(10.08,10.11,10.09,10.10,10.12)
> mean(titrant2)
[1] 10.1
> sd(titrant2)
[1] 0.01581139

จะเห็นว่าสามารถทำได้อย่างง่ายๆ ด้วยการกรอกข้อมูลลงไปพร้อมกับสร้างตัวแปรขึ้นมาเพื่อเก็บค่าไว้ ในที่นี่คือ titrant โดยใช้รูปแบบการกรอกดังตัวอย่างโดยใช้ฟังก์ชั่น c() หลังจากนั้นจึงใชัคำสั่ง built-in สองอันคือ mean() กับ sd() เพื่อคำนวณค่าเฉลี่ย (arithmatic mean) และ ค่าเบี่ยงเบนมาตรฐาน (standard deviation) ในตัวอย่างด้านล่างเป็นการกรอกข้อมูลในอีกรูปแบบหนึ่งให้เก็บอยู่ในชื่อ titrant2 และรูปแบบการสั่งโดยเปลี่ยนชื่อตัวแปร ซึ่งให้ผลเหมือนกัน

การแก้ไขข้อมูลที่กรอกลงไปแล้วสามารถจะทำได้ด้วยคำสั่ง edit() ตามรูปแบบในตัวอย่างต่อไปนี้

> x<-c(55,57,59,56,56,59)
> x
[1] 55 57 59 56 56 59
> x=edit(x)
> x
[1] 55 57 59 56 56 59 70

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

การเก็บข้อมูลที่ทำการกรอกไปแล้วเพื่อคำนวณในภายหลัง

เมื่อสั่ง exit โปรแกรม จะมีการถามให้ save workspace image ให้ตอบ yes เพื่อจะได้สามารถใช้ข้อมูลที่กรอกลงไปทั้งหมดในภายหลังได้ หรือจากเมนู File >Save workspace เพื่อเก็บไฟล์ไว้ในชื่ออื่นหรือเก็บลงแผ่นดิสค์เพื่อนำไปใช้ในคอมพิวเตอร์เครื่องอื่น ไฟล์ที่เซฟจะมีนามสกุล *.RData

หรือหากต้องการเก็บเพียงข้อมูลบางอันไว้ อาจจะทำได้โดยใช้คำสั่ง save() เช่นทำการเซฟออปเจ็ค titrant ไว้ในไฟล์ชื่อ x.RData ด้วยคำสั่ง

> save(titrant, file="titrant.RData")

หรือเซฟหลายๆออปเจ็คไว้ในไฟล์เดียวกัน

> x<- c(1,2,3)
> y<- c(4,5,6)
> save(x,y,file="xy.RData")

เมื่อทำการเปิดโปรแกรมขึ้นมาใหม่ โดยที่เราไม่ได้ทำการเซฟ workspace หรือทำงานบนเครื่องอื่นก็สามารถโหลดข้อมูลเก่ามาใช้ได้อีกด้วยคำสั่ง load()

> load("xy.RData")

ไฟล์ที่ทำการเซฟจะถูกบันทึกลงใน working directory ซึ่งสามารถเปลี่ยนได้โดยเรียกเมนู File > Chang dir...

สร้างฮิสโตแกรม

เราต้องการสร้างฮิสโตแกรมสำหรับความเข้มข้นของ nitrate จำนวน 50 ค่า ซึ่งการกรอกข้อมูลจำนวนมากนี้สามารถทำได้สะดวกมากขึ้นด้วยการกรอกผ่าน spreadsheet ด้วยการสั่ง

> data.entry(y=c(NA))

ก็จะเห็นหน้าต่าง spreadsheet ปรากฏขึ้นมา โดยที่หัวของคอลัมน์แรกทางซ้ายนั้นจะชื่อ c(NA) ให้คลิกแล้วทำการเปลี่ยนชื่อของเวคเตอร์ข้อมูลเช่นในกรณีนี้ใช้ nitrate หลังจากนั้นทำการกรอกตัวเลขทั้ง 50 ค่าลงในคอลัมน์ที่หนึ่งไล่ลงไปทีละแถว

เมื่อกรอกครบแล้วให้คลิกขวา เลือก close เพื่อปิด spreadsheet และ R จะทำการเก็บข้อมูลให้โดยอัตโนมัติ เมื่อกลับมาสู่หน้าต่างหลักแล้วให้ลองเรียกดูข้อมูล ก็จะได้ดังนี้

> nitrate
 [1] 0.51 0.51 0.49 0.51 0.51 0.51 0.52 0.48 0.51 0.50 0.51 0.53 0.46 0.51 0.50
[16] 0.50 0.48 0.49 0.48 0.53 0.51 0.49 0.49 0.50 0.52 0.49 0.50 0.48 0.47 0.52
[31] 0.52 0.52 0.49 0.50 0.50 0.53 0.49 0.49 0.51 0.50 0.50 0.49 0.51 0.49 0.51
[46] 0.47 0.50 0.47 0.48 0.51

ทำการสร้างฮิสโตแกรมด้วยคำสั่ง

> hist(nitrate)

หรือทำการเพิ่มอ็อปชั่นให้กับคำสั่ง hist() เช่นเพิ่มเลเบลให้กับแกน x

> hist(nitrate, xlab="Nitrate ion concentration, ug/mL")

ก็จะเห็นหน้าต่างกราฟฟิกแสดงขึ้นมาดังภาพ

การกระจายแบบปรกติ

ผลการไตเตรตซ้ำมีการกระจายแบบปรกติ โดยมีค่าเฉลี่ย 10.15 mL และค่าเบี่ยงเบนมาตรฐาน 0.02 mL จากหาสัดส่วนของผลการวัดที่อยู่ในช่วง 10.12 ถึง 10.29 mL

> pnorm(10.20,10.15,.02)
[1] 0.9937903
> pnorm(10.12,10.15,.02)
[1] 0.0668072
> .9937903-.0668072
[1] 0.9269831

pnorm() ใช้หา cumulative distribution function, P(X<=x) เมื่อ X มีการกระจายแบบปรกติ สัดส่วนของค่าที่อยู่ในช่วง 10.12-10.20 คือ 0.927

คำนวณช่วงความเชื่อมั่นของค่าเฉลี่ยเมื่อมีจำนวนตัวอย่างขนาดใหญ่

คำนวณช่วงความเชื่อมั่นของค่าเฉลี่ยของความเข้มข้นไนเตรตจากตัวอย่างข้างบน (ขนาดตัวอย่างคือ 50) จากสูตร ช่วงความเชื่อมั่นคือ mean ± z*sigma/sqrt(n) โดยค่า z คำนวณได้สำหรับที่ความเชื่อมั่นต่างๆ ด้วย qnorm() นั่นคือค่า inverse ของ cumulative distribution function

#z for 95%, 99% confidence limit, respectively 
> qnorm(.975)
[1] 1.959964
> qnorm(.995)
[1] 2.575829

> mean(nitrate)
[1] 0.4998
> qnorm(.975)*sd(nitrate)/sqrt(50)
[1] 0.004566235
> qnorm(.995)*sd(nitrate)/sqrt(50)
[1] 0.006001049

ช่วงระดับความเชื่อมั่นจะเป็น 0.500 ±0.0046 ug/mLและ 0.500 ±0.0060 ug/mLที่ 95 และ 99% ตามลำดับ (เมื่อทำการปัดค่าตัวเลข)

คำนวณช่วงความเชื่อมั่นของค่าเฉลี่ยเมื่อมีจำนวนตัวอย่างขนาดเล็ก

ในกรณีที่มีขนาดตัวอย่างมีขนาดเล็กเราจะใช้การคำนวณด้วย t-test แทน z-test โดย R มีฟังก์ชั่น t.test() มาให้แล้ว ตัวอย่างเช่นค่าโซเดียมในตัวอย่างปัสสะวะตรวจได้เป็น 102, 97, 99, 98, 101, 106 mM คำนวณช่วงความเชื่อมั่นของความเข้มข้นของโซเดียมที่ 95 และ 99% ได้ดังนี้

> Na=c(102,97,99,98,101,106)
> t.test(Na)

        One Sample t-test

data:  Na 
t = 75.2575, df = 5, p-value = 7.848e-09
alternative hypothesis: true mean is not equal to 0 
95 percent confidence interval:
  97.0672 103.9328 
sample estimates:
mean of x 
    100.5 

นั่นคือช่วงความเชื่อมั่นที่ 95% คือ 100.5±3.4

> t.test(Na,conf.level=.99)

        One Sample t-test

data:  Na 
t = 75.2575, df = 5, p-value = 7.848e-09
alternative hypothesis: true mean is not equal to 0 
99 percent confidence interval:
  95.11542 105.88458 
sample estimates:
mean of x 
    100.5 

นั่นคือช่วงความเชื่อมั่นที่ 99% คือ 100.5±5.4 โปรดสังเกตุว่า ค่า conf.level คือระดับความเชื่อมั่นซึ่งค่า defualt คือ 0.95 (95%)

คำนวน t-test เปรียบเทียบค่าเฉลี่ยของตัวอย่างกับค่าจริง

ในการตรวจสอบหา selenourea ในน้ำ ได้ค่าจากน้ำประปาที่ทำการเติม 50 ng/mL selenourea ลงไป ด้งนี้ 50.4, 50.7, 49.1, 49.0, 51.1 ng/mL จงหาว่าผลที่ได้แสดงให้เห็นถึง systematic error หรือไม่

> x<-c(50.4,50.7,49.1,49.0,51.1)
> t.test(x,mu=50.0,conf.level=0.95)

        One Sample t-test

data:  x 
t = 0.1404, df = 4, p-value = 0.8951
alternative hypothesis: true mean is not equal to 50 
95 percent confidence interval:
 48.87358 51.24642 
sample estimates:
mean of x 
    50.06 

ค่า t ที่คำนวณได้คือ 0.1404 โดย P(|t|>0.1404) =0.891 ซึ่งค่า P นี้มีค่ามากกว่า 0.05 นั่นคือผลนั้นไม่มีนัยสำคัญที่ P=0.05 และยอมรับ null hypothesis ได้ นั่นคือไม่แสดงให้เห็นว่ามี systematic error เกิดขึ้น นอกจากนี้การคำนวนยังบอกค่าช่วงความเชื่อมั่นของค่าเฉลี่ยที่ 95% ให้เป็น 48.87-51.25 ng/mL

คำนวน t-test เปรียบเทียบค่าเฉลี่ยจากการทดลองสองอันเมื่อข้อมูลมีการกระจายเท่ากัน

ในการตรวจหาปริมาณตะกั่วในอาหาร ได้ทดลองต้มรีฟลักตัวอย่างในกรดไฮโดรคลอริกด้วยเวลาต่างกัน และผลที่ได้เป็นดังนี้
เวลา 30 นาที ได้ผลความเข้มข้นของตะกั่วเป็น 55, 57, 59, 56, 56, 59 mg/kg
เวลา 75 นาที ได้ผลความเข้มข้นของตะกั่วเป็น 57, 55, 58, 59, 59, 59 mg/kg
ค่าเฉลี่ยของปริมาณตะกั่วจากสองอุณหภูมิมีความแตกต่างกันอย่างมีนัยสำคัญหรือไม่

> x<-c(55,57,59,56,56,59)
> y<-c(57,55,58,59,59,59)
> var(x);var(y)
[1] 2.8
[1] 2.566667

ทำการคำนวณโดย assume ว่าสองตัวอย่างมีการกระจายข้อมูลเท่ากัน

> t.test(x,y,var.equal=TRUE)

        Two Sample t-test

data:  x and y 
t = -0.8811, df = 10, p-value = 0.3989
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -2.940597  1.273931 
sample estimates:
mean of x mean of y 
 57.00000  57.83333

P(|t|>0.88)=0.399 เนื่องจากค่าความน่าจะเป็นนี้สูงกว่า 0.05 ดังนั้นผลไม่มีนัยสำคัญที่ระดับ 5% หรือนั่นคือไม่มีหลักฐานยืนยันได้ว่าระยะเวลาในการต้มมีผลต่อปริมาณของตะกั่วที่ตรวจวัดได้

ในกรณีที่เราไม่มีข้อมูลดิบสำหรับแต่ละตัวอย่าง แต่มีเพียงค่าสรุปเช่นตัวอย่างเปรียบเทียบสองวิธีการในการตรวจหาโครเมียมในตัวอย่างได้ผลดังนี้ วิธีแรกได้ค่าเฉลี่ย =1.48; sd = 0.28 ส่วนวิธีที่สองได้ค่าเฉลี่ย=2.33; sd = 0.31 ทั้งสองมีขนาดตัวอย่าง= 5 ให้คำนวณว่าทั้งสองให้ค่าเฉลี่ยต่างกันอย่างมีนัยสำคัญหรือไม่

ในกรณีนี้มีค่าการกระจายไม่ต่างกันมากและเราจะใช้การคำนวณด้วย t-test โดยสูตรการคำนวณสามารถดูได้จากตำราสถิติทั่วไป ใน R คำนวณได้ดังนี้

> pooled.s=sqrt(((4*.28^2)+(4*.31^2))/(5+5-2))
> pooled.s
[1] 0.2953811
> tstat<-(2.33-1.48)/(.295*sqrt(1/5+1/5))
> tstat
[1] 4.555824

เมื่อเราคำนวณค่า statistic t แล้วเราก็จะทำการเปรียบเทียบค่านี้กับค่า critical value t (df=8) ซึ่งปรกติสามารถเปิดดูตาราง t-distribution ได้หรือสามารถคำนวณด้วย R ได้เช่นกัน เช่นที่ (P=0.05) ใน R ด้วยคำสั่ง qt() จะให้ค่า t สำหรับ one-tailed test ดังนั้นในกรณีนี้เราต้องการสำหรับ two-taied test นั่นคือที่ P=0.05 เราต้องคำนวณโดยการใช้ค่า P=0.025 แทน

> qt(.975,8)
[1] 2.306004
> qt(.995,8)
[1] 3.355387

ดังนั้นค่า t จากการคำนวณ (4.56) มีค่ามากกว่า critical value (2.31) ดังนั้นค่าเฉลี่ยจากสองวิธีมีความแตกต่างกันอย่างมีนัยสำคัญ ที่ 5% level และจริงแล้วค่า critical value สำหรับ t ที่ 1% จะมีค่าประมาณ 3.36 นั่นคือความแตกต่างมีนัยสำคัญที่ 1% level

หรือเราสามารถใช้ t.test() โดยตรงได้โดยการสร้างข้อมูลดิบหลอกขึ้นมาโดยใช้คำสั่ง rnorm() ซึ่งเป็นการสร้างตัวเลขสุ่มสำหรับการกระจายแบบ Normal โดยหากต้องการสร้างตัวอย่างแบบเทียม (psudo random sample) จำนวน n ค่าที่มี่ค่าเฉลี่ย mu และ ค่าเบี่ยงเบนมาตรฐานเป็น sigma เราสามารถสั่งใน R ได้ดังนี้ x<-mu+sigma*scale(rnorm(n))(คำแนะนำจาก R-help mailing list)

> a<-1.48+0.28*scale(rnorm(5))
> b<-2.33+0.31*scale(rnorm(5))
> t.test(a,b,var.equal=T)

        Two Sample t-test

data:  a and b 
t = -4.5499, df = 8, p-value = 0.001875
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -1.2807971 -0.4192029 
sample estimates:
mean of x mean of y 
     1.48      2.33

พบว่าค่า P(|t| > 4.55)=0.002 ซึ่งต่ำว่า 0.01 นั่นคือความแตกต่างของทั้งสองวิธีมีนัยสำคัญที่ระดับ 1% level

คำนวน t-test เปรียบเทียบค่าเฉลี่ยจากการทดลองสองอันเมื่อข้อมูลมีการกระจายไม่เท่ากัน

ในกรณีที่สองตัวอย่างมีค่าการกระจายข้อมูลต่างกันมากเราไม่สามารถจะ assume ดังตัวอย่างด้านบน การคำนวณใน R จะทำได้ดังตัวอย่างต่อไปนี้

ทำการตรวจปริมาณไธออล (mM) ตัวอย่างเลือดจากคนสองกลุ่มคือกลุ่ม normal และกลุ่มที่เป็น rheumatoid ได้ผลดังนี้
Normal: 1.84, 1.92, 1.94, 1.92, 1.85, 1.91, 2.07
Rheumatoid: 2.81, 4.06, 3.62, 3.27, 3.27, 3.76

> normal<-c(1.84,1.92,1.94,1.92,1.85,1.91,2.07)
> rheumatoid<-c(2.81,4.06,3.62,3.27,3.27,3.76)
> sd(normal);sd(rheumatoid)
[1] 0.0755929
[1] 0.4404884
> t.test(normal,rheumatoid)

        Welch Two Sample t-test

data:  normal and rheumatoid 
t = -8.4772, df = 5.253, p-value = 0.0002937
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -2.004938 -1.082205 
sample estimates:
mean of x mean of y 
 1.921429  3.465000 

ในกรณีที่ค่าการกระจายของสองตัวอย่างมีค่าไม่เท่ากัน R จะเลือกใช้ Welch's t-test ในกรณีนี้สรุปได้ว่า P(|t| >8.48 = 0.0003 ซึ่งมีค่าต่ำมากๆ นั่นคือแสดงให้เห็นว่าค่าเฉลี่ยของความเข้มข้นของไธออลมีความแตกต่างกันในสองกลุ่มอย่างมีนัยสำคัญถึงที่ P=0.001

Paired t-test

ทำการตรวจสอบปริมาณพาราเซตามอล(%)ในเม็ดยาด้วยสองวิธี (a,b) โดยเก็บตัวอย่างเม็ดยา 10 เม็ดจาก batch ต่างๆ ทำการวิเคราะห์เพื่อดูว่าวิธีการทั้งสองให้ผลต่างกันหรือไม่
วิธี a ให้ผลสำหรับ Batch ที่ 1-10 ดังนี้ 84.63,84.38,84.08,84.41,83.82,83.55,83.92,83.69,84.06,84.03
วิธี b ให้ผลสำหรับ Batch ที่ 1-10 ดังนี้ 83.15,83.72,83.84,84.20,83.92,84.16,84.02,83.60,84.13,84.24

ในกรณีนี้เราจะใช้การเปรียบเทียบด้วยวิธี paired t-test ได้ดังนี้

> a<-c(84.63,84.38,84.08,84.41,83.82,83.55,83.92,83.69,84.06,84.03)
> b<-c(83.15,83.72,83.84,84.20,83.92,84.16,84.02,83.60,84.13,84.24)
> t.test(a,b,paired=T)

        Paired t-test

data:  a and b 
t = 0.8821, df = 9, p-value = 0.4007
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.2487527  0.5667527 
sample estimates:
mean of the differences 
                  0.159 

P(|t|>0.882)=0.40 ซึ่งค่าความน่าจะเป็นนี้สูงว่า 0.05 ดังนั้นสองวิธีไม่มีความแตกต่างกันอย่างมีนัยสำคัญที่ P=0.05

One-sided test

มีความสงสัยว่าในการไตเตรตมีความผิดพลาดของตัวอินดิเคเตอร์และทำให้ได้ผลมากกว่าเป็นจริง นั่นคือมี positive bias จึงทำการตรวจสอบด้วยการใช้สารละลานกรดเข้มข้น 0.1 M พอดีไตเตรตสารละลายด่างเข้มข้น 0.1 M พอดีปริมาตร 25 mL ได้ผลดังนี้ (mL): 25.06, 25.18, 24.87, 25.51, 25.34, 25.41 เราสามารถตรวจสอบ positive bias ของผลนี้ได้ดังนี้

> titrant<-c(25.06,25.18,24.87,25.51,25.34,25.41)

> t.test(titrant,mu=25.00,var.equal=T,alternative=c("greater"))

        One Sample t-test

data:  titrant 
t = 2.3473, df = 5, p-value = 0.03289
alternative hypothesis: true mean is greater than 25 
95 percent confidence interval:
 25.03232      Inf 
sample estimates:
mean of x 
 25.22833 

ค่า P(t > 2.35)=0.03 ซึ่งน้อยกว่า 0.05 ดังนั้น null hypothesis ถูกปฏิเสธนั่นคือมี positive bias เกิดขึ้น

F-test สำหรับเปรียบเทียบค่าการเบี่ยงเบนมาตรฐาน

การหาค่า critical value ที่ significant level=alpha สำหรับ F distribution โดย df1 เป็นดีกรีอิสระของตัวเศษและ df2 เห็นดีกรีอิสระของตัวหาร

เปรียบเทียบวิธีการตรวจสอบ COD ในน้ำทิ้งแบบ standard กับวิธีใหม่ ได้ผลดังนี้
standard method : mean=72, sd=3.31 mg/L
proposed method : mean=72, sd = 1.51 mg/L
ทั้งสองวิธีทำการวัด 8 ครั้ง ต้องการทราบว่าความแม่นยำ (precision) ของวิธีที่เสนอใหม่นี้ดีกว่าวิธีมาตรฐานอย่างมีนัยสำคัญหรือไม่ ในกรณีนี้เราจะใช้การทดสอบสำหรับ One-sided

> F.stat<-(3.31^2)/(1.51^2)
> F.stat
[1] 4.805096
> qf(.05,7,7,lower.tail=FALSE)
[1] 3.787044

ค่าคำนวณ F (4.80) มีค่ามากกว่า critical value =3.787 (P=0.05) นั่นคือค่าการกระจายของวิธีมาตรฐานมีค่ามากกว่าวิธีเสนอใหม่อย่างมีนัยสำคัญที่ 5% level

ในการทดสอบ t-test ด้านบนเรา assume ว่าค่าการกระขายของทั้งสองวิธีในการตรวจหาโครเมียมนั้นไม่ต่างกัน เราสามารถทดสอบ assumption นี้ได้ ค่า sd ของสองวิธีซึ่งได้จากการวัด 5 ครั้งเป็น 0.28 และ 0.31 ค่าเฉลี่ยเป็น 1.48 และ 2.33

> F=(.31^2)/(.28^2)
> F
[1] 1.225765
> qf(c(.025,1-.025),4,4)
[1] 0.1041175 9.6045299

ค่า F (1.23)คำนวณมีค่าน้อยกว่ากว่า critical value (9.60) นั่นคือไม่มีความแตกต่างของค่าการกระจายจากสองวิธีอย่างมีนัยสำคัญที่ 5% level

การเปรียบเทียบค่าเฉลี่ยหลายค่า และการทำ box plot

ผลการเก็บสารฟลูออเรสเซนต์ที่สภาวะต่าง a, b, c, และ d ให้ค่า signal 0tในการการวัดแต่ละตัวอย่าง 3 ซ้ำดังนี้
A: 102, 100, 101
B: 101, 101, 104
C: 97, 95, 99
D: 90, 92, 94
เราสามารถใช้ ANOVA test ทดสอบว่าความแตกต่างของค่าเฉลี่ยที่ได้จากแต่ละสภาวะการเก็บมีความแตกต่างกันมากเกินกว่าจะเป็นผลจาก random error หรือไม่

> a<-c(102,100,101)
> b<-c(101,101,104)
> c<-c(97,95,99)
> d<-c(90,92,94)

ทำการวาด box-plot เพื่อดูข้อมูลแบบคร่าวๆ

> signal=data.frame(a,b,c,d)
> signal
    a   b  c  d
1 102 101 97 90
2 100 101 95 92
3 101 104 99 94
> boxplot(signal, xlab="condition",ylab="Signal")

data.frame() เป็นคำสั่งในการสร้าง data frame ซึ่งเป็นรูปแบบการเก็บข้อมูลในรูปแถวและคอลัมน์อย่างที่แสดงในตัวอย่างด้านบน แต่ละคอลัมน์คือข้อมูลของสภาวะที่ต่างกัน และแต่ละแถวเป็นข้อมูลแต่ละซ้ำของการวัด หลังจากนั้นเราใช้ boxplot() พร้อมกับอ็อปชั่นในการกำหนด label แกน x และ y ตามลำดับเพื่อสร้างกราฟ

ในการทำการวิเคราะห์ ANOVA จะใช้ฟังก์ชั่น oneway.test() ซึ่งรับข้อมูลที่อยู่ในลักษณะที่มีหนึ่งตัวแปรเก็บค่า signal และอีกหนึ่งตัวแปรที่จะระบุกลุ่ม เราสามารถสร้างข้อมูลแบบนี้ได้ด้วยคำสั่ง stack()

> signal.stack=stack(signal)
> signal.stack
   values ind
1     102   a
2     100   a
3     101   a
4     101   b
5     101   b
6     104   b
7      97   c
8      95   c
9      99   c
10     90   d
11     92   d
12     94   d

นั่นคือเราได้ข้อมูลที่ตัวแปรชื่อ values เก็บค่า และ ind บอกหมู่ของค่า

> oneway.test(values ~ ind, data=signal.stack,var.equal=T)

        One-way analysis of means

data:  values and ind 
F = 20.6667, num df = 3, denom df = 8, p-value = 0.0004002

ค่า p-value 0.0004 หมายความว่า null hypothesis ที่ว่าค่าเฉลี่ยมีค่าเท่ากันถูกปฏิเสธ นอกจากนี้เรายังสามารถใช้ฟังก์ชั่น anova() เรียกผลของฟังก์ชั่น lm() ซึ่งในกรณีจะได้รายละเอียดของการคำนวณมากกว่า oneway.test()

> anova(lm(values ~ ind,data=signal.stack))
Analysis of Variance Table

Response: values
          Df Sum Sq Mean Sq F value    Pr(>F)    
ind        3    186      62  20.667 0.0004002 ***
Residuals  8     24       3                      
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

 

การจำแนกและประมาณค่า variance ด้วย ANOVA

ในการทดสอบความบริสุทธิ์ของเกลือจากตัวอย่าง a-e ซึ่งเก็บจากถัง ทำการวัดแต่ละตัวอย่าง 4 ซ้ำ

     a    b    c    d
1 98.8 99.3 98.3 98.0
2 98.7 98.7 98.5 97.7
3 98.9 98.8 98.8 97.4
4 98.8 99.2 98.8 97.3
ทำการทดสอบด้วย ANOVA ดังนี้
> a<-c(98.8,98.7,98.9,98.8)
> b<-c(99.3,98.7,98.8,99.2)
> c<-c(98.3,98.5,98.8,98.8)
> d<-c(98.0,97.7,97.4,97.3)
> e<-c(99.3,99.4,99.9,99.4)
> purity=data.frame(a,b,c,d)
> purity.stack=stack(purity)
> anova(lm(values~ind,data=purity.stack))
Analysis of Variance Table

Response: values
          Df Sum Sq Mean Sq F value    Pr(>F)    
ind        4 7.8400  1.9600      30 5.343e-07 ***
Residuals 15 0.9800  0.0653                      
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

ผลการทดสอบ between-sample Mean sq มีค่ามากกว่า winthin-sample Mean sq และผล F-test พบว่ามีความแตกต่างอย่างมีนัยสำคัญ นั่นคือ sampling variance แตกต่างจาก 0 อย่างมีนัยสำคัญ และเราสามารถประมาณค่าได้จาก var=[(between-sample Mean sq) - (within-sample Mean sq)]/n = (1.96-0.0653)/4 = 0.47

การทดสอบความปรกติของการกระจาย (Normality of distribution)

ใน R มีฟังก์ชั่นสำหรับสร้าง normal probability plot ง่ายๆ คือ qqnorm() และ qqline() ซึ่งลากเส้นอ้างอิง (เส้นที่ลากผ่านจุด q1 และ q3)

> measurement<-c(109,89,99,99,107,111,86,74,115,107,134,113,110,88,104)
> qqnorm(measurement);qqline(measurement)

การสร้างกราฟมาตรฐาน

การคำนวณ linear regression เพื่อหาเส้นตรงความสัมพันธ์ระหว่างข้อมูลสองชุดด้วยวิธี least squares method อย่างเช่นการเตรียม calibration curve ในการวิเคราะห์ทางเคมี สามารถทำได้ดังตัวอย่างการวัดความเข้มของสี (intensity) กับความเข้มข้น(conc)ของสารดังนี้

> conc=c(0,2,4,6,8,10,12)
> intensity=c(2.1,5.0,9.0,12.6,17.3,21,24.7)
> model=lm(intensity~conc)
> model

Call:
lm(formula = intensity ~ conc)

Coefficients:
(Intercept)         conc  
      1.518        1.930  

นั่นคือเราได้สัมประสิทธิ์ตัวแรกพจน์แรกซึ่งเป็นจุดตัดแกน Y เป็น 1.518 และสัมประสิทธิ์ตัวที่สองคือความชันของเส้นกราฟเป็น 1.930 จากนั้นเราสามารถ Plot กราฟเส้นตรงนี้ได้ด้วยคำสั่ง plot ()เพื่อแสดงจุดของข้อมูล แล้วใช้ abline() เพื่อแสดงเส้นตรงจากการคำนวณ จากนั้นแสดงสมการ y=a+bx ได้โดยการสร้างตัวแปรเช่น txt แล้วทำการแปะข้อความลงกราฟด้วยคำสั่ง legend()

> plot(conc,intensity)
> abline(model$coef)
> txt<-paste("Y= ", format(model$coef[1],digits=3)," + ",format(model$coef[2],digits=3),"x")
> legend(5,10,txt)

การคำนวณหาสัมประสิทธิ์สำหรับ first degree polynomial

ในกรณีนี้เราสนใจจะ fit ข้อมูลระหว่าง y กับตัวแปร 3 ตัวคือ x1,x2,x3 ซึ่งอยู่ในรูปรหัส ด้วยสมการ polynomial อันดับหนึ่ง นั่นคือในรูป y=b0 + b1x1 + b2x2 + b3x3

> x1<-c(-1,1,-1,1,-1,1,-1,1)
> x2<-c(-1,-1,1,1,-1,-1,1,1)
> x3<-c(-1,-1,-1,-1,1,1,1,1)
> y<-c(2.83,3.56,2.23,3.06,2.47,3.30,1.95,2.56)
> lm(y~x1+x2+x3)

Call:
lm(formula = y ~ x1 + x2 + x3)

Coefficients:
(Intercept)           x1           x2           x3  
      2.745        0.375       -0.295       -0.175  

นั่นคือจะได้ว่าสัมประสิทธิ์ b0 = 2.745, b1 = 0.375, b2 = -0.295, b3 = -0.175

การคำนวณหาสัมประสิทธิ์สำหรับ polynomial ดีกรีสูงขึ้นไป

ในกรณีที่เราต้องการ fit ข้อมูลกับสมการ polynomial ที่มีอันดับมากกว่าหนึ่ง เช่น quadratic ในรูป y = b0 + b1x + b11x^2 เราจะใช้ฟังก์ชั่น I() เพื่อป้องกันการแปลความหมายของเครื่องหมายเช่น ^ (ยกกำลัง) หรือ * (คูณ) ในสูตร ว่าเป็น formular operator

> temp<-c(-4,0,1,3,5)
> yield<-c(45.9,79.8,78.9,77.1,62.5)
> model<-lm(yield~temp+I(temp^2))
> summary(model)
Call:
lm(formula = yield ~ temp + I(temp^2))

Residuals:
       1        2        3        4        5 
-0.08442  1.18547 -1.55796  0.53213 -0.07522 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 78.61453    0.94795   82.93 0.000145 ***
temp         3.10624    0.21977   14.13 0.004968 ** 
I(temp^2)   -1.26282    0.07083  -17.83 0.003132 ** 
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Residual standard error: 1.437 on 2 degrees of freedom
Multiple R-Squared: 0.9952,     Adjusted R-squared: 0.9904 
F-statistic: 206.3 on 2 and 2 DF,  p-value: 0.004823 

เราสามารถ Plot เส้นโค้งของสมการนี้ด้วยฟังก์ชั่น lines() ดังนี้

> plot(temp,yield)
> x<-seq(-4,5,by=0.01)
> lines(x,model$coef[1]+model$coef[2]*x+model$coef[3]*x^2)

Stem and leaf diagram

stem เป็นฟังก์ชั่นอยู่ในแพ็คเกจ graphics ต้องทำการ Load เข้ามาเสียก่อนจากเมนู Packages > Load Packages...

> x=c(.03,.05,.08,.08,.10,.11,.18,.19,.20,.20,.22,.22,.23,.29,.30,.32,.34,.40,.47,
.48,.55,.56,.58,.64,.66,.78,.78,.86,.89,.96)
> x
 [1] 0.03 0.05 0.08 0.08 0.10 0.11 0.18 0.19 0.20 0.20 0.22 0.22 0.23 0.29 0.30
[16] 0.32 0.34 0.40 0.47 0.48 0.55 0.56 0.58 0.64 0.66 0.78 0.78 0.86 0.89 0.96

> stem(x)

  The decimal point is 1 digit(s) to the left of the |

  0 | 3588
  1 | 0189
  2 | 002239
  3 | 024
  4 | 078
  5 | 568
  6 | 46
  7 | 88
  8 | 69
  9 | 6

สร้าง random number

สมมติว่าต้องการจัดเรียง 3 treamentsให้กับ 30 หน่วยการทดลอง เราสามารถสร้างตัวเลขสุ่มขึ้นมาได้ง่ายๆ ด้วยฟังก์ชั่น sample ()

> sample(1:30)
 [1] 19 28 27 30 20 14 18 12  5  2 26 25 23  4  7 17  3 16 11 21 10  8 13 15  9
[26] 22  1 24 29  6
> sample(1:30)
 [1] 12 27 24 16  8 22  6 23  5  3  7 17 11 15 10 30 28 26  2  1 21 25 29 13 14
[26] 18  4 20 19  9
จัด treatment ที่ 1 ให้กับ 10 หน่วยทดลองแรก และtreatement ที่ 2, 3 ให้กับอีก 10 หน่วยทดลองต่อไปตามหมายเลขที่ได้จากเลขสุ่ม

Principal Component Analysis

จากตัวอย่างความเข้มของ fluorescence spectrum ที่ 4 ความยาวคลื่นแสงของสาร 12 ตัว

> x
   L300 L350 L400 L450
1    16   62   67   27
2    15   60   69   31
3    14   59   68   31
4    15   61   71   31
5    14   60   70   30
6    14   59   69   30
7    17   63   68   29
8    16   62   69   28
9    15   60   72   30
10   17   63   69   27
11   18   62   68   28
12   18   64   67   29
> summary(x.pc<-princomp(x,cor=T))
Importance of components:
                          Comp.1    Comp.2     Comp.3     Comp.4
Standard deviation     1.6972580 0.8032755 0.62425085 0.29047284
Proportion of Variance 0.7201712 0.1613129 0.09742228 0.02109362
Cumulative Proportion  0.7201712 0.8814841 0.97890638 1.00000000
> Eigenvalue<-x.pc$sdev^2
> Eigenvalue
    Comp.1     Comp.2     Comp.3     Comp.4 
2.88068482 0.64525159 0.38968912 0.08437447 
> x.pc$loadings

Loadings:
     Comp.1 Comp.2 Comp.3 Comp.4
L300  0.547 -0.238 -0.395  0.699
L350  0.546 -0.299 -0.324 -0.712
L400 -0.400 -0.913              
L450 -0.493  0.145 -0.856       

               Comp.1 Comp.2 Comp.3 Comp.4
SS loadings      1.00   1.00   1.00   1.00
Proportion Var   0.25   0.25   0.25   0.25
Cumulative Var   0.25   0.50   0.75   1.00
> x.pc$scores
    
        Comp.1      Comp.2      Comp.3      Comp.4
  1   1.666158  0.80083759  1.03603742 -0.19405319
  2  -1.348447  0.48694126 -0.58608198  0.13446699
  3  -1.799595  1.47573673 -0.15416733  0.06136844
  4  -1.559386 -0.96809150 -0.68963100 -0.25383367
  5  -1.663676 -0.08162762  0.34467898 -0.29254615
  6  -1.730128  0.73999495  0.49862767  0.12592699
  7   1.423589  0.01670531 -0.60066790 -0.18972679
  8   0.764515 -0.36365119  0.53444140 -0.16809432
  9  -1.833764 -1.51562089  0.16751234  0.25937049
  10  1.839914 -0.82136807  0.65452244 -0.09078217
  11  1.811299 -0.06458679 -0.07189018  0.78487657
  12  2.429521  0.29473021 -1.13338186 -0.17697319
> plot(x.pc$scores[,1:2],type="n")
> text(x.pc$scores[,1:2],rownames(x))
> abline(h=0,v=0)

Back to home page