BÀI TIỂU LUẬN TÌM HIỂU VỀ NLME
Phạm Thị Thu Thảo - Huỳnh Trọng Phúc
2023-06-28
1 CHƯƠNG 1: Giới thiệu
1.1 Đặt vấn đề
Trong quá trình phân tích và xử lý dữ liệu thì chắc ai cũng có thể gặp những khó khăn trong quá trình làm như xử lý dữ liệu thiếu làm ảnh hướng đến quá trình thực hành và đánh giá dữ liệu, ví dụ như “missForest” hay “mice” để có thể ước lượng các giá trị bị thiếu. Tuy nhiên, việc xử lý dữ liệu còn thiếu là một thử thách và đòi hỏi phân tích dữ liệu, xác định mô hình phù hợp hay không? Để lựa chọn mô hình có phù hợp cũng là một thách thức đối với người làm. Một điều quan trọng của mô hình tuyến tính là biến phụ thuộc và các biến độc lập có phân phối chuẩn không. Nó ảnh hưởng đến ước lượng hệ số hồi quy và giả thuyết thống kê, sai số dư trong phân phối chuẩn cũng là một vấn đề để tâm, vì nó có thể ảnh hưởng đến việc kiểm định giả thuyết, độ tin cậy của kết quả.
Để giải quyết các vấn đề trên thì nhóm chúng em chọn package nlme bởi vì nó cung cấp những công cụ mạnh mẽ để xây dựng và phân tích mô hình tuyến tính hỗn hợp, giúp mô hình hóa sự biến đổi không đồng nhất và mối quan hệ không tuyến tính giữa các biến độc lập và biến phụ thuộc. Xử lý được dữ liệu lặp lại và phụ thuộc: Package nlme được thiết kế để xử lý dữ liệu lặp lại và phụ thuộc, các quan sát không độc lập với nhau. Điều này rất hữu ích trong nghiên cứu lâm sàng, thí nghiệm đo lường lặp lại và các mô hình dữ liệu dạng bảng. Ngoài ra, package nlme cho phép mô hình hóa cấu trúc phân tán của các quan sát không cùng phương sai hoặc tương quan thành mô hình tạo ra ước lượng chính xác hơn và điều chỉnh cho sự biến đổi không đồng nhất dữ liệu. Nó còn giúp chúng ta có những công cụ đánh giá kiểm định mô hình, mô hình có phù hợp không bằng các kiểm định thống kê và các chỉ số đánh giá hiệu suất mô hình và biểu đồ sai số.
1.2 Lý do chọn vấn đề
Hiện nay, ngôn ngữ R được mọi người sử dụng rất rộng rãi bởi vì nó giúp chúng ta có thể tối ưu cho thống kê, phân tích, vẽ biểu đồ và trình bày dữ liệu, nếu như bạn là một người lập trình viên thì ngôn ngữ lập trình R không thể nào thiếu được khi bạn muốn làm việc với doanh nghiệp, code R có rất nhiều tính năng hữu dụng giúp cho chúng ta có thể ứng dụng trong phân tích dữ liệu của một công ty nào đó. Ngoài ra, thì cộng đồng R rất đông và tích cực ở đó có rất nhiều nhà làm thống kê, phân tích dữ liệu, nhà làm chính sách, phát triển phần mềm, họ có thể giúp mình trả lời các câu hỏi thắc mắc, đây cũng là nơi cung cấp những tài liệu vô cùng quý giá để chúng ta có thể học tập và trao đổi kinh nghiệm. Mặc khác, thì R có những điểm mạnh như sau là chạy code không cần compiler, thực hiện được các phép tính trên vectors, và còn một điểm mạnh đáng chú ý là Statistical-Language được ứng dụng trong thống kê dữ liệu, sinh học và cả di truyền học giúp chúng ta đáp ứng việc hoàn thành tất cả các thuật toán vì code R là loại ngôn ngữ turing-complete.
Để có thể chứng minh được tính ứng dụng trong thống kê dữ liệu thì hôm nay nhóm chúng tôi sử dụng package “nlme” (Linear and Nonlinear Mixed Effects Models) cung cấp các công cụ và phương pháp để phân tích mô hình tuyến tính, phi tuyến tính sử dụng các hiệu ứng kết hợp để phân tích dữ liệu có cấu trúc phân cụm hoặc lặp lại. Package này có thể hỗ trợ loạt các mô hình như hồi quy tuyến tính, mô hình tuyến tính tổng hợp, mô hình tuyến tính đa cấp và nhiều mô hình khác. Package này giúp cho chúng ta có thể phân tích dữ liệu lặp lại như dữ liệu chuỗi thời gian, dữ liệu thu thập. Ngoài ra, còn cung cấp phương pháp ước lượng tham số hiệu quả như ước lượng Generalized Estimating Equations(Phương trình ước lượng tổng quát) và ước lượng tỷ lệ suy giảm trong mô hình tuyến tính hỗn hợp. Package này hỗ trợ mô hình hóa dữ liệu có cấu trúc phân cụm, trong đó dữ liệu được phân loại vào nhóm không gian(dữ liệu trung tâm nghiêm cứu). Nó cung cấp các phương trình tuyến tính đa cấp và mô hình tuyến tính tổng quát để mô hình hóa các hiệu ứng nhóm và cá nhân. Package nlme cung cấp các mô hình ước lượng tham số cho dữ liệu thiếu thông qua phương pháp tối đa khả năng (maximum likelihood) hoặc phương pháp không tham số. Package nlme được phát triển để đáp ứng các yêu cầu phân tích phức tạp trong phân tích dữ liệu có cấu trúc phân cụm hoặc dữ liệu lặp lại.
1.3 Giới thiệu về Package nlme
Package nlme trong R được phát triển bởi Pinheiro và Bates và được xuất bản lần đầu tiên vào năm 1993. “nlme” là viết tắt của “Nonlinear Mixed-Effects Models” (Mô hình tuyến tính hỗn hợp phi tuyến). Cả hai đều là các nhà khoa học và nhà phân tích dữ liệu hàng đầu, có đóng góp đáng kể vào lĩnh vực mô hình hỗn hợp phi tuyến và phân tích dữ liệu lặp lại.
Pinheiro và Bates phát triển package nlme nhằm cung cấp một phương pháp mô hình hóa và phân tích dữ liệu lặp lại và phụ thuộc trong môi trường R. Package này tập trung vào mô hình tuyến tính hỗn hợp, trong đó các mô hình tuyến tính được áp dụng cho dữ liệu có cấu trúc phân tán và các quan sát không độc lập với nhau.
Package nlme đã trở thành một công cụ phổ biến và được sử dụng rộng rãi trong lĩnh vực phân tích dữ liệu lặp lại, nghiên cứu lâm sàng, và các lĩnh vực khác liên quan đến phân tích dữ liệu có cấu trúc phân tán. Đây là một trong những thành tựu nổi bật trong sự nghiên cứu và ứng dụng của cả hai tác giả trong lĩnh vực này
Package nlme được tạo ra dùng để xây dựng và phân tích các mô hình hỗn hợp tuyến tính. Một số hàm thông dụng trong package nlme:
Lme() Hàm này được sử dụng để xây dựng mô hình tuyến tính hỗn hợp. Nó cho phép bạn xác định các thành phần ngẫu nhiên và cấu trúc phụ thuộc đa cấp trong mô hình.
gnls(): Hàm này được sử dụng để xây dựng mô hình phi tuyến tính hỗn hợp. Nó cho phép bạn xác định các hàm liên kết phi tuyến và cấu trúc phụ thuộc đa cấp.
getVarCov() và vcov(): Hai hàm này được sử dụng để tính toán ma trận hiệp phương sai của các tham số ước lượng trong mô hình. getVarCov() tính toán ma trận hiệp phương sai trực tiếp, trong khi vcov() trích xuất ma trận hiệp phương sai từ kết quả mô hình.
ranef(): Hàm này được sử dụng để trích xuất các ước lượng cốt lỗi cho từng nhóm trong mô hình tuyến tính hỗn hợp.
fixed.effects(): Hàm này được sử dụng để trích xuất các ước lượng hiệu chỉnh cho các hiệu ứng cố định trong mô hình tuyến tính hỗn hợp.
residuals(): Hàm này được sử dụng để tính toán các giá trị dư của mô hình tuyến tính hỗn hợp.
anova(): Hàm này được sử dụng để thực hiện phân tích phương sai (ANOVA) cho mô hình tuyến tính hỗn hợp.
1.3.1 Tổng quan về các gói trong package nlme
Package nlme trong R được tạo ra bởi José Pinheiro, Douglas Bates, Saikat DebRoy, Deepayan Sarkar và John C. Team. Tuy nhiên, có một số gói khác đã đóng góp vào sự phát triển của nlme hoặc có liên quan trực tiếp đến nó. Dưới đây là một số gói quan trọng liên quan đến package nlme:
1.3.1.1 Các gói cốt lõi trong nlme
1.effects: cung cấp các công cụ cho phân tích hiệu ứng trong mô hình tuyến tính và trực quan hóa kết quả của mô hình thống kê.
Một số hàm quan trong package effects
effect(): Hàm này được sử dụng để tính toán và trực quan hóa hiệu ứng ước lượng từ mô hình. Nó cho phép người dùng xem sự ảnh hưởng của một hoặc nhiều biến độc lập đối với biến phụ thuộc đã được điều chỉnh cho các biến khác.
allEffects(): Hàm này tính toán tất cả các hiệu ứng ước lượng từ mô hình và tổng hợp chúng thành một đối tượng có thể được sử dụng để trực quan hóa các hiệu ứng.
plot() (trong ngữ cảnh của package effects): Hàm này được sử dụng để vẽ biểu đồ của các hiệu ứng ước lượng. Nó cho phép người dùng tùy chỉnh kiểu biểu đồ, nhãn trục, tiêu đề và các yếu tố khác của biểu đồ.
termplot(): Hàm này được sử dụng để vẽ biểu đồ của các hiệu ứng ước lượng cho một thuật ngữ (term) cụ thể trong mô hình. Nó cho phép người dùng xem sự ảnh hưởng của một biến độc lập cụ thể lên biến phụ thuộc.
eff() và effpoly(): Hai hàm này được sử dụng để tính toán hiệu ứng ước lượng cho các mô hình phi tuyến và cho phép người dùng trực quan hóa chúng.
mod2way(): Hàm này được sử dụng để tính toán và hiển thị hiệu ứng ước lượng trong các mô hình hai chiều, bao gồm hiệu ứng chéo (interaction effects) giữa các biến độc lập.
setContrasts(): Hàm này được sử dụng để thiết lập phương pháp so sánh và kiểm tra giả thiết trong mô hình.
2.glmnet: Cung cấp thuật toán Elastic Net để huấn luyện mô hình hồi quy tuyến tính với điều kiện sự lựa chọn biến tự động.Dưới đây là một số chức năng và lợi ích chính của gói “glmnet”:
Hồi quy Ridge: Gói cho phép bạn xây dựng mô hình hồi quy Ridge, một phương pháp giảm thiểu quá khớp bằng cách thêm hàm mất mát L2 regularization vào mô hình. Điều này giúp giảm hiện tượng quá khớp và cải thiện khả năng tổng quát hóa của mô hình.
Hồi quy Lasso: Gói cung cấp các công cụ để xây dựng mô hình hồi quy Lasso, một phương pháp lựa chọn biến đặc trưng thông qua hàm mất mát L1 regularization. Phương pháp này có thể giúp loại bỏ các biến không quan trọng và tạo ra một mô hình đơn giản chỉ với các biến quan trọng nhất.
Elastic Net: Gói “glmnet” cung cấp cả hai phương pháp ridge và lasso trong một mô hình duy nhất gọi là elastic net. Phương pháp này kết hợp cả hai hàm mất mát L1 và L2 regularization để cân bằng giữa lựa chọn biến đặc trưng và giảm thiểu quá khớp.
Đánh giá và lựa chọn mô hình: Gói cung cấp các công cụ để đánh giá hiệu suất của mô hình dựa trên các độ đo như sai số bình phương trung bình (MSE), hệ số xác định (R-squared) và các độ đo khác. Bạn có thể sử dụng các công cụ này để lựa chọn mô hình tốt nhất dựa trên hiệu suất trên tập dữ liệu kiểm tra.
Xử lý dữ liệu thiếu: Gói “glmnet” hỗ trợ xử lý dữ liệu thiếu trong quá trình xây dựng mô hình. Bạn có thể sử dụng các phương pháp như điền giá trị trung bình, điền giá trị phổ biến nhất hoặc sử dụng các phương pháp khác để giải quyết vấn đề dữ liệu thiếu, giúp tạo ra các mô hình có khả năng tổng quát hóa tốt và hiệu suất dự đoán cao.
3.glm: Cung cấp các phương pháp hồi quy tuyến tính và logistic. Đây là một gói là một gói cung cấp chức năng liên quan đến mô hình tuyến tính tổng quát (Generalized Linear Models - GLMs). GLMs là mô hình thống kê mở rộng của mô hình tuyến tính truyền thống, cho phép mô hình hóa dữ liệu với phân phối không chuẩn và mối quan hệ phi tuyến. Gói “glm” cung cấp các chức năng chính sau:
Hàm “glm()”: Đây là hàm chính để xây dựng mô hình tuyến tính tổng quát. Hàm này cho phép bạn tùy chỉnh các thông số quan trọng như hàm liên kết, phân phối lỗi, các biến giải thích, v.v.
Hàm “predict()”: Hàm này được sử dụng để dự đoán giá trị của biến phụ thuộc dựa trên mô hình GLM đã được xây dựng.
Hàm “summary()”: Hàm này cho phép bạn xem các thông tin thống kê về mô hình GLM đã được xây dựng, bao gồm các hệ số ước lượng, giá trị p, đánh giá chất lượng mô hình, v.v.
Gói “glm” rất hữu ích trong việc mô hình hóa dữ liệu với các biến phụ thuộc không tuân theo phân phối chuẩn hoặc khi quan hệ giữa biến phụ thuộc và biến giải thích không phải là một hàm tuyến tính. Nó cung cấp các phương pháp ước lượng và kiểm định cho các mô hình tuyến tính tổng quát và là một công cụ mạnh mẽ cho phân tích thống kê trong R.
4.gnls: được sử dụng để ước lượng các mô hình tuyến tính tổng quát (generalized linear models - GLMs) không đồng nhất. Gói này cung cấp các hàm để thực hiện việc ước lượng và kiểm định các mô hình GNLS.
các hàm chính trong package gnls:
gnls(): Hàm này được sử dụng để ước lượng mô hình GNLS từ dữ liệu. Nó cho phép bạn chỉ định công thức của mô hình, dữ liệu và loại phân phối của biến phụ thuộc.
summary.gnls(): Hàm này tính toán bảng tóm tắt cho kết quả của một đối tượng gnls đã được ước lượng. Bảng tóm tắt bao gồm thông tin về các tham số ước lượng, sai số tiêu chuẩn, giá trị p và giá trị AIC (Akaike Information Criterion).
predict.gnls(): Hàm này được sử dụng để dự đoán giá trị biến phụ thuộc từ một model gnls đã được ước lượng.
anova.gnls(): Hàm này tính toán bảng ANOVA (Analysis of Variance) cho các mô hình gnls khác nhau. Nó giúp so sánh sự khác biệt giữa các mô hình và xem xét sự ảnh hưởng của các biến độc lập.
plot.gnls(): Hàm này được sử dụng để vẽ đồ thị phân tích cho model gnls đã được ước lượng, bao gồm đồ thị phân tán, đường cong phù hợp và sai số dự báo.
5.MASS cung cấp các phương pháp thống kê cho mô hình tuyến tính hỗn hợp và ước lượng các tham số trong mô hình tuyến tính.
Mô hình hồi quy:
lm.ridge(): Mô hình hồi quy ridge.
stepAIC(): Chọn mô hình tốt nhất bằng tiêu chuẩn AIC.
Phân tích thành phần chính (PCA): princomp(): Phân tích thành phần chính.
Phân tích hỗn hợp (Mixture analysis): Mclust(): Phân tích hỗn hợp Gaussian.
Phân tích phân cụm (Cluster analysis): kmeans(): Phân tích phân cụm K-means.
Phân tích chuỗi thời gian (Time series analysis):
arima(): Mô hình chuỗi thời gian ARIMA.
acf(): Hàm tương quan tự của chuỗi thời gian
Mô hình tuyến tính tổng quát (Generalized linear models): glm(): Mô hình tuyến tính tổng quát.
- Phân tích thích nghi (Adaptive analysis): mvrnorm(): Tạo mẫu ngẫu nhiên từ phân phối đa biến.
6.lattice dùng để cung cấp các công cụ để tạo biểu đồ thống kê, giúp hiển thị và trực quan hóa dữ liệu sau khi xử lý và sử dụng trong package nlme
xyplot(): Tạo đồ thị hai biến (biểu đồ tương quan) với các điểm dữ liệu được hiển thị trên một lưới.
bwplot(): Tạo biểu đồ hình hộp cho một biến phân loại và một biến liên tục.
densityplot(): Tạo đồ thị mật độ của một biến liên tục.
histogram(): Tạo biểu đồ histogram của một biến liên tục.
dotplot(): Tạo biểu đồ chấm cho một biến phân loại.
stripplot(): Tạo biểu đồ dọc (strip plot) cho một biến phân loại và một biến liên tục.
levelplot(): Tạo đồ thị mức (level plot) cho một biến phân loại và hai biến liên tục.
splom(): Tạo biểu đồ ma trận sương sương (scatterplot matrix) cho một tập hợp các biến.
7.Matrix cung cấp các công cụ và hàm xử lí dữ liệu dưới dạng ma trận và ma trận hiệu quả.
Matrix(): Hàm này được sử dụng để tạo và chuyển đổi ma trận. Nó cho phép bạn tạo ma trận từ các đối tượng khác nhau như vector, danh sách, ma trận thưa, v.v.
rbind2(), cbind2(): Các hàm này thực hiện việc gộp (rbind) hoặc ghép (cbind) các ma trận cùng chiều.
dim(), nrow(), ncol(): Các hàm này trả về số hàng (nrow) và số cột (ncol) của ma trận, hoặc cung cấp kích thước (dim) của ma trận.
diag(): Hàm này tạo một ma trận đường chéo từ một vector.
solve(): Hàm này giải hệ phương trình tuyến tính hoặc tìm ma trận nghịch đảo của một ma trận.
t(): Hàm này chuyển đổi ma trận bằng cách hoán đổi hàng và cột.
crossprod(), tcrossprod(): Các hàm này tính tích chéo (cross-product) và tích chéo chuyển vị (transpose cross-product) của hai ma trận.
det(): Hàm này tính định thức của một ma trận.
eigen(): Hàm này tính các giá trị riêng và vector riêng của một ma trận.
svd(): Hàm này thực hiện phân rã giá trị đơn (singular value decomposition) của một ma trận.
1.3.1.2 Các gói phát triển trong package nlme
1.Knitr cung cấp các hàm để tạo và định dạng báo cáo, tài liệu hoặc trang web có chứa mã R và kết quả tính toán
knit(): Hàm này được sử dụng để biên dịch và thực thi mã R trong một tệp markdown. Nó nhận vào một tệp markdown và tạo ra một tệp markdown đã được xử lý với mã R thực thi và kết quả tính toán.
knit2html(): Hàm này dùng để biên dịch và thực thi mã R trong một tệp markdown và tạo ra một tệp HTML chứa kết quả của mã R. Nó tự động chuyển đổi tệp markdown thành HTML và thêm kết quả tính toán vào tệp HTML.
knit2pdf(): Hàm này được sử dụng để biên dịch và thực thi mã R trong một tệp markdown và tạo ra một tệp PDF chứa kết quả của mã R. Nó tự động chuyển đổi tệp markdown thành PDF và thêm kết quả tính toán vào tệp PDF.
kable(): Hàm này được sử dụng để tạo bảng từ dữ liệu trong R. Nó cho phép người dùng định dạng bảng bằng cách chỉ định các tham số như tiêu đề, chú thích và kiểu định dạng.
include_graphics(): Hàm này được sử dụng để chèn hình ảnh vào tài liệu. Nó cho phép người dùng chèn hình ảnh từ các tệp ảnh trong R Markdown và xác định kích thước, vị trí và các thuộc tính khác của hình ảnh.
fig.cap(): Hàm này được sử dụng để thêm chú thích (caption) cho hình ảnh được tạo ra trong tài liệu. Người dùng có thể chỉ định nội dung của chú thích và vị trí hiển thị.
opts_chunk\(set(): Hàm này được sử dụng để đặt các tùy chọn môi trường chunk. Với opts_chunk\)set(), người dùng có thể thiết lập các tùy chọn như định dạng mã R, kích thước hình ảnh, đầu ra hiển thị hay ẩn, vv.
2.Lme4 cung cấp các phương pháp mô hình hóa tuyến tính hỗn hợp khác, được sử dụng để phân tích dữ liệu lặp lại và phân tán và các chức năng để điều chỉnh và phân tích các mô hình hỗn hợp: tuyến tính (lmer), tuyến tính tổng quát (glmer) và phi tuyến tính (nlmer.)
Gói này sử dụng các phương pháp đại số tuyến tính hiệu quả, hiện đại như được triển khai trong gói Eigen và sử dụng các lớp tham chiếu để tránh sao chép quá mức các đối tượng lớn. lme4 bao gồm các khả năng của mô hình hỗn hợp tuyến tính tổng quát (GLMM), thông qua chức năng glmer và không triển khai các tính năng của nlme để lập mô hình phương sai thay đổi và tương quan của phần dư. Thực hiện các hiệu ứng ngẫu nhiên chéo theo cách dễ dàng hơn cho người dùng và nhanh hơn nhiều.
lme4 cung cấp các tiện ích tích hợp để lập hồ sơ khả năng và khởi động tham số và cho phép linh hoạt hơn trong việc chỉ định các chức năng khác nhau để tối ưu hóa các tham số phương sai-hiệp phương sai của hiệu ứng ngẫu nhiên.
lme4 cung cấp mô hình tuyến tính đa cấp (linear mixed-effects model) và là một phiên bản nâng cấp của nlme. lme4 đã đạt được sự phát triển mạnh mẽ và sử dụng cấu trúc dữ liệu được lưu trữ theo cách khác, gọi là “phân vùng hiệu suất” (efficient partitioned data).
Một số hàm quan trong trong lme4
lmer(): Hàm này được sử dụng để xây dựng mô hình tuyến tính hỗn hợp (mixed-effects linear models). Nó cho phép mô hình hóa hiệu ứng cố định và ngẫu nhiên lên kết quả quan sát.
glmer(): Hàm này được sử dụng để xây dựng mô hình logistic hỗn hợp (mixed-effects logistic regression models). Nó thích ứng cho các mô hình khi kết quả là một biến nhị phân.
nlmer(): Hàm này được sử dụng để xây dựng mô hình phi tuyến tính hỗn hợp (mixed-effects nonlinear models). Nó cho phép mô hình hóa mối quan hệ phi tuyến tính giữa biến phụ thuộc và các biến đầu vào.
glmer.nb(): Hàm này dùng để xây dựng mô hình hỗn hợp với phân phối Negative Binomial. Nó thích ứng cho các mô hình khi kết quả là biến đếm với phân phối Negative Binomial.
lmerControl(): Hàm này được sử dụng để tạo ra một đối tượng kiểm soát (control object) cho mô hình tuyến tính hỗn hợp. Nó cho phép người dùng tinh chỉnh các thiết lập và yếu tố liên quan đến ước lượng tham số và tốc độ hội tụ của mô hình.
ranef(): Hàm này được sử dụng để trích xuất hiệu ứng của yếu tố ngẫu nhiên từ mô hình tuyến tính hỗn hợp. Nó cho phép người dùng xem hiệu ứng của yếu tố ngẫu nhiên lên kết quả quan sát.
3.nlmeU: Đây là một gói phát triển mở rộng cho nlme, được sử dụng để thực hiện mô hình tuyến tính đa cấp. Gói này cung cấp các công cụ mở rộng và tính năng bổ sung cho nlme.
4.lmeSplines: Gói này kết hợp các tính năng của nlme và splines để thực hiện mô hình tuyến tính đa cấp sử dụng các spline functions. Điều này cho phép mô hình hồi quy đa cấp được ứng dụng cho dữ liệu không tuyến tính.
1.3.1.3 Các gói hỗ trợ cho package nlme
1.Dplyr: Package này cung cấp các công cụ mạnh mẽ để thực hiện các hoạt động xử lý, truy vấn và biến đổi dữ liệu. Dplyr cho phép bạn thực hiện các thao tác như chọn cột, lọc dữ liệu, sắp xếp, nhóm và tóm tắt dữ liệu.
• Mutate(): mutate() là một hàm cực kì tiện ích trong khi phân tích dữ liệu, vì nó cho phép bạn tạo ra biến số mới tùy thích trong data frame. Bạn có thể dùng mutate cho 2 việc: hoán chuyển biến số hiện có và tạo ra biến số mới hoàn toàn.
• Select(): chọn các biến dựa trên tên của chúng.
• Filter(): chọn các trường hợp dựa trên giá trị của chúng.
• Summarise(): dùng để tóm tắt dữ liệu
• Arrange(): thay đổi thứ tự của các hàng
2.Ggplot2: là một package hỗ trợ visualization rất mạnh trong R. Dựa trên package này ta có thể vẽ được các đồ thị dạng barchart, line, plot, density, candle chart,pie,… và rất nhiều các đồ thị khác.Ngoài ra ggplot2 còn cho phép người dùng tùy chỉnh màu sắc, kích cỡ, theme, title, … để đồ thị được đẹp hơn và cho phép bạn tạo các biểu đồ dựa trên dữ liệu và áp dụng các thay đổi mô phỏng phức tạp và tùy chỉnh .Cấu trúc của ggplot2 được chia rõ ràng làm 2 phần chính.
ggplot(): phần này qui định đồ thị sẽ sử dụng data nào làm đầu vào. Lưu ý data phải có dạng data.frame. Dạng vector sẽ không được support.
geom_(aes(x,y)): Phần này qui định kiểu đồ thị và các trục tọa độ từ dữ liệu đầu vào. Nếu chỉ có ggplot() mà không thêm geom_() thì chúng ta chỉ nhận được background mà không có đồ thị mặc dù data đã được khai báo. Trong geom_() chúng ta phải khai báo thêm trục tọa độ vào các arguments x và y của aes() chẳng hạn như geom_point(aes(x=bienx,y=bieny)).
• Histograms: geom_histogram()
• Biểu đồ cột:* geom_bar()* hoặc geom_col()
• Box plots: geom_boxplot()
• Điểm (vd: biểu đồ phân tán): geom_point()
• Biểu đồ đường: geom_line() hoặc geom_path()
• Đường xu hướng: geom_smooth()
3.Hmisc chứa nhiều chức năng hữu ích cho phân tích dữ liệu, đồ họa cấp cao, hoạt động tiện ích, chức năng tính toán kích thước và sức mạnh mẫu, mô phỏng, nhập và chú thích bộ dữ liệu, gán giá trị bị thiếu, tạo bảng nâng cao, phân cụm biến, thao tác chuỗi ký tự, chuyển đổi đối tượng R sang mã LaTeX và html, mã hóa lại các biến, bộ nhớ đệm, tính toán song song được đơn giản hóa, ước tính thống kê chung về cửa sổ di chuyển và hỗ trợ diễn giải phân tích thành phần chính.
4.Readr: Package này hỗ trợ đọc và ghi dữ liệu từ và đến các định dạng dữ liệu phổ biến như CSV, Excel và các tập tin văn bản. Readr giúp đọc dữ liệu nhanh chóng và chính xác, đồng thời cung cấp các công cụ để xử lý và chuyển đổi dữ liệu . Những hàm quan trọng trong package này gồm:
• Đọc tệp tin dữ liệu có các trường phân tách bằng dấu phẩy: read_csv()
• Đọc tệp tin dữ liệu có các trường phân tách bằng tab (TSV): read_csv()
• Đọc tệp tin dữ liệu có các trường phân cách bằng dấu chấm phẩy (;): read_csv2()
• Đọc tệp tin dữ liệu có các trường phân cách bằng kí tự khác: read_delim()
• Đọc tệp tin dữ liệu có kích thước của trường được quy định trước: read_table()
5.splines cung cấp các công cụ và hàm liên quan để xử lý đường cong spline, từ một hàm phức tạp vẽ một đường cong đơn giản hơn.
1.3.2 Mục tiêu và lợi ích của package nlme
Mục tiêu chính của package nlme (Nonlinear Mixed Effects) là cung cấp các công cụ và phương pháp để mô hình hóa phân tích dữ liệu không tuyến tính và có yếu tố ngẫu nhiên, đặc biệt tập trung vào mô hình tuyến tính hỗn hợp (mixed-effects models) trong đó có sự kết hợp giữa hiệu ứng cố định (fixed effects) và hiệu ứng ngẫu nhiên (random effects). Ngoài ra, còn cung cấp một bộ khung mạnh mẽ và linh hoạt để xử lý mô hình tuyến tính hỗn hợp phúc tạp. Package này cho phép người dùng mô hình hóa các mối quan hệ không tuyến tính và mô hình tác động của các yếu tố ngẫu nhiên trong dữ liệu. Nó cung cấp các công cụ ước lượng tham số, kiểm định giả thuyết, đánh giá mô hình và khảo sát tương quan giữa các biến. Hơn nữa, còn có thể áp dụng vào một số lĩnh vực như y học, lâm nghiệp, công nghiệp ….
mục tiêu chung: package nlme là cung cấp phương pháp và công cụ để phân tích dữ liệu lặp lại và dư liệu chuỗi thời gian
mục tiêu cụ thể: gói nlme dùng để mô hình tuyến tính và phi tuyến tính, xử lý dữ liệu lặp lại, mô hình chuỗi thời gian, đánh giá và so sánh mô hình
Package nlme (Nonlinear Mixed-Effects Models) giúp chúng ta rất nhiều trong lúc làm việc với các mô hình tuyến tính và các gói có liên quan. Một số lợi ích mà gói này có thể mang lại đó là: ước lượng mô hình hỗn hợp phi tuyến tính trong đó có mối quan hệ phi tuyến tính của biến phụ thuộc và biến giải thích thì mô hình này cho biết mối quan hệ phức tạp và khám phá các hiệu ứng phi tuyến tính trong dữ liệu. Hơn nữa, nó còn cho phép mô hình hóa cốt lỗi,nghĩa là ccs mô hình ngẫu nhiên không hteer giairi tihs được mô hình, cốt lỗi giúp mô hình hóa sự biến đổi không đồng nhất trong dữ liệu và có thể đánh giá mức độ ảnh hưởng của các yếu tố không quan sát được. Ngoài ra, thì gói package anyf giúp chúng ta mô hình hóa cấu trúc phụ thuộc đa cấp độ và xem xét biến động không đồng nhất, không những, nó còn cung cấp phương pháp đánh giá mô hình hỗn hợp bao gồm tính toán giá trị hàm hợp lý log-likelihood, kiểm định giả thuyết, và phân tích phương sai (ANOVA) giữa các mô hình, Phân tích dữ liệu dạng chuỗi thời gian mô hình hóa các xu hướng thay đổi theo thời gian và sự phụ thuộc không đồng nhất trong chuỗi thời gian. Cuối cùng, cung cấp các công cụ để tính toán độ tin cậy thống kê cho các ước lượng mô hình, bao gồm sai số chuẩn, khoảng tin cậy…
1.3.3 Chức năng của package nlme
Xây dựng mô hình tuyến tính và phi tuyến tính: Package nlme cung cấp các hàm để xây dựng mô hình tuyến tính và phi tuyến tính, bao gồm mô hình tuyến tính hỗn hợp (mixed-effects linear models) và mô hình phi tuyến tính hỗn hợp (mixed-effects nonlinear models).
Xử lý dữ liệu có cấu trúc theo thời gian: Đối với dữ liệu dạng chuỗi thời gian, package nlme giúp xử lý sự phụ thuộc giữa các quan sát liên tiếp trong thời gian. Nó cho phép mô hình hóa hiệu ứng thay đổi theo thời gian và đánh giá ảnh hưởng của các biến đầu vào.
Mô hình hỗn hợp: Package nlme cho phép xây dựng mô hình hỗn hợp (mixed-effects models) để mô hình hóa sự ảnh hưởng của các yếu tố ngẫu nhiên và cố định lên kết quả quan sát. Điều này rất hữu ích trong việc mô hình hóa dữ liệu có cấu trúc phân nhóm hoặc dữ liệu lặp lại.
Kiểm định mô hình: Package nlme cung cấp các công cụ để kiểm tra tính phù hợp của mô hình đã được xây dựng, bao gồm kiểm tra giả thiết, đánh giá độ tin cậy của các ước lượng tham số, và so sánh mô hình khác nhau.
Với package nlme, người dùng có thể phân tích các mô hình có cấu trúc phức tạp và đối tượng nghiên cứu có sự phụ thuộc không đồng nhất, giúp tăng khả năng hiểu được quá trình diễn biến dữ liệu theo thời gian và cung cấp các thông tin hữu ích cho việc ra quyết định và dự báo.
1.4 Đối tượng nghiên cứu của gói nlme
Nghiên cứu dược lý: Package nlme thường được sử dụng để mô hình hóa và phân tích dữ liệu về sự hấp thụ, phân phối và loại trừ của các loại thuốc trong cơ thể con người hoặc động vật. Các nghiên cứu này thường sử dụng dữ liệu lặp lại từ nhiều người hoặc đối tượng trong khoảng thời gian nhất định.
Nghiên cứu y tế: Package nlme cũng có thể được áp dụng trong các nghiên cứu y tế, nơi dữ liệu được thu thập từ các bệnh nhân hoặc đối tượng y tế khác qua nhiều lần quan sát. Ví dụ, nghiên cứu về hiệu quả điều trị, tiến triển bệnh, hoặc ảnh hưởng của các yếu tố y tế đến kết quả sức khỏe có thể sử dụng package nlme để mô hình hóa và phân tích dữ liệu.
Nghiên cứu sinh thái học: Package nlme có thể được sử dụng để nghiên cứu các hệ sinh thái và mô hình hóa các yếu tố ảnh hưởng đến sự biến đổi trong quần thể động vật hoặc thực vật qua thời gian. Các nghiên cứu này có thể sử dụng dữ liệu lặp lại hoặc theo dõi từ các khu vực khác nhau hoặc từ cùng một khu vực trong thời gian dài.
Nghiên cứu tâm lý học và xã hội học: Package nlme cung cấp công cụ cho việc mô hình hóa và phân tích dữ liệu theo dõi trong các lĩnh vực như tâm lý học, xã hội học và giáo dục. Ví dụ, nghiên cứu về sự phát triển cá nhân, hiệu quả của phương pháp giảng dạy hoặc tác động của các yếu tố xã hội có thể sử dụng package nlme để phân tích dữ liệu.
1.5 Phương pháp nghiên cứu của gói nlme
Xác định mô hình tuyến tính hỗn hợp: Đầu tiên, bạn cần xác định mô hình tuyến tính hỗn hợp phù hợp cho dữ liệu của mình. Mô hình này bao gồm các thành phần tuyến tính và phi tuyến tính, cùng với các yếu tố ngẫu nhiên để mô hình hóa sự biến đổi giữa các quan sát lặp lại hoặc theo dõi theo thời gian. Bạn có thể sử dụng các công cụ trong package nlme để xây dựng và ước lượng mô hình tuyến tính hỗn hợp.
Xác định hiệp phương sai/covariance structure: Package nlme cung cấp các loại hiệp phương sai/covariance structure khác nhau, như AR(1), unstructured, compound symmetry, etc., để mô hình hóa mối quan hệ giữa các quan sát lặp lại hoặc theo dõi theo thời gian. Bạn cần xác định hiệp phương sai/covariance structure phù hợp cho dữ liệu của mình dựa trên kiến thức chuyên môn và kiểm định thống kê.
Estimation and inference: Sau khi xác định mô hình và hiệp phương sai/covariance structure, bạn có thể sử dụng các phương pháp ước lượng như Maximum Likelihood (ML) hoặc Restricted Maximum Likelihood (REML) để ước lượng các tham số của mô hình. Package nlme cung cấp các chức năng để tiến hành việc ước lượng này và tính toán các giá trị p-value, khoảng tin cậy và các chỉ số khác để kiểm tra giả thiết và đánh giá mô hình.
Kiểm định mô hình: Bạn cần kiểm tra tính phù hợp của mô hình bằng cách sử dụng các phân phối chuẩn hóa của các hệ số ước lượng, các đồ thị kiểm tra tuân theo phân phối, kiểm định tương quan và các phương pháp khác. Package nlme cung cấp các công cụ để thực hiện các kiểm định này và đánh giá tính phù hợp của mô hình.
Phân tích và diễn giải kết quả: Cuối cùng, sau khi ước lượng và kiểm định mô hình, bạn có thể phân tích và diễn giải kết quả từ mô hình. Package nlme cung cấp các phương pháp để trực quan hóa dữ liệu, tính toán các giá trị dự báo, kiểm tra hiệu quả của mô hình và tiến hành phân tích nhóm.
1.6 Đóng góp của gói nlme
Mô hình hóa dữ liệu lặp lại: Package nlme cung cấp các công cụ mạnh mẽ để xử lý và mô hình hóa dữ liệu lặp lại từ các quan sát được thực hiện trên cùng một đối tượng hoặc nhóm đối tượng. Điều này rất hữu ích trong việc nghiên cứu sự ảnh hưởng của các yếu tố ngẫu nhiên và kiểm soát các yếu tố confounding.
Mô hình hóa dữ liệu theo dõi theo thời gian: Package nlme hỗ trợ mô hình hóa dữ liệu theo dõi theo thời gian, nơi các quan sát được thu thập liên tục qua thời gian. Điều này cho phép các nhà nghiên cứu nghiên cứu sự thay đổi và tiến triển của các biến quan trọng qua thời gian.
Xử lý hiệp phương sai/Covariance structure: Package nlme cung cấp một loạt các hiệp phương sai/covariance structure để mô hình hóa mối quan hệ giữa các quan sát lặp lại hoặc theo dõi theo thời gian. Điều này cho phép người dùng tùy chỉnh và lựa chọn hiệp phương sai phù hợp với dữ liệu của họ.
Phân tích và diễn giải kết quả: Package nlme cung cấp các công cụ để phân tích và diễn giải kết quả từ mô hình, bao gồm việc tính toán giá trị dự báo, kiểm tra hiệu quả của mô hình và phân tích nhóm. Điều này giúp người dùng có cái nhìn sâu hơn về mối tương quan và tác động của các yếu tố nghiên cứu.
Mở rộng và linh hoạt: Package nlme có khả năng linh hoạt trong việc xử lý nhiều dạng dữ liệu và hỗ trợ cho các phân tích phức tạp. Nó cũng có khả năng mở rộng với sự phát triển của các phương pháp mới và nghiên cứu tiên tiến.
2 CHƯƠNG 2: Xây dựng mô hình tuyến tính hỗn hợp
2.1 Thử nghiệm
Nghiên cứu có 19 chuột, được chia thành 2 nhóm: nhóm thứ nhất được cho uống thuốc (n = 9 chuột) và nhóm thứ hai là nhóm chứng (không uống thuốc, gồm 10 chuột). Ở mỗi chuột, nồng độ đường trong máu được đo 4 thời điểm: trước khi uống thuốc (T0), 2 giờ, 3 giờ, và 4 giờ sau khi uống thuốc (tạm kí hiệu T2, T3 và T4). Kết quả của thí nghiệm như sau:
Trong đó:
i là mã số định danh của chuột 1 đến 19 (i = 1, 2, 3, …, 19);
yi là nồng độ glucose đo lường được cho chuột i;
ai là nồng độ glucose trước khi uống thuốc của chuột i;
bi là tốc độ giảm glucose của chuột i.
Tốc độ giảm glucose tùy thuộc vào thời gian và thời gian có thể tạm kí hiệu bằng T. Ở đây T = 0, 2, 3, và 4.
Đây là mô hình hồi qui tuyến tính (linear regression model) mà có lẽ các bạn đã từng biết qua. Đối với các con chuột có id từ 1 đến 19 chúng ta có thể ước tính thông số a và b trên bằng lệnh R như sau:
Tạo dữ liệu
library(nlme)
data <- data.frame(
Treatment = c(rep("Test", 9), rep("Control", 10)),
Id = 1:19,
T0 = c(5.9, 5.3, 4.6, 6.2, 6.0, 6.4, 7.6, 5.9, 7.5, 6.2, 6.9, 5.6, 5.1, 5.7, 5.0, 5.2, 7.7, 8.0, 7.7),
T2 = c(3.9, 4.7, 3.7, 4.6, 5.4, 4.7, 4.1, 3.1, 6.1, 5.3, 5.6, 4.7, 3.9, 4.7, 4.0, 4.2, 6.2, 5.8, 5.0),
T3 = c(3.9, 3.5, 3.3, 4.3, 5.2, 4.8, 3.8, 3.6, 5.4, 4.9, 5.9, 4.6, 2.9, 4.3, 3.5, 4.0, 6.1, 6.5, 6.3),
T4 = c(3.6, 3.2, 3.2, 3.9, 4.8, 4.3, 4.1, 3.3, 4.6, 4.5, 5.9, 4.0, 2.9, 4.6, 3.3, 3.8, 5.7, 6.0, 6.2)
)Lệnh Treatment giúp chúng ta có thể chia nhóm: Nhóm điều trị, có hai giá trị “Test” (nhóm được thử nghiệm, những con chuột uống thuốc) và “Control” (nhóm chứng những con chuột không uống thuốc).
Id: Mã số định danh cho mỗi chuột.
T0, T2, T3, T4: Các giá trị nồng độ đường trong máu đo được tại các thời điểm T0, T2, T3 và T4 tương ứng.
Để thực hiện các phân tích thống kê, ví dụ như kiểm tra sự khác biệt giữa nhóm điều trị và nhóm chứng, tìm hiểu biểu đồ thay đổi nồng độ đường trong máu theo thời gian, và đánh giá tác động của thuốc (nếu có) lên nồng độ đường trong máu. Điều này có thể được thực hiện bằng cách sử dụng hàm lm() để lập mô hình hồi quy tuyến tính.
Dựa vào thời gian và nồng độ glucose thì chúng ta có thể biết .
2.1.1 Con chuột có id = 1
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 5.6171 -0.5743Do đó, a1 = 5.6171 và b1 = -0.5743. Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 5.6171 mmol/L, và giảm khoảng 0.5743 mmol/L mỗi giờ (60 phút) sau khi uống thuốc. Phân tích tương tự chúng ta sẽ có ước số cho các con chuột còn lại.
Tính toán tốc độ thay đổi của nồng độ glucose:
## [1] -1.0 0.0 -0.3Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.5, tức là nồng độ glucose giảm 0.25 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ không đổi , tức là nồng độ glucose không đổi trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm đi 0.3 , tức là nồng độ glucose giảm 0.3 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 1 giảm dần theo thời gian sau 2 giờ thì con chuột thứ nhất giảm từ 5.9 xuống 3.9 sau đó thì nồng độ glucose giữ nguyên cho tới 3 giờ, 1 tiếng sau thì nồng độ glucose từ 3.9 giảm còn 3.6, cho thấy nồng độ glucose giảm dần theo thời gian.
2.1.2 Con chuột có id = 2
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 5.4286 -0.5571Sau khi hồi quy thì con chuột thứ 2 có a1 = 5.4286 và b1 = -0.557. Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 5.4286 mmol/L, và giảm khoảng 0.557 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose:
## [1] -0.3 -1.2 -0.3Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.3, tức là nồng độ glucose giảm 0.15 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ giảm đi 1.2 , tức là nồng độ glucose giảm đi 1.2 trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm đi 0.3 , tức là nồng độ glucose giảm 0.3 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu:
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét nồng độ glucose của con chuột số 2 trước khi uống thuốc là 5.3 sau 2 giờ thì nồng độ glucose đã giảm đi còn 4.7 tức giảm đi 0.6, sau 1 giờ tiếp theo thì đã giảm đi còn 3.5 ( đã giảm đi 1,2 so với lúc giờ thứ 2), sau giờ thứ 3 thì nồng độ glucose giảm còn 3.2
2.1.3 Con chuột id = 3
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 4.5229 -0.3657Con số 3 có a1= 4.5229 và b1 = -0.3657. Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 4.5229 mmol/L, và giảm khoảng 0.3657 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -0.45 -0.40 -0.10Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.45, tức là nồng độ glucose giảm 0.225 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ giảm đi 0.4 , tức là nồng độ glucose giảm đi 0.4 trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm đi 0.1 , tức là nồng độ glucose giảm 0.1 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét nồng độ glucose của con chuột số 3 trước khi uống thuốc là 4.6 sau 2 giờ thì nồng độ glucose đã giảm đi còn 3.7 (tức giảm bớt đi 0.9), sau 1 giờ tiếp theo thì đã giảm đi còn 3.3 (giảm đi 0.4), sau giờ thứ 3 thì nồng độ glucose giảm còn 3.2 (giảm đi 0,1)
2.1.4 Con chuột id = 4
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 6.0486 -0.5771Con số 4 có a1 = 6.0486 và b1 = - 0.5771, Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 6.0486 mmol/L, và giảm khoảng 0.5771 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -0.8 -0.3 -0.4Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.8, tức là nồng độ glucose giảm 0.4 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ giảm đi 0.3 , tức là nồng độ glucose giảm đi 0.3 trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm đi 0.4 , tức là nồng độ glucose giảm 0.4 đơn vị trong khoảng thời gian từ T3 đến T4.
-Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 4 giảm dần theo thời gian sau 2 giờ thì con chuột thứ nhất giảm từ 6.2 xuống 4.6 (giảm 1.6 ), 1 tiếng sau thì nồng độ glucose từ 4.6 giảm còn 4.3 (giảm đi 0.3), 1 tiếng sau nữa (4h) thì giảm xuống còn 3.9 (giảm 0.4)
2.1.5 Con chuột id = 5
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 6.0057 -0.2914Con số 5 có a1 = 6.0057 và b1 = -0.2914 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 6.0057 mmol/L, và giảm khoảng 0.2914 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose:
## [1] -0.3 -0.2 -0.4Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.3, tức là nồng độ glucose giảm 0.15 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ giảm đi 0.2 , tức là nồng độ glucose giảm đi 0.2 trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm đi 0.4 , tức là nồng độ glucose giảm 0.4 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 5 giảm dần theo thời gian sau 2 giờ thì con chuột thứ nhất giảm từ 6.0 xuống 5.4 (giảm 0.6 ), 1 tiếng sau thì nồng độ glucose từ 5.4 giảm còn 5.2 (giảm đi 0.2), 1 tiếng sau nữa (4h) thì giảm xuống còn 4.8 (giảm 0.4).
2.1.6 Con chuột id = 6
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 6.1943 -0.5086Con số 6 có a1 = 6.1943 và b1 = -0.5086 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 6.1943 mmol/L, và giảm khoảng - 0.5086 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -0.85 0.10 -0.50Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.85, tức là nồng độ glucose giảm 0.4225 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ tăng lên 0.1 , tức là nồng độ glucose tăng lên 0.1 trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm đi 0.5 , tức là nồng độ glucose giảm 0.5 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 6 giảm dần theo thời gian sau 2 giờ thì con chuột thứ nhất giảm từ 6.4 xuống 4.7 (giảm 1.7 ), 1 tiếng sau thì nồng độ glucose từ 4.7 tăng lên 4.8 (tăng lên 0.1), 1 tiếng sau nữa (4h) thì giảm xuống còn 4.3 (giảm 0.5).
2.1.7 Con chuột id = 7
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 6.9829 -0.9257Con số 7 có a1 = 6.9829 và b1 = -0.9257 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 6.9829 mmol/L, và giảm khoảng 0.9257 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -1.75 -0.30 0.30Từ T0 đến T2: Có một tốc độ thay đổi giảm 1.75, tức là nồng độ glucose giảm 0.875 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ giảm đi 0.3 , tức là nồng độ glucose giảm đi 0.3 trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi tăng lên 0.3 , tức là nồng độ glucose tăng 0.3 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 7 giảm dần theo thời gian sau 2 giờ thì con chuột thứ nhất giảm từ 6.4 xuống 4.7 (giảm 1.7), 1 tiếng sau thì nồng độ glucose từ 4.7 tăng lên 4.8 (tăng lên 0.1), 1 tiếng sau nữa (4h) thì giảm xuống còn 4.3 (giảm 0.5).
2.1.8 Con chuột id = 8
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 5.4086 -0.6371Con số 8 có a1 = 5.4086 và b1 = -0.6371 .Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 6.1943 mmol/L, và giảm khoảng 0.5086 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -1.4 0.5 -0.3Từ T0 đến T2: Có một tốc độ thay đổi giảm 1.4, tức là nồng độ glucose giảm 0.7 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ tăng lên 0.5 , tức là nồng độ glucose tăng 0.5 trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm 0.3 , tức là nồng độ glucose giảm 0.3 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 8 giảm dần theo thời gian sau 2 giờ thì con chuột giảm từ 5.9 xuống 3.1 (giảm 2.8), 1 tiếng sau thì nồng độ glucose từ 3.1 tăng lên 3.6 (tăng lên 0.5), 1 tiếng sau nữa (4h) thì giảm xuống còn 3.3 (giảm 0.3).
2.1.9 Con chuột id = 9
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 7.52 -0.72Con số 9 có a1 = 7.52 và b1 = -0.72 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 7.52 mmol/L, và giảm khoảng 0.72 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -0.7 -0.7 -0.8Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.7, tức là nồng độ glucose giảm 0.35 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ giảm 0.7 , tức là nồng độ glucose giảm 0.7 trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm 0.8, tức là nồng độ glucose giảm 0.8 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 9 giảm dần theo thời gian sau 2 giờ thì con chuột giảm từ 7.5 xuống 6.1 (giảm 1.4), 1 tiếng sau thì nồng độ glucose từ 6.1 giảm xuống 3.6 (giảm 0.7), 1 tiếng sau nữa (4h) thì giảm xuống còn 4.6 (giảm 0.8).
2.1.10 Con chuột id = 10
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 6.1829 -0.4257Con số 10 có a1 = 6.1829 và b1 = -0.4257 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 6.1829 mmol/L, và giảm khoảng 0.4257 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -0.45 -0.40 -0.40Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.45, tức là nồng độ glucose giảm 0.225 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ giảm 0.4 , tức là nồng độ glucose giảm 0.4 trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm 0.4, tức là nồng độ glucose giảm 0.4 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 10 giảm dần theo thời gian sau 2 giờ thì con chuột giảm từ 6 xuống 3 (giảm 3), 1 tiếng sau thì nồng độ glucose từ 3 tăng lên 3.6 (tăng 0.6), 1 tiếng sau nữa (4h) thì giảm xuống còn 3.4 (giảm 0.2).
2.1.11 Con chuột id = 11
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 6.6343 -0.2486Con số 11 có a1 = 6.6343 và b1 = -0.2486 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 6.6343 mmol/L, và giảm khoảng 0.2486 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -0.65 0.30 0.00Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.65, tức là nồng độ glucose giảm 0.325 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ tăng 0.3 , tức là nồng độ glucose tăng 0.3 trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi là 0, tức là nồng độ glucose không thay đổi đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 11 giảm dần theo thời gian sau 2 giờ thì con chuột giảm từ 6.9 xuống 5.6 (giảm 1.3), 1 tiếng sau thì nồng độ glucose từ 5.6 tăng lên 5.9 (tăng 0.3), 1 tiếng sau nữa (4h) thì vẫn giữ nguyên mức 5.9.
2.1.12 Con chuột id = 12
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 5.58 -0.38Con số 12 có a1 = 5.58 và b1 = -0.38 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 5.58 mmol/L, và giảm khoảng 0.38 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -0.45 -0.10 -0.60Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.45, tức là nồng độ glucose giảm 0.225 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ giảm 0.1, tức là nồng độ glucose giảm 0.1 trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm 0.6, tức là nồng độ glucose giảm 0.6 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 12 giảm dần theo thời gian sau 2 giờ thì con chuột giảm từ 5.6 xuống 4.7 (giảm 0.9), 1 tiếng sau thì nồng độ glucose từ 4.7 giảm còn 4.5 (giảm 0.2), 1 tiếng sau nữa (4h) thì giảm xuống còn 4 (giảm 0.5).
2.1.13 Con chuột id = 13
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 5.0371 -0.5943Con số 13 có a1 = 5.0371 và b1 = -0.5943 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 5.0371 mmol/L, và giảm khoảng 0.5943 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -0.6 -1.0 0.0Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.6, tức là nồng độ glucose giảm 0.3 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ giảm 1, tức là nồng độ glucose giảm 1 đơn vị trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi là 0, tức là nồng độ glucose không thay đổi đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 13 giảm dần theo thời gian sau 2 giờ thì con chuột giảm từ 5.1 xuống 3.8 (giảm 1.3), 1 tiếng sau thì nồng độ glucose từ 3.8 giảm còn 3 (giảm 0.8), 1 tiếng sau nữa (4h) thì vẫn giữ nguyên mức 3.
2.1.14 Con chuột id = 14
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 5.5257 -0.3114Con số 14 có a1 = 5.5257 và b1 = -0.3114 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 5.5257 mmol/L, và giảm khoảng 0.3114 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -0.5 -0.4 0.3Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.5, tức là nồng độ glucose giảm 0.25 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ giảm 0.4, tức là nồng độ glucose giảm 0.4 đơn vị trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi tăng 0.3, tức là nồng độ glucose tăng 0.3 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 14 giảm dần theo thời gian sau 2 giờ thì con chuột giảm từ 5.8 xuống 4.7 (giảm 1.1), 1 tiếng sau thì nồng độ glucose từ 4.7 giảm còn 4.2 (giảm 0.5), 1 tiếng sau nữa (4h) thì tăng từ 4.2 lên 4.6 (tăng 0.4).
2.1.15 Con chuột id = 15
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 4.94 -0.44Con số 15 có a1 = 4.94 và b1 = -0.44 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 4.94 mmol/L, và giảm khoảng 0.44 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -0.5 -0.5 -0.2Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.5, tức là nồng độ glucose giảm 0.25 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ giảm 0.5, tức là nồng độ glucose giảm 0.5 đơn vị trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm 0.2, tức là nồng độ glucose giảm 0.2 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 15 giảm dần theo thời gian sau 2 giờ thì con chuột giảm từ 5.0 xuống 4.0 (giảm 1.0), 1 tiếng sau thì nồng độ glucose từ 4.0 giảm còn 3.5 (giảm 0.5), 1 tiếng sau nữa (4h) thì giảm từ 3.5 xuống còn 3.2 (giảm 0.3). Chứng tỏ con chuột có id = 15 sẽ giảm dần theo thời gian.
2.1.16 Con chuột id = 16
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 5.0971 -0.3543Con số 16 có a1 = 5.0971 và b1 = -0.3543 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 5.0971 mmol/L, và giảm khoảng 0.3543 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -0.5 -0.2 -0.2Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.5, tức là nồng độ glucose giảm 0.25 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2.
Từ T2 đến T3: Có một tốc độ giảm 0.2, tức là nồng độ glucose giảm 0.2 đơn vị trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm 0.2, tức là nồng độ glucose giảm 0.2 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 16 giảm dần theo thời gian sau 4 giờ thử nghiệm thì nó giảm từ 5.2 xuống còn 3.8 (giảm đi 1.4 đơn vị).
2.1.17 Con chuột id = 17
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 7.5371 -0.4943Con số 17 có a1 = 7.5371 và b1 = -0.4943 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 7.5371 mmol/L, và giảm khoảng 0.4943 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -0.75 -0.10 -0.40Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.75, tức là nồng độ glucose giảm 0.375 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2
Từ T2 đến T3: Có một tốc độ giảm 0.1, tức là nồng độ glucose giảm 0.1 đơn vị trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm 0.4, tức là nồng độ glucose giảm 0.4 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 17 giảm dần theo thời gian sau 4 giờ thử nghiệm thì nó giảm từ 7.7 xuống còn 5.5 (giảm đi 2.2 đơn vị).
2.1.18 Con chuột id = 18
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 7.6229 -0.4657Con số 18 có a1 = 7.5371 và b1 = -0.4943 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 7.5371 mmol/L, và giảm khoảng 0.4943 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -1.1 0.7 -0.5Từ T0 đến T2: Có một tốc độ thay đổi giảm 0.75, tức là nồng độ glucose giảm 0.375 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2
Từ T2 đến T3: Có một tốc độ giảm 0.1, tức là nồng độ glucose giảm 0.1 đơn vị trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm 0.4, tức là nồng độ glucose giảm 0.4 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 18 giảm dần theo thời gian sau 4 giờ thử nghiệm thì nó giảm từ 7.7 xuống còn 5.5 (giảm đi 2.2 đơn vị).
2.1.19 Con chuột id = 19
##
## Call:
## lm(formula = y ~ T)
##
## Coefficients:
## (Intercept) T
## 7.0714 -0.3429Con số 19 có a1 = 7.0714 và b1 = -0.3429 , Nói cách khác ước số nồng độ glucose lúc ban đầu của chuột 1 là 7.0714 mmol/L, và giảm khoảng 0.3429 mmol/L mỗi giờ (60 phút) sau khi uống thuốc.
Tính toán tốc độ thay đổi của nồng độ glucose
## [1] -1.35 1.30 -0.10Từ T0 đến T2: Có một tốc độ thay đổi giảm 1.35, tức là nồng độ glucose giảm 0.675 đơn vị glucose/giờ trong khoảng thời gian từ T0 đến T2
Từ T2 đến T3: Có một tốc độ tăng 1.3, tức là nồng độ glucose tăng 1.3 đơn vị trong khoảng thời gian từ T2 đến T3.
Từ T3 đến T4: tốc độ thay đổi giảm 0.1, tức là nồng độ glucose giảm 0.1 đơn vị trong khoảng thời gian từ T3 đến T4.
Vẽ biểu đồ sự dao động của nồng độ glucose lúc ban đầu
plot(T, y, type = "o", pch = 16, col = "blue", xlab = "Thời gian", ylab = "Nồng độ glucose",
main = "Biểu đồ sự dao động của nồng độ glucose")
- Nhận xét: Con chuột id số 19 giảm dần theo thời gian sau 2 giờ thì con chuột giảm từ 7.7 xuống 5.0 (giảm 2.7), 1 tiếng sau thì nồng độ glucose từ 5.0 tăng lên 6.3 (tăng 1.3), 1 tiếng sau nữa (4h) thì giảm từ 6.3 xuống còn 6.0 (giảm 0.3).
2.2 Mô hình tuyến tính hỗn hợp
- Dữ liệu ban đầu
## mixtools package, version 2.0.0, Released 2022-12-04
## This package is based upon work supported by the National Science Foundation under Grant No. SES-0518772 and the Chan Zuckerberg Initiative: Essential Open Source Software for Science (Grant No. 2020-255193).T0 <- c(5.9, 5.3, 4.6, 6.2, 6.0, 6.4, 7.6, 5.9, 7.5, 6.2, 6.9, 5.6, 5.1, 5.7, 5.0, 5.2, 7.7, 8.0, 7.7)
T2 <- c(3.9, 4.7, 3.7, 4.6, 5.4, 4.7, 4.1, 3.1, 6.1, 5.3, 5.6, 4.7, 3.9, 4.7, 4.0, 4.2, 6.2, 5.8, 5.0)
T3 <- c(3.9, 3.5, 3.3, 4.3, 5.2, 4.8, 3.8, 3.6, 5.4, 4.9, 5.9, 4.6, 2.9, 4.3, 3.5, 4.0, 6.1, 6.5, 6.3)
T4 <- c(3.6, 3.2, 3.2, 3.9, 4.8, 4.3, 4.1, 3.3, 4.6, 4.5, 5.9, 4.0, 2.9, 4.6, 3.3, 3.8, 5.7, 6.0, 6.2)- Tạo dữ liệu dạng long format
##
## Attaching package: 'dplyr'## The following object is masked from 'package:nlme':
##
## collapse## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
data_long <- data.frame(
Treatment = rep(c("Test", "Control"), each = 19),
T = c(rep(0, 19), rep(2, 19), rep(3, 19), rep(4, 19)),
Glucose = c(T0, T2, T3, T4)
)- Thêm cột id để nhận dạng từng chuột trong dữ liệu
- Xem dữ liệu đã được sắp xếp dạng long format
## Treatment T Glucose id
## 1 Test 0 5.9 1
## 2 Test 0 5.3 2
## 3 Test 0 4.6 3
## 4 Test 0 6.2 4
## 5 Test 0 6.0 5
## 6 Test 0 6.4 6
## 7 Test 0 7.6 7
## 8 Test 0 5.9 8
## 9 Test 0 7.5 9
## 10 Test 0 6.2 10
## 11 Test 0 6.9 11
## 12 Test 0 5.6 12
## 13 Test 0 5.1 13
## 14 Test 0 5.7 14
## 15 Test 0 5.0 15
## 16 Test 0 5.2 16
## 17 Test 0 7.7 17
## 18 Test 0 8.0 18
## 19 Test 0 7.7 19
## 20 Control 2 3.9 1
## 21 Control 2 4.7 2
## 22 Control 2 3.7 3
## 23 Control 2 4.6 4
## 24 Control 2 5.4 5
## 25 Control 2 4.7 6
## 26 Control 2 4.1 7
## 27 Control 2 3.1 8
## 28 Control 2 6.1 9
## 29 Control 2 5.3 10
## 30 Control 2 5.6 11
## 31 Control 2 4.7 12
## 32 Control 2 3.9 13
## 33 Control 2 4.7 14
## 34 Control 2 4.0 15
## 35 Control 2 4.2 16
## 36 Control 2 6.2 17
## 37 Control 2 5.8 18
## 38 Control 2 5.0 19
## 39 Test 3 3.9 1
## 40 Test 3 3.5 2
## 41 Test 3 3.3 3
## 42 Test 3 4.3 4
## 43 Test 3 5.2 5
## 44 Test 3 4.8 6
## 45 Test 3 3.8 7
## 46 Test 3 3.6 8
## 47 Test 3 5.4 9
## 48 Test 3 4.9 10
## 49 Test 3 5.9 11
## 50 Test 3 4.6 12
## 51 Test 3 2.9 13
## 52 Test 3 4.3 14
## 53 Test 3 3.5 15
## 54 Test 3 4.0 16
## 55 Test 3 6.1 17
## 56 Test 3 6.5 18
## 57 Test 3 6.3 19
## 58 Control 4 3.6 1
## 59 Control 4 3.2 2
## 60 Control 4 3.2 3
## 61 Control 4 3.9 4
## 62 Control 4 4.8 5
## 63 Control 4 4.3 6
## 64 Control 4 4.1 7
## 65 Control 4 3.3 8
## 66 Control 4 4.6 9
## 67 Control 4 4.5 10
## 68 Control 4 5.9 11
## 69 Control 4 4.0 12
## 70 Control 4 2.9 13
## 71 Control 4 4.6 14
## 72 Control 4 3.3 15
## 73 Control 4 3.8 16
## 74 Control 4 5.7 17
## 75 Control 4 6.0 18
## 76 Control 4 6.2 19- Hàm arrange() trong package dplyr chúng ta có thể sắp xếp data_long trước tiên theo cột id theo thứ tự tăng dần và sau đó theo cột T theo thứ tự tăng dần và in kết quả ra.
## Treatment T Glucose id
## 1 Test 0 5.9 1
## 2 Control 2 3.9 1
## 3 Test 3 3.9 1
## 4 Control 4 3.6 1
## 5 Test 0 5.3 2
## 6 Control 2 4.7 2
## 7 Test 3 3.5 2
## 8 Control 4 3.2 2
## 9 Test 0 4.6 3
## 10 Control 2 3.7 3
## 11 Test 3 3.3 3
## 12 Control 4 3.2 3
## 13 Test 0 6.2 4
## 14 Control 2 4.6 4
## 15 Test 3 4.3 4
## 16 Control 4 3.9 4
## 17 Test 0 6.0 5
## 18 Control 2 5.4 5
## 19 Test 3 5.2 5
## 20 Control 4 4.8 5
## 21 Test 0 6.4 6
## 22 Control 2 4.7 6
## 23 Test 3 4.8 6
## 24 Control 4 4.3 6
## 25 Test 0 7.6 7
## 26 Control 2 4.1 7
## 27 Test 3 3.8 7
## 28 Control 4 4.1 7
## 29 Test 0 5.9 8
## 30 Control 2 3.1 8
## 31 Test 3 3.6 8
## 32 Control 4 3.3 8
## 33 Test 0 7.5 9
## 34 Control 2 6.1 9
## 35 Test 3 5.4 9
## 36 Control 4 4.6 9
## 37 Test 0 6.2 10
## 38 Control 2 5.3 10
## 39 Test 3 4.9 10
## 40 Control 4 4.5 10
## 41 Test 0 6.9 11
## 42 Control 2 5.6 11
## 43 Test 3 5.9 11
## 44 Control 4 5.9 11
## 45 Test 0 5.6 12
## 46 Control 2 4.7 12
## 47 Test 3 4.6 12
## 48 Control 4 4.0 12
## 49 Test 0 5.1 13
## 50 Control 2 3.9 13
## 51 Test 3 2.9 13
## 52 Control 4 2.9 13
## 53 Test 0 5.7 14
## 54 Control 2 4.7 14
## 55 Test 3 4.3 14
## 56 Control 4 4.6 14
## 57 Test 0 5.0 15
## 58 Control 2 4.0 15
## 59 Test 3 3.5 15
## 60 Control 4 3.3 15
## 61 Test 0 5.2 16
## 62 Control 2 4.2 16
## 63 Test 3 4.0 16
## 64 Control 4 3.8 16
## 65 Test 0 7.7 17
## 66 Control 2 6.2 17
## 67 Test 3 6.1 17
## 68 Control 4 5.7 17
## 69 Test 0 8.0 18
## 70 Control 2 5.8 18
## 71 Test 3 6.5 18
## 72 Control 4 6.0 18
## 73 Test 0 7.7 19
## 74 Control 2 5.0 19
## 75 Test 3 6.3 19
## 76 Control 4 6.2 19Chúng ta có thể thấy các giá trị được sắp xếp theo id tăng dần và thời gian tăng dần của từng con chuột.
Sử dụng hàm lm() để hồi quy tuyến tính với mối quan hệ giữa Glucose và T.
## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmList##
## Call:
## lm(formula = Glucose ~ T, data = data_long)
##
## Coefficients:
## (Intercept) T
## 6.050 -0.485Kết quả cho thấy:
Hệ số chặn là giá trị của biến phụ thuộc ( Glucose) khi biến độc lập (T) bằng 0. Trong trường hợp này, khi T bằng 0, giá trị ước tính Glucose là 6,050.
Hệ số của T (-0,485) biểu thị sự thay đổi của biến phụ thuộc ( Glucose) khi biến độc lập tăng một đơn vị T. Vì vậy, cứ tăng một đơn vị trong T, Glucose giá trị ước tính giảm 0,485 đơn vị.
Mô hình hồi quy tuyến tính có thể được biểu diễn bằng phương trình: Glucose = 6.050 - 0.485 * T.
Sử dụng hàm lmer() để xây dựng mô hình hiệu ứng hỗn hợp, trong đó
Glucoselà biến phụ thuộc
T là biến dự đoán
(1 | id) xác định id Biến đổi.
##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer- Chúng ta có thể viết được mô hình hồi quy dựa vào:
Glucose là biến phụ thuộc
T là biến dự đoán
(1 | id) xác định
id Biến đổi
- Trích xuất thông tin chi tiết về mô hình
## Linear mixed model fit by REML ['lmerMod']
## Formula: data_long$Glucose ~ data_long$T + (1 | id)
## Data: data_long
##
## REML criterion at convergence: 165.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6211 -0.4514 -0.0399 0.5125 3.1788
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.8057 0.8976
## Residual 0.2548 0.5048
## Number of obs: 76, groups: id, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.05038 0.23133 26.15
## data_long$T -0.48496 0.03915 -12.39
##
## Correlation of Fixed Effects:
## (Intr)
## data_long$T -0.381Từ kết quả cho thấy
Phương sai ước tính của phần chặn ngẫu nhiên là 0,8057. Giá trị này thể hiện sự thay đổi của mức glucose cơ bản giữa các con chuột hoặc nhóm khác nhau.
Độ lệch chuẩn ước tính của chặn ngẫu nhiên là 0,8976. Độ lệch chuẩn cho biết mức độ đường cơ bản khác nhau giữa các con chuột khác nhau trong nhóm “id”.
Phương sai ước tính của phần dư là 0,2548. Giá trị này thể hiện tính biến thiên của các điểm dữ liệu xung quanh các giá trị phù hợp (giá trị dự đoán) của mô hình.
Độ lệch chuẩn ước tính của phần dư là 0,5048. Độ lệch chuẩn cho biết các phép đo glucose thực tế sai lệch bao nhiêu so với các giá trị phù hợp trong mỗi nhóm.
In thông tin Estimate, Std. Err., và 95% Conf. Interval cho các hệ số
## Estimate Std. Error t value
## (Intercept) 6.0503759 0.23133426 26.15426
## data_long$T -0.4849624 0.03914693 -12.38826Hệ số ước tính cho data_long$T là -0,4849624. Điều này cho thấy sự thay đổi trung bình về mức glucose (mg/l) khi tăng một đơn vị trong thời gian.
Sai số chuẩn của ước tính hệ số cho là data_long$T 0,03914693. Nó thể hiện sự thay đổi của hệ số ước tính trên các mẫu khác nhau
Giá trị t là -12.33826. Nó đo lường có bao nhiêu lỗi tiêu chuẩn mà ước tính khác không. Trong trường hợp này, hệ số cho khác data_long$T không đáng kể
Xem lại thông tin kết quả mô hình
## Linear mixed model fit by REML ['lmerMod']
## Formula: data_long$Glucose ~ data_long$T + (1 | id)
## Data: data_long
## REML criterion at convergence: 165.3081
## Random effects:
## Groups Name Std.Dev.
## id (Intercept) 0.8976
## Residual 0.5048
## Number of obs: 76, groups: id, 19
## Fixed Effects:
## (Intercept) data_long$T
## 6.050 -0.485- Trích xuất ma trận hiệp phương sai của các hiệu ứng ngẫu nhiên
- Xem lại ma trận hiệp phương sai
## Groups Name Std.Dev.
## id (Intercept) 0.89761
## Residual 0.50475- Vẽ biểu đồ mối quan hệ giữa thời gian và nồng độ glucose.
data_plot <- ggplot(data_long, aes(x = T, y = Glucose, color = id, group = id)) +
geom_point(size = 2) +
geom_line(color = "black") +
xlab("time (h)") + ylab("glucose (mg/l)")
print(data_plot)
Trong nghiên cứu này, hai thông số quan trọng cần quan tâm là ai và bi, nhưng chúng dao động giữa các chuột, do đó chúng ta sẽ sử dụng ký hiệu i để nhắc nhở điều này.
Nếu chuột trong nghiên cứu được phân nhóm một cách ngẫu nhiên, chúng ta kỳ vọng rằng nồng độ trung bình ban đầu (trung bình ai) giữa hai nhóm sẽ không khác nhau. Đồng thời, ảnh hưởng của thuốc chủ yếu được phản ánh qua tốc độ tăng glucose giữa hai nhóm (trung bình bi).
Mục tiêu của nghiên cứu là đánh giá sự ảnh hưởng của thuốc lên sự biến đổi của hai thông số này, và viết lại các thông số dưới dạng ký hiệu:
Ai: Nồng độ trung bình ban đầu
Bi: Tốc độ tăng glucose giữa hai nhóm
Chúng ta sử dụng ký hiệu i để chỉ rằng hai thông số này thay đổi giữa các chuột trong nghiên cứu, và việc phân nhóm một cách ngẫu nhiên giúp xác định tác động chủ yếu của thuốc lên biến đổi của chúng.
data_plot <- ggplot(data_long, aes(x = T, y = Glucose, group = id)) +
geom_point(color = "black", size = 2) +
geom_line(color = "black") +
facet_wrap(~ id) +
xlab("time (h)") + ylab("glucose (mg/l)")
print(data_plot)
- Thứ nhất, nồng độ glucose lúc ban đầu (baseline khác nhau giữa các chuột; và thứ hai, tốc độ (rate) giảm glucose cũng khác nhau giữa các chuột. Tại các thời điểm khác nhau thì các con chuột này có hàm lượng glucose khác nhau.
3 CHƯƠNG 3: Xây dựng mô hình phi tuyến tính hỗn hợp
4 CHƯƠNG 4: Các mô hình tuyến tính và phi tuyến tính trong mô hình hóa rừng (FOREST MODELLING)
4.1 Mô hình tuyến tính
4.1.1 Mô hình tuyến tính đơn
Mô hình quan hệ giữa một biến ảnh hưởng, độc lập và một biến phụ thuộc, bị ảnh hưởng, dự báo trong đó yi và xi được đo đạc từ các mẫu i, i = 1, … n. Mô hình tuyến tính đơn biến cho cá thể. Mẫu thứ i ñược viết như sau (Laar và Akca, 2007; Mehtatalo, 2013):
\(y_i= b_0 + b_1x_i+ ԑ_i\)
Trong mô hình trên biến yi gồm hai bộ phận: Một là phụ thuộc vào biến xi theo mô hình tuyến tính và hai là sai số ngẫu nhiên của mô hình \(ԑ_i\). \(b_0\) và \(b_1\) là hai tham số của mô hình tuyến tính đơn.Phương pháp ước lượng các tham số của mô hình chủ yếu là bình phương tối thiểu.
Việc thiết lập mô hình tuyến tính một biến số có thể áp dụng chương trình lập sẵn trong Statgraphics hoặc lập Codes để chạy trong R. Nếu áp dụng Statgraphics thì các chi tiết thống kê, đồ thị chỉ sử dụng theo chương trình lập sẵn, trong khi ñó, nếu sử dụng R sẽ có ñiều kiện tính toán các chỉ tiêu thống kê, đồ thị khác nhau theo mong muốn.
Một ví dụ cho sử dụng mô hình tuyến tính đơn (một biến độc lập) đó là thiếp lập mối quan hệ giữa trữ lượng rừng (M, \(m^3\)/ha) theo tổng tiết diện ngang (BA, \(m^2\)/ha). Quan hệ M = f(BA) có khả năng biểu diễn tốt ở dạng tuyến tính M = a + b × BA. Từ 120 ô mẫu 1000 \(m^2\) rừng trồng tếch ở Tây Nguyên (Bảo Huy và cộng sự, 1998) đã tính toán được các biến số lâm phần tếch gồm M, BA, chiều cao và đường trung bình (Hg (m), Dg (cm)), mật độ (N, cây/ha),… trong Dữ liệu 10 ở phần phụ lục. Sử dụng bộ dữ liệu này để thiết lập mô hình M = a + b × BA theo chương trình thống kê là mã nguồn mở R.
THIẾT LẬP MÔ HÌNH TUYẾN TÍNH MỘT BIẾN:
Sử dụng dữ liệu từ file excel ggg để lập quan hệ tuyến tính M = a + b × BA từ dữ liệu 110 ô mẫu của rừng trồng tếch. Sau đây là các codes để lập và đánh giá mô hình tuyến tính một biến theo phương pháp bình phương tối thiểu trong R, theo chương trình “lm” (Chambers, 1992):
- Gọi các gói package cần thiết
##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine- Đọc dữ liệu từ file excel
- Dùng lệnh view để xem bảng dữ liệu vừa được đọc
Lập mô hình theo lm
Giá trị dự đoán và sai số của mô hình
Tóm tắt kết quả mô hình
Kết quả trên, mục Coeficients cung cấp các giá trị tham số mô hình và kiểm tra bằng tiêu chuẩn t với Pr được cung cấp. Với kết quả này có được mô hình: M = 6.3207 - 0.2815 x s
Với \(R^2adj\) = 0.0205; các tham số đều có giá trị Pr < 0.0000000 , có nghĩa là không tồn tại các tham số.
Bảng kết quả ANOVA cho thấy giá Pr kiểm tra sự tồn tại của R là 0.05381 > 0.00000, có nghĩa là R # 0 rõ rệt, hay R tồn tại, hay tồn lại mối quan hệ giữa M và BA ở mức 1.506%.
Codes trong R tính toán AIC và các sai số
Cacsc chỉ tiêu đánh giá mô hình:
Kết quả tính toán các chỉ tiêu đánh giá, sai số mô hình là
AIC = 398.8678
Bias = -4.208456\(e^-16\); Bias gần bằng 0 có nghĩa là sai số âm và dương của giá trị dự báo là xấp xỉ nhau, hay nói khác mô hình không dự báo quá cao hay quá thấp so với quan sát.
RMSE = 0.7120327 \(m^3\)/ha; Sai số trung phương cho thấy ước tính theo mô hình này sẽ cho sai số trung bình về M là 0.7 \(m^3\)/ha
MAPE = 10.30375 %; cho thấy sai số của mô hình ước tính M là khoảng 10%.
Đồ thị đánh giá mô hình tuyến tính 1 biến:
Đồ thị sai số và Q-Q
Đồ thị của hình trước cho thấy sai số khá phân tán khi giá trị dự báo tăng lên và đồ thị Q-Q thì phân bố có phần nằm xa đường chéo, chứng tỏ mô hình chưa thực sự phù hợp.
Kết quả hình sau cho thấy giá trị dự báo và quan sát khá bám sát nhau, như vậy mô hình khá tốt, tuy nhiên ñồ thị sai số chỉ ra biến động sai số tăng khi giá trị dự báo tăng lên, kết quả này phù hợp với đồ thị Q-Q, cho thấy mô hình tuyến tính mô tả chưa thực sự tốt quan hệ giữa r và s. Vì vậy cần thử nghiệm mô phỏng theo các mô hình phi tuyến tính khác.
Đồ thị mô hình và quan sát
4.1.2 Mô hình tuyến tính đa biến
Trong thực tế, một biến số thuộc y thông thường không chỉ bị ảnh hưởng, chi phối hoặc dự báo thông qua một biến số x mà có thể là nhiều biến số xi. Lúc này cần xem xét thiết lập mối quan hệ theo mô hình tuyến tính ña biến.
Mô hình tuyến tính ña biến ñược biểu diễn như sau (Mehtatalo, 2013): \(Y_i = B_0 + B_1X_{i1} + B_2X_{i2} +...+ B_pX_{ip} + ԑ_i\)
Trong mô hình trên với n dữ liệu quan sát có một biến phụ thuộc và p biến độc lập, ảnh hưởng. \(Y_i\) là biến phụ thuộc với giá trị quan sát thứ i. \(X_{ij}\) là giá trị quan sát thứ i của biến độc lập thứ j, với j đi từ 1 đến p. \(β_j\) là các tham số được ước lượng qua mô hình tuyến tính nhiều lớp, và \(ԑ_i\) sai số của biến phụ thuộc thứ i, có phân bố chuẩn.
4.1.2.1 Lựa chọn các biến số độc lập ảnh hưởng trong mô hình tuyến tính đa biến
Trong thực tế, đôi khi chưa rõ biến phụ thuộc \(Y_i\) bị ảnh hưởng bởi những biến độc lập \(X_{ij}\) nào. Vì vậy sử dụng chỉ số Cp của Mallows (1973) sẽ hỗ trợ cho việc xác định số biến số độc lập ảnh hưởng.
Chỉ số Mallow’ Cp (1973) được sử dụng để lựa chọn số biến số tham gia mô hình tốt nhất trong trường hợp có nhiều biến độc lập nhưng chưa rõ có ảnh hưởng đến \(Y_i\) hay không. Chỉ số Cp càng bé và càng gần với số biến số độc lập p thì đó là các biến độc lập có ảnh hưởng rõ rệt đến biến phụ thuộc; dựa vào đây để xác định p biến số \(X_{ij}\) tham gia mô hình khi có quá nhiều biến số được giả định là có ảnh hưởng đến \(Y_i\).
Nếu một mô hình có p biến số độc lập được lựa chọn từ một tập hợp K> p, chỉ tiêu thống kê Cp được tính toán:
\(C_p = \frac{SSE_p}{S^2} - N + 2P\)
\(SSE_p =\sum\limits_{i=1}^{N}(Y_i - Y_{pi})^2\)
Trong đó: \(Y_{pi}\) là giá trị dự đoán từ giá trị quan sát thứ i là của \(Y_i\) của mô hình có p biến số độc lập; \(S^2\) bình phương trung bình phần dư (residual mean square) sau khi mô hình quan hệ hoàn thành với K biến số độc lập và được ước tính từ sai số trung bình bình phương (mean square error – MSE); N là dung lượng mẫu quan sát.
4.2 Mô hình phi tuyến tính
Mô hình phi tuyến tính có dạng tổng quát (Linton và Hardle, 1998; Mehtatalo, 2013; Picardet al., 2015):
\(y_i= f(x_{1i}, ….., x_{qi}, β_1, ….β_p) + ԑ_i\)
Trong đó \(y_i\) là biến phụ thuộc theo mẫu i, f là dạng hàm phi tuyến, \(x_{qi}\) là biến độc lập thứ q của mẫu i, \(β_p\) là tham số thứ p của mô hình và \(ԑ_i\) là sai số dự báo \(y_i\) của mẫu thứ i. Dạng hàm f phi tuyến tính rất đa dạng, tùy thuộc rất lớn vào mối quan hệ giữa \(y_i\) và \(x_{qi}\) nghiên cứu.
Mối quan hệ giữa các nhân tố cây rừng, lâm phần, các thành phần của hệ sinh thái rừng thường rất phức tạp, khó diễn tả đầy đủ và chính xác theo mô hình tuyến tính. Vì vậy các mô hình phi tuyến từ một đến nhiều biến, hoặc tổ hợp biến, từ đơn giản ở dạng mũ đến phức tạp như mũ, exp nhiều tầng được áp dụng. Mô hình phi tuyến tính được sử dụng rất đa dạng trong mô hình hóa các mối quan hệ trong hệ sinh thái rừng, mô hình hóa quá trình sinh trưởng, sản lượng, sinh khối, carbon rừng.
Đối với quan hệ chiều cao (H) với đường kính ngang ngực (DBH), các quan hệ phi tuyến tính sau thường được sử dụng (Sola et al., 2014):
\(H = 1.3 + a×DBH^b\)
\(H = 1.3 + a×exp(\frac{-b}{DBH+c})\) (Ratkowsky and Giles, 1990)
\(H = 1.3 + a×(1 – exp(-b×DBH^c))\) (Weibull in Yang et al., 1978)
Đối với quan hệ chiều cao bình quân tầng trội (\(H_o\)) theo tuổi (A) để lập biểu cấp năng suất rừng trồng:
- Với rừng trồng tếch (Tectona grandis L.f.), Bảo Huy (1995), Bảo Huy và cộng sự (1998) đã sử dụng mô hình phi tuyến tính Schumacher:
\(H_o = a_i×exp(-b_iA^{-m})\)
- Với rừng trồng trám trắng (Canarium album Raeusch) Bảo Huy và Đào Công Khanh (2008) cũng sử dụng hàm Schumacher:
\(Ho = 1617.456×exp(- 7.72465 ×A^{-0.15})\)
Đối với mô hình sinh trưởng thể tích cây rừng (V) theo tuổi (A), một số mô hình phi tuyến mũ nhiều tầng thường được sử dụng như hàm Schumacher, Korf (Nguyễn Ngọc Lung, 1989). Ví dụ sinh trưởng thể tích loài bằng lăng (Lagerstroemia calyculata Kurz) theo hàm Korf (Bảo Huy, 1993):
\(V = 8.73218×exp( - 25.27369 /A^{0.55211})\)
Đối với phương trình thể tích cây rừng (V) theo các nhân tố điều tra như đường kính ngang ngực (DBH), chiều cao (H), một số dạng phi tuyến power sau được sử dụng; trong đó các biến số có thể là biến đơn hay tổ hợp (Vũ Tiến Hinh, 2012):
\(V = b_o×DBH^{b1}×H^{b2}\)
\(V = b_o×(DBH^2×H)^{b1}\)
Đối với mô hình ước tính sinh khối – carbon cây rừng tự nhiên trên mặt đất (AGB) theo một đến nhiều biến số như DBH, H và khối lượng thể tích gỗ (WD) ở vùng nhiệt đới, Việt Nam, một số mô hình phi tuyến sau thường được sử dụng:
Brown (1997): AGB = exp(-2.134 + 2.530×ln(DBH))
IPCC (2003): AGB = exp(-2.289 + 2.649×ln(DBH)-0.021×(ln(DBH)^2))
Chave et al. (2014): \(AGB = 0.0673×(WD×DBH^2×H)^{0.976}\)
Huy et al., 2016b: AGB cho rừng lá rộng thường xanh vùng Nam Trung Bộ, trong đó CA là diện tích tán lá cây rừng:
\(AGB = 0.61345× (DBH^2×H×WD)^{0.86983} × CA^{0.18834}\)
Huy et al. 2016c thiết lập mô hình AGB cho rừng khộp Việt Nam
\(AGB = 0.06203×DBH^{2.26430}×H^{0.51415}×WD^{0.79456}\)
Việc ước lượng các mô hình phi tuyến tính có thể áp dụng một trong các phương pháp: tuyến tính hóa và bình phương tối thiểu (chương trình “lm” trong R, Chambers, 1992) hoặc phi tuyến bình phương tối thiểu (chương trình “nls” trong R: Nonlinear Least Squares) (Bates và Watts, 1988) hoặc phương pháp phi tuyến của Marquardt (StatPoint-Inc., 2005) hoặc phi tuyến ảnh hưởng phức hợp hợp lý tối da (chương trình “nlme” trong R: Nonlinear Mixed-Effects Models - Maximum Likelihood) (Davidian và Giltinan, 1995; Pinheiro et al., 2014).
4.2.1 Lựa chọn các biến số ảnh hưởng trong mô hình phi tuyến
Cũng như mô hình tuyến tính đa biến, mối quan hệ phi tuyến cũng có dạng từ một biến số ảnh hưởng cho đến nhiều biến. Ví dụ đối với mô hình ước tính sinh khối cây rừng trên mặt đất (AGB), trong khi Chave et al. (2014) lập quan hệ với tối đa 3 biến số là DBH, H và WD,trong khi đó Huy et al. (2016b) đã phát hiện thêm biến CA đã làm tăng độ tin cậy của mô hình.
Chỉ tiêu Mallow’s Cp (1973) được sử dụng để xác định số biến số độc lập \(x_i\) ảnh hưởng có ý nghĩa tới biến phụ thuộc y; trong đó, khi \(C_p\) tiến dến giá trị p gần bằng nhất với số biến số độc lập thì đó là các biến tối ưu ảnh hưởng đến y.
4.2.2 Ước lượng mô hình phi tuyến theo phương pháp tuyến tính hóa
Phương pháp ước lượng mô hình phi tuyến tính phổ biến là tuyến tính hóa nó, thường là logarit hóa nếu mô hình dạng hàm mũ, từ đó ước lượng các tham số của mô hình tuyến tính từ một đến nhiều biến theo phương pháp bình phương tối thiểu.
Một số hàm phi tuyến và được tuyến tính hóa thông qua logarit phổ biến:
\(y = (b_0)×(x_1^{b1})×(x2^{b2})……×(x_n^{bn})×ԑ\), tuyến tính hóa: \(log(y) = log(b_0) + b_1×log(x_1) + b_2×log(x_2) + ……. b_n×log(x_n) + log(ԑ)\)
\(y = (b_0)×exp((x_1^{b1})×(x_2^{b2})……×(x_n^{bn})) ×ԑ\), tuyến tính hóa: \(log(y) = log(b_o) + b_1×x_1 + b_2×x_2 +……. b_n×x_n + log(ԑ)\)
Trong đó y là biến phụ thuộc, \(x_i\) là biến độc lập thứ i đến n, \(b_i\) là tham số theo biến số thứ i đến n của mô hình và ԑ là sai số dự báo y.
Đối với mô hình tuyến tính hóa dạng logarit, khi chuyển về lại dạng nguyên thủy là phi tuyến để ước tính y, cần sử dụng một hệ số điều chỉnh. Sau khi tuyến tính hóa dạng logarit để lượng các tham số của mô hình theo phương pháp bình phương tối thiểu, cần tính giá trị CF; sau đó mô hình được chuyển về dạng nguyên thủy phi tuyến và khi dự đoán y cần nhân thêm hệ số CF. Ví dụ với mô hình power sau khi tính toán các tham số, dự đoán y sẽ là:
\(y =CF× b_o×(x_1^{b1})×(x_2^{b2})……×x_n^{b_n}\)
Trong đó CF được tính:
\(CF = exp(\frac{RSE^2}{2})\)
CF luôn lớn hơn 1. Trong đó RSE (Residual standard error) là sai tiêu chuẩn của phần dư (sai lệch giữa quan sát và dự đoán). Khi RSE càng lớn thì CF càng lớn, có nghĩa mô hình càng có độ tin cậy thấp. Mô hình tốt khi CF càng tiến đến dần 1.
- Thiết lập mô hình phi tuyến được tuyến tính hóa theo logarit sử dụng codes “lm” trong chương trình R:
Code “lm” (Chambers, 1992) thực hiện trong R giúp cho việc ước lượng các mô hình tuyến tính hoặc được tuyến tính hóa từ một đến nhiều biến. Chương trình này ước lượng các tham số của mô hình theo phương pháp bình phương tối thiểu. Đồng thời R còn cho phép chúng ta rất linh hoạt trong việc tính toán các chỉ tiêu thống kê và lập các đồ thị để đánh giá, so sánh các mô hình với nhau. Sau đây giới thiệu lần lượt sử dụng codes “lm” chạy trong R để ước lượng các mô hình phi tuyến được logarit hóa từ một đến nhiều biến độc lập hoặc tổ hợp biến và tính toán các chỉ tiêu sai số, đồ thị.
Tiếp tục sử dụng dữ liệu AGB theo các nhân tố điều tra (DBH, H, WD và CA) của 110 cây rừng mẫu của rừng lá rộng thường xanh vùng Nam Trung Bộ, ước lượng mô hình AGB theo một đến nhiều nhân tố, tổ hợp nhân tố dạng hàm power được logarit hóa theo code “lm” trong R.
Tiến hành ước lượng mô hình phi tuyến power một biến được logarit hóa trong chương trình “lm” của R theo phương pháp bình phương tối thiểu:
\(AGB = b_0 × DBH^{b1} × ԑ\), logarit hóa: \(log(AGB) = log(b_0) + b_1×log(DBH)+log(ԑ)\)
Lập mô hình logarit: \(log(AGB) = a + b × log(DBH)\)
Mô hình hỗn hợp phi tuyến tính hóa theo logarit
Khai báo các gói package cần sử dụng
- Đọc dữ liệu từ file excel
- Dùng lệnh view để xem bảng dữ liệu vừa được đọc
Mô hình tuyến tính logarit hóa (lm) có công thức: log(AGB) = a + b×log(DBH)
Xây dựng mô hình tuyến tính logarit hóa với biến phụ thuộc là “log(AGB)” (logarit của biến AGB) và biến độc lập là “log(DBH)” (logarit của biến DBH). Dữ liệu được sử dụng để xây dựng mô hình là từ bảng dữ liệu “ggg” và dùng hàm summary() in ra thông tin tóm tắt của mô hình tuyến tính logarit hóa vừa được tạo ra
##
## Call:
## lm(formula = log(AGB) ~ log(DBH), data = ggg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.07384 -0.13351 -0.00078 0.18722 0.50357
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.6247 0.1696 -33.17 <2e-16 ***
## log(DBH) 2.4715 0.0321 77.01 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2732 on 108 degrees of freedom
## Multiple R-squared: 0.9821, Adjusted R-squared: 0.9819
## F-statistic: 5930 on 1 and 108 DF, p-value: < 2.2e-16- Phân tích phương sai (ANOVA) cho mô hình tuyến tính
## Analysis of Variance Table
##
## Response: log(AGB)
## Df Sum Sq Mean Sq F value Pr(>F)
## log(DBH) 1 442.57 442.57 5929.9 < 2.2e-16 ***
## Residuals 108 8.06 0.07
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Tính các giá trị dự đoán (fitted values) và sai số (residuals) của mô hình hồi quy tuyến tính
- Tính giá trị của hệ số góc của mô hình lmt1, tức là log(DBH), bằng cách lấy số mũ (exponential) của giá trị hệ số tương ứng.
## (Intercept)
## 0.003607577- kết quả của giá trị exponential của phần tử đầu tiên trong danh sách
các hệ số của mô hình tuyến tính
lmt1là 0.003607577
## log(DBH)
## 2.471502Bạn đã lấy giá trị của hệ số log(DBH) (hay hệ số góc) của mô hình lmt1. Giá trị này là 2.471502, và tên của nó là log(DBH).
Kết quả trên lập được mô hình: \(AGB = 0.106836×DBH^{2.471502}\) với các tham số tồn tại rõ rệt với Pr = <2e-16 < 0.0001
4.2.2.1 Tính các chỉ tiêu thống kê:
- Sử dụng hàm summary() để trích xuất các thông số tóm tắt của mô hình tuyến tính lmt1, bao gồm các thông số về phương sai của các residual trong mô hình. summary(lmt1) trả về một đối tượng chứa các thông số tóm tắt của mô hình và \(\sigma^2\) là giá trị phương sai của residual, tức là sai số dư cuối cùng của mô hình.
## [1] 0.07463452Giá trị phương sai \(\sigma^2\) = 0.07463452 sau khi chạy mô hình lmt1
Viết cách khác cũng có thể tính giá trị phương sai
## [1] 1.038022Kết quả trên cho thấy giá trị phương sai là 1.038022
CF = 1.038022. CF gần bằng 1 cho thấy mô hình được ước lượng khá tốt
Đưa mô hình về dạng nguyên thủy phi tuyến:
lmt1.CF: Hệ số Confidence Factor đã tính trước đó từ mô hình tuyến tính lmt1.
a: Giá trị của hệ số hồi quy tuyến tính (Intercept) đã được tính từ mô hình lmt1.
ggg$DBH: Cột dữ liệu chứa giá trị của biến DBH trong dataframe ggg.
b: Giá trị của hệ số hồi quy tuyến tính log(DBH) đã được tính từ mô hình lmt1.
Lệnh này tính giá trị dự đoán (backtr1.fit) của mô hình tuyến tính logarit hóa (lmt1) cho biến phụ thuộc AGB (lượng gỗ trồng) dựa trên giá trị của biến DBH (đường kính gỗ). Công thức tính giá trị dự đoán là lmt1.CF * a * ggg$DBH^b.
Lệnh này tính sai số dự đoán (backtr1.res) của mô hình tuyến tính logarit hóa (lmt1). Sai số dự đoán là sự khác biệt giữa giá trị thực tế của biến phụ thuộc AGB và giá trị dự đoán backtr1.fit từ mô hình.
Đây là một cột mới được thêm vào trong dataframe ggg, chứa giá trị dự đoán của mô hình tuyến tính logarit hóa (lmt1) cho biến phụ thuộc AGB.
lmt1.CF: Đây là hệ số Confidence Factor đã được tính trước đó từ mô hình tuyến tính lmt1. Hệ số Confidence được sử dụng để tính khoảng tin cậy cho giá trị dự đoán.
a: Giá trị của hệ số hồi quy tuyến tính (Intercept) đã được tính từ mô hình lmt1. Hệ số này là một phần của mô hình tuyến tính logarit hóa và được sử dụng trong công thức tính giá trị dự đoán.
+ggg$DBH: Cột dữ liệu chứa giá trị của biến DBH (đường kính gỗ) trong dataframe ggg. Đây là biến độc lập được sử dụng trong mô hình lmt1.
+b: Giá trị của hệ số hồi quy tuyến tính log(DBH) đã được tính từ mô hình lmt1. Hệ số này là một phần của mô hình tuyến tính logarit hóa và được sử dụng trong công thức tính giá trị dự đoán.
- ggg$backtr1.res: Đây là một cột mới được thêm vào trong dataframe ggg, chứa giá trị sai số dự đoán của mô hình tuyến tính logarit hóa (lmt1). Sai số dự đoán là sự khác biệt giữa giá trị thực tế của biến phụ thuộc AGB và giá trị dự đoán backtr1.fit từ mô hình. Sai số dự đoán là một thước đo cho hiệu suất của mô hình trong việc dự đoán lượng gỗ trồng AGB dựa trên đường kính gỗ DBH.
- Kiểm định sự phù hợp của mô hình (các chỉ tiêu thống kê đánh giá mô hình)
## [1] 0.9353272- R-squared (R hiệu điều chỉnh) là 0.9353272, mô hình phù hợp, có nghĩa là mô hình giải thích được khoảng 93.53% phương sai của biến AGB trong dữ liệu. Và tồn tại với Pr = 2.2e-16
## [1] 0.9334968- R-squared adjusted là 0.9334968, gần giá trị R2, cho thấy mô hình vẫn duy trì mức độ phù hợp tốt sau khi điều chỉnh cho số lượng biến độc lập trong mô hình.
Tính \(AIC = -2L + 2*p\)
trong đó p là số tham số của mô hình
ML.backtr1: Đây là biểu thức tính Maximum Likelihood cho mô hình hồi quy tuyến tính đã fit (backtr1). Công thức tính Log hợp lý cực đại như sau:
Tổng các phần tử trong ngoặc đơn: ggg\(backtr1.res^2/var : Đây là bình phương sai số giữa giá trị dự đoán ggg\)backtr1.fit và giá trị thực tế ggg$AGB, chia cho phương sai của sai số .
log(2*pi): Đây là logarit cơ số e của 2 nhân với số pi (3.141592653589793). Đây là một hằng số trong công thức tính Log hợp lý cực đại.
log(var(ggg$backtr1.res)): Đây là logarit cơ số e của phương sai của sai số cũng là một hằng số trong công thức.
Tổng của tất cả các phần tử trên sẽ là giá trị Log hợp lý cực đại.
AIC.backtr1: Đây là tiêu chuẩn thông tin Akaike (Akaike Information Criterion) cho mô hình đã fit (backtr1). Công thức tính AIC như sau:
-2*ML.backtr1: Đây là giá trị Log hợp lý cực đại nhân với -2.
2*2: Đây là 2 lần số lượng thông số mô hình (trong trường hợp này có 2 thông số: hệ số chặn và hệ số hồi quy).
Tổng của các phần tử trên sẽ là giá trị AIC.
Bias: Đây là giá trị độ lệch trung bình của các giá trị dự đoán ggg\(backtr1.fit so với giá trị thực tế ggg\)AGB. Được tính bằng trung bình của sai số.
RMSE: Đây là giá trị giai đoạn trung bình bình phương của sai số dự đoán giữa mô hình và dữ liệu thực tế. Công thức tính RMSE là căn bậc hai của trung bình của bình phương sai số.
MAPE: Đây là giá trị tỷ lệ phần trăm trung bình của sai số tuyệt đối giữa mô hình và dữ liệu thực tế. Công thức tính MAPE là trung bình của tỷ lệ phần trăm giữa sai số tuyệt đối và giá trị thực tế, nhân với 100 để tính tỷ lệ phần trăm.
ML.backtr1 <- -1/2*sum(ggg$backtr1.res^2/var(ggg$backtr1.res)+log(2*pi)
+log(var(ggg$backtr1.res)))
AIC.backtr1 <- -2*ML.backtr1 +2*2
AIC.backtr1## [1] 2126.817Bias <- mean(ggg$backtr1.res)
RMSE <- sqrt(mean((ggg$backtr1.res)^2))
MAPE <- 100*mean(abs(ggg$backtr1.res)/ggg$AGB)
Bias## [1] 55.53353## [1] 3752.613## [1] 22.61417Nhận xét:
AIC = 1620.2
Giá trị độ lệch trung bình của các giá trị dự đoán so với giá trị thực tế Bias: 55.53353 kg
Giá trị giai đoạn trung bình bình phương của sai số dự đoán giữa mô hình và dữ liệu thực tế RMSE là 3752.613 kg
Giá trị tỷ lệ phần trăm trung bình của sai số tuyệt đối giữa mô hình và dữ liệu thực tế: 22.61417 %
p1 <- ggplot(ggg, aes(x=ggg$backtr1.fit , y=AGB))
p1 <- p1 + geom_point(cex=2)
p1 <- p1 + geom_abline(cex = 1.5, intercept = 0, slope = 1, col="black")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.p1 <- p1 + xlab("Dự đoán AGB (hg)") + ylab("Quan sát AGB (hg)") + theme_bw()+
theme_bw()
p1 = p1 + theme(axis.title.y = element_text(size = rel(1.5)))
p1 = p1 + theme(axis.title.x = element_text(size = rel(1.5)))
p1 <- p1 + theme(plot.title = element_text(size = rel(1.7)))
p1 = p1 + theme(axis.text.x = element_text(size=15))
p1 = p1 + theme(axis.text.y = element_text(size=15))
p1 ## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.
Kết quả ở trên cho thấy giá trị dự đoán và quan sát khá bám sát nhau. Tuy nhiên, mô hình một biến số DBH có sai số rộng và phân tán khi giá trị dự đoán lớn (hay DBH tăng). Điều này có thể được cải thiện nếu áp dụng mô hình có trọng số
p2 <- ggplot(ggg, aes(x=ggg$backtr1.fit, y=ggg$backtr1.res ))
p2 <- p2 + geom_point(cex = 2)
p2 <- p2 + geom_line(cex = 1.5, aes(x=ggg$backtr1.fit, y=0))
p2 <- p2 + xlab("Dự đoán AGB (hg)") + ylab("Sai số (hg)") + theme_bw()+ theme_bw()
p2 = p2 + theme(axis.title.y = element_text(size = rel(1.5)))
p2 = p2 + theme(axis.title.x = element_text(size = rel(1.5)))
p2 <- p2 + theme(plot.title = element_text(size = rel(1.7)))
p2 = p2 + theme(axis.text.x = element_text(size=15))
p2 = p2 + theme(axis.text.y = element_text(size=15))
p2## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.## Warning: Use of `ggg$backtr1.res` is discouraged.
## ℹ Use `backtr1.res` instead.## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.
## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.
## Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.## Warning: Use of `ggg$backtr1.res` is discouraged.
## ℹ Use `backtr1.res` instead.## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.
Để nâng cao độ tin cậy của ước lượng AGB, tiếp tục thử nghiệm tiến hành ước lượng mô hình phi tuyến power mở rộng gồm hai biến DBH và H được logarit hóa trong chương trình “lm” của R theo phương pháp bình phương tối thiểu:
\(AGB = b_o×DBH^{b1}×H^{b2}×ԑ\),
logarit hóa:
\(log(AGB) = log(b_o) + b_1×log(DBH) +b_2×log(H) + log(ԑ)\) Lập mô hình tuyến tính hóa logarit nhiều biến:
##
## Call:
## lm(formula = log(AGB) ~ log(DBH) + log(H), data = ggg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.04288 -0.13952 0.01563 0.18718 0.48347
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.81757 0.30558 -22.31 < 2e-16 ***
## log(DBH) 2.14010 0.07875 27.18 < 2e-16 ***
## log(H) 0.57932 0.12761 4.54 1.48e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2513 on 107 degrees of freedom
## Multiple R-squared: 0.985, Adjusted R-squared: 0.9847
## F-statistic: 3514 on 2 and 107 DF, p-value: < 2.2e-16## Analysis of Variance Table
##
## Response: log(AGB)
## Df Sum Sq Mean Sq F value Pr(>F)
## log(DBH) 1 442.57 442.57 7006.553 < 2.2e-16 ***
## log(H) 1 1.30 1.30 20.609 1.482e-05 ***
## Residuals 107 6.76 0.06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Tính toán sai số và dự đoán
- Tính toán các tham số của mo hình
## (Intercept)
## 0.001094376## log(DBH)
## 2.140101## log(H)
## 0.5793182Ta lập được mô hình sau: \(AGB = 0.057356×DBH^{2.140101}×H^{0.579318}\) với các tham số và \(R^2\) tồn tại rõ rệt với Pr < 0.0001
Tính các chỉ tiêu thống kê của mô hình đã lập ở trên
Mô hình trở lại chuyển đổi: \(Y = exp(a)*X^b*CF\)
Tính hệ số hiệu chỉnh: \(CF = exp(\frac{RSE^2}{2})\)
## [1] 0.0631657## [1] 1.032087- Tính sai số và dự đoán
- Các chỉ tiêu thống kê
## [1] 0.9420377- R-squared (R hiệu điều chỉnh) là 0.9420377, mô hình phù hợp, có nghĩa là mô hình giải thích được khoảng 94.2% phương sai của biến AGB trong dữ liệu. Và tồn tại với Pr = 2.2e-16
## [1] 0.9398296- R-squared adjusted là 0.9398296, gần giá trị \(R^2\), cho thấy mô hình vẫn duy trì mức độ phù hợp tốt sau khi điều chỉnh cho số lượng biến độc lập trong mô hình.
Tính \(AIC = -2L + 2*p\)
trong đó p là số tham số của mô hình
ML.backtr1 <- -1/2*sum(ggg$backtr1.res^2/var(ggg$backtr1.res)+log(2*pi)
+log(var(ggg$backtr1.res)))
AIC.backtr1 <- -2*ML.backtr1 +2*3
AIC.backtr1## [1] 2116.764## [1] 180.1414## [1] 3552.594## [1] 20.93618\(R^2adj = 0.939829\); CF = 1.032087; AIC = 2116.764
Các sai số: Bias = 180.1414 kg; RMSE = 3552.594 kg và MAPE = 20.93%
Vẽ các đồ thị của mô hình: \(AGB = 0.057356×DBH^{2.140101}×H^{0.579318}\)
- Giá trị quan sát và dự đoán
p1 <- ggplot(ggg, aes(x=ggg$backtr1.fit , y=AGB))
p1 <- p1 + geom_point(cex=2)
p1 <- p1 + geom_abline(cex = 1.5, intercept = 0, slope = 1, col="black")
p1 <- p1 + xlab("Dự đoán AGB (kg)") + ylab("Quan sát AGB (kg)") + theme_bw()+
theme_bw()
p1 = p1 + theme(axis.title.y = element_text(size = rel(1.5)))
p1 = p1 + theme(axis.title.x = element_text(size = rel(1.5)))
p1 <- p1 + theme(plot.title = element_text(size = rel(1.7)))
p1 = p1 + theme(axis.text.x = element_text(size=15))
p1 = p1 + theme(axis.text.y = element_text(size=15))
p1## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.
- Kết quả và giá trị dự đoán
p2 <- ggplot(ggg, aes(x=ggg$backtr1.fit, y=ggg$backtr1.res ))
p2 <- p2 + geom_point(cex = 2)
p2 <- p2 + geom_line(cex = 1.5, aes(x=ggg$backtr1.fit, y=0))
p2 <- p2 + xlab("Dự đoán AGB (kg)") + ylab("Sai số (kg)") + theme_bw()+ theme_bw()
p2 = p2 + theme(axis.title.y = element_text(size = rel(1.5)))
p2 = p2 + theme(axis.title.x = element_text(size = rel(1.5)))
p2 <- p2 + theme(plot.title = element_text(size = rel(1.7)))
p2 = p2 + theme(axis.text.x = element_text(size=15))
p2 = p2 + theme(axis.text.y = element_text(size=15))
p2## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.## Warning: Use of `ggg$backtr1.res` is discouraged.
## ℹ Use `backtr1.res` instead.## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.
## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.
## Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.## Warning: Use of `ggg$backtr1.res` is discouraged.
## ℹ Use `backtr1.res` instead.## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.
Nhận xét:
Kết quả trên đồ thị cho thấy, khi tăng thêm một biến số H, mô hình ước tính AGB bám sát hơn giá trị quan sát trên đường chéo (trái) và biến động sai số có giảm, tuy nhiên, cũng rất phân tán khi AGB lớn vì mô hình thiếu trọng số.
Tiếp tục thử nghiệm tiến hành ước lượng mô hình phi tuyến power mở rộng gồm bốn biến DBH, H, WD và CA, trong đó ba biến DBH, H và WD tạo thành một tổ hợp biến đại diện cho \(AGB (DBH^2HWD = (DBH/100)^2×H×WD×1000 kg)\) và được logarit hóa trong chương trình “lm” của R theo phương pháp bình phương tối thiểu:
\(AGB = bo×(DBH^2HWD)^{b1}×CA^{b2}×ԑ\), logarit hóa: \(log(AGB) = log(b_o) +b_1×log(DBH^2HWD) +b_2×log(CA) + log(ԑ)\)
Thiết lập mô hình logrit hóa mô hình đa biến/tổ hợp biến \(log(AGB) = log(b_o) +b_1×log(DBH^2HWD) +b_2×log(CA)\)
- Tổ hợp biến: \(DBH^2HWD\) xấp xỉ AGB (Tạo tổ hợp biến)
- Lập mô hình
##
## Call:
## lm(formula = log(AGB) ~ log(DBH2HWD) + log(CA), data = ggg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.94405 -0.09722 0.02236 0.14768 0.49252
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11.10416 0.20230 -54.891 < 2e-16 ***
## log(DBH2HWD) 0.87591 0.01977 44.295 < 2e-16 ***
## log(CA) 0.17360 0.03880 4.474 1.92e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2194 on 107 degrees of freedom
## Multiple R-squared: 0.9886, Adjusted R-squared: 0.9884
## F-statistic: 4625 on 2 and 107 DF, p-value: < 2.2e-16## Analysis of Variance Table
##
## Response: log(AGB)
## Df Sum Sq Mean Sq F value Pr(>F)
## log(DBH2HWD) 1 444.52 444.52 9230.842 < 2.2e-16 ***
## log(CA) 1 0.96 0.96 20.021 1.919e-05 ***
## Residuals 107 5.15 0.05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Đầu ra của mô hình (Tính dự đoán và sai số mô hình)
- Tính toán các hệ số mô hình ( Tính toán các tham số mô hình)
## (Intercept)
## 1.504965e-05## log(DBH2HWD)
## 0.8759102## log(CA)
## 0.1735963-Mô hình ước lượng: \(AGB = 0.603378×(DBH^2HWD)^{0.875814}×CA^{0.173520}\) với các tham số và \(R^2\) tồn tại rõ rệt với Pr < 0.0001
Tính toán các chỉ tiêu thống kê của mô hình: \(AGB = 0.603378×(DBH^2HWD)^{0.875814}×CA^{0.173520}\)
Mô hình được chuyển đổi trở lại: \(Y = exp(a)*X^b*CF\)
Tính hệ số hiệu chỉnh: \(CF = exp(RSE^2/2)\):
## [1] 0.04815569## [1] 1.02437- Tính giá trị dự báo và sai số
ggg$backtr1.fit <- lmt1.CF * a * ggg$DBH2HWD^b*ggg$CA^c
ggg$backtr1.res <- ggg$AGB - ggg$backtr1.fit- Tính các chỉ tiêu thống kê, sai số
## [1] 0.9617249## [1] 0.9602668- Tính AIC = -2L + 2*p trong đó p là số tham số của mô hình
ML.backtr1 <- -1/2*sum(ggg$backtr1.res^2/var(ggg$backtr1.res)+log(2*pi)
+log(var(ggg$backtr1.res)))
AIC.backtr1 <- -2*ML.backtr1 +2*3
AIC.backtr1## [1] 2071.117Bias <- mean(ggg$backtr1.res)
RMSE <- sqrt(mean((ggg$backtr1.res)^2))
MAPE <- 100*mean(abs(ggg$backtr1.res)/ggg$AGB)
Bias## [1] 30.58959## [1] 2886.896## [1] 17.27431Kết quả mô hình: \(AGB = 0.603378×(DBH^2HWD)^{0.875814}×CA^{0.173520}\) có các chỉ tiêu thống kê, sai số như sau : CF = 1.02437; \(R^2adj\) = 0.960269; AIC = 2071.117 Các sai số: Bias = 30.58959 kg; RMSE = 2886.896 kg và MAPE = 17.27%
Vẽ đồ thị của mô hình \(AGB = 0.603378×(DBH^2HWD)^{0.875814}×CA^{0.173520}\)
Các giá trị được quan sát và dự đoán:
p1 <- ggplot(ggg, aes(x=ggg$backtr1.fit , y=AGB))
p1 <- p1 + geom_point(cex=2)
p1 <- p1 + geom_abline(cex = 1.5, intercept = 0, slope = 1, col="black")
p1 <- p1 + xlab("Dự đoán AGB (kg)") + ylab("Quan sát AGB (kg)") + theme_bw()+
theme_bw()
p1 = p1 + theme(axis.title.y = element_text(size = rel(1.5)))
p1 = p1 + theme(axis.title.x = element_text(size = rel(1.5)))
p1 <- p1 + theme(plot.title = element_text(size = rel(1.7)))
p1 = p1 + theme(axis.text.x = element_text(size=15))
p1 = p1 + theme(axis.text.y = element_text(size=15))
p1 ## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.
- Kết quả và giá trị dự đoán
p2 <- ggplot(ggg, aes(x=ggg$backtr1.fit, y=ggg$backtr1.res ))
p2 <- p2 + geom_point(cex = 2)
p2 <- p2 + geom_line(cex = 1.5, aes(x=ggg$backtr1.fit, y=0))
p2 <- p2 + xlab("Dự đoán AGB (kg)") + ylab("Sai số (kg)") + theme_bw()+ theme_bw()
p2 = p2 + theme(axis.title.y = element_text(size = rel(1.5)))
p2 = p2 + theme(axis.title.x = element_text(size = rel(1.5)))
p2 <- p2 + theme(plot.title = element_text(size = rel(1.7)))
141## [1] 141p2 = p2 + theme(axis.text.x = element_text(size=15))
p2 = p2 + theme(axis.text.y = element_text(size=15))
p2## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.## Warning: Use of `ggg$backtr1.res` is discouraged.
## ℹ Use `backtr1.res` instead.## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.
## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.
## Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.## Warning: Use of `ggg$backtr1.res` is discouraged.
## ℹ Use `backtr1.res` instead.## Warning: Use of `ggg$backtr1.fit` is discouraged.
## ℹ Use `backtr1.fit` instead.
- Từ đồ thị trên cho thấy mô hình ước tính AGB với bốn biến số DBH, H, WD và CA có giá trị dự đoán bám sát quan sát; tuy nhiên, với mô hình không có trọng số thì sai số vẫn biến động mạnh khi AGB tăng lên. Vì vậy, cần sử dụng mô hình có trọng số.
Từ các kết quả trên cho thấy, khi gia tăng số biến số đầu vào từ một biến DBH lên hai biến (thêm H) và bốn biến (DBH, H, WD, CA) thì mô hình có độ tin cậy càng cao và sai số giảm. Mô hình tốt nhất có bốn biến số đầu vào ở dạng tổ hợp biến với \(R^2adj\) cao nhất, CF bé nhất và gần bằng 1, AIC bé nhất, các sai số Bias, RMSE, MAPE bé nhất. Sử dụng mô hình 4 biến số giảm sai số tương đối MAPE 5% so với mô hình một biến số.
4.2.3 Ước lượng mô hình theo phương pháp phi tuyến tính
Ngoài việc ước lượng các mô hình phi tuyến tính bằng cách tuyến tính hóa mô hình theo hàm logarit và áp dụng phương pháp bình phương tối thiểu, thì có hàng loạt các phương pháp ước lượng trực tiếp mô hình phi tuyến tính mà không phải tuyến tính hóa như là phương pháp phi tuyến Marquardt (StatPoint-Inc., 2005), phi tuyến bình phương tối thiểu (chương trình “nls” trong R: Nonlinear Least Squares) (Bates và Watts, 1988), phi tuyến ảnh hưởng phức hợp hợp lý tối đa (chương trình “nlme” trong R: Nonlinear Mixed-Effects - Maximum Likelihood) (Davidian và Giltinan, 1995; Pinheiro et al., 2014). Các phương pháp phi tuyến ước lượng trực tiếp các tham số mô hình mà không phải thông qua hàm trung gian để tuyến tính hóa, do đó, về mặt nguyên tắc chung nó đã làm giảm được sai số ước lượng trung gian, cụ thể là không cần sử dụng hệ số điều chỉnh CF như mô hình tuyến tính hóa, do đó, tăng độ tin cậy của mô hình; ngoài ra còn có kỹ thuật trọng số để giảm biến động sai số của mô hình khi ước lượng bằng phương pháp phi tuyến tính.
Sau đây lần lượt giới thiệu các phương pháp ước lượng các hàm phi tuyến khác nhau và cách thức đánh giá để lựa chọn phương pháp thích hợp trong từng trường hợp cụ thể.
4.2.3.1 Phương pháp phi tuyến tính Marquardt
Phương pháp phi tuyến Marquardt (StatPoint-Inc., 2005) ước lượng trực tiếp các tham số của mô hình phi tuyến với thông tin đầu vào là các tham số khởi đầu của mô hình.
4.2.3.2 Phương pháp phi tuyến bình phương tối thiểu (Code nls trong R)
Phương pháp phi tuyến tính bình phương tối thiểu (nls - Nonlinear Least Squares) ước lượng trực tiếp các tham số của mô hình phi tuyến với thông tin đầu vào là các tham số khởi đầu của mô hình. Phương pháp này ñược áp dụng thuận lợi trong phần mềm mã nguồn mở R theo chương trình nls (Bates et al., 1988).
Sử dụng Dữ liệu 12 của Bảo Huy và Đào Công Khanh (2008) trong lập biểu sản lượng rừng trồng trám trắng (Canarium album Raeusch) tại các tỉnh Lạng Sơn, Bắc Giang, Quảng Ninh. Từ dữ liệu 73 lâm phần trám trắng được thu thập, minh họa phương pháp để thiết lập mô hình ước tính thể tích bình quân (V, m3) theo đường kính (Dg, cm) và chiều cao bình quân lâm phần (Hg, m): V = f(Dg, Hg) theo mô hình power tổ hợp biến:
\(V = a×(Dg^2Hg)^b + ԑ\)
Trong đó \(Dg^2Hg\) là tổ hợp biến đại diện cho thể tích cây (m3) \(Dg^2Hg = (\frac{Dg}{100})^2×Hg\)
Tiến hành thiết lập mô hình phi tuyến tính theo codes nls như sau:
- Gọi các package cần thiết:
- Đọc dữ liệu từ file excel
- Dùng lệnh view để xem bảng dữ liệu vừa được đọc
Lập mô hình \(V = a*(Dg^2*Hg)^b\) theo code nls
- Tạo tổ hợp biến đại diện thể tích
- Tham số khởi đầu của mô hình
- Áp dụng chương trình lập mô hình phi tuyến nls
- Tóm tắt kết quả mô hình
##
## Formula: V ~ a * Dg2Hg^b
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 0.33842 0.01051 32.19 <2e-16 ***
## b 0.88971 0.01145 77.73 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.001051 on 71 degrees of freedom
##
## Number of iterations to convergence: 5
## Achieved convergence tolerance: 9.393e-07- Tính giá trị dự đoán và sai số của mô hình
ggg1$nls_least_square.fit <- fitted.values(nls_least_square)
ggg1$nls_least_square.res <- residuals(nls_least_square) Mô hình được thiết lập: \(V = 0.33842×(Dg^2Hg)^{0.88971}\)
Các tham số của mô hình có các Pr < 2e-16, chứng tỏ các tham số tồn tại rõ rệt.
Tính các chỉ tiêu thống kê, và sai số của mô hình phi tuyến lập theo codes nls:
- Tính toán AIC, R2 và lỗi
## [1] -790.0709## [1] 0.9935933## [1] 0.9933147## [1] 6.956402e-05## [1] 0.001036942## [1] 7.850645AIC = -790.0709; \(R^2adj\) = 0.9933147 (Quan hệ là chặt chẽ)
Bias = 6.95e-05\(m^3\); RMSE = 0.001036\(m^3\) và MAPE = 7.85%. Các sai số bé cho thấy dạng mô hình tổ hợp biến mô phỏng tốt quan hệ V với Dg và Hg của loài trám trắng.
4.2.4 Mô hình phi tuyến có hay không có trọng số (Weight)
Biến số trọng số được tính: Weight = 1/X^a trong đó X là biến số độc lập ảnh hưởng rõ rệt nhất đến biến phụ thuộc và khi nó tăng lên thì biến phụ thuộc có sự phân hóa cao; a thường biến động từ -20 đến +20; thay đổi a để mô hình có sai số phân bố đều theo các giá trị dự đoán trên đồ thị.
4.2.4.1 So sánh mô hình theo phương pháp phi tuyến bình phương tối thiểu (nls) có hay không có trọng số
Mô hình phi tuyến tính có hay không có trọng số có thể được thiết lập theo chương trình phi tuyến bình phương tối thiểu nls (nonlinear least square) (Bate et al., 1988).
Dưới đây minh họa cách thiết lập các mô hình ước tính sinh khối cây rừng trên mặt đất (AGB, kg) theo các biến số DBH, H, WD được tổ hợp biến là DBH2HWD dạng hàm power (Huy et al.2016b) trên cơ sở dữ liệu 11 của 110 cây mẫu ở vùng Nam Trung Bộ, so sánh mô hình có hay không có trọng số thực hiện theo code nls.
\(AGB = a× (DBH^2HWD)^b+ ԑ\)
Lập mô hình phi tuyến không có trọng số
- Tổ hợp biến đại diện sinh khối AGB
start <- coefficients(lm(log(AGB)~log(DBH2HWD), data=ggg))
names(start) <- c("a","b")
start[1] = exp(start[1])
nls_least_square <- nls(AGB~a*DBH2HWD^b, data=ggg,
start=start)- Tóm tắt kết quả mô hình
##
## Formula: AGB ~ a * DBH2HWD^b
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 6.946e-06 4.718e-06 1.472 0.144
## b 9.778e-01 2.943e-02 33.224 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3148 on 108 degrees of freedom
##
## Number of iterations to convergence: 5
## Achieved convergence tolerance: 1.735e-074.2.5 Phương pháp phi tuyến ảnh hưởng tổng hợp (Nonlinear MixedEffects - nlme) Maximum Likelihood có trọng số để ước lượng mô hình phi tuyến tính
Ngoài phương pháp phi tuyến bình phương tối thiểu có trọng số sử dụng codes nls trong phần mềm mã nguồn mở R, còn có phương pháp phi tuyến ảnh hưởng phức hợp hợp lý tối đa (nlme: nonlinear mixed effects - Maximum Likelihood) có trọng số (weights), sau đây viết tắt là phương pháp phi tuyến Maximum Likelihood có trọng số (Davidian et al., 1995; Pinheiro et al., 2014).
Trong một số quan hệ phức tạp như, mô hình bị ảnh hưởng của nhiều nhân tố khác ngoài các biến số độc lập, thì mô hình theo phương pháp nlme Maximum Likelohood có trọng số tỏ ra có hiệu quả để nâng cao độ tin cậy, vì vậy, cần có thử nghiệm để áp dụng phương pháp này so với nls thông thường và hay áp dụng.
Để minh họa áp dụng phương pháp phi tuyến Maximum Likelihood có trọng số, sử dụng mô hình power. Phần mềm mã nguồn mở R được áp dụng theo chương trình nlme và chẩn đoán qua sơ đồ sử dụng ggplot2.
Phân tích ban đầu cho thấy, biến động của sai số có xu hướng gia tăng khi gia tăng \(X_i\) trong các mô hình. Vì vậy, một hàm phương sai theo trọng số đã được áp dụng để điều chỉnh các tham số của mô hình nhằm giảm biến động sai số này.
4.3 Chỉ Số Furnival’s Index để lựa chọn dạng phương trình khác nhau hoặc phương pháp ước lượng mô hình phi tuyến: Tuyến tính hóa hay phi tuyến Maximum Likelihood
Các mô hình có thể khác nhau về biến phụ thuộc y, một mô hình không ñổi biến số vẫn là y, mô hình khác thì đổi biến ví dụ như là log(y), 1/y, sqrt(y). Khi các mô hình khác nhau về đổi biến số thì cũng áp dụng phương pháp ước lượng khác nhau. Phổ biến là mô hình phi tuyến có thể ước lượng các tham số của mô hình theo một trong hai phương pháp chính là, tuyến tính hóa bình phương tối thiểu (sử dụng codes lm trong R) hoặc phi tuyến bao gồm: phi tuyến bình phương tối thiểu (sử dụng codes nls trong R) hoặc phi tuyến Marquardt hoặc phi tuyến Maximum Likelihood (code nlme trong R).
Để so sánh các mô hình khác nhau về biến y (y và log(y), 1.y,…) được ước lượng theo phương pháp khác nhau cơ bản như là, tuyến tính hoặc phi tuyến tính; lúc này cần sử dụng chỉ số Furnival’s Index (FI) (1961) (Jayaraman, 1999)). Chỉ số Furnival’s Index dùng để so sánh các mô hình không giống nhau về biến số phụ thuộc (ví dụ y và ln(y)), vì lúc này các chỉ tiêu thống kê như R2, AIC, các loại sai số sẽ không thể dùng để so sánh do khác nhau về giá trị của biến phụ thuộc. Các mô hình theo phương pháp ước lượng áp dụng có chỉ số Furnival Index (FI) thấp hơn là tốt hơn.
Công thức tính Furnival’s Index (FI) như sau:
\[FI = \frac{(RMSE * 1)}{Geometric mean (y')}\]
Trong đó: RMSE: Root Mean Squared Error: Sai số trung phương; y’ là đạo hàm bậc nhất của biến phụ thuộc y và bằng 1, nếu là biến phụ thuộc được đổi biến số là ln(y) thì sẽ bằng 1/y.
Công thức tính trung bình hình học (Geometric mean) (Huy et al., 2016b):
\(Geometric Mean = (\prod\limits_{i=1}^{N}x_i)^{1/n} = \sqrt[n]{x_1 x x_2 x ... x_n}\)
Tính FI cho từng dạng mô hình và phương pháp lập mô hình có thể được thực hiện trong R như sau: Để minh họa, sử dụng codes nêu trên để tính FI cho các mô hình ước tính Earthquake theo một biến số s hoặc tổ hợp biến số s + q theo hai dạng mô hình tuyến tính hóa logarit và power ứng với hai phương pháp tuyến tính bình phương tối thiểu (lm trong R) và phi tuyến Maximum likelihood có trọng số (nlme trong R)
5 CHƯƠNG 5: KẾT LUẬN
5.1 Ý nghĩa
Gói “nlme” (non-linear mixed-effects model) là một công cụ mạnh mẽ được sử dụng để xây dựng và phân tích mô hình tuyến tính. Công cụ này hỗ trợ nhiều bài toán phân tích dữ liệu thực tế, đặc biệt là trong lĩnh vực nghiên cứu y học, sinh học, kinh tế học và nhiều lĩnh vực khác. Trong nghiên cứu y học và thử nghiệm lâm sàng, gói “nlme” được sử dụng để mô hình hóa sự biến đổi giữa các cá thể và giữa các thời điểm trong các nhóm điều trị khác nhau. Phân tích dữ liệu lặp lại: Trong các nghiên cứu theo dõi dài hạn, dữ liệu có thể được thu thập từ cùng một cá thể theo thời gian hoặc từ cùng một vùng địa lý. Gói “nlme” hỗ trợ xử lý dữ liệu lặp lại này bằng cách giúp xác định và mô hình hoá các hiệu ứng kết hợp từ dữ liệu lặp lại này. Đánh giá hiệu quả điều trị: Gói “nlme” cung cấp các công cụ để đánh giá hiệu quả của các phương pháp điều trị, thuốc hoặc liệu pháp trong các thử nghiệm lâm sàng.
5.2 Ưu và Nhược điểm của package nlme
5.2.1 Ưu điểm của nlme
Package nlme trong R có một số ưu điểm vượt trội khi sử dụng để mô hình hóa và phân tích dữ liệu có cấu trúc theo thời gian:
Mô hình hỗn hợp tuyến tính không có giả định về phân phối: Package nlme cho phép xây dựng mô hình tuyến tính hỗn hợp không yêu cầu các giả định về phân phối của dữ liệu. Điều này rất hữu ích khi làm việc với dữ liệu có tính chất phi-Gaussian.
Xử lý dữ liệu lặp lại: Package nlme được thiết kế để xử lý các dữ liệu có tính chất lặp lại, như dữ liệu thu thập từ nhiều quan sát trên cùng một đối tượng hoặc dữ liệu theo dõi theo thời gian. Nó cung cấp các công cụ mạnh mẽ để mô hình hóa và phân tích dữ liệu lặp lại.
Mô hình phân tích hiệu quả: Package nlme cung cấp các công cụ phân tích hiệu quả để xác định các thông số mô hình, đánh giá độ tin cậy của ước lượng và kiểm tra giả thiết. Nó cung cấp các phương pháp ước lượng tham số hiệu quả và các chỉ số đánh giá mô hình như AIC (Akaike’s Information Criterion) và BIC (Bayesian Information Criterion).
Đa dạng mô hình: Package nlme cung cấp các công cụ để xây dựng nhiều loại mô hình phức tạp, bao gồm mô hình tuyến tính, phi tuyến tính, logistic và các mô hình dự báo. Điều này cho phép người dùng tìm hiểu và mô phỏng nhiều khía cạnh của dữ liệu.
Tích hợp sử dụng với các packages khác: Package nlme có thể được kết hợp với các package khác trong R như ggplot2, dplyr và tidyr để trực quan hóa dữ liệu, xử lý dữ liệu và thực hiện các phân tích thêm.
5.2.2 Nhược điểm của nlme
Các câu lệnh trong package nlme khá là khó học đòi người dùng phải sử dụng các kiến thức về thống kê và mô hình tuyến tính đa cấp. Điều này làm cho mấy người mới bắt đầu học cảm thấy khó khăn trong quá trình làm bài
Thời gian chạy chậm bởi vì số lượng quan sát lớn hoặc phải thực hiện các phân tích phức tạp, việc chạy mô hình tuyến tính đa cấp trong nlme có thể mất thời gian rất lâu, giảm hiệu suất và làm cho việc phân tích dữ liệu trở nên chậm chạp.
Gói package nlme không cung cấp nhiều tính năng linh hoạt để mô hình hóa các mối quan hệ phi tuyến tính hoặc mô hình hồi quy phi tuyến, mô hình này chủ yếu trong phạm vi tuyến tính và đa cấp.
Thiếu dữ liệu package này gặp nhiều khó khăn không thể xử lí được việc đó
Hạn chế trong việc trực quan hóa: nlme cung cấp một số tính năng trực quan hóa kết quả, nhưng không phải là một công cụ trực quan hóa chuyên sâu. Điều này có thể làm cho việc trực quan hóa và hiển thị kết quả mô hình trở nên khó khăn.
Phức tạp trong việc sử dụng: Package nlme được thiết kế cho việc xử lý và phân tích dữ liệu lặp lại, nhưng cách sử dụng của nó cũng đòi hỏi người dùng có kiến thức và kinh nghiệm về thống kê cao. Điều này có thể làm cho nó khó khăn và phức tạp đối với những người mới bắt đầu hoặc người không có nền tảng vững chắc về các phương pháp phân tích dữ liệu lặp lại.
Khó khăn trong việc tìm hiểu: Vì package nlme có rất nhiều chức năng và cú pháp phức tạp, việc tìm hiểu và làm quen với nó có thể đòi hỏi thời gian và nỗ lực. Điều này có thể làm cho việc sử dụng package trở nên khó khăn và mất thời gian đối với người mới bắt đầu hoặc người không có kinh nghiệm trong việc sử dụng R.
Hạn chế trong việc xử lý dữ liệu phi-Gaussian: Mặc dù package nlme được thiết kế để xử lý dữ liệu không yêu cầu các giả định phân phối, nhưng nó không hỗ trợ tốt cho việc xử lý dữ liệu phi-Gaussian. Điều này có thể khiến nó hạn chế trong việc mô hình hóa và phân tích dữ liệu phi-Gaussian.
6 CHƯƠNG 6: Tài liệu tham thảo
https://www.ttn.edu.vn/attachments/article/2875/Luan%20an%20NTT.pdf
https://rpubs.com/lengockhanhi/396521
https://jchiquet.github.io/MAP566/docs/mixed-models/map566-lecture-nonlinear-mixed-model.html
https://r-universe.dev/manuals/nlme.html#residuals.nlmeStruct
https://quantdev.ssri.psu.edu/sites/qdev/files/GCM_Chp3_Tutorial_2.html
https://asancpt.github.io/nlme/chapter-1.html
https://epirhandbook.com/vn/importing.html
https://www.ttn.edu.vn/attachments/article/2875/Luan%20an%20NTT.pdf
https://rstudio-pubs-static.s3.amazonaws.com/168431_ba8dd35e0de3468aaff301b9f08d2862.html#22_readxl