การสร้าง response surface จากแผนการทดลอง Box-Behnken Design ด้วยโปรแกรม R

1. เลือกเมนู Design >> Analyze design >> Response surface plots ...

2. เลือกชนิดของแผนภาพว่าต้องการแบบ Contour plot, 3D plot หรือ Image plot
    กรณีนี้จะเลือกแบบ contour plot โดยมี image plot เป็นฉากหลัง
    ในช่อง Type of specification ให้เลือก Select manually
    จากนั้นเลือกคู่ปัจจัยที่ต้องการ แล้วคลิกลูกศร > ให้ไปอยู่ช่องด้านขวา
  

3. เลือกแถบ Modify slice positions ด้านบน
    กำหนดค่าของปัจจัยที่สาม (ที่ไม่ใช่คู่ปัจจัยที่เลือกก่อนหน้านี้)
   จากนั้นคลิก OK

4. โปรแกรมจะเปิดหน้าต่างใหม่เป็นภาพ contour plot โดยมีฉากหลังเป็น image plot

5. กรณีที่ต้องการสร้าง 3D response surface ให้เลืิอกเมนูเดิม
    คลิกที่ช่อง 3D perspective plot
    หากต้องการให้ surface มีสีให้คลิกที่ช่อง Colour the surface

6. โปรแกรมจะสร้างกราฟ response surface ในหน้าต่างกราฟเดิม (กรณีที่ไม่ได้ปิดหน้าต่างกราฟเดิม) ดังตัวอย่างต่อไปนี้

การวิเคราะห์ผลการทดลอง Box-Behnken Design ด้วยโปรแกรม R

1. ในไฟล์ csv ที่สร้างไว้ สามารถนำผลการทดลองมาเติมลงในคอลัมน์ถัดไป

    ตัวอย่างในภาพด้านล่างเป็นผลผลิตเจลาติน (Yield) จากการสกัดที่อุณหภูมิ (Temp), pH และเวลา (Time) ในระดับต่าง ๆ

    เมื่อนำข้อมูลมาลงแล้ว ให้บันทึก (save) ไว้ในชื่อเดิม

2. นำเข้าข้อมูล โดยเปิดโปรแกรม R Commander แล้วเลือกเมนู Design >> Import... >> Re-Import experiment from csv and rda files...

3. เลือกไดเรกทอรี (โฟลเดอร์) ที่บันทึกข้อมูลไว้
    เลือกไฟล์ rda
    เลือกแผนการทดลองในช่อง Select design
    เช็คถูกในช่อง Replace responses, if they exist already
    เลือกไฟล์ csv
    ตั้งชื่อแผนการทดลองใหม่ หรือใช้ชื่อแผนการทดลองเดิม
    คลิก OK

4. กดปุ่ม View data set เพื่อดูข้อมูลใหม่

5. วิเคราะห์ผล โดยเลือกเมนู Design >>  Analyze design >> Response surface model

6. เลือกค่าตอบสนอง (Response) และปัจจัย (Factors)
    - หากต้องการให้โปรแกรมใช้ค่ารหัส (coded values) ในการคำนวณ ให้เช็คถูกในช่อง Use coded values for calculation ไว้
    - หากต้องการใช้ค่าจริงในการคำนวณ ให้เช็คช่อง Use coded values for calculation ออก
    จากนั้นคลิก OK

7. ตั้งชื่อแบบจำลองในช่อง Enter name for model
    กำหนดเทอมต่าง ๆ ของแบบจำลอง
    จากนั้นคลิก OK

8. จะปรากฏผลการคำนวณในช่อง Output Window ดังนี้


Call:
rsm(formula = Yield ~ FO(Temp, pH, Time) + TWI(Temp, pH, Time) + 
    PQ(Temp, pH, Time), data = Gelatin.withresp)


Residuals:
     Min       1Q   Median       3Q      Max 
