Hồi quy tuyến tính Bootstrap trong R

Có vẻ như một trong những yêu cầu lớn là có thể vẽ các kết quả đã khởi động. Đây là một lựa chọn

Đầu tiên, làm cho dữ liệu

df <- structure(list(replicate = c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2), level = c("low", "low", "mid", "mid", "low", "low", "mid", "mid", "up", "up", "up"), high.density = c(14, 54, 82, 24, 12, 11, 12, NA, 40, NA, 20), low.density = c(36, 31, 10, 
                                                                                                                                                                                                                                        NA, 28, 45, 17, 24, 10, 5, 2)), class = c("spec_tbl_df","tbl_df","tbl", "data.frame"), row.names = c(NA, -11L), spec = structure(list(cols = list(replicate = structure(list(), class = c("collector_double", "collector")), level = structure(list(), class = c("collector_character","collector")), high.density = structure(list(), class = c("collector_double","collector")), low.density = structure(list(), class = c("collector_double", 
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     "collector"))), default = structure(list(), class = c("collector_guess", "collector")), skip = 1L), class = "col_spec"))

df$replicate <- as.factor(as.numeric(df$replicate))
df$level <- as.factor(as.character(df$level))

Tiếp theo, tạo một số dữ liệu giả định di chuyển biến quan tâm low.density trong phạm vi của nó. Ở đây chúng tôi đang làm điều đó cho tất cả các kết hợp có thể có của replicate và level, nhưng bạn cũng có thể chỉ chọn một giá trị cho mỗi

hyp <- expand.grid(replicate = as.factor(c(1,2)), 
                   level=factor(1:3, labels=c("low", "mid", "up")), 
                   low.density = seq(min(df$low.density, na.rm=TRUE), 
                                     max(df$low.density, na.rm=TRUE), 
                                     length=25))

Sau đó, chúng tôi thực hiện bootstrapping. Trong chức năng bên dưới, chúng tôi vẽ dữ liệu, ước tính mô hình và sau đó tạo dự đoán. Nếu mô hình hoặc các dự đoán không thành công, thì mô hình cụ thể đó sẽ bị loại bỏ và một mô hình khác sẽ được rút cho đến khi bạn nhận được 1000 lượt rút hợp lệ

res <- NULL
i <- 1
while(i <= 1000){
  tmp <- df_shuffle(df)
  mod <- try(lm(high.density ~ low.density + replicate + level, data=tmp))
  if(!inherits(mod, "try-error")){
    pred <- try(predict(mod, newdata=hyp))
    if(!inherits(pred, "try-error")){
      res <- cbind(res, pred )
      i <- i+1
    }
  }
}

Ước tính mô hình ban đầu để có được các giá trị dự đoán

orig <- lm(high.density ~ low.density + level + replicate, data=df)
hyp$fit <- predict(orig, newdata=hyp)

Tính khoảng tin cậy lượng tử cho từng dự đoán được khởi động và thêm chúng vào tập dữ liệu

