Fitting Langmuir and Freundlich curve in R
Duc Nguyen | Chuyên đào tạo kỹ năng R |
www.tuhocr.com
June 25, 2023
Bài này mình hướng dẫn cách fitting đường đẳng nhiệt Langmuir và Freundlich trong R phục vụ nghiên cứu về hấp phụ vật liệu nhé. Cơ bản là trước đây mình có viết một bài báo về khảo sát khả năng hấp phụ ammonium của than sinh học từ trấu, lúc đó mình sử dụng Origin Pro để fitting đường cong hấp phụ để thu nhận các tham số của phương trình Langmuir. Nay mình modify function từ package PUPAIM chuyên phân tích hóa lý trong R để tìm tham số cho phương trình Langmuir và Freundlich.
Kết quả này mình so sánh với Origin Pro thì thấy cũng không khác biệt nhiều (đặc biệt ở tham số Qsat, K hay n đều chính xác đến 4 con số sau dấu thập phân). Các bạn có thể tham khảo function này và cải tiến nếu cần thiết nhé.
Data từ bài báo Adsorption of ammonium, nitrite, and nitrate onto rice husk biochar for nitrogen removal
https://www.researchgate.net/publication/350495918_Adsorption_of_ammonium_nitrite_and_nitrate_onto_rice_husk_biochar_for_nitrogen_removal
ammonium_rhb <- read.delim("ammonium-rhb-16h.txt")
print(ammonium_rhb)## Ce Qe
## 1 0.34 0.0036
## 2 0.36 0.0032
## 3 0.38 0.0028
## 4 0.77 0.0092
## 5 0.83 0.0080
## 6 0.78 0.0090
## 7 1.48 0.0108
## 8 1.53 0.0098
## 9 1.52 0.0100
## 10 3.03 0.0196
## 11 3.02 0.0198
## 12 2.94 0.0214
## 13 4.54 0.0304
## 14 4.63 0.0286
## 15 4.66 0.0280Function để fitting curve và thu nhận các tham số đường Langmuir isotherm.
library(ggplot2)
library(gridExtra)
langmuir_curve_v1 <- function (Ce, Qe, ...)
{
x <- Ce
y <- Qe
data <- data.frame(x, y)
fit1 <- (y ~ Qsat * K * x/(1 + K * x))
start1 <- list(Qsat = 1, K = 1)
fit2 <- stats::nls(formula = fit1,
data = data,
start = start1,
algorith = "port",
...)
print("Langmuir Isotherm Nonlinear Analysis")
print(summary(fit2))
rsqq <<- lm(Qe ~ predict(fit2))
print(summary(rsqq))
parshill <- as.vector(coefficients(fit2))
Qsat_ok <- parshill[1L]
K_ok <- parshill[2L]
unclass(summary(rsqq))$r.squared -> r.squared_1
round(r.squared_1, digits = 4) -> r.squared_1
unclass(summary(rsqq))$adj.r.squared -> adj.r.squared_1
round(adj.r.squared_1, digits = 4) -> adj.r.squared_1
round(Qsat_ok, digits = 4) -> Qsat_ok
round(K_ok, digits = 4) -> K_ok
unclass(summary(fit2))$coefficients -> matrix_1
matrix_1[, 2] -> thamso_a
round(thamso_a, digits = 4) -> thamso_a
Qsat_ok_data <- paste(Qsat_ok, "±", thamso_a[1])
K_ok_data <- paste(K_ok, "±", thamso_a[2])
#############
table_langmuir <<- data.frame("Qsat" = Qsat_ok_data,
"K" = K_ok_data,
"R-Squared" = r.squared_1,
"Adj.R-Squared" = adj.r.squared_1)
as.data.frame(t(table_langmuir)) ->> table_langmuir_ok
names(table_langmuir_ok)[1] <<- "Giá trị"
table_langmuir_ok$"Tham số" <<- row.names(table_langmuir_ok)
table_langmuir_ok <<- table_langmuir_ok[, c(2, 1)]
#############
rhs <- function(x) {
(Qsat_ok * K_ok * x/(1 + K_ok * x))
}
ggplot2::theme_set(ggplot2::theme_bw(10))
ggplot2::ggplot(data, ggplot2::aes(x = x, y = y)) + ggplot2::geom_point(color = "blue") +
ggplot2::geom_function(color = "#D35400", fun = rhs) +
ggplot2::labs(x = expression(bold(C[e]~(mg/L))),
y = expression(bold(Q[e]~(mg/L))),
title = "Phương trình đường đẳng nhiệt Langmuir\ncho quá trình hấp phụ ammonium của than sinh học từ trấu",
subtitle = "Qe = (Qsat * K * Ce) / (1 + K * Ce)",
caption = "Duc Nguyen | tuhocr.com modified from PUPAIM package") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
theme_bw() +
theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 1, "cm"))
}Kết quả tính toán đường Langmuir isotherm.
langmuir_curve_v1(ammonium_rhb$Ce, ammonium_rhb$Qe, lower = c(0, 0.0001)) -> plot_1## [1] "Langmuir Isotherm Nonlinear Analysis"
##
## Formula: y ~ Qsat * K * x/(1 + K * x)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Qsat 0.10032 0.02918 3.438 0.00441 **
## K 0.08665 0.03295 2.630 0.02080 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.001607 on 13 degrees of freedom
##
## Algorithm "port", convergence message: relative convergence (4)
##
##
## Call:
## lm(formula = Qe ~ predict(fit2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0022154 -0.0009668 -0.0003820 0.0009352 0.0024802
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0006431 0.0007258 0.886 0.392
## predict(fit2) 0.9684109 0.0428681 22.590 8.13e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.001561 on 13 degrees of freedom
## Multiple R-squared: 0.9752, Adjusted R-squared: 0.9732
## F-statistic: 510.3 on 1 and 13 DF, p-value: 8.127e-12plot_1 + theme(axis.text = element_text(size = 10, face = "bold"),
axis.title = element_text(size = 12, colour = "darkgreen"),
plot.title = element_text(size = 12, colour = "darkgreen", face = "bold")) +
annotation_custom(tableGrob(table_langmuir_ok, rows=NULL),
xmin = 3.43, xmax = 3.43, ymin = 0.00752, ymax = 0.00752)
Nếu ta thay đổi giá trị xlim và
ylim thì sẽ nhìn rõ hơn đường cong hấp phụ Langmuir với giá
trị bão hòa Qe gần bằng 0.1.
plot_1 + xlim(0, 1000) + ylim(0, 0.1) + theme(axis.text = element_text(size = 10, face = "bold"),
axis.title = element_text(size = 12, colour = "darkgreen"),
plot.title = element_text(size = 12, colour = "darkgreen", face = "bold")) +
annotation_custom(tableGrob(table_langmuir_ok, rows=NULL),
xmin = 710, xmax = 710, ymin = 0.018, ymax = 0.018)
Function để fitting curve và thu nhận các tham số đường Freundlich isotherm.
Tương tự như trên, ta modify function này theo phương trình Freundlich isotherm.
library(ggplot2)
library(gridExtra)
freundlich_curve_v1 <- function (Ce, Qe, ...)
{
x <- Ce
y <- Qe
data <- data.frame(x, y)
fit1 <- (y ~ Kf * (x^(1/nf)))
start1 <- list(Kf = 1, nf = 1)
fit2 <- stats::nls(formula = fit1,
data = data,
start = start1,
algorith = "port",
...)
print("Freundlich Isotherm Nonlinear Analysis")
print(summary(fit2))
rsqq <<- lm(Qe ~ predict(fit2))
print(summary(rsqq))
parshill <- as.vector(coefficients(fit2))
Kf_ok <- parshill[1L]
nf_ok <- parshill[2L]
Kf_ok -> Qsat_ok
nf_ok ->K_ok
unclass(summary(rsqq))$r.squared -> r.squared_1
round(r.squared_1, digits = 4) -> r.squared_1
unclass(summary(rsqq))$adj.r.squared -> adj.r.squared_1
round(adj.r.squared_1, digits = 4) -> adj.r.squared_1
round(Qsat_ok, digits = 4) -> Qsat_ok
round(K_ok, digits = 4) -> K_ok
unclass(summary(fit2))$coefficients -> matrix_1
matrix_1[, 2] -> thamso_a
round(thamso_a, digits = 4) -> thamso_a
Qsat_ok_data <- paste(Qsat_ok, "±", thamso_a[1])
K_ok_data <- paste(K_ok, "±", thamso_a[2])
############
table_langmuir <<- data.frame("K" = Qsat_ok_data,
"n" = K_ok_data,
"R-Squared" = r.squared_1,
"Adj.R-Squared" = adj.r.squared_1)
as.data.frame(t(table_langmuir)) ->> table_langmuir_ok
names(table_langmuir_ok)[1] <<- "Giá trị"
table_langmuir_ok$"Tham số" <<- row.names(table_langmuir_ok)
table_langmuir_ok <<- table_langmuir_ok[, c(2, 1)]
############
rhs <- function(x) {
(Kf_ok * (x^(1/nf_ok)))
}
ggplot2::theme_set(ggplot2::theme_bw(10))
ggplot2::ggplot(data, ggplot2::aes(x = x, y = y)) + ggplot2::geom_point(color = "blue") +
ggplot2::geom_function(color = "#D35400", fun = rhs) +
ggplot2::labs(x = expression(bold(C[e]~(mg/L))),
y = expression(bold(Q[e]~(mg/L))),
title = "Phương trình đường đẳng nhiệt Freundlich\ncho quá trình hấp phụ ammonium của than sinh học từ trấu",
subtitle = "Qe = K * (Ce^(1/n))",
caption = "Duc Nguyen | tuhocr.com modified from PUPAIM package") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
theme_bw() +
theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 1, "cm"))
}freundlich_curve_v1(ammonium_rhb$Ce, ammonium_rhb$Qe, lower = c(0, 0.0001)) -> plot_2## [1] "Freundlich Isotherm Nonlinear Analysis"
##
## Formula: y ~ Kf * (x^(1/nf))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Kf 0.0083249 0.0005031 16.55 4.09e-10 ***
## nf 1.2324467 0.0699890 17.61 1.88e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.001471 on 13 degrees of freedom
##
## Algorithm "port", convergence message: relative convergence (4)
##
##
## Call:
## lm(formula = Qe ~ predict(fit2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0020089 -0.0009137 -0.0005391 0.0010942 0.0023804
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0001284 0.0007005 0.183 0.857
## predict(fit2) 0.9936304 0.0413479 24.031 3.7e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.001469 on 13 degrees of freedom
## Multiple R-squared: 0.978, Adjusted R-squared: 0.9763
## F-statistic: 577.5 on 1 and 13 DF, p-value: 3.703e-12plot_2 + theme(axis.text = element_text(size = 10, face = "bold"),
axis.title = element_text(size = 12, colour = "darkgreen"),
plot.title = element_text(size = 12, colour = "darkgreen", face = "bold")) +
annotation_custom(tableGrob(table_langmuir_ok, rows=NULL),
xmin = 3.43, xmax = 3.43, ymin = 0.00752, ymax = 0.00752)
Nếu ta thay đổi giá trị xlim và
ylim thì phương trình đường Freundlich thu được không đạt
được điểm bão hòa như phương trình đường Langmuir ở trên.
plot_2 + xlim(0, 1000) + ylim(0, 0.1) + theme(axis.text = element_text(size = 10, face = "bold"),
axis.title = element_text(size = 12, colour = "darkgreen"),
plot.title = element_text(size = 12, colour = "darkgreen", face = "bold")) +
annotation_custom(tableGrob(table_langmuir_ok, rows=NULL),
xmin = 710, xmax = 710, ymin = 0.018, ymax = 0.018)
So sánh với kết quả fitting từ Origin Pro cho thấy việc fitting trên R là chính xác. Right-click để phóng lớn hình ảnh.
Sơ kết
Trên đây là hướng dẫn fitting đường đẳng nhiệt Langmuir và Freundlich trong R phục vụ nghiên cứu về hấp phụ vật liệu. Để học R bài bản từ A đến Z, thân mời Bạn tham gia khóa học “HDSD R để xử lý dữ liệu” để có nền tảng vững chắc về R nhằm tự tay làm các câu chuyện dữ liệu của riêng mình!
ĐĂNG KÝ NGAY:
https://www.tuhocr.com/register
