Bootstrap для оценки дов.интервалов

[marratt]

Здравствуйте,
Передо мной стоит задача определить доверительные интервалы для оценки отношения двух выборочных переменных: R = x/y,
где х - число животных в данном квадрате с признаком L;
y - число животных в данном квадрате.
Распределение сильно ассиметрично.
Решил использовать bootstrap для определения доверительных интервалов полученной оценки.
Для этого использую пакет 'boot'. Исходные данные в виде data frame - lv, котороя состоит из 2 столбцов - х и у, соответственно.
Код следующий:

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

lv.boot <- boot(lv, sum(lv[,1])/sum(lv[,2]), R = 1000, stype = "f")
boot.ci(lv.boot, conf = c(0.90, 0.95),
        type = c("norm", "basic", "perc", "bca"))

Однако на этапе команды boot выходит ошибка:
Error in statistic(data, original, ...) :
could not find function "statistic"
Пробовал разные варианты, писать предварительные функции, но ничего не помогает. Не могу понять, где у меня ошибка.
Как сделать так, чтобы он сделал ресэмплинг исходных данных, посчитал их отношение и определил доверительные интервалы?
Спасибо

[gamm]

То, что нужно считать, нужно оформить в виде функции (это про sum(...) - у вас выражение). См. пример в документации.

А что мешает просто использовать glm(...,family="binomial")?

Получите среднее и доверительный интервал в предположении, что все квадраты однородны. Либо нужно добавлять в модель факторы (свойства квадратов)

P.S. А если уж очень хочется бутстрапа, кто мешает сделать его руками в цикле, выбирая индексы квадратов с возвратом. Это просто, и все будет под полным контролем.

[marratt]

Здравствуйте, спасибо за советы,

То, что нужно считать, нужно оформить в виде функции (это про sum(...) - у вас выражение). См. пример в документации.

Да, я попробовал оформить в виде функции:

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

> f3 <- function(lv, id) {
+   sum(lv$l)/sum(lv$m)}
> boot.out <- boot(lv,f3,1000)

После этого процесс пошёл, но я что-то неверное указал в функции, потому что он мне выдал:

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

ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = lv, statistic = f3, R = 1000)
Bootstrap Statistics :
     original  bias    std. error
t1* 0.2272727       0           0

Если он и делал ресэмплинг, то все 1000 проб дали одно и то же отношение - что вряд ли.
Видимо, у меня в функции есть ошибка.

А что мешает просто использовать glm(...,family="binomial")?

Я не настолько искушён в регрессионном анализе, даже и не подумал его использовать для оценки интервалов. У меня просто данные не дихотомические. У меня есть квадраты, в которых было несколько животных M и несколько животных с признаком L. То есть я не могу разделить всю выборку на два класса. А в биномиальной регрессии, насколько я понимаю, используется дихотомия признаков.
Квадраты у меня неоднородны - выборка была стратифицированная (4 страты).

А если уж очень хочется бутстрапа, кто мешает сделать его руками в цикле, выбирая индексы квадратов с возвратом. Это просто, и все будет под полным контролем.

Скорее всего так и сделаю, это, конечно, не очень изящно и эргономично, зато, действительно, всё на виду, всё под контролем.
Спасибо

[gamm]

Я не настолько искушён в регрессионном анализе, даже и не подумал его использовать для оценки интервалов. У меня просто данные не дихотомические. У меня есть квадраты, в которых было несколько животных M и несколько животных с признаком L.

данные у вас как раз дихотомические, каждое из М животных либо имеет признак, либо нет
Это совершенно классическая биномиальная модель, данные для нее готовятся в виде двух столбцов: two-column matrix with the columns giving the numbers of successes and failures, т.е. в первом столбце число животных с признаком, во втором - без него

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

M<-c(10,5,7,9)
L<-c( 2,4,2,3)
Y<-cbind(L,M-L)
cur.lm<-glm(Y~1,family="binomial")
summary(cur.lm)
# prediction on logit-scale
fitted.prob<-predict(cur.lm,newdata=data.frame(1),se.fit=TRUE)
# get fitted proportion and 95% CI
fitted.prop<-1/(1+exp(-fitted.prob$fit))
fitted.prop.lo<-1/(1+exp(-(fitted.prob$fit-qnorm(0.975)*fitted.prob$se.fit)))
fitted.prop.hi<-1/(1+exp(-(fitted.prob$fit+qnorm(0.975)*fitted.prob$se.fit)))
c(fitted.prop.lo,fitted.prop,fitted.prop.hi)
#get proportion
sum(L)/sum(M)

А теперь по существу задачи. У вас данные наверняка пространственно коррелированы между соседними квадратами, и лучше бы сочинить модель с пространственной корреляцией, и посмотреть, не меняется ли доля по территории или в зависимости от разных факторов вроде ландшафта или растительности ....

[marratt]

Здравствуйте, большое спасибо за столь ценные советы и за код для glm.
Вы знаете, я когда планировал это исследование, я думал просто расчитать долю животных с признаком L, а затем доверительные интервалы по Clopper-Pearson.
Но здесь есть один нюанс, который меня смутил: моя выборка включает 123 квадрата. Но животных мы нашли всего 44.
Если я расчитываю обычную долю и Clopper-Pearson CI, то дов.интервалы будут считаться исходя из n=44.
Я почитал У.Кокрена (Методы выборочного исследования, 1976), у него есть описания оценки отношения и стандартной ошибки для него. При таком подходе n=123. Вот только мои данные не нормально распределены и ассимитричны, поэтому эту оценку стандартной ошибки я не могу использовать для обычного расчёта дов.интервалов. Вот и решил попробовать бутстреп.
Но ваш совет о пространственной корреляции всё меняет, действительно, у меня есть все основания предполагать, что распределение признака L имеет пространственную приуроченность (по сути L - это наличие определённой инфекции в животном). Буду читать литературу по пространственному анализу (если регрессионный анализ я хоть чуть-чуть да понимаю, то с пространственным анализом у меня всё плохо - ничего о нём не знаю).
Спасибо ещё раз!

[gamm]

Но ваш совет о пространственной корреляции всё меняет, действительно, у меня есть все основания предполагать, что распределение признака L имеет пространственную приуроченность (по сути L - это наличие определённой инфекции в животном).

если можно, опишите, как собирались данные, и про кого.

[marratt]

если можно, опишите, как собирались данные, и про кого.

Отправил личным сообщением