พล็อตกราฟ 2 แกน Y

   กราฟสองแกน Y เป็นที่นิยมในการนำเสนอค่า Y สองค่าที่มีสเกลแตกต่างกันมาก เช่นในตัวอย่าง

pH	Protein	1/Protein
2	10	0.1
4	20	0.05
6	30	0.033333
8	40	0.025
10	50	0.02

    จะเห็นได้ว่าค่าปริมาณ Protein กับค่าส่วนกลับของค่าเหล่านี้มีค่าต่างกันมาก หากเอามาพล็อตบนแกน Y เดียวกันเราแทบจะมองไม่เห็นค่า 1/Protein เลย

    เราสามารถพล็อตกราฟ 2 แกน Y บน R ได้หลายวิธี ซึ่งวิธีหนึ่งที่ง่ายก็คือการใช้แพ็กเกจ plotrix

> install.packages("plotrix")
> library(plotrix)

> pH<-c(2,4,6,8,10)
> protein<-c(10,20,30,40,50)
> reciprocal<-1/protein

จากนั้นเราจะใช้ฟังก์ชั่น twoord.plot() ดังนี้

> twoord.plot(pH,protein,pH,reciprocal,lcol=1,rcol=1,xlab="pH",ylab="Protein (mg/L)",rylab="1/Protein (L/mg)")

 นอกจากนี้เรายังสามารถกำหนด option ต่างๆ ของการพล็อตได้ เช่น lyim และ rylim เป็นตัวกำหนดช่วงแกนของการพล็อต หรือใช้ option ของ graphic ปรกติ ได้เช่น type = "o" เป็นการพล็อตแบบที่เส้นทับบนจุด ไม่ใช่เว้นเหมือนด้านบน

> twoord.plot(pH,protein,pH,reciprocal,type="o",lcol=1,rcol=1,xlab="pH",ylab="Protein (mg/L)",rylab="1/Protein (L/mg)",lylim=c(0,55),rylim=c(0,0.11))