Nhớ lại rằng trong hồi quy tuyến tính, chúng ta khớp một mô hình có dạng sau $\( \begin{aligned} f_\hat{\theta} (x) = \hat{\theta_0} + \hat{\theta_1} x_1 + \ldots

Chúng tôi muốn suy ra các hệ số thực của mô hình. Vì \(\hat{\theta_0}\), \(\hat{\theta_1}\), \(\ldots\) \(\hat{\theta_p}\) là các công cụ ước tính khác nhau dựa trên dữ liệu đào tạo của chúng tôi/ . Bootstrapping là một cách tiếp cận phi tham số đối với suy luận thống kê cung cấp cho chúng tôi các lỗi tiêu chuẩn và khoảng tin cậy cho các tham số của chúng tôi

Hãy cùng xem một ví dụ về cách chúng tôi sử dụng các phương pháp bootstrapping trong hồi quy tuyến tính

18. 3. 1. Dữ liệu

Otis Dudley Duncan là một nhà xã hội học định lượng quan tâm đến việc đo lường mức độ uy tín của các nghề nghiệp khác nhau. Chỉ có 90 nghề nghiệp được đánh giá về mức độ uy tín trong cuộc khảo sát của Trung tâm Nghiên cứu Ý kiến ​​Quốc gia (NORC) năm 1947. Duncan muốn “điền” điểm uy tín cho những nghề nghiệp chưa được xếp hạng bằng cách sử dụng dữ liệu về thu nhập và giáo dục về từng nghề nghiệp được ghi lại trong cuộc điều tra dân số năm 1950. Khi kết hợp dữ liệu NORC với dữ liệu điều tra dân số năm 1950, chỉ có 45 nghề nghiệp có thể khớp. Cuối cùng, mục tiêu của Duncan là tạo ra một mô hình để giải thích uy tín bằng cách sử dụng các đặc điểm khác nhau;

Bộ dữ liệu Duncan là một mẫu ngẫu nhiên chứa thông tin về uy tín và các đặc điểm khác của 45 U. S. nghề nghiệp năm 1950. Các biến là

________số 8_______

occupationincomeeducationprestige0accountant6286821pilot7276832architect759290. 42người gác cổng720843cảnh sát34474144người phục vụ83210

45 hàng × 4 cột

Thông thường, bạn nên khám phá dữ liệu thông qua trực quan hóa để hiểu được mối quan hệ giữa các biến của mình. Dưới đây, chúng ta sẽ hình dung mối tương quan giữa thu nhập, học vấn và uy tín

plt.scatter(x=duncan["education"], y=duncan["prestige"])

plt.scatter(x=duncan["income"], y=duncan["prestige"])

plt.scatter(x=duncan["income"], y=duncan["education"])

Từ các biểu đồ trên, chúng ta thấy rằng cả giáo dục và thu nhập đều có mối tương quan thuận với uy tín; . Hãy điều chỉnh một mô hình tuyến tính bằng cách sử dụng các biến giải thích này để giải thích uy tín

18. 3. 2. lắp mô hình

Chúng ta sẽ điều chỉnh mô hình sau, giải thích uy tín của một nghề nghiệp như là một hàm tuyến tính của thu nhập và giáo dục

\[ \begin{aligned} \texttt{prestige}_i = \theta_0^* + \theta_\texttt{ income}^*\cdot\texttt{ income}_i + \theta_\texttt{education}^*\cdot\texttt

Để phù hợp với mô hình này, chúng tôi sẽ xác định ma trận thiết kế (X) và biến phản hồi của chúng tôi (y)

X = duncan.loc[:, ["income", "education"]]
X.head()

thu nhậpgiáo dục0628617276275923559046486

y = duncan.loc[:, "prestige"]
y.head()

0    82
1    83
2    90
3    76
4    90
Name: prestige, dtype: int64

Dưới đây, chúng tôi khớp với mô hình tuyến tính của mình và in tất cả các hệ số \(\widehat \theta\) của mô hình (từ phương trình trên) sau khi mô hình đã khớp với dữ liệu. Lưu ý rằng đây là các hệ số mẫu của chúng tôi

plt.scatter(x=duncan["education"], y=duncan["prestige"])

plt.scatter(x=duncan["education"], y=duncan["prestige"])

Các hệ số trên cho chúng ta ước tính của các hệ số thực. Nhưng nếu dữ liệu mẫu của chúng tôi khác, chúng tôi sẽ điều chỉnh mô hình của mình cho phù hợp với dữ liệu khác, khiến các hệ số này khác nhau. Chúng tôi muốn khám phá những gì các hệ số của chúng tôi có thể đã sử dụng các phương pháp bootstrapping

Trong các phương pháp và phân tích bootstrapping của chúng tôi, chúng tôi sẽ tập trung vào hệ số giáo dục. Chúng tôi muốn khám phá mối quan hệ một phần giữa uy tín và giáo dục khi thu nhập không đổi (chứ không phải là mối quan hệ cận biên giữa uy tín và giáo dục bỏ qua thu nhập). Hệ số hồi quy từng phần \(\widehat\theta_\texttt{education}\) minh họa mối quan hệ từng phần giữa uy tín và học vấn trong dữ liệu của chúng tôi

18. 3. 3. Khởi động các quan sát

Trong phương pháp này, chúng tôi coi các cặp \((X_i, y_i)\) là mẫu của mình, vì vậy chúng tôi xây dựng mẫu khởi động lại bằng cách lấy mẫu có thay thế từ các cặp này

\[ \begin{aligned} (X_i^*, y_i^*) = (X_I, y_I), \text{ trong đó } I=1,\dots,n \text{ được lấy mẫu thống nhất một cách ngẫu nhiên. } \end{aligned} \]

Nói cách khác, chúng tôi lấy mẫu n quan sát với sự thay thế từ các điểm dữ liệu của chúng tôi; . Sau đó, chúng tôi sẽ điều chỉnh mô hình hồi quy tuyến tính mới cho dữ liệu được lấy mẫu này và ghi lại hệ số giáo dục \(\tilde\theta_\texttt{education}\);

plt.scatter(x=duncan["education"], y=duncan["prestige"])

plt.scatter(x=duncan["education"], y=duncan["prestige"])

plt.scatter(x=duncan["education"], y=duncan["prestige"])

Lưu ý cách phân phối lấy mẫu ở trên hơi lệch về bên trái

18. 3. 3. 1. Ước tính các hệ số thực

Mặc dù chúng ta không thể đo lường trực tiếp \(\theta^*_\texttt{education}\) nhưng chúng ta có thể sử dụng khoảng tin cậy bootstrap để tính đến sự thay đổi trong hệ số hồi quy mẫu \(\widehat\theta_{\texttt{education}}\). Bên dưới, Chúng tôi xây dựng khoảng tin cậy xấp xỉ 95% cho hệ số thực \(\theta^*_\texttt{education}\), sử dụng phương pháp phân vị bootstrap. Khoảng tin cậy kéo dài từ 2. phần trăm thứ 5 đến 97. Phân vị thứ 5 trong số 10.000 hệ số khởi động

plt.scatter(x=duncan["education"], y=duncan["prestige"])

plt.scatter(x=duncan["education"], y=duncan["prestige"])

Từ khoảng tin cậy ở trên, chúng ta khá chắc chắn rằng hệ số thực nằm giữa 0. 236 và 0. 775

18. 3. 4. Khoảng tin cậy sử dụng lý thuyết bình thường

Chúng ta cũng có thể tạo khoảng tin cậy dựa trên lý thuyết bình thường. Vì các giá trị \(\widehat \theta_{educ}\) có vẻ như được phân phối chuẩn, nên chúng ta có thể xây dựng khoảng tin cậy bằng cách tính toán như sau

\[ \begin{aligned} [\widehat \theta - z_{\frac{\alpha}{2}}*SE(\theta^*), \widehat \theta + z_{\frac{\alpha}{2}

trong đó \(SE(\theta^*)\) là sai số chuẩn của các hệ số khởi động của chúng tôi, \(z\) là một hằng số và \(\widehat \theta\) là hệ số mẫu của chúng tôi. Lưu ý rằng \(z\) thay đổi tùy thuộc vào mức độ tin cậy của khoảng mà chúng tôi đang xây dựng. Vì chúng tôi đang tạo khoảng tin cậy 95%, chúng tôi sẽ sử dụng 1. 96

plt.scatter(x=duncan["education"], y=duncan["prestige"])

plt.scatter(x=duncan["education"], y=duncan["prestige"])

plt.scatter(x=duncan["education"], y=duncan["prestige"])

** Quan sát. ** Lưu ý khoảng tin cậy sử dụng lý thuyết bình thường hẹp hơn khoảng tin cậy sử dụng phương pháp phân vị, đặc biệt là về phía bên trái của khoảng

Chúng ta sẽ không đi sâu vào khoảng tin cậy của lý thuyết chuẩn tắc, nhưng nếu bạn muốn tìm hiểu thêm, hãy tham khảo X

18. 3. 4. 1. Hệ số thực sự có thể là 0?

Mặc dù chúng tôi đã quan sát thấy mối quan hệ một phần tích cực giữa giáo dục và uy tín (từ 0. 55), điều gì sẽ xảy ra nếu hệ số thực sự là 0 và không có mối quan hệ một phần nào giữa học vấn và uy tín?

Để chính thức kiểm tra xem mối quan hệ từng phần giữa giáo dục và uy tín có thật hay không, chúng tôi muốn kiểm tra các giả thuyết sau

giả thuyết không. Hệ số riêng phần thực sự là 0

Giả thuyết thay thế. Hệ số riêng phần thực sự không phải là 0

Vì chúng tôi đã xây dựng khoảng tin cậy 95% cho hệ số thực, chúng tôi chỉ cần xem liệu 0 có nằm trong khoảng này hay không. Lưu ý rằng 0 không nằm trong khoảng tin cậy của chúng tôi ở trên;

Nếu khoảng tin cậy cho hệ số thực chứa 0, thì chúng ta sẽ không có đủ bằng chứng để bác bỏ giả thuyết không. Trong trường hợp này, hệ số quan sát được \(\widehat\theta_{\texttt{education}}\) có thể là giả

18. 3. 4. 2. Phương pháp 1 Bootstrapping Reflection

Để xây dựng phân phối lấy mẫu của hệ số \(\widehat\theta_{\texttt{education}}\) và xây dựng khoảng tin cậy cho hệ số thực, chúng tôi đã lấy mẫu lại trực tiếp các quan sát và trang bị các mô hình hồi quy mới trên các mẫu bootstrap của chúng tôi. Phương pháp này hoàn toàn coi các biến hồi quy \(X_i\) là ngẫu nhiên thay vì cố định

Trong một số trường hợp, chúng tôi có thể coi \(X_i\) là cố định (ví dụ: nếu dữ liệu được lấy từ một thiết kế thử nghiệm). Trong trường hợp các biến giải thích được kiểm soát hoặc giá trị của các biến giải thích được đặt bởi người thực nghiệm, thì chúng ta có thể xem xét phương pháp khởi động thay thế sau

18. 3. 5. Thay thế. Bootstrapping phần còn lại

Một cách tiếp cận khác để kiểm tra giả thuyết trong hồi quy tuyến tính là khởi động phần dư. Cách tiếp cận này đưa ra nhiều giả định cơ bản và ít được sử dụng hơn trong thực tế. Trong phương pháp này, chúng tôi xem xét phần dư \(e_i. = y_i - X_i\widehat\beta \) làm mẫu của chúng tôi, vì vậy chúng tôi xây dựng mẫu lại bootstrap bằng cách lấy mẫu có thay thế từ các phần dư này. Khi chúng tôi xây dựng từng mẫu bootstrap, chúng tôi có thể tính toán các giá trị được trang bị mới bằng cách sử dụng các phần dư này. Sau đó, chúng tôi hồi quy các giá trị Y mới này thành các giá trị X cố định để thu được các hệ số hồi quy bootstrap

Để rõ ràng hơn, hãy chia phương pháp này thành các bước

1. Ước tính các hệ số hồi quy cho mẫu ban đầu và tính giá trị phù hợp \(\widehat y\) và phần dư \(e_i\) cho mỗi quan sát

2. Chọn các mẫu bootstrap của phần dư; . Sau đó, tính toán các giá trị \(\tilde y_i\) đã khởi động bằng cách tính toán \(\widehat y + \tilde e_i\) trong đó các giá trị phù hợp \(\widehat y_i = X_i\widehat\beta\) thu được từ hồi quy ban đầu

3. Hồi quy các giá trị \(\tilde y_i\) đã khởi động trên các giá trị \(X\) cố định để thu được hệ số hồi quy bootstrap \(\tilde \theta\)

4. Lặp lại các bước hai và ba nhiều lần để có được một số hệ số hồi quy bootstrap \(\tilde \theta_1, \tilde \theta_2, \dots \tilde \theta_n\). Chúng có thể được sử dụng để tính toán các lỗi tiêu chuẩn bootstrap và khoảng tin cậy

18. 3. 5. 1. Ước tính các hệ số thực

Bây giờ chúng ta có các hệ số hồi quy bootstrapped, chúng ta có thể xây dựng khoảng tin cậy bằng cách sử dụng các kỹ thuật tương tự như trước đây. Chúng tôi sẽ để điều này như một bài tập

18. 3. 5. 2. Bootstrapping Phản xạ dư

Hãy suy nghĩ về phương pháp này. Bằng cách gắn lại ngẫu nhiên các phần dư đã lấy mẫu lại với các giá trị đã khớp, quy trình này ngầm giả định rằng các lỗi được phân phối giống hệt nhau. Cụ thể hơn, phương pháp này giả định rằng sự phân bố dao động xung quanh đường cong hồi quy là giống nhau đối với tất cả các giá trị của đầu vào \(x_i\). Đây là một bất lợi vì các lỗi thực sự có thể có phương sai không cố định;

Mặc dù phương pháp này không đưa ra bất kỳ giả định nào về hình dạng của phân bố sai số, nhưng nó hoàn toàn giả định rằng dạng chức năng của mô hình là đúng. Bằng cách dựa vào mô hình để tạo từng mẫu bootstrap, chúng tôi giả định rằng cấu trúc mô hình là phù hợp

18. 3. 6. Tóm lược

Trong phần này, chúng tôi nêu bật các kỹ thuật bootstrapping được sử dụng trong cài đặt hồi quy tuyến tính

Nói chung, bootstrapping các quan sát được sử dụng phổ biến hơn cho bootstrapping. Phương pháp này thường hiệu quả hơn các kỹ thuật khác vì nó đưa ra ít giả định cơ bản hơn;

Chúng tôi cũng nhấn mạnh một phương pháp thay thế, có một số nhược điểm. Bootstrapping phần dư có thể được sử dụng khi chúng tôi muốn xử lý các quan sát của mình \(X\) là cố định. Lưu ý rằng phương pháp này nên được sử dụng một cách thận trọng vì nó đưa ra các giả định bổ sung về sai số và hình thức của mô hình

Bootstrap trong hồi quy tuyến tính là gì?

Làm cách nào để khởi động một mô hình trong R?

Làm cách nào để khởi động dữ liệu trong R?

Phương pháp bootstrap trong R là gì?