-1.01125 -0.56375 -0.01125  0.56375  1.01125 


Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -20.385625   9.966808  -2.045 0.080078 .  
Temp          1.732583   0.262960   6.589 0.000308 ***
pH           -7.547523   1.909353  -3.953 0.005511 ** 
Time          5.634375   1.396219   4.035 0.004961 ** 
Temp:pH       0.164955   0.017931   9.199 3.70e-05 ***
Temp:Time    -0.061750   0.016586  -3.723 0.007427 ** 
pH:Time       0.636486   0.134482   4.733 0.002125 ** 
Temp^2       -0.018633   0.002155  -8.645 5.54e-05 ***
pH^2         -0.518627   0.141705  -3.660 0.008071 ** 
Time^2       -0.771250   0.121246  -6.361 0.000381 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 


# ส่วนที่ระบายสีเหลืองด้านบนนี้เป็นสัมประสิทธิ์หน้าเทอมต่าง ๆ ของสมการ
# ส่วนค่า Pr (>|t|) ด้านหลังแสดงระดับนัยสำคัญของเทอมต่าง ๆ


Residual standard error: 0.9952 on 7 degrees of freedom
Multiple R-squared: 0.9909, Adjusted R-squared: 0.9793 
F-statistic:  84.9 on 9 and 7 DF,  p-value: 2.59e-06 


# ค่า Adjusted R-squared และ p-value ของสมการแสดงในแถบสีฟ้า


Analysis of Variance Table


Response: Yield
                    Df Sum Sq Mean Sq  F value    Pr(>F)
FO(Temp, pH, Time)   3 496.75 165.583 167.1949 7.084e-07
TWI(Temp, pH, Time)  3 119.73  39.908  40.2968 8.739e-05
PQ(Temp, pH, Time)   3 140.26  46.754  47.2089 5.181e-05
Residuals            7   6.93   0.990                   
Lack of fit          3   4.12   1.375   1.9573    0.2624
Pure error           4   2.81   0.702                   


# ระดับนัยสำคัญของ lack of fit แสดงในแถบสีเขียว


Stationary point of response surface:
     Temp        pH      Time 
58.531674  3.796844  2.876294 


Eigenanalysis:
$values
[1] -0.004681207 -0.311172729 -0.992656132


$vectors
           [,1]       [,2]       [,3]
[1,] 0.98316426 -0.1671253  0.0738727
[2,] 0.17936671  0.8055700 -0.5646987
[3,] 0.03486582  0.5684419  0.8219843

การออกแบบแผนการทดลอง Box-Behnken Design ด้วยโปรแกรม R

ตัวอย่าง
ต้องการศึกษาปัจจัยในขั้นตอนการสกัดเจลาตินที่มีผลต่อปริมาณเจลาตินที่สกัดได้ (yield) 3 ปัจจัย ได้แก่
(1) อุณหภูมิในช่วง 40-70 C
(2) pH ในช่วง 4.7 ถึง 7.4
(3) ระยะเวลาการสกัด 1-5 ชม.
การวางแผนการทดลองแบบ Box-Behnken ด้วยโปรแกรม R มีขั้นตอนดังนี้

1. หลังจากเปิด plugin ใน R Commander ที่มีชื่อว่า Rcmdrplugin.DoE แล้ว
    ให้เลือกเมนู Design >> Create design >> Box-Behnken...

2. ในแถบ Base Settings
    กำหนดชื่อของแผนการทดลองในช่อง Name of new design
    กำหนดจำนวนปัจจัยในช่อง Number of factors
    และกำหนดจำนวนซ้ำของจุดกึ่งกลางในช่อง Number of center points

3. ในแถบ Factor Details
    3.1 หากต้องการกำหนดค่าสูงต่ำของแต่ละปัจจัยให้เช็คเครื่องหมายถูกในช่อง Common factor levels ออก แล้วกำหนดชื่อ ค่าต่ำสุด และค่าสูงสุดของแต่ละปัจจัยในช่องด้านล่าง
    3.2 หากต้องการให้ค่าสูงต่ำของแต่ละปัจจัยเป็นค่ารหัส (coded value) ให้ละเครื่องหมายถูกหน้า Common factor levels ไว้
  

