2013-01-15 101 views
6

Tôi có một câu hỏi có lẽ là một truy vấn thống kê hơn một truy vấn liên quan đến r trực tiếp, tuy nhiên có thể là tôi chỉ gọi một gói r không chính xác vì vậy tôi sẽ đăng câu hỏi ở đây. Tôi có số liệu sau:Hồi quy tuyến tính và phi tuyến tính một cách khôn ngoan trong R

x<-c(1e-08, 1.1e-08, 1.2e-08, 1.3e-08, 1.4e-08, 1.6e-08, 1.7e-08, 
1.9e-08, 2.1e-08, 2.3e-08, 2.6e-08, 2.8e-08, 3.1e-08, 3.5e-08, 
4.2e-08, 4.7e-08, 5.2e-08, 5.8e-08, 6.4e-08, 7.1e-08, 7.9e-08, 
8.8e-08, 9.8e-08, 1.1e-07, 1.23e-07, 1.38e-07, 1.55e-07, 1.76e-07, 
1.98e-07, 2.26e-07, 2.58e-07, 2.95e-07, 3.25e-07, 3.75e-07, 4.25e-07, 
4.75e-07, 5.4e-07, 6.15e-07, 6.75e-07, 7.5e-07, 9e-07, 1.15e-06, 
1.45e-06, 1.8e-06, 2.25e-06, 2.75e-06, 3.25e-06, 3.75e-06, 4.5e-06, 
5.75e-06, 7e-06, 8e-06, 9.25e-06, 1.125e-05, 1.375e-05, 1.625e-05, 
1.875e-05, 2.25e-05, 2.75e-05, 3.1e-05) 

y2<-c(-0.169718017273307, 7.28508517630734, 71.6802510299446, 164.637259265704, 
322.02901173786, 522.719633360006, 631.977073772459, 792.321270345847, 
971.810607095548, 1132.27551798986, 1321.01923840546, 1445.33152600664, 
1568.14204073109, 1724.30089942149, 1866.79717333592, 1960.12465709003, 
2028.46548012508, 2103.16027631327, 2184.10965255236, 2297.53360080873, 
2406.98288043262, 2502.95194879366, 2565.31085776325, 2542.7485752473, 
2499.42610084412, 2257.31567571328, 2150.92120390084, 1998.13356362596, 
1990.25434682546, 2101.21333152526, 2211.08405955931, 1335.27559108724, 
381.326449703455, 430.9020598199, 291.370887491989, 219.580548355043, 
238.708972427248, 175.583544448326, 106.057481792519, 59.8876372379487, 
26.965143266819, 10.2965349811467, 5.07812046132922, 3.19125838983254, 
0.788251933518549, 1.67980552001939, 1.97695007279929, 0.770663673279958, 
0.209216903989619, 0.0117903221723813, 0.000974437796492681, 
0.000668823762763647, 0.000545308757270207, 0.000490042305650751, 
0.000468780182460397, 0.000322977916070751, 0.000195423690538495, 
0.000175847622407421, 0.000135771259866332, 9.15607623591363e-05) 

mà khi âm mưu trông như thế này: Segmentation test http://i48.tinypic.com/25pltoy.png

Tôi đã sau đó cố gắng sử dụng các gói Phân khúc để tạo ra ba hồi quy tuyến tính (rắn đường màu đen) ở ba miền (10^⁻8--10^⁻7,10^⁻7--10^⁻6 và> 10^-6) vì tôi có cơ sở lý thuyết để tìm các mối quan hệ khác nhau ở các vùng khác nhau này. Tuy nhiên, rõ ràng nỗ lực của tôi bằng cách sử dụng mã sau không thành công:

lin.mod <- lm(y2~x) 
segmented.mod <- segmented(lin.mod, seg.Z = ~x, psi=c(0.0000001,0.000001)) 

Vì vậy, câu hỏi đầu tiên của tôi - có thêm thông số nào của phân khúc mà tôi có thể tinh chỉnh khác với điểm ngắt không? Cho đến nay, tôi hiểu rằng tôi có các lần lặp được đặt ở mức tối đa như mặc định tại đây.