ภาษาไทยใน R

    ต้องยอมรับว่าการใช้ภาษาไทยบน RGui ซึ่งเป็น GUI เบื้องต้นที่มากับ R นั้นมีความยากลำบากมากถึงเป็นไปไม่ได้เลยทีเดียว ทางออกที่ดูจะดีที่สุดน่าจะเป็นการไปใช้ RStudio แทน ซึ่งก็สามารถดาวโหลดได้ฟรีเช่นกัน
   เราจะลองสร้างไฟล์ข้อมูลที่มีภาษาไทยใน Excel ดังตัวอย่าง

 

    ก่อนอื่นให้กำหนดค่าต่อไปนี้ (ผมได้จาก https://th-th.facebook.com/BusinessAnalyticsNIDA/videos/1541641359470302/)

> Sys.setlocale("LC_CTYPE", "thai")
[1] "Thai_Thailand.874"
> options(encoding="UTF-8")

    เราจะใช้การ Import ข้อมูลด้วย File > Import Dataset > From Excel  ซึ่งช่วยให้การนำข้อมูลเข้าจาก Excel โดยตรงเป็นไปด้วยความสะดวกอย่างยิ่ง

   ถ้าลองพล็อตกราฟดู

> attach(thai)
> pie(ปริมาณ,label=ประเทศ)

  ก็ได้กราฟที่พอดูได้นะครับ แต่เรายังพอจะปรับให้ดูดีกว่านี้ได้ เช่น ปรับ font ให้เป็นแบบที่เราต้องการใช้ เช่น THSarabunPSK และขนาดตัวหนังสือด้วยพารามีเตอร์ cex ของการแสดงกราฟฟิค

> windowsFonts(Sarabun=windowsFont("TT TH SarabunPSK"))
> pie(ปริมาณ,label=ประเทศ,cex=2.0,family="Sarabun")

 

mean comparison ด้วยแพกเกจ agricolae

      ก่อนหน้านี้ผมได้โพสเกี่ยวกับการวิเคราะห์ ANOVA และการเปรียบเทียบทรีตเมนต์เมื่อพบอิทธิพลของทรีนเมนต์ที่ศึกษามีความสำคัญอย่างมีนัยสำคัญ โดยใช้แพกเกจ multcomp แต่ดูเหมือนจะยังไม่ตอบโจทย์ความต้องการในการวิเคราะห์เมื่อเทียบกับการใช้ซอฟท์แวร์สถิติอื่นๆ โดยเฉพาะการให้สัญลักษณ์ว่าทรีตเมนต์ใดต่างทรีตเมนต์ใดไม่ต่างกัน ในครั้งนี้ผมจะใช้แพกเกจ agricolae ซึ่งช่วยให้การวิเคราะห์ ANOVA และการเปรียบเทียบทรีตเมนต์ทำได้อย่างสะดวกมากทีเดียว

     เริ่มจากที่ข้อมูลที่ใช้เป็นตัวอย่างในไฟล์ Excel ดังรูป

     
     บน RStudio ใช้ File > Import Dataset > From Excel  เข้ามาในชื่อ packing

     ทำการติดตั้งแพกเกจ agricolae และเรียกใช้

> install.packages("agricolae")

> library(agricolae)

  สั่ง attach ชุดข้อมูลที่ต้องการใช้

> attach(packing)  

  
    วิเคราะห์ ANOVA ด้วย aov() ตามปรกติ

> aov.out<-aov(count~treatment)

> summary(aov.out)
            Df Sum Sq Mean Sq F value   Pr(>F)    
treatment    3  32.87  10.958   94.58 1.38e-06 ***
Residuals    8   0.93   0.116                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

  พบว่าทรีตเมนต์มีอิทธิพลอย่างมีนัยสำคัญ จึงทำการเปรียบเทียบค่าเฉลี่ยของแต่ละทรีตเมนต์ เช่น การใช้ LSD จะใช้ฟังก์ชั่น LSD.test()

> LSD.test(aov.out,"treatment",console=TRUE)

Study: aov.out ~ "treatment"

LSD t Test for count 

Mean Square Error:  0.11585 

treatment,  means and individual ( 95 %) CI

           count       std r      LCL      UCL  Min  Max
CO2         3.36 0.3968627 3 2.906844 3.813156 2.91 3.66
commercial  7.48 0.4386342 3 7.026844 7.933156 6.98 7.80
mixed gas   7.26 0.1946792 3 6.806844 7.713156 7.04 7.41
vacuum      5.50 0.2749545 3 5.046844 5.953156 5.26 5.80

alpha: 0.05 ; Df Error: 8
Critical Value of t: 2.306004 

t-Student: 2.306004
Alpha    : 0.05
Least Significant Difference 0.640859
Means with the same letter are not significantly different


Groups, Treatments and means
a   commercial   7.48 
a   mixed gas   7.26 
b   vacuum   5.5 
c   CO2   3.36 

จะเห็นได้ว่าฟังก์ชั่น LSD.test แสดงผลทั้งค่าเฉลี่ยของแต่ละทรีตเมนต์และผลการเปรียบเทียบด้วย LSD ในรูปของตารางในด้านล่างซึ่งช่วยให้ผู้วิเคราะห์สามารถอ่านผลได้ง่ายเช่นเดียวกับซอฟท์แวร์อื่นๆ

เช่นเดียวกัน เราสามารถใช้เทคนิค Duncan ในการเปรียบเทียบทรีตเมนต์ได้ด้วยฟังก์ชั่น duncan.test() ซึ่งในกรณีนี้ให้ผลไม่แตกต่างกัน

> duncan.test(aov.out,"treatment",console=TRUE)

Study: aov.out ~ "treatment"

Duncan's new multiple range test
for count 

Mean Square Error:  0.11585 

treatment,  means

           count       std r  Min  Max
CO2         3.36 0.3968627 3 2.91 3.66
commercial  7.48 0.4386342 3 6.98 7.80
mixed gas   7.26 0.1946792 3 7.04 7.41
vacuum      5.50 0.2749545 3 5.26 5.80

alpha: 0.05 ; Df Error: 8 

Critical Range
        2         3         4 
0.6408590 0.6678356 0.6829140 

Means with the same letter are not significantly different.

Groups, Treatments and means
a   commercial   7.48 
a   mixed gas    7.26 
b   vacuum       5.5 
c   CO2          3.36  

ในแพกเกจ agricolae ยังมีฟังก์ชั่นในการเปรียบเทียบทรีตเมนต์ที่ใช้กันอื่นๆ ไม่ว่าจะเป็น Tukey, SNK, bonferroni,Sche ffe

อบรมการใช้ R

        ภาควิชาเทคโนโลยีอาหาร คณะวิศวกรรมศาสตร์และเทคโนโลยีอุตสาหกรรม มหาวิทยาลัยศิลปากรขอเชิญเจ้าหน้าที่ที่เกี่ยวข้อง และผู้สนใจทั่วไปเข้าอบรมเชิงปฏิบัติการเรื่อง "การใช้ซอฟท์แวร์ฟรี R ในการวิเคราะห์สถิติเบื้องต้นสำหรับอุตสาหกรรมอาหาร"


วันที่:  31 มีนาคม 2560 เวลา 8.30 - 16.30 น.
สถานที่: ห้องคอมพิวเตอร์ ภาควิชาเทคโนโลยีอาหาร ชั้น 3 อาคารคณะวิศวกรรมศาสตร์และเทคโนโลยีอุตสาหกรรม มหาวิทยาลัยศิลปากร อ.เมือง จ.นครปฐม

เปิดรับสมัครถึง 28 มีนาคม 2560 ค่าลงทะเบียนการอบรมท่านละ 1,500 บาท

รายละเอียดกำหนดการและสมัครออนไลน์: http://bit.ly/2lwHMoP

 

แพกเกจ FADA สำหรับเลือกตัวแปรและทำ LDA

       ในกรณีข้อมูลที่มีตัวแปรจำนวนมากอย่างในกรณีของข้อมูล NIR spectrum นั้น มักมีความสนใจในการเลือกเฉพาะตัวแปรที่มีความสำคัญในการสร้างโมเดลในการทำนาย ซึ่งวิธีการเลือกตัวแปรนั้นมีหลายวิธี แพกเกจ FADA เป็นอีกแพกเกจที่ออกแบบมาเพื่อใช้สำหรับการเลือกตัวแปรและวิเคราะห์ discriminant analysis โดยเฉพาะ โดยขั้นตอนของการวิเคราะห์มีอยู่ 3 ขั้นตอนคือ

1. decorrelation ข้อมูลใน training set ด้วยฟังก์ชั่น decorrelate.train
2. decorrelation ข้อมูล test set ด้วยฟังก์ชั่น decorrelate.test
3. สร้าง supervised classification model ด้วยฟังก์ชั่น FADA

บน  R console

> install.packages("FADA")
> library(FADA)

สร้างข้อมูล training set โดยจะต้องมีรูปแบบของข้อมูลเป็น list ที่มี 2 ตัวแปรคือ x เป็น matrix ที่ประกอบด้วยด้วยข้อมูลตัวแปร (ในกรณีนี้คือค่า absorbance ที่แต่ละความยาวคลื่น ของแต่ละตัวอย่าง) และ y เป็นตัวแปรที่ใช้ระบุกลุ่มของตัวอย่าง (ต้องเป็นตัวเลขเท่านั้น)

> source("C:\\Users\\Windows\\Documents\\csv2matrix.R")
> coffee<-csv2matrix()
> groups<-rep(c(1,2,3),each=50)
> even<-seq(2,300,2)
> odd<-seq(1,300,2)

ชุด traning set

> coffee.list.train <-list(x=coffee[odd,],y=groups)
> str(coffee.list.train)
List of 2
 $ x: num [1:150, 1:558] 0.609 0.638 0.619 0.614 0.631 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:558] "12481.8432" "12466.41447" "12450.98574" "12435.55701" ...
 $ y: num [1:150] 1 1 1 1 1 1 1 1 1 1 ...


