Воскрешение коэффициентов из смоделированных данных в регрессии Пуассона - PullRequest
0 голосов
/ 30 мая 2020

Я пытаюсь понять, как воскресить оценки модели из смоделированных данных в регрессиях Пуассона. Есть и другие похожие сообщения об интерпретации коэффициентов на StackExchange / CrossValidated (https://stats.stackexchange.com/questions/11096/how-to-interpret-coefficients-in-a-poisson-regression, https://stats.stackexchange.com/questions/128926/how-to-interpret-parameter-estimates-in-poisson-glm-results), но я думаю, что мой вопрос другой (хотя, по общему признанию, связанный). Я пытаюсь воскресить известные отношения, чтобы понять, что происходит с моделью. Я публикую здесь вместо CrossValidated, потому что я думаю, что это не столько статистическая интерпретация, сколько то, как я мог бы получить известную / смоделированную связь с помощью кода.

Вот некоторые смоделированные данные y и x с известными отношениями к некоторому ответу resp

set.seed(707)
x<-rnorm(10000,mean=5,sd=1)
y<-rnorm(10000,mean=5,sd=1)        
resp<-(0.5*x+0.7*y-0.1*x*y) # where I define some relationships

С линейной регрессией это очень просто:

summary(lm(resp~y+x+y:x))

Вывод показывает точную линейную зависимость между x, y и взаимодействием.

Coefficients:
              Estimate Std. Error    t value Pr(>|t|)    
(Intercept)  1.592e-14  1.927e-15  8.260e+00   <2e-16 ***
y            7.000e-01  3.795e-16  1.845e+15   <2e-16 ***
x            5.000e-01  3.800e-16  1.316e+15   <2e-16 ***
y:x         -1.000e-01  7.489e-17 -1.335e+15   <2e-16 ***

Теперь, если меня интересует регрессия Пуассона, мне нужны целые числа, я просто округляю, но сохраните связь между предикторами и ответом:

resp<-round((0.5*x+0.7*y-0.1*x*y),0)
glm1<-glm(resp~y+x+y:x,family=poisson())
summary(glm1)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.419925   0.138906   3.023   0.0025 ** 
y            0.163919   0.026646   6.152 7.66e-10 ***
x            0.056689   0.027375   2.071   0.0384 *  
y:x         -0.011020   0.005261  -2.095   0.0362 *  

Насколько я понимаю, из-за функции связи нужно возвести в степень результаты, чтобы понять их. Но здесь ни экспоненциальная оценка, ни перехват + оценка не вернут меня к исходным значениям.

> exp(0.419925+0.163919)
[1] 1.792917
> exp(0.163919)
[1] 1.178119

Как мне интерпретировать эти значения как связанные с исходным соотношением 0.7*y?

Теперь, если я помещу то же линейное уравнение в экспоненциальную функцию, я получу значения напрямую - нет необходимости использовать exp():

resp<-round(exp(0.5*x+0.7*y-0.1*x*y),0)
summary(glm(resp~y+x+y:x,family=poisson()))
Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.002970   0.045422   0.065    0.948    
y            0.699539   0.008542  81.894   <2e-16 ***
x            0.499476   0.008912  56.047   <2e-16 ***
y:x         -0.099922   0.001690 -59.121   <2e-16 ***

Может ли кто-нибудь объяснить мне, что я здесь неправильно интерпретирую и как я может найти исходные значения известной взаимосвязи без предварительного использования функции exp(), как указано выше?

1 Ответ

3 голосов
/ 30 мая 2020

Вы пренебрегаете тем фактом, что GLM Пуассона по умолчанию использует ссылку журнала (экспоненциальную обратную ссылку) (или, скорее, вы не используете эту информацию постоянно). Вы должны либо сгенерировать свои «данные» с экспоненциальной обратной ссылкой:

resp <- round(exp(0.5*x+0.7*y-0.1*x*y))

, либо соответствовать модели с помощью идентификационной ссылки (family=poisson(link="identity")). (Я бы не рекомендовал последнее, поскольку это редко бывает разумной моделью.)

Как бы то ни было, сложнее смоделировать пуассоновские данные, которые точно будут соответствовать указанному набору параметров , потому что (в отличие от гауссова, где вы можете уменьшить дисперсию до произвольно малых значений), вы не можете генерировать реальные пуассоновские данные с произвольно малым шумом. (Ваш оператор round() дает целые числа, но не результаты, распределенные по Пуассону.)

...