Câu hỏi thứ hai của tôi là: tôi có thể thử phân đoạn bằng gói nls không? Dường như hai khu vực đầu tiên trên cốt truyện (10^⁻8--10^⁻7 và 10^-7--10^-6) là tiếp tục từ tuyến tính sau đó phần cuối cùng vì vậy có lẽ một hàm đa thức sẽ tốt hơn đây?

Ví dụ về kết quả tôi thấy có thể chấp nhận được Tôi đã chú thích cốt truyện ban đầu bằng tay: Annotated segmentation example http://i45.tinypic.com/zjb439.jpg.

Chỉnh sửa: Lý do sử dụng khớp tuyến tính là sự đơn giản mà chúng cung cấp, với mắt chưa được đào tạo của tôi, nó sẽ yêu cầu hàm phi tuyến khá phức tạp để khôi phục tập dữ liệu dưới dạng đơn vị duy nhất. Một suy nghĩ đã vượt qua tâm trí của tôi là để phù hợp với một mô hình lognormal dữ liệu vì điều này có thể làm việc cho skew cùng một trục x đăng nhập. Tôi không có đủ năng lực trong R để làm điều này tuy nhiên như kiến ​​thức của tôi chỉ kéo dài đến fitdistr mà cho đến nay như tôi hiểu sẽ không làm việc ở đây.

Bất kỳ trợ giúp hoặc hướng dẫn nào có liên quan sẽ được đánh giá cao nhất.

Trả lời

4

Nếu bạn không hài lòng với gói segmented, bạn có thể thử gói earth với thuật toán mars. Nhưng ở đây, tôi thấy rằng kết quả của mô hình phân đoạn là rất chấp nhận được. xem R-Squared bên dưới.

lin.mod <- lm(y2~x) 
segmented.mod <- segmented(lin.mod, seg.Z = ~x, psi=c(0.0000001,0.000001)) 
summary(segmented.mod) 

Meaningful coefficients of the linear terms: 
       Estimate Std. Error t value Pr(>|t|)  
(Intercept) -2.163e+02 1.143e+02 -1.893 0.0637 . 
x   4.743e+10 3.799e+09 12.485 <2e-16 *** 
U1.x  -5.360e+10 3.824e+09 -14.017  NA  
U2.x   6.175e+09 4.414e+08 13.990  NA  

Residual standard error: 232.9 on 54 degrees of freedom 
Multiple R-Squared: 0.9468, Adjusted R-squared: 0.9419 

Convergence attained in 5 iterations with relative change 3.593324e-14 

Bạn có thể kiểm tra kết quả bằng cách vẽ các mô hình:

plot(segmented.mod) 

enter image description here

Để có được hệ số lô, bạn có thể làm điều này:

 intercept(segmented.mod) 
$x 
              Est. 
intercept1 -216.30 
intercept2 3061.00 
intercept3   46.93 

> slope(segmented.mod) 
$x 
             Est.   St.Err.  t value  CI(95%).l  CI(95%).u 
slope1  4.743e+10 3.799e+09  12.4800  3.981e+10  5.504e+10 
slope2 -6.177e+09 4.414e+08 -14.0000 -7.062e+09 -5.293e+09 
slope3 -2.534e+06 5.396e+06  -0.4695 -1.335e+07  8.285e+06 
+0

Ah lạ tôi không thể làm việc đó lần đầu tiên. Cảm ơn những lời khuyên và các thông tin thêm về việc thu thập các hệ số vv. Bạn quan tâm đến việc tôi có thể đi về một thứ gì đó giống như một điều bất thường đối với dữ liệu đó? Điều này có thể đảm bảo một câu hỏi riêng biệt nhưng trong trường hợp có một cách giải quyết nhanh chóng, tôi nghĩ rằng tôi sẽ thêm nó như là một bình luận ở đây. Cảm ơn một lần nữa. – user1912925

Các vấn đề liên quan