จากข้อมูลตัวอย่างการสลายตัวของสาร Betanin ของ van Boekel (2009) ดังรูป
บน R console
> file <- file.choose()
> data <- read.csv(file,header=TRUE)
> data
time Betanin
1 0 4.29
2 10 3.65
3 20 3.04
4 25 2.83
5 35 2.30
6 40 2.20
7 50 1.83
8 60 1.49
9 90 0.96
10 100 0.70
ลองดูลักษณะความสัมพันธ์ของข้อมูล
> plot(data)
สร้างแบบจำลองด้วยสมการอันดับหนึ่ง C = C0 exp(-kt) โดยใช้ non linear regression ซึ่งใน R คือฟังก์ชั่น nls () โดยเราจะใส่ค่าเริ่มต้นให้กับพารามีเตอร์ดังนี้
ความเข้มข้นเริ่มต้น Betanin0 ประมาณได้จากจุดตัดแกน Y ของกราฟประมาณ 4.5 สำหรับค่า k ประมาณค่าโดยแทนค่า C และ C0 ลงในสมการ เช่น ที่เวลา 25 นาที 2.83 = 4.5 exp(k*25) หา k โดย take ln ได้เป็น k= - ln (2.83/4.5)/25 หรือประมาณ 0.018 และใช้ trace =TRUE เพื่อแสดงค่าพารามีเตอร์ที่โปรแกรมคำนวณได้ในแต่ละรอบของการหาคำตอบ
> model <- nls(Betanin~Betanin0*exp(-k*time), start=list(Betanin0=4.5,k=0.018), data=data,trace=TRUE)
0.08501524 : 4.500 0.018
0.0150839 : 4.31203118 0.01733814
0.01507111 : 4.31285522 0.01732814
0.01507111 : 4.3128601 0.0173282
ใช้ฟังก์ชัน summary() เพื่อดูค่าพารามีเตอร์
> summary(model)
Formula: Betanin ~ Betanin0 * exp(-k * time)
Parameters:
Estimate Std. Error t value Pr(>|t|)
Betanin0 4.3128601 0.0326408 132.13 1.2e-14 ***
k 0.0173282 0.0002667 64.97 3.5e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.0434 on 8 degrees of freedom
Number of iterations to convergence: 3
Achieved convergence tolerance: 4.296e-07
นั่นคือเราได้ค่า k ของสมการเป็น 0.0173 โดยมีค่า standard error เป็น 0.00027 นอกจากนี้อาจจะคำนวณค่าช่วงความเชื่อมั่น (confidence interval) ได้ด้วยฟังก์ชั่น confint()
> confint(model)
Waiting for profiling to be done...
.
.
.
2.5% 97.5%
Betanin0 4.23791259 4.38824654
k 0.01671977 0.01794663
นั่นคือ 95% CI ของ k มีค่า (0.0167 - 0.0179)
จากนั้นเราสามารถวาดเส้นที่เป็น fitted line จากสมการดังกล่าวได้โดยการสร้างตัวแปร x ขึ้นมาใหม่ให้ครอบคลุมช่วงเวลาที่สนใจ เช่น 0 - 100
> x <- 0:100
> x
[1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
[19] 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
[37] 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
[55] 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
[73] 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
[91] 90 91 92 93 94 95 96 97 98 99 100
วาดเส้นด้วย lines() ร่วมกับฟังก์ชั่น predict() ซึ่งใช้คำนวณค่า y จาก model
> lines(x,predict(model,list(time=x)))
van Boekel, M. A. J. S. 2009. Kinetic Modeling of Reactions in Foods. Boca Raton, FL: CRC Press. 767 p.
ไม่มีความคิดเห็น:
แสดงความคิดเห็น