ชุด test set

> coffee.list.test <-list(x=coffee[even,],y=groups)
> str(coffee.list.test)
List of 2
 $ x: num [1:150, 1:558] 0.618 0.625 0.64 0.65 0.623 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:558] "12481.8432" "12466.41447" "12450.98574" "12435.55701" ...
 $ y: num [1:150] 1 1 1 1 1 1 1 1 1 1 ...

 decorrelation ข้อมูลใน training set ด้วยฟังก์ชั่น decorrelate.train

> fa <- decorrelate.train(coffee.list.train)
[1] "Number of factors: 10 factors"
[1] "Objective criterion: "
[1] 270.0801
[1] 140.3937
[1] 31.69263
[1] 1.058545
[1] 0.02515622
[1] 0.001721782
[1] 1.541182e-05

decorrelation ข้อมูล test set ด้วยฟังก์ชั่น decorrelate.test

> fa2 <- decorrelate.test(fa,coffee.list.test)


สร้าง supervised classification model ด้วยฟังก์ชั่น FADA

> class <- FADA(fa2,method="sda",sda.method="lfdr")

ตัวแปรที่ใช้ทั้งหมด 22 ตัว (ลำดับของตัวแปร)

> class$selected
 [1] 528 348 347 536 527 349 543 470 542 471 469 532 529 472 468 537  54 544 541 192 351 350

ซึ่งตรงกับค่าเลขคลื่นดังต่อไปนี้

 [1]  4350.902  7128.074  7143.502  4227.472  4366.331  7112.645  4119.471
 [8]  5245.768  4134.900  5230.340  5261.197  4289.187  4335.473  5214.911
[15]  5276.626  4212.043 11664.120  4104.042  4150.329  9534.956  7081.787
[22]  7097.216

ผลการจัดกลุ่มในตัวอย่าง test set จะเห็นว่าถูกต้อง 100%

