Если коротко, то известен оптимум u=14.0 и нужно посчитать для него 95% confidence interval. Формулы для расчета приведены во вложенной статье (выделены маркером в Appendix на с. 10). Набор данных в формате csv прилагаю.
Ниже приведен код по которому это все считается в R:
Код: Выделить всё
opt_int_func <- function(path2file = "ex_3_2_data.csv", alpha = 0.05, u = 14.0){
data.df <- read.csv(path2file)
n <- nrow(data.df)
print(paste0("n = ", n))
t <- abs(qt(alpha/2, n-2))
print(paste0("t = ", t))
quad.regression <- lm(formula = log(abundance) ~ poly(temp, 2, raw = TRUE), data = data.df)
v_11 <- vcov(quad.regression)[2,2]
print(paste0("v_11 = ", v_11))
v_22 <- vcov(quad.regression)[3,3]
print(paste0("v_22 = ", v_22))
v_12 <- vcov(quad.regression)[2,3]
print(paste0("v_12 = ", v_12))
b_2 <- coef(summary(quad.regression))["poly(temp, 2, raw = TRUE)2","Estimate"]
print(paste0("b_2 = ", b_2))
g <- (t^2*v_22)/b_2^2
print(paste0("g = ", g))
varu <- (v_11 + 4*u*v_12 + 4*u^2*v_22)/4*b_2^2
print(paste0("varu = ", varu))
D <- 4*b_2^2*varu - g*(v_11 - v_12^2/v_22)
print(paste0("D = ", D))
u_low <- (u + 0.5*g*v_12/v_22 + 0.5*t*sqrt(abs(D))/b_2)/(1 - g)
u_up <- (u + 0.5*g*v_12/v_22 - 0.5*t*sqrt(abs(D))/b_2)/(1 - g)
print(paste0("u_low = ", u_low))
print(paste0("u_up = ", u_up))
}
opt_int_func()
Код: Выделить всё
[1] "n = 34"
[1] "t = 2.0369333434601"
[1] "v_11 = 0.00126917422648957"
[1] "v_22 = 2.67714463650853e-06"
[1] "v_12 = -5.70434434576893e-05"
[1] "b_2 = -0.00893575905415818"
[1] "g = 0.139111351618348"
[1] "varu = 3.46584897102688e-09"
[1] "D = -7.47266172586016e-06"
[1] "u_low = 14.1788018032742"
[1] "u_up = 14.902630917845"
Подскажите пожалуйста, где у меня может быть ошибка?