Phương pháp spline khối cho dữ liệu chuỗi dọc?

Question

Tôi có một dữ liệu nối tiếp được định dạng như sau:

time milk Animal_ID 
30  25.6 1 
31  27.2 1 
32  24.4 1 
33  17.4 1 
34  33.6 1 
35  25.4 1 
33  29.4 2 
34  25.4 2 
35  24.7 2 
36  27.4 2 
37  22.4 2 
80  24.6 3 
81  24.5 3 
82  23.5 3 
83  25.5 3 
84  24.4 3 
85  23.4 3 
. . . 

Nói chung, 300 loài động vật có hồ sơ của sữa trong thời điểm khác nhau của thời gian ngắn. Tuy nhiên, nếu chúng ta kết hợp dữ liệu của chúng với nhau và không quan tâm đến animal_ID khác nhau, chúng ta sẽ có một đường cong giữa sữa ~ thời gian như thế này, dòng bên dưới: Ngoài ra, trong hình trên, chúng ta có dữ liệu cho 1 ví dụ động vật, chúng ngắn và rất biến đổi. Mục đích của tôi là để làm mịn từng dữ liệu động vật nhưng sẽ là nếu mô hình cho phép học tập tổng quát từ toàn bộ dữ liệu được đưa vào. Tôi đã sử dụng mô hình mịn khác nhau (ns, bs, smooth.spline) với định dạng sau nhưng nó chỉ không làm việc:

mod <- lme(milk ~ bs(time, df=3), data=dat, random = ~1|Animal_ID) 

tôi hy vọng nếu ai đó đã giải quyết vấn đề này sẽ cung cấp cho tôi một lời khuyên. Cảm ơn Bạn có thể truy cập toàn bộ tập dữ liệu tại đây:

Answer 1

Tôi khuyên bạn nên sử dụng gói mgcv. Đây là một trong các gói R được khuyến nghị, thực hiện một lớp mô hình được gọi là mô hình hỗn hợp phụ gia tổng quát. Bạn chỉ có thể tải nó bằng cách library(mgcv). Đây là một thư viện rất mạnh, có thể xử lý từ mô hình hồi quy tuyến tính đơn giản nhất, đến các mô hình tuyến tính tổng quát, các mô hình phụ, cho các mô hình phụ tổng quát, cũng như các mô hình có hiệu ứng hỗn hợp (hiệu ứng cố định + hiệu ứng ngẫu nhiên). Bạn có thể liệt kê tất cả các chức năng (xuất khẩu) của mgcv qua

ls("package:mgcv") 

Và bạn có thể thấy rất nhiều người trong chúng.

Đối với dữ liệu và vấn đề cụ thể của bạn, bạn có thể sử dụng một mô hình với công thức:

model <- milk ~ s(time, bs = 'cr', k = 100) + s(Animal_ID, bs = 're') 

Trong mgcv, s() là một thiết lập cho các chức năng mượt mà, thể hiện bằng cơ sở spline ngụ ý bởi bs. "cr" là cơ sở spline khối, chính xác là những gì bạn muốn. k là số nút. Nó sẽ được chọn tùy thuộc vào số lượng giá trị duy nhất của biến số time trong tập dữ liệu của bạn. Nếu bạn đặt k cho chính xác con số này, bạn sẽ kết thúc với một spline làm mịn; trong khi bất kỳ giá trị nhỏ hơn có nghĩa là một spline hồi quy. Tuy nhiên, cả hai sẽ bị phạt (nếu bạn biết nghĩa phạt nào). Tôi đã đọc dữ liệu của bạn trong:

dat <- na.omit(read.csv("data.txt", header = TRUE)) ## I saved you data into file "data.txt" 
dat$Animal_ID <- factor(dat$Animal_ID) 
nrow(dat) ## 12624 observations 
length(unique(dat$time)) ## 157 unique time points 
length(ID <- levels(dat$Animal_ID)) ## 355 cows 

Có 157 giá trị duy nhất, vì vậy tôi cho rằng k = 100 có thể phù hợp.

Đối với Animal_ID (bị ép buộc làm yếu tố), chúng tôi cần thuật ngữ mô hình cho hiệu ứng ngẫu nhiên. "re" là một lớp đặc biệt cho hiệu ứng ngẫu nhiên i.i.d. Nó được chuyển đến bs đối với một số lý do xây dựng ma trận nội bộ (vì vậy đây không phải là một chức năng trơn tru!).

Bây giờ để phù hợp với mô hình GAM, bạn có thể gọi di sản gam hoặc liên tục phát triển bam (gam cho dữ liệu lớn). Tôi nghĩ rằng bạn sẽ sử dụng sau này. Họ có cùng quy ước gọi điện thoại tương tự như lm và glm.Ví dụ, bạn có thể làm:

fit <- bam(model, data = dat, family = "gaussian", discrete = TRUE, nthreads = 2) 

Như bạn thấy, bam cho phép tính toán song song đa lõi qua nthreads. Trong khi discrete là một tính năng mới được phát triển có thể tăng tốc độ hình thành ma trận.

Vì bạn đang xử lý dữ liệu chuỗi thời gian, cuối cùng bạn có thể xem xét một số tương quan tự thời gian. mgcv cho phép cấu hình tương quan AR1, có hệ số tương quan được chuyển bởi bam đối số rho. Tuy nhiên, bạn cần thêm chỉ mục AR_start để thông báo cho mgcv cách chuỗi thời gian chia nhỏ thành nhiều phần. Ví dụ: khi đạt đến một số khác nhau Animal_ID, AR_start hãy lấy số TRUE để chỉ ra một phân đoạn mới của chuỗi thời gian. Xem ?bam để biết chi tiết.

mgcv cũng cung cấp

