Я пытаюсь запустить регрессию с нулевым раздувом для переменной с непрерывным откликом в R. Я знаю о реализации gamlss, но я действительно хотел бы попробовать этот алгоритм Дейла Маклеррана, который концептуально немного более прост. К сожалению, код находится в SAS, и я не уверен, как переписать его для чего-то вроде nlme.
Код выглядит следующим образом:
proc nlmixed data=mydata;
parms b0_f=0 b1_f=0
b0_h=0 b1_h=0
log_theta=0;
eta_f = b0_f + b1_f*x1 ;
p_yEQ0 = 1 / (1 + exp(-eta_f));
eta_h = b0_h + b1_h*x1;
mu = exp(eta_h);
theta = exp(log_theta);
r = mu/theta;
if y=0 then
ll = log(p_yEQ0);
else
ll = log(1 - p_yEQ0)
- lgamma(theta) + (theta-1)*log(y) - theta*log(r) - y/r;
model y ~ general(ll);
predict (1 - p_yEQ0)*mu out=expect_zig;
predict r out=shape;
estimate "scale" theta;
run;
От: http://listserv.uga.edu/cgi-bin/wa?A2=ind0805A&L=sas-l&P=R20779
ДОБАВИТЬ:
Примечание: здесь нет смешанных эффектов - только исправлено.
Преимущество этой подгонки состоит в том, что (хотя коэффициенты такие же, как если бы вы отдельно подгоняли логистическую регрессию к P (y = 0) и регрессию гамма-ошибок с лог-ссылкой на E (y | y> 0)), вы можете оценить объединенную функцию E (y), которая включает в себя нули. Можно предсказать это значение в SAS (с CI), используя линию predict (1 - p_yEQ0)*mu
.
Кроме того, можно написать собственные контрастные операторы, чтобы проверить значимость предикторных переменных на E (y). Например, вот еще одна версия кода SAS, которую я использовал:
proc nlmixed data=TestZIG;
parms b0_f=0 b1_f=0 b2_f=0 b3_f=0
b0_h=0 b1_h=0 b2_h=0 b3_h=0
log_theta=0;
if gifts = 1 then x1=1; else x1 =0;
if gifts = 2 then x2=1; else x2 =0;
if gifts = 3 then x3=1; else x3 =0;
eta_f = b0_f + b1_f*x1 + b2_f*x2 + b3_f*x3;
p_yEQ0 = 1 / (1 + exp(-eta_f));
eta_h = b0_h + b1_h*x1 + b2_h*x2 + b3_h*x3;
mu = exp(eta_h);
theta = exp(log_theta);
r = mu/theta;
if amount=0 then
ll = log(p_yEQ0);
else
ll = log(1 - p_yEQ0)
- lgamma(theta) + (theta-1)*log(amount) - theta*log(r) - amount/r;
model amount ~ general(ll);
predict (1 - p_yEQ0)*mu out=expect_zig;
estimate "scale" theta;
run;
Затем, чтобы оценить «подарок1» против «подарок2» (b1 против b2), мы можем написать это утверждение:
estimate "gift1 versus gift 2"
(1-(1 / (1 + exp(-b0_f -b1_f))))*(exp(b0_h + b1_h)) - (1-(1 / (1 + exp(-b0_f -b2_f))))*(exp(b0_h + b2_h)) ;
Может ли R сделать это?
Ответы:
Потратив некоторое время на этот код, мне кажется, что он в основном:
1) выполняет ли логистическую регрессию с правой стороны
b0_f + b1_f*x1
иy > 0
в качестве целевой переменной,2) Для тех наблюдений, для которых y> 0, выполняется регрессия с правой стороны
b0_h + b1_h*x1
, гамма-вероятность иlink=log
,3) Также оценивается параметр формы гамма-распределения.
Это максимизирует вероятность совместно, что приятно, потому что вам нужно сделать только один вызов функции. Однако вероятность в любом случае разделяется, поэтому в результате вы не получите улучшенных оценок параметров.
Вот некоторый R-код, который использует
glm
функцию для экономии усилий при программировании. Это может быть не то, что вам нужно, поскольку это затеняет сам алгоритм. Код, конечно, не так чист, как мог / должен быть.Параметр формы для гамма-распределения равен 1 / параметр дисперсии для гамма-семейства. Коэффициенты и другие элементы, к которым вы могли бы получить программный доступ, могут быть доступны в отдельных элементах списка возвращаемых значений:
Прогнозирование может быть выполнено с использованием результатов процедуры. Вот еще немного кода R, который показывает, как генерировать ожидаемые значения, и некоторую другую информацию:
И пример прогона:
Теперь для извлечения коэффициентов и контрастов:
источник
foo.pred$fit
дает точечную оценку E (y), но компонентfoo.pred$pred.ygt0$pred
даст вам E (y | y> 0). Я добавил в стандартный расчет ошибок для y, кстати, вернулся как se.fit. Коэффициенты могут быть получены из компонентов коэффициентами (foo.pred$pred.ygt0
) и коэффициентами (foo.pred$pred.p.ygt0
); Я напишу процедуру извлечения и процедуру контрастирования через некоторое время.