4. ในแถบ Export
    ให้กำหนดไฟล์ที่บันทึกเป็น csv and rda
    เลือกไดเรกทอรี (โฟลเดอร์) ที่ต้องการบันทึกไฟล์
    และตั้งชื่อไฟล์ที่ต้องการบันทึก
    จากนั้นคลิก OK

5. คลิกปุ่ม View data set เพื่อดูแผนการทดลองที่สร้างขึ้น

6. แผนการทดลองจะถูกบันทึกไว้ในไฟล์นามสกุล csv ที่สร้างขึ้น ไฟล์ดังกล่าวสามารถเปิดด้วยโปรแกรม spreadsheet เช่น Excel, OpenOffice.org Calc เป็นต้น

การวางแผนการทดลอง Central composite design (CCD) ด้วย RcmdrPlugin.DoE

1. โหลด plug-in ที่ติดตั้งไว้ โดยเลือกเมนู Tools >> Load Rcmdr plugin(s)...
   

   
2. เลือก RcmdrPlugin.DoE แล้วคลิก OK
   

   

   
   
3. จะปรากฏข้อความให้เริ่ม R Commander ใหม่อีกครั้ง คลิก Yes
   

   

   
   
4. R Commander ที่เปิดขึ้นมาใหม่จะมีเมนู Design ปรากฏขึ้น
   ให้เลือกเมนู Design >> Create Design >> Central Composite
   

   

   
   
5. ตั้งชื่อแผนการทดลองในช่อง Name of new design
   ในช่อง Position of star points (alpha) ให้พิมพ์คำว่า “orthogonal” หรือ “rotatable” แล้วแต่แผนการทดลองที่ต้องการ ในที่นี้จะพิมพ์ “rotatable”
   จากนั้นให้คลิกปุ่ม Create new 2-level design
   

   

   
   
6. กำหนดจำนวนปัจจัยในช่อง Number of factors ในที่นี้จะกำหนดเป็น 2 ปัจจัย
   ในช่อง Minimum resolution ให้เลือก III
   จากนั้นคลิกแถบ Factor Details ด้านบน
   

   

     
7. โปรแกรมจะกำหนดระดับของแต่ละปัจจัยเป็น -1 ถึง 1
   หากต้องการกำหนดระดับของปัจจัยเอง ให้เช็คเครื่องหมายถูกหน้า Common factor levels ออก
   

   

   
   
8. เลือกปัจจัย A เพื่อตั้งชื่อปัจจัย (Factor name) กำหนดระดับต่ำสุด (First level) และสูงสุด (Second level)
   จากนั้นเลือกปัจจัย B แล้วกำหนดค่าต่าง ๆ เช่นเดียวกัน
   คลิก OK
   

   

   
   
9. จะปรากฏข้อความให้กำหนดส่วน star ของแผนการทดลอง คลิก OK
   

   

   
   
10. ระบุจำนวน center point ในช่อง Number of center points for star portion
    จำนวนซ้ำของ center point ในแผนการทดลองจะเป็น 2 เท่าของค่าที่ใส่ในช่องนี้
    จากนั้นคลิก OK

11. คลิกปุ่ม View data set ด้านบนเพื่อดูแผนการทดลองที่สร้างขึ้น

การติดตั้งและเรียกใช้งาน R Commander

          R Commander หรือ Rcmdr เป็นแพคเกจหนึ่งของ R ที่แสดงผลเป็นแบบ graphic user interface (GUI) มีลักษณะเป็นเมนู ซึ่งช่วยให้เรียกใช้คำสั่งของ R ได้ง่ายขึ้น

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

> install.packages("Rcmdr", dependencies=TRUE)

          หลังจากติดตั้งแล้ว สามารถเรียกใช้ R commander ด้วยคำสั่ง

> library(Rcmdr)

     จะปรากฏหน้าต่างของ R commander ดังภาพ

 หมายเหตุ:

          กรณีที่ต้องการให้เปิด R commander ทันทีเมื่อเปิดโปรแกรม R ขึ้นมา ให้เปิดไฟล์ rprofile.site ซึ่งอยู่ในไดเรกทอรี C:\Program Files\R\R-…\etc แล้วเพิ่มคำสั่งต่อไปนี้

