Sử dụng tính không chắc chắn của thử nghiệm khi thực hiện hồi quy tuyến tính trong R

Question

Tôi có bộ dữ liệu thử nghiệm mà tôi muốn phù hợp với đa thức. Các dữ liệu bao gồm một biến độc lập, biến phụ thuộc và một sự không chắc chắn trong các phép đo sau này, ví dụ

2000 0.2084272 0.002067834 
2500 0.207078125 0.001037248 
3000 0.2054202 0.001959138 
3500 0.203488075 0.000328942 
4000 0.2013152 0.000646088 
4500 0.198933825 0.001375657 
5000 0.196375 0.000908696 
5500 0.193668575 0.00014721 
6000 0.1908432 0.000526976 
6500 0.187926325 0.001217318 
7000 0.1849442 0.000556495 
7500 0.181921875 0.000401883 
8000 0.1788832 0.001446992 
8500 0.175850825 0.0
9000 0.1728462 0.001676249 
9500 0.169889575 0.001011735 
10000 0.167  0.000326678 

(cột x, y, + -y).

tôi có thể thực hiện một sự phù hợp đa thức bằng cách sử dụng trên với ví dụ

mydata = read.table("example.txt") 
model <- lm(V2~V1+I(V1^2)+I(V1^3)+I(V1^4), data = mydata) 

nhưng điều này không tận dụng các giá trị không chắc chắn. Làm thế nào để tôi thông báo cho R rằng cột thứ ba của tập dữ liệu là một sự không chắc chắn và do đó nó nên được sử dụng trong phân tích hồi quy?

Answer 1

Với lỗi đo lường trong biến phụ thuộc không tương quan với các biến độc lập, hệ số ước tính không thiên vị nhưng lỗi chuẩn quá nhỏ. Đây là tài liệu tham khảo tôi đã sử dụng (trang 1 & 2): http://people.stfx.ca/tleo/econ370term2lec4.pdf

Tôi nghĩ bạn chỉ cần điều chỉnh các lỗi chuẩn từ những tính toán bởi lm(). Và đó là những gì tôi đã cố gắng thực hiện trong mã bên dưới. Tôi không phải là người thống kê, do đó bạn có thể muốn đăng bài để vượt qua xác nhận và yêu cầu trực giác tốt hơn.

Đối với ví dụ bên dưới, tôi giả định cột "không chắc chắn" là độ lệch chuẩn (hoặc lỗi chuẩn). Để đơn giản, tôi đã thay đổi mô hình thành: y ~ x.

# initialize dataset 
df <- data.frame(
     x = c(2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000,7500,8000,8500,9000,9500,10000), 
     y = c(0.2084272,0.207078125,0.2054202,0.203488075,0.2013152,0.198933825, 0.196375,0.193668575, 0.1908432, 0.187926325, 0.1849442, 0.181921875, 0.1788832, 0.175850825, 0.1728462,0.169889575, 0.167), 
     y.err = c(0.002067834, 0.001037248, 0.001959138, 0.000328942, 0.000646088, 0.001375657, 0.000908696, 0.00014721, 0.000526976, 0.001217318, 0.000556495, 0.000401883, 0.001446992, 0.0, 0.001676249, 0.001011735, 0.000326678) 
) 

df 

# model regression 
model <- lm(y~x, data = df) 
summary(model) 

# get the variance of the measurement error 
# thanks to: http://schools-wikipedia.org/wp/v/Variance.htm 
# law of total variance 
pooled.var.y.err <- mean((df$y.err)^2) + var((df$y.err)^2) 

# get variance of beta from model 
# thanks to: http://stats.stackexchange.com/questions/44838/how-are-the-standard-errors-of-coefficients-calculated-in-a-regression 
X <- cbind(1, df$x) 
#  (if you add more variables to the model you need to modify the following line) 
var.betaHat <- anova(model)[[3]][2] * solve(t(X) %*% X) 

# add betaHat variance to measurement error variance 
var.betaHat.adj <- var.betaHat + pooled.var.y.err 

# calculate adjusted standard errors 
sqrt(diag(var.betaHat.adj)) 

# compare to un-adjusted standard errors 
sqrt(diag(var.betaHat))