> class$proba.test
       1 2 3
  [1,] 1 0 0
  [2,] 1 0 0
  [3,] 1 0 0
  [4,] 1 0 0
  [5,] 1 0 0
  [6,] 1 0 0
  [7,] 1 0 0
  [8,] 1 0 0
  [9,] 1 0 0
 [10,] 1 0 0
 [11,] 1 0 0
 [12,] 1 0 0
 [13,] 1 0 0
 [14,] 1 0 0
 [15,] 1 0 0
 [16,] 1 0 0
 [17,] 1 0 0
 [18,] 1 0 0
 [19,] 1 0 0
 [20,] 1 0 0
 [21,] 1 0 0
 [22,] 1 0 0
 [23,] 1 0 0
 [24,] 1 0 0
 [25,] 1 0 0
 [26,] 1 0 0
 [27,] 1 0 0
 [28,] 1 0 0
 [29,] 1 0 0
 [30,] 1 0 0
 [31,] 1 0 0
 [32,] 1 0 0
 [33,] 1 0 0
 [34,] 1 0 0
 [35,] 1 0 0
 [36,] 1 0 0
 [37,] 1 0 0
 [38,] 1 0 0
 [39,] 1 0 0
 [40,] 1 0 0
 [41,] 1 0 0
 [42,] 1 0 0
 [43,] 1 0 0
 [44,] 1 0 0
 [45,] 1 0 0
 [46,] 1 0 0
 [47,] 1 0 0
 [48,] 1 0 0
 [49,] 1 0 0
 [50,] 1 0 0
 [51,] 0 1 0
 [52,] 0 1 0
 [53,] 0 1 0
 [54,] 0 1 0
 [55,] 0 1 0
 [56,] 0 1 0
 [57,] 0 1 0
 [58,] 0 1 0
 [59,] 0 1 0
 [60,] 0 1 0
 [61,] 0 1 0
 [62,] 0 1 0
 [63,] 0 1 0
 [64,] 0 1 0
 [65,] 0 1 0
 [66,] 0 1 0
 [67,] 0 1 0
 [68,] 0 1 0
 [69,] 0 1 0
 [70,] 0 1 0
 [71,] 0 1 0
 [72,] 0 1 0
 [73,] 0 1 0
 [74,] 0 1 0
 [75,] 0 1 0
 [76,] 0 1 0
 [77,] 0 1 0
 [78,] 0 1 0
 [79,] 0 1 0
 [80,] 0 1 0
 [81,] 0 1 0
 [82,] 0 1 0
 [83,] 0 1 0
 [84,] 0 1 0
 [85,] 0 1 0
 [86,] 0 1 0
 [87,] 0 1 0
 [88,] 0 1 0
 [89,] 0 1 0
 [90,] 0 1 0
 [91,] 0 1 0
 [92,] 0 1 0
 [93,] 0 1 0
 [94,] 0 1 0
 [95,] 0 1 0
 [96,] 0 1 0
 [97,] 0 1 0
 [98,] 0 1 0
 [99,] 0 1 0
[100,] 0 1 0
[101,] 0 0 1
[102,] 0 0 1
[103,] 0 0 1
[104,] 0 0 1
[105,] 0 0 1
[106,] 0 0 1
[107,] 0 0 1
[108,] 0 0 1
[109,] 0 0 1
[110,] 0 0 1
[111,] 0 0 1
[112,] 0 0 1
[113,] 0 0 1
[114,] 0 0 1
[115,] 0 0 1
[116,] 0 0 1
[117,] 0 0 1
[118,] 0 0 1
[119,] 0 0 1
[120,] 0 0 1
[121,] 0 0 1
[122,] 0 0 1
[123,] 0 0 1
[124,] 0 0 1
[125,] 0 0 1
[126,] 0 0 1
[127,] 0 0 1
[128,] 0 0 1
[129,] 0 0 1
[130,] 0 0 1
[131,] 0 0 1
[132,] 0 0 1
[133,] 0 0 1
[134,] 0 0 1
[135,] 0 0 1
[136,] 0 0 1
[137,] 0 0 1
[138,] 0 0 1
[139,] 0 0 1
[140,] 0 0 1
[141,] 0 0 1
[142,] 0 0 1
[143,] 0 0 1
[144,] 0 0 1
[145,] 0 0 1
[146,] 0 0 1
[147,] 0 0 1
[148,] 0 0 1
[149,] 0 0 1
[150,] 0 0 1

ถ้าเราใช้ตัวแปรเหล่านี้ไปใช้วิเคราะห์ lda ใน MASS เหมือนก่อนหน้านี้

> grouping<-rep(c("black","unpeel","good"),each=100)
> coffee.df<-as.data.frame(coffee)
> coffee.trn <-coffee.df[odd,class$selected]
> str(coffee.trn)
'data.frame':   150 obs. of  22 variables:
 .
 .
 .




> coffee.test <-coffee.df[even,class$selected]> coffee.lda <- lda(coffee.trn,grouping =grouping[odd],prior = rep(1,3)/3)
> plot(coffee.lda)