local({
old <- getOption("defaultPackages")
options(defaultPackages = c(old, "Rcmdr"))
})

การวิเคราะห์ความแปรปรวนสำหรับแผนการทดลองแบบสุ่มสมบูรณ์ (completely randomized design, CRD) ด้วยโปรแกรม R

การวิเคราะห์ความแปรปรวนสำหรับแผนการทดลองแบบสุ่มสมบูรณ์

(completely randomized design, CRD) ด้วยโปรแกรม R

          การวิเคราะห์ความแปรปรวนของผลการทดลองที่ได้จากการวางแผนแบบ CRD สามารถใช้คำสั่งพื้นฐานของ R โดยตรง ไม่ต้องใช้แพคเกจเสริม ขั้นตอนการวิเคราะห์มี 2 ขั้นตอน

          1) วิเคราะห์ ANOVA ด้วยคำสั่ง aov

> ชื่อแบบจำลอง = aov(ค่าตอบสนอง~ปัจจัย, data=ข้อมูล)

          2) แสดงตาราง ANOVA ด้วยคำสั่ง summary

> summary(ชื่อแบบจำลอง)

ตัวอย่าง:

         ปริมาณกรดอะมิโน hydroxyproline ในคอลลาเจนจากปลาที่อาศัยในภูมิอากาศ 4 เขตเป็นดังตาราง

Hydroxyproline (% w/w)
Polar	Temperate	Subtropical	Tropical
7.24	6.79	7.74	9.00
7.24	6.61	10.87	10.28
7.24	6.63	7.80	14.59
6.45	7.81	10.13	14.96
7.79	7.28	10.25	9.38
7.25	7.32	10.01	12.03
6.42	7.28	8.41	10.86
6.03	7.27	9.65	10.86
	7.98	7.6	10.75
	7.85	8.9	10.92
		9.05
		9.41
		11.05
		8.14

          การวิเคราะห์ความแปรปรวนด้วยโปรแกรม R ทำได้ดังนี้

# นำเข้าข้อมูลจาก Excel ที่บันทึกไว้ในชื่อ collagen.csv ไว้ในเขตข้อมูลที่ตั้งชื่อว่า collagen > collagen=read.table('collagen.csv',header=TRUE,sep=',') # สร้างแบบจำลอง ANOVA ที่มีชื่อว่า hyp จากข้อมูล collagen โดยให้ hydroxyproline เป็นค่าตอบสนอง และ climate เป็นปัจจัย > hyp=aov(hydroxyproline~climate, data=collagen) # แสดงตาราง ANOVA ของแบบจำลอง hyp > summary(hyp)             Df  Sum Sq Mean Sq F value    Pr(>F)    climate      3 118.248  39.416  25.919 2.676e-09 *** Residuals   38  57.788   1.521                      --- Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

          จากผลที่ได้แสดงว่าค่าเฉลี่ยปริมาณกรดอะมิโน hydroxyproline ในคอลลาเจนจากปลาทั้ง 4 เขตภูมิอากาศนั้น มีอย่างน้อย 1 เขตที่แตกต่างจากเขตอื่น ที่ระดับนัยสำคัญ 0.001 ดังนั้น จึงต้องทดสอบความแตกต่างของค่าเฉลี่ยด้วยวิธี Duncan โดยใช้แพคเกจ laercio

# เรียกใช้แพคเกจ laercio > library(laercio) # เปรียบเทียบค่าเฉลี่ยจากแบบจำลอง hyp ด้วยวิธี Duncan > LDuncan(hyp)  DUNCAN TEST TO COMPARE MEANS  Confidence Level:  0.95  Dependent Variable:  hydroxyproline  Variation Coefficient:  13.95678 %  Independent Variable:  climate   Factors     Means                 tropical    11.365225089     a    subtropical 9.21465036892857  b   temperate   7.2829372959       c   polar       6.95821078         c