1. summary.gam chức năng cho bản tóm tắt mô hình
2. gam.check cho mô hình cơ bản kiểm tra
3. plot.gam chức năng cho âm mưu điều khoản cá nhân
4. predict.gam (hoặc predict.bam) cho dự đoán trên dữ liệu mới.

Ví dụ, bản tóm tắt của mô hình đề nghị trên là:

> summary(fit) 

Family: gaussian 
Link function: identity 

Formula: 
milk ~ s(time, bs = "cr", k = 100) + s(Animal_ID, bs = "re") 

Parametric coefficients: 
     Estimate Std. Error t value Pr(>|t|)  
(Intercept) 26.1950  0.2704 96.89 <2e-16 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Approximate significance of smooth terms: 
       edf Ref.df  F p-value  
s(time)  10.81 13.67 5.908 1.99e-11 *** 
s(Animal_ID) 351.43 354.00 136.449 < 2e-16 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

R-sq.(adj) = 0.805 Deviance explained = 81.1% 
fREML = 29643 Scale est. = 5.5681 n = 12624 

Các edf (mức độ hiệu quả của tự do) có thể được coi như là một thước đo của mức độ phi tuyến tính. Vì vậy, chúng tôi đặt vào k = 100, trong khi kết thúc bằng edf = 10.81. Điều này cho thấy rằng spline s(time) đã bị phạt nặng. Bạn có thể xem những gì trông giống như s(time) bởi:

plot.gam(fit, page = 1) 

Lưu ý rằng tác động ngẫu nhiên s(Animal_ID) cũng có một "trơn tru", đó là một hằng số bò cụ thể. Đối với các hiệu ứng ngẫu nhiên, một âm mưu Gaussian QQ sẽ được trả về.

Các số liệu chẩn đoán được trả về bởi

invisible(gam.check(fit)) 

vẻ OK, vì vậy tôi nghĩ rằng mô hình này là chấp nhận được (tôi không cung cấp cho bạn lựa chọn mô hình, vì vậy nghĩ ra một mô hình tốt hơn nếu bạn nghĩ có).

Nếu bạn muốn chắc dự đoán cho Animal_ID = 26, bạn có thể làm

newd <- data.frame(time = 1:150, Animal_ID = 26) 
oo <- predict.gam(fit, newd, type = `link`, se.fit = TRUE) 

Lưu ý rằng

• Bạn cần phải bao gồm cả các biến trong newd (nếu không mgcv than phiền biến mất)
• từ bạn chỉ có một spline mịn s(time) và cụm từ hiệu ứng ngẫu nhiên s(Animal_ID) là hằng số trên Animal_ID. do đó, bạn có thể sử dụng type = 'link' để dự đoán riêng lẻ.Nhân tiện, type = 'terms' chậm hơn type = 'link'.

Nếu bạn muốn chắc dự đoán cho nhiều hơn một con bò, hãy thử một cái gì đó như thế này:

pred.ID <- ID[1:10] ## predict first 10 cows 
newd <- data.frame (time = rep (1:150, times = n), Animal_ID = factor (rep (pred.ID, each = 150))) 
oo <- predict.bam (fit, newd, type = "link", se.fit = TRUE) 

Lưu ý rằng

• Tôi đã sử dụng predict.bam đây, như bây giờ chúng ta có 150 * 10 = 1500 dữ liệu điểm để dự đoán. Ngoài ra: chúng tôi yêu cầu se.fit = TRUE. Điều này khá tốn kém, vì vậy hãy sử dụng predict.bam nhanh hơn predict.gam. Đặc biệt, nếu bạn đã lắp mô hình của mình bằng cách sử dụng bam(..., discrete = TRUE), bạn có thể có predict.bam(..., discrete = TRUE). quá trình dự đoán đi qua các bước tương tự hình thành ma trận như trong mô hình phù hợp (xem ?smoothCon sử dụng trong mô hình phù hợp và ?PredictMat sử dụng trong dự đoán, nếu bạn đang quan tâm để biết cấu trúc bên trong hơn mgcv.)
• tôi đã chỉ định Animal_ID như các yếu tố, bởi vì đây là một hiệu ứng ngẫu nhiên.

Để biết thêm chi tiết về mgcv, bạn có thể tham khảo hướng dẫn sử dụng thư viện. Kiểm tra đặc biệt ?mgcv, ?gam, ?bam?s.

---

cập nhật cuối cùng

Mặc dù tôi nói rằng tôi sẽ không giúp bạn với phần mô hình, nhưng tôi nghĩ rằng mô hình này là tốt hơn (nó mang lại cao hơn adj-Rsquared) và cũng là hợp lý hơn theo nghĩa:

model <- milk ~ s(time, bs = 'cr', k = 20) + s(Animal_ID, bs = 're') + s(Animal_ID, time, bs = 're') 

Thuật ngữ cuối cùng áp đặt slop ngẫu nhiên. Điều này ngụ ý rằng chúng tôi giả định rằng mỗi con bò riêng lẻ có mô hình tăng/giảm sản lượng sữa khác nhau. Đây là giả định hợp lý hơn trong vấn đề của bạn. Mô hình trước đó với chỉ chặn ngẫu nhiên là không đủ. Sau khi thêm slop ngẫu nhiên này, thuật ngữ trơn tru s(time) trông mượt mà hơn. Đây là một dấu hiệu tốt không phải là một dấu hiệu xấu, bởi vì chúng tôi muốn một số lời giải thích đơn giản cho s(time), phải không? So sánh số s(time) bạn nhận được từ cả hai mô hình và xem những gì bạn khám phá.

Tôi cũng đã giảm k = 100 thành k = 20. Như chúng ta đã thấy trong lần khớp trước, số edf cho cụm từ này là khoảng 10, vì vậy k = 20 là khá đủ.