> coffee.predict <-predict(coffee.lda,new=coffee.test) 
> table(grouping[even],coffee.predict$class)
       
         black good unpeel
  black     50    0      0
  good       0   50      0
  unpeel     0    0     50

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


ถ้าสร้างตัวแปร wavenum ซึ่งมีข้อมูล wavenumber ที่ใช้

> wave<-read.csv(file.choose(),header=FALSE)
> wavenum<-wave$V1

> plot(wavenum,coffee[1,],type="l")

> plot(wavenum[class$selected],coffee[1,class$selected])

วิเคราะห์ LDA (Linear Discriminant Analysis) ด้วย R

    LDA เป็นเทคนิคในการจัดกลุ่ม (classification) โดยการสร้างโมเดลจาก training set เพื่อใช้ทดสอบการจัดกลุ่มในตัวอย่างในอนาคต ในตัวอย่างก่อนหน้านี้ซึ่งเราต้องการจัดกลุ่มตัวอย่างออกเป็น 3 กลุ่ม เราจะแบ่งข้อมูลตัวอย่างจำนวน 300 ตัวอย่างเป็น Training 150 ตัวอย่าง และอีก 150 ตัวอย่างไว้ทดสอบ

> coffee<-csv2matrix()
> coffee.df<-as.data.frame(coffee)
> groups<-rep(c("good","unpeel","black"),each=100)

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

> even<-seq(2,300,2)
> odd<-seq(1,300,2) 

แบ่งข้อมูลเป็น 2 ชุด
> coffee.trn <-coffee.df[odd,]
> coffee.test <-coffee.df[even,]

ใช้ฟังก์ชั่น lda ในแพกเกจ MASS

> install.packages("MASS")
> library(MASS)

> coffee.lda <- lda(coffee.trn,grouping =groups[odd],prior = rep(1,3)/3)
Warning message:
In lda.default(x, grouping, ...) : variables are collinear

เกิดความข้อความเตือนเนื่องจากตัวแปร x ใน NIR มีลักษณะ collinear  คือมีความสัมพันธ์กันเอง เราจะแก้ปัญหานี้ในภายหลัง


> plot(coffee.lda)

จากการพล็อตข้อมูลของ Training set จะเห็นว่าสามารถแยกกลุ่มออกได้ชัดเจนมาก

ทำนายการจัดกลุ่มตัวอย่างอีก 150 ตัวอย่างด้วยฟังก์ชั่น predict

> coffee.predict <-predict(coffee.lda,new=coffee.test) 

แสดงผลการทำนายกลุ่ม

> table(groups[even],coffee.predict$class)
       
         black good unpeel
  black     50    0      0
  good       0   50      0
  unpeel     0    0     50

จะพบว่าโมเดลที่สร้างขึ้นสามารถจัดกลุ่มตัวอย่างใน test set ได้ถูกต้อง 100 %


เนื่องจากตัวแปรของ NIR มีจำนวนมากและมักจะมีลักษณะสัมพันธ์กัน เราสามารถใช้ตัวแปรใหม่จาก PCA นั่นคือ PC มาใช้ในการทำ LDA ได้



> coffee.pca.trn <-prcomp(coffee.trn)
> summary(coffee.pca.trn)


predict ค่า PC ใน test set
> coffee.pca.test<-predict(coffee.pca.trn,new=coffee.test)


จาก summary จะพบว่า PC1-10 อธิบายความแปรปรวนได้เกิน 99% ดังนั้นเราจะใช้เพียง PC1-10 เป็นตัวแปรสำหรบ LDA

> coffee.pca.lda <- lda(coffee.pca.trn$x[,1:10],grouping =groups[odd],prior = rep(1,3)/3) 
> plot(coffee.pca.lda)

> coffee.pca.predict <-predict(coffee.pca.lda,new=coffee.pca.test[,1:10])
> table(groups[even],coffee.pca.predict$class)
       
         black good unpeel
  black     50    0      0
  good       0   50      0
  unpeel     0    0     50

วิเคราะห์ PCA ด้วย R

    เราสามารถวิเคราะห์ PCA (Principal Component Analysis) บน R ได้หลายวิธี แต่วิธีที่นิยมคือการใช้ ฟังก์ชั่น prcomp

   ในตัวอย่างนี้ผมจะลองวิเคราะห์ PCA เพื่อแบ่งกลุ่ม NIR spectra ของตัวอย่าง 3 กลุ่ม กลุ่มละ 100 ตัวอย่าง ในแต่เส้น NIR spectrum นั้นได้มากจากการวัดการดูดกลืนแสงในช่วง12500 - 4000 cm-1 จึงมีข้อมูลจำนวน 558 ค่า ข้อมูลซึ่งเดิมมาจากโปรแกรม OPUS ถูกแปลงเป็น csv ตามโพสต์ก่อนหน้านี้ จากนั้นแปลง csv เป็น matrix อีกที
    จาก matrix เราจะแปลงข้อมูลเป็น data.frame เพื่อให้ง่ายต่อการวิเคราะห์


