GLM о полиции останавливает данные в книге Гельмана / Хилла

Question

Кто-нибудь работал с полицией Нью-Йорка, останавливает данные, упомянутые в книге Гелмана, Хилла Анализ данных с использованием Reg.и многоуровневое моделирование (ARM).Данные находятся в

http://www.stat.columbia.edu/~gelman/arm/examples/police/

файл frisk_with_noise.dat.Я удалил часть описания этих данных, переименовал в past.arrests как аресты, сохранил их как frisk.dat.Вызывается glm из R следующим образом:

library ("foreign")
frisk <- read.table ("frisk.dat", header=TRUE)
attach (frisk)
glm(formula = stops ~ 1, family=poisson, offset=log(arrests))

Вызов glm прямо из книги ARM.В любом случае я получаю ошибку:

Error: NA/NaN/Inf in foreign function call (arg 4)

Есть идеи?У Гельмана есть фрагмент кода в том же каталоге с именем Police_setup.R, который должен иметь некоторый код очистки, но он тоже не работает.

Answer 1

Я не вернулся, чтобы посмотреть, что именно делает Гельман в этой главе (моя копия книги находится в хранилище ...), но особая проблема этого примера в том, что в некоторых случаях «аресты» равны нулю случаев, поэтому использование журнала (арестов) в качестве смещения вызывает проблемы. (Вам не нужна библиотека (иностранная), и использование аргумента данных для glm обычно безопаснее / лучше, чем использование attach ().)

X <- read.table("frisk_with_noise.dat",skip=6,header=TRUE)
names(X)[3] <- "arrests"
glm(stops~1,family=poisson,offset=log(arrests),data=X,
    subset=arrests>0)

Answer 2

Приведенный выше код работает, но результаты анализа отличаются от книги.Согласно этому посту авторам пришлось вручную изменять данные из-за проблем с конфиденциальностью.

Answer 3

Для примера из главы 6 Гельмана я считаю, что сначала вам нужно объединиться по преступности. Файл frisk_with_noise.dat содержит 900 наблюдений, по одной записи на этническую группу, по каждому участку, за преступление (75 * 3 * 4 = 900). Но пример выходных данных в главе 6 показывает n = 225 (75 * 3). Таким образом, расширение кода Бена чем-то вроде этого немного приблизит вас к репликации вывода:

library(arm) # for display() function
X <- read.table("frisk_with_noise.dat",skip=6,header=TRUE)
names(X)[3] <- "arrests"
X <- aggregate(cbind(stops, arrests) ~ precinct + eth, data=X, sum)
fit.1 <- glm(stops~1,family=poisson,offset=log(arrests),data=X,
             subset=arrests>0)
display(fit.1)

Но, как отмечено в примечании вверху файла frisk_with_noise.dat, добавлен шум, поэтому результаты не могут быть точно воспроизведены.