Блог

Я нашел функции onls::onls() и pracma::odregress() которые вычисляют модели ортогональной регрессии. Я хотел бы построить такие модели в том же стиле geom_smooth() , что и линия регрессии, окруженная 95% доверительным интервалом.

Пример:

 example <- structure(list(y = c(-28.9143374863044, -28.5783512160246, -29.1751498307569, 
-28.5613677412358, -29.2441600709021, -29.1848482932202, -29.469712350617, 
-29.1212786695474, -29.3338385227209, -29.0582324840251, -29.1159002526588, 
-29.1384485361936, -29.4743426548081, -29.242305699462, -29.5517891592378, 
-29.1701701877517, -29.2337122509592, -29.150317639976, -29.139526754614, 
-29.05974643127, -29.0540797909476, -29.0859798970361, -29.27517072563, 
-29.1907525452561, -30.0965246973573, -28.9734662257987, -29.6953578711591, 
-28.2014460687026, -30.0621997994278, -27.9399550295493, -29.8886842413551, 
-29.6609659140518, -29.6920474706673, -30.2418230320867, -29.8334571372628, 
-29.8626462112615, -29.9051818751105, -29.6518825347484, -29.5380886463871, 
-29.7500527026688, -29.6095990506199, -29.6049957701729, -29.5368579894466, 
-29.5861340837645, -29.5737037489314, -29.5773848425703, -28.0265409956043, 
-28.0899954900073, -28.265152586989, -28.0062832808179, -27.7205565228848, 
-27.4041257575861, -28.1113851658386, -26.914663492446, -27.877772497213, 
-27.0684956870887, -27.9276723508022, -27.7588907638397, -27.3710663654935, 
-27.3623535825255, -27.7783142763593, -28.5132310123219, -28.5193067297636, 
-28.5283974320574, -28.6153706663899, -28.6816032262091, -29.1043640141426, 
-28.44589108955, -28.6614098552091, -28.7403207700811), x = c(33.1158714294, 
18.6527993810972, 17.0276514703819, 22.3627925702962, 18.170924813473, 
32.0677953809724, 46.5216445923, 34.9911138888596, 25.0910229505442, 
13.9473438263, 17.381641499988, 17.014380035215, 40.9107205320526, 
52.2695803285185, 58.9499627404227, 40.5894751586832, 23.496896254444, 
33.6412616569372, 14.7548102820616, 46.3057677573658, 14.280050708175, 
31.2877073530984, 18.8534870545271, 16.5168182808868, 63.9908365598676, 
33.7277991683148, 35.4163778417314, 32.1050571361531, 51.3240160147292, 
24.4237814340378, 39.3334452128324, 53.8079129732769, 43.26844558712, 
58.3003234863, 43.934151887875, 76.8046441618721, 64.8779439305438, 
46.8684772359235, 66.4989547729, 41.9780584414396, 50.2248225396345, 
58.8492643072032, 64.5647735596, 48.3225469025232, 60.4074024077677, 
57.3789336302925, 11.2785320282, 11.3491302769043, 7.59091310831495, 
18.4789668943737, 5.84773873549871, 10.6156844347299, 15.7432512138035, 
11.4885938379565, 7.74754936760848, 12.1071624756, 14.9075944237136, 
20.9201573163328, 30.2789412366595, 33.8582180028129, 15.4269225956373, 
8.53801707561128, 10.1814249853966, 7.33018941782735, 8.42749268077253, 
9.74786459733547, 10.5363144200841, 10.7873065304121, 16.7602893825786, 
12.7551904319156)), class = "data.frame", row.names = c(1L, 2L, 
3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 
17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 
30L, 31L, 32L, 33L, 34L, 35L, 36L, 39L, 40L, 41L, 42L, 43L, 44L, 
45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 
58L, 59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L, 69L, 70L, 
71L, 72L))

mod <- onls(y ~ a*x b, data = example, start = list(a = 0.03, b = -28))
newData <- data.frame(x = seq(min(example$x), max(example$x), 0.1))
newData$y <- predict(mod, newdata = newData)

plot(y ~ x, data = newData, type = "l", col = "red")
points(y ~ x, data = example)

# for a regular lm() model the subsequent steps would be
conf <- predict(mod, newdata = newData, interval = "confidence", level = 0.95)
lines(newData$x, conf[,2])
lines(newData$x, conf[,3])
 

Однако последние шаги не дадут никакого полезного результата при применении к onls модели. Существуют ли какие-либо методы для вычисления или оценки этих доверительных интервалов?

Редактировать: как упоминал Дэни, onls пакет содержит confint.onls() функцию. Это даст мне верхнюю и нижнюю оценки для каждого параметра регрессии при заданном уровне достоверности:

 confint(mod, level = 0.95)
 

Я мог бы сделать что-то вроде

 conf_a <- confint(mod, param = "a", level = 0.95)
conf_b <- confint(mod, param = "b", level = 0.95)
 

и вычислите экстремумы для каждого x

 x <- seq(min(example$x), max(example$x), 0.1)
test <- cbind(
   conf_a[1]*x conf_b[1],
   conf_a[1]*x conf_b[2],
   conf_a[2]*x conf_b[1],
   conf_a[2]*x conf_b[2]
)

maxima <- vector()
for(i in 1:length(x)){
   maxima[i] <- max(test[i,])
}
 

но это не совсем похоже на то, что я ожидал, и я не совсем уверен, что это правильный подход.