> source("D:\\csv2matrix.R")
> data<-csv2matrix()
> data.df<-as.data.frame(data)

สร้างข้อมูลจัดกลุ่มตัวอย่าง 3 กลุ่ม คือ good, black, unpeel

> groups<-rep(c("good","black","unpeel"),each=100)

วิเคราะห์ PCA ด้วยฟังก์ชั่น prcomp()

> data.pca<-prcomp(data.df)

แสดง PC ต่างๆ 
> summary(data.pca)


จากนั้นเราจะแสดง biplot ด้วย ggbiplot ซึ่งจะต้องติดตั้งจาก github

> library(devtools)

> install_github("vqv/ggbiplot")

> library(ggbiplot)
> ggbiplot(data.pca,obs.scale=0,groups=groups,ellipse=TRUE,var.axes=FALSE) 

confidence interval สำหรับ nls model

    มีนักศึกษามาคุยกับผมเรื่องการสร้าง confidence interval สำหรับ nonlinear model ลองๆ หาดูเลยพบว่ามีวิธีการที่ค่อนข้างสะดวกอยู่คือการใช้ฟังก์ชั้่น plotFit ในแพกเกจชื่อ investr
    จากข้อมูลตามโพสต์ก่อนหน้านี้


บน R console

> install.packages("investr")
> library(investr)
> plotFit(model,interval="confidence",pch=19,shade=FALSE)

mean comparison สำหรับ RCBD

     โพสต์ก่อนนี้ได้กล่าวถึงการทำ mean comparison สำหรับการทดลองแบบ CRD คราวนี้ลองมาดูกรณีของการทดลองแบบ RCBD บ้าง จากตัวอย่างก่อนหน้านี้ เราพบว่า treatment มีความแตกต่างกันอย่างนัยสำคัญเราจะใช้แพกเกจ  multcomp ในการเปรียบเทียบค่าเฉลี่ยของ treatment ดูว่า Treatment ใดแตกต่างกันบ้าง

   บน R console
เพื่อให้การใช้การเปรียบเทียบ treatment กับ control ทำได้สะดวกเราจะเปลี่ยนวิธีการกรอกข้อมูลเล็กน้อย โดยใช้ 0 แทน control

> data
   trt block nitrogen
1    0     1    34.98
2    0     2    41.22
3    0     3    36.94
4    0     4    39.97
5    2     1    40.89
6    2     2    46.69
7    2     3    46.65
8    2     4    41.90
9    3     1    42.07
10   3     2    49.42
11   3     3    52.68
12   3     4    42.91
13   4     1    37.18
14   4     2    45.85
15   4     3    40.23
16   4     4    39.20
17   5     1    37.99
18   5     2    41.99
19   5     3    37.61
20   5     4    40.45
21   6     1    34.89
22   6     2    50.15
23   6     3    44.57
24   6     4    43.29

 ทำหให้ทั้ง trt และ block เป็น factor
> data$block<-as.factor(data$block)
> data$trt<-as.factor(data$trt)


 > str(data)
'data.frame':   24 obs. of  3 variables:
 $ trt     : Factor w/ 6 levels "0","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
 $ block   : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
 $ nitrogen: num  35 41.2 36.9 40 40.9 ...

> model <-aov(nitrogen~trt+block,data=data)
> summary(model)
            Df Sum Sq Mean Sq F value  Pr(>F) 
trt          5  201.3   40.26   5.592 0.00419 **
block        3  197.0   65.67   9.120 0.00112 **
Residuals   15  108.0    7.20                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> library(multcomp)
Loading required package: mvtnorm
Loading required package: survival
Loading required package: TH.data

คำส่งจะคล้ายการใช้ Tukey เพียงแต่ใช้ "Dunnett"

> compare<-glht(model,linfct=mcp(trt="Dunnett"))
> summary(compare)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = nitrogen ~ trt + block, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)  
2 - 0 == 0    5.755      1.897   3.033  0.03331 *
3 - 0 == 0    8.492      1.897   4.476  0.00181 **
4 - 0 == 0    2.337      1.897   1.232  0.62267  
5 - 0 == 0    1.232      1.897   0.650  0.94550  
6 - 0 == 0    4.947      1.897   2.607  0.07415 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)


ในกรณีที่ 0 ไม่ใช่ตัว control ที่เราอยากเปรียบเทียบด้วย เช่นเป็น treatment 4 ตามตามตัวอย่างใน Kuehl (2000) เราสามารถกำหนดให้เปรียบเทียบ treatment ต่างๆ กับ 4 ได้โดยการสร้าง contrast matrix ดังนี้

