Optimum confidence interval

Вопросы по статистическому пакету R. Не обязательно гео.
Ответить
Аватара пользователя
darsvid
Активный участник
Сообщения: 135
Зарегистрирован: 29 июн 2012, 12:40
Статьи: 5
Проекты: 1
Репутация: 87
Откуда: Kyїv, Ukraine
Контактная информация:

Optimum confidence interval

Сообщение darsvid » 03 июл 2017, 12:37

Пытаюсь воспроизвести упражнение 3.2.6 из Data analysis in community and landscape ecology.

Если коротко, то известен оптимум 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"
По ответу в книжке интервал [12.8; 16.2] (что тоже меня смущает, т.к. он должен быть симметричный по отношению к оптимуму (скорее всего, должно быть 12.8;15.2 или 11.8;16.2). Но у меня, как видите, получаются совсем не те значения.
Подскажите пожалуйста, где у меня может быть ошибка?
Вложения
ex_3_2_data.csv
данные
(387 байт) 10 скачиваний
ter Braak, Looman_1986_Weighted averaging, logistic regression and the Gaussian response model.pdf
см. Appendix, c. 10
(658.47 КБ) 20 скачиваний
Последний раз редактировалось darsvid 04 июл 2017, 14:56, всего редактировалось 1 раз.

gamm
Гуру
Сообщения: 2534
Зарегистрирован: 15 окт 2010, 08:33
Репутация: 517
Ваше звание: программист
Откуда: Казань

Re: Optimum confidence interval

Сообщение gamm » 03 июл 2017, 15:04

а зачем вообще вся эта имитация ручного счета, если есть метод максимального правдоподобия для Гауссовского распределения? студентов мучить?

сделал в лоб по формулам, вроде попал. Асимметрия заложена в формуле 3.31, там к u_hat добавлена добавка, и только потом плюс-минус

Код: Выделить всё

p<-read.table("ex_3_2_data.csv",header=TRUE,sep=",")
# [1] "site" "abundance" "temp"

p.s <-sum(p$temp*p$abundance)  /sum(p$abundance)
p.ss<-sum(p$temp^2*p$abundance)/sum(p$abundance)
p.mean<-p.s
p.var <-p.ss-p.s^2
p.sd  <-sqrt(p.var)

p.lm<-lm(abundance~1+temp+I(temp^2),data=p)
b0<-coef(p.lm)[1]
b1<-coef(p.lm)[2]
b2<-coef(p.lm)[3]

cur.u<- (-b1/(2*b2))
cur.t<-1/sqrt(-2*b2)

p.vcov<-vcov(p.lm)

cur.v11<-p.vcov[2,2]
cur.v22<-p.vcov[3,3]
cur.v12<-p.vcov[2,3]

u.var<-(cur.v11+4*cur.u*cur.v12+4*cur.u^2*cur.v22)/(4*b2^2)
t.var<-cur.v22/(-8*b2^3)

t.alpha<-qt(0.975,nrow(p)-3)
cur.g<-t.alpha^2*cur.v22/(b2^2)
cur.D<-4*b2^2*u.var-cur.g*(cur.v11-cur.v12^2/cur.v22)
p.delta<-0.5*t.alpha*sqrt(cur.D)/b2
cur.CI<-(cur.u+0.5*cur.g*cur.v12/cur.v22 + c(1,-1)*p.delta)/(1-cur.g)
#  12.81210 16.21133
я бы не так делал ... у нас переменная count, это Пуассон. Разница во втором-третьем знаке, но все же ...

Код: Выделить всё

p.glm<-glm(abundance~1+temp+I(temp^2),data=p,family="poisson")
summary(p.glm)

mtr<-rbind(coef(p.lm),coef(p.glm))
colnames(mtr)<-c("Intercept","temp","temp^2")
rownames(mtr)<-c("LM","GLM")
print(round(mtr,3))

#    Intercept    temp   temp^2
#LM    2.11916 0.24973 -0.00894
#GLM   2.11201 0.26524 -0.00966

 vcov(p.glm)
 #             (Intercept)          temp     I(temp^2)
#(Intercept)  1.124954e-02 -2.040539e-03  7.721716e-05
#temp        -2.040539e-03  4.907720e-04 -2.071235e-05
#I(temp^2)    7.721716e-05 -2.071235e-05  9.183323e-07

Аватара пользователя
darsvid
Активный участник
Сообщения: 135
Зарегистрирован: 29 июн 2012, 12:40
Статьи: 5
Проекты: 1
Репутация: 87
Откуда: Kyїv, Ukraine
Контактная информация:

Re: Optimum confidence interval

Сообщение darsvid » 04 июл 2017, 14:55

gamm,

большое спасибо, особенно за альтернативный вариант с Пуассоном.

По вашему коду еще раз перепроверила у себя формулы - проблема была в скобках знаменателя varu.

Код: Выделить всё

# неправильный вариант
varu <- (v_11 + 4*u*v_12 + 4*u^2*v_22) / 4*b_2^2
# правильный вариант
varu <- (v_11 + 4*u*v_12 + 4*u^2*v_22) / (4*b_2^2)

Ответить

Вернуться в «R»