Предсказание бессмыслицы с использованием пакета, сегментированного в R

Сначала я подобрал пуассоновский glm в R следующим образом:

> Y<-c(13,21,12,11,16,9,7,5,8,8)
> X<-c(74,81,80,79,89,96,69,88,53,72)
> age<-c(50.45194,54.89382,46.52569,44.84934,53.25541,60.16029,50.33870,
+ 51.44643,38.20279,59.76469)
> dat=data.frame(Y=Y,off.set.term=log(X),age=age)
> fit.1=glm(Y~age+offset(off.set.term),data=dat,family=poisson)

Затем я попытался получить прогнозы ответа (в масштабе журнала) для нового набора данных с помощью функции predict. Обратите внимание, что я установил смещение равным нулю.

> newdat=data.frame(age=c(52.09374,50.89329,50.61472,39.13358,44.79453),off.set.term=rep(0,5))
> predict(fit.1,newdata =newdat,type="link")
        1         2         3         4         5 
-1.964381 -1.956234 -1.954343 -1.876416 -1.914839

Затем я попробовал пакет сегментированный (версия 0.3-0.0) в R и подогнал сегментированный glm следующим образом. (Последняя версия сегментированного пакета (т.е. 0.3–1.0), похоже, не поддерживает термин смещения при использовании функции прогнозирования.)

> library(segmented)
> fit.2=segmented(fit.1,seg.Z=~age,psi=list(age=mean(age)),
+ offs=off.set.term,data=newdat)

Затем я использовал функцию прогнозирования с fit.2, чтобы получить прогнозируемые значения:

> library(segmented)
> fit.2=segmented(fit.1,seg.Z=~age,psi=list(age=mean(age)),offs=off.set.term,data=newdat)
> 
> predict(fit.2,newdata =newdat,type="link")
        1         2         3         4         5 
-26.62968 -26.08611 -25.95997 -20.76125 -23.32456

Эти предсказанные значения значительно отличаются от того, которое я получил с помощью fit.1.

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

> fit.3=glm(Y~age,data=dat,family=poisson)
> newdat.2=data.frame(age=c(52.09374,50.89329,50.61472,39.13358,44.79453))
> predict(fit.3,newdata =newdat.2,type="link")
       1        2        3        4        5 
2.406016 2.395531 2.393098 2.292816 2.342261 
> fit.4=segmented(fit.3,seg.Z=~age,psi=list(age=mean(age)),data=newdat.2)
> predict(fit.4,newdata =newdat.2,type="link")
       1        2        3        4        5 
2.577669 2.524503 2.512165 2.003679 2.254396

r package prediction glm
person Stat    schedule 24.04.2014    source источник
comment
Можно ожидать разных прогнозов от разных моделей. Непонятно, почему вы должны рассматривать numberOfDrugs как смещение, если только healthValue не обязательно изменяется с количеством лекарств - а затем установка смещения на ноль означает, что все новые пациенты получают по одному лекарству? Определить лучшую модель будет ролью человека со знанием предметной области. Если они знают, что есть внезапное изменение реакции (значение для здоровья или значение для здоровья на количество лекарств, в зависимости от исключенного смещения или нет) из-за возраста, то можно использовать сегментированный термин - в противном случае я бы предпочел квадратичный член.   -  person Gavin Kelly    schedule 24.04.2014
comment
Это просто числовой пример с поддельными данными, и имя переменной не означает ничего, кроме ответа, независимой переменной и смещения. Я не определял переменные, поэтому действительно странно, что вы делаете выводы на основе своих интерпретаций! Как бы то ни было, я изменил имена, чтобы было понятнее. Мы устанавливаем смещение равным нулю, чтобы получить значения в пуассоновских glm. Здесь вместо ставок меня интересует ссылка, так что это нормально. А вы видели последние написанные мной коды? Прогнозируемые значения должны быть близки друг к другу.   -  person Stat    schedule 24.04.2014
comment
Я не уверен, как ваш комментарий может помочь ... так что, возможно, вы захотите его удалить.   -  person Stat    schedule 24.04.2014

Ответы (1)


arrow_upward
1
arrow_downward

Поскольку я получил ответ от сопровождающего сегментированного пакета, я решил поделиться им здесь. Сначала обновите пакет до версии 0.3-1.0 с помощью

install.packages("segmented",type="source")

После обновления выполнение тех же команд приводит к:

> Y<-c(13,21,12,11,16,9,7,5,8,8)
> X<-c(74,81,80,79,89,96,69,88,53,72)
> age<-c(50.45194,54.89382,46.52569,44.84934,53.25541,60.16029,50.33870,
+ 51.44643,38.20279,59.76469)
> dat=data.frame(Y=Y,off.set.term=log(X),age=age)
> fit.1=glm(Y~age+offset(off.set.term),data=dat,family=poisson)
> 
> newdat=data.frame(age=c(52.09374,50.89329,50.61472,39.13358,44.79453),off.set.term=rep(0,5))
> predict(fit.1,newdata =newdat,type="link")
        1         2         3         4         5 
-1.964381 -1.956234 -1.954343 -1.876416 -1.914839 
> 
> library(segmented)
> fit.2=segmented(fit.1,seg.Z=~age,psi=list(age=mean(age)),offs=off.set.term,data=newdat)
> predict(fit.2,newdata =newdat,type="link")
Error in offset(off.set.term) : object 'off.set.term' not found

Таким образом, срок смещения не может быть найден. Теперь уловка (на данный момент) состоит в том, чтобы сначала прикрепить newdat, а затем предсказать следующее:

> attach(newdat)
The following object is masked _by_ .GlobalEnv:

    age
> predict(fit.2,newdata =newdat,type="link")
        1         2         3         4         5 
-1.825831 -1.853842 -1.860342 -2.128237 -1.996147

Теперь результаты действительно имеют смысл. Ваше здоровье!

person Stat    schedule 24.04.2014
comment
Проголосовал за правильный ответ ... Нет предела глупости;) - person Stat; 25.04.2014
comment
из любопытства, объяснил ли сопровождающий, почему смещение перестало работать в более поздней версии? Намереваются ли они восстановить эту функциональность (или она явно работала, но фактически не работала в предыдущей версии)? - person Ben Bolker; 18.06.2014
comment
Нет, он не объяснил, почему перестал работать смещение. Он показал мне, как это сделать. Он собирается исправить это в следующей версии, если это еще не сделано. - person Stat; 18.06.2014
comment
PS на самом деле вполне приемлемо ответить на свой вопрос, например meta.stackexchange.com/questions/147697/ - person Ben Bolker; 18.06.2014