> contrast <-rbind("0-4"=c(1,0,0,-1,0,0),
                   "2-4"=c(0,1,0,-1,0,0),
                   "3-4"=c(0,0,1,-1,0,0),
                   "5-4"=c(0,0,0,-1,1,0),
                   "6-4"=c(0,0,0,-1,0,1))
> contrast
    [,1] [,2] [,3] [,4] [,5] [,6]
0-4    1    0    0   -1    0    0
2-4    0    1    0   -1    0    0
3-4    0    0    1   -1    0    0
5-4    0    0    0   -1    1    0
6-4    0    0    0   -1    0    1

> summary(glht(model,linfct=mcp(trt=contrast)))

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: aov(formula = nitrogen ~ trt + block, data = data)

Linear Hypotheses:
         Estimate Std. Error t value Pr(>|t|) 
0-4 == 0   -2.337      1.897  -1.232   0.6226 
2-4 == 0    3.417      1.897   1.801   0.2941 
3-4 == 0    6.155      1.897   3.244   0.0218 *
5-4 == 0   -1.105      1.897  -0.582   0.9643 
6-4 == 0    2.610      1.897   1.376   0.5284 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

ทำให้สรุปได้ว่ามีเพียง treatment 3 เท่านั้นที่ให้ค่าแตกต่างกับ treatment 4

mean comparison ด้วย R

    หลังจากการวิเคราะห์ ANOVA  เมื่อพบความแตกต่างอย่างมีนัยสำคัญของ treatment เราก็อยากทราบว่า Treatment ใดบ้างที่แตกต่างกัน โดยการวิเคราะห์ที่เรียกว่า post-hoc mean comparison ก่อนหน้านี้ผมได้เขียนเกี่ยวกับการใช้วิธี Tukey สำหรับการเปรียบเทียบค่าเฉลี่ยไปแล้ว ด้วยฟังก์ชั่น TukeyHSD

    ใน R มีแพกเกจหนึ่งชื่อ multcomp ซึ่งใช้ในการเปรียบเทียบค่าเฉลี่ยได้หลายรูปแบบ เช่นจากตัวอย่างเดิม

บน R console

> install.packages("multcomp")
> library(multcomp)

ส่วนนี้เป็นการวิเคราะห์ ANOVA

> file<-file.choose()
> data<-read.csv(file,header=TRUE)
> data
    treatment count
1  commercial  7.66
2  commercial  6.98
3  commercial  7.80
4      vacuum  5.26
5      vacuum  5.44
6      vacuum  5.80
7   mixed gas  7.41
8   mixed gas  7.33
9   mixed gas  7.04
10        CO2  3.51
11        CO2  2.91
12        CO2  3.66
> model <- aov(count~treatment,data=data)
> summary(model)
            Df Sum Sq Mean Sq F value   Pr(>F)   
treatment    3  32.87  10.958   94.58 1.38e-06 ***
Residuals    8   0.93   0.116                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

การเปรียบเทียบแบบ pairwise comparison ด้วยฟังก์ชั่น glht
#ใน  treatment="Tukey" นั้น treatment คือชื่อของส่วนที่เราจะเปรียบเทียบซึ่งในที่นี่ชื่อ treatment ตามข้อมูลด้านบน

> compare<-glht(model,linfct=mcp(treatment="Tukey"))
> summary(compare)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = count ~ treatment, data = data)

Linear Hypotheses:
                            Estimate Std. Error t value Pr(>|t|)   
commercial - CO2 == 0         4.1200     0.2779  14.825   <0.001 ***
mixed gas - CO2 == 0          3.9000     0.2779  14.033   <0.001 ***
vacuum - CO2 == 0             2.1400     0.2779   7.700   <0.001 ***
mixed gas - commercial == 0  -0.2200     0.2779  -0.792    0.856   
vacuum - commercial == 0     -1.9800     0.2779  -7.125   <0.001 ***
vacuum - mixed gas == 0      -1.7600     0.2779  -6.333   <0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

จะเห็นว่าได้ผลลักษณะเดียวกันกับที่ใช้ TukeyHSD 


        เราสามารถทดสอบโดยใช้ Fisher's LSD โดยใช้ summary และกำหนด test=univariate() จะเห็นได้ว่าค่า p ที่ได้จะมีค่าน้อยลงเนื่องจากค่า critical difference ใน LSD จะมีค่าน้อย

> summary(compare,test=univariate())

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = count ~ treatment, data = data)

Linear Hypotheses:
                            Estimate Std. Error t value Pr(>|t|)   
commercial - CO2 == 0         4.1200     0.2779  14.825 4.22e-07 ***
mixed gas - CO2 == 0          3.9000     0.2779  14.033 6.45e-07 ***
vacuum - CO2 == 0             2.1400     0.2779   7.700 5.74e-05 ***
mixed gas - commercial == 0  -0.2200     0.2779  -0.792 0.451410   
vacuum - commercial == 0     -1.9800     0.2779  -7.125 9.95e-05 ***
vacuum - mixed gas == 0      -1.7600     0.2779  -6.333 0.000225 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Univariate p values reported)








หากต้องการใช้ Bonferroni test

> summary(compare,test=adjusted(type="bonferroni"))

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = count ~ treatment, data = data)

Linear Hypotheses:
                            Estimate Std. Error t value Pr(>|t|)   
commercial - CO2 == 0         4.1200     0.2779  14.825 2.53e-06 ***
mixed gas - CO2 == 0          3.9000     0.2779  14.033 3.87e-06 ***
vacuum - CO2 == 0             2.1400     0.2779   7.700 0.000345 ***
mixed gas - commercial == 0  -0.2200     0.2779  -0.792 1.000000   
vacuum - commercial == 0     -1.9800     0.2779  -7.125 0.000597 ***
vacuum - mixed gas == 0      -1.7600     0.2779  -6.333 0.001348 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- bonferroni method)

วิเคราะห์ ANOVA สำหรับ RCBD ด้วย R

     ในโพสต์ก่อนหน้านี้ผมได้เขียนถึงการวิเคราะห์ ANOVA สำหรับแผนการทดลองแบบ CRD ไปแล้ว คราวนี้จะได้พูดถึงแผนการทดลองอีกแบบที่พบบ่อยคือ RCBD (Randomized Complete Block Design) โดยผมใช้ตัวอย่างจาก Kuehl (2001) เช่นกัน โดยการทดลองนี้มี 6 treatment และจัดเป็น 4  block

    บน R console
> file<-file.choose()   #เลือกไฟล์ csv ที่บันทึกข้อมูลไว้
> data <-read.csv(file,header=TRUE)    #อ่านไฟล์ที่เลือก
> data
       trt block nitrogen
1  control     1    34.98
2  control     2    41.22
3  control     3    36.94
4  control     4    39.97
5        2     1    40.89
6        2     2    46.69
7        2     3    46.65
8        2     4    41.90
9        3     1    42.07
10       3     2    49.42
11       3     3    52.68
12       3     4    42.91
13       4     1    37.18
14       4     2    45.85
15       4     3    40.23
16       4     4    39.20
17       5     1    37.99
18       5     2    41.99
19       5     3    37.61
20       5     4    40.45
21       6     1    34.89
22       6     2    50.15
23       6     3    44.57
24       6     4    43.29

ดูโครงสร้างข้อมูล

> str(data)
'data.frame':   24 obs. of  3 variables:
 $ trt     : Factor w/ 6 levels "2","3","4","5",..: 6 6 6 6 1 1 1 1 2 2 ...
 $ block   : int  1 2 3 4 1 2 3 4 1 2 ...
 $ nitrogen: num  35 41.2 36.9 40 40.9 ...


จะเห็นว่าข้อมูลที่อ่านเข้ามาอยู่ในรูปของ data.frame โดย trt (หมายถึง treatment ในการทดลองนี้) เป็นข้อมูลแบบจัดกลุ่ม (factor) ที่มี 6 ระดับอยู่แล้ว แต่ว่าข้อมูลของ block ยังไม่เป็น เราจะแปลงข้อมูลนี้เป็น factor

> data$block<-as.factor(data$block)
> str(data)
'data.frame':   24 obs. of  3 variables:
 $ trt     : Factor w/ 6 levels "2","3","4","5",..: 6 6 6 6 1 1 1 1 2 2 ...
 $ block   : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
 $ nitrogen: num  35 41.2 36.9 40 40.9 ...

จะเห็นว่าข้อมูล block กลายเป็น factor ที่มี 4 ระดับ

วิเคราะห์ ANOVA โดยใช้ฟังก์ชั่น aov คล้ายกับในกรณี CRD แต่เราเพิ่มเทอมของ block เข้าในโมเดลดังนี้

> model <-aov(nitrogen~trt+block,data=data)

แสดงผล
> summary(model)
            Df Sum Sq Mean Sq F value  Pr(>F)  
trt          5  201.3   40.26   5.592 0.00419 **
block        3  197.0   65.67   9.120 0.00112 **
Residuals   15  108.0    7.20                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
    จากผลการวิเคราะห์แสดงให้เห็นว่า treatment มีความแตกต่างกันอย่างมีนัยสำคัญ (p = 0.00419)


เอกสารอ้างอิง
Kuehl RO. 2000. Design of experiment: statistical principles of research design and analysis. Pacific Grove: Duxbury Press. 666 p.