Próbuję uruchomić regresję z zerowym napełnieniem dla zmiennej ciągłej odpowiedzi w R. Jestem świadomy implementacji gamlss, ale naprawdę chciałbym wypróbować ten algorytm Dale'a McLerrana, który jest koncepcyjnie nieco prostszy. Niestety kod znajduje się w SAS i nie jestem pewien, jak go ponownie napisać dla czegoś takiego jak nlme.
Kod jest następujący:
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;
Od: http://listserv.uga.edu/cgi-bin/wa?A2=ind0805A&L=sas-l&P=R20779
DODAJ:
Uwaga: Nie ma tu żadnych efektów mieszanych - tylko naprawione.
Zaletą tego dopasowania jest to, że (mimo że współczynniki są takie same, jakbyś osobno dopasował regresję logistyczną do P (y = 0) i regresję błędu gamma z logarytmicznym łączem do E (y | y> 0)), możesz oszacuj połączoną funkcję E (y), która zawiera zera. Można przewidzieć tę wartość w SAS (z CI) za pomocą wiersza predict (1 - p_yEQ0)*mu
.
Ponadto można napisać niestandardowe instrukcje kontrastu, aby przetestować istotność zmiennych predykcyjnych na E (y). Na przykład, oto inna wersja kodu SAS, którego użyłem:
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;
Następnie, aby oszacować „prezent1” w porównaniu z „prezentem2” (b1 kontra b2), możemy napisać to wyrażenie szacunkowe:
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)) ;
Czy R może to zrobić?
Odpowiedzi:
Po spędzeniu trochę czasu na tym kodzie wydaje mi się, że zasadniczo:
1) Dokonuje regresji logistycznej z prawą stroną
b0_f + b1_f*x1
iy > 0
jako zmienną docelową,2) W przypadku obserwacji, dla których y> 0, wykonuje regresję z prawą stroną
b0_h + b1_h*x1
, prawdopodobieństwo gamma ilink=log
,3) Ocenia także parametr kształtu rozkładu gamma.
Maksymalnie zwiększa to prawdopodobieństwo, co jest miłe, ponieważ wystarczy wykonać tylko jedno wywołanie funkcji. Jednak prawdopodobieństwo i tak się rozdziela, więc nie otrzymujesz lepszych oszacowań parametrów.
Oto kod R, który wykorzystuje tę
glm
funkcję, aby zaoszczędzić na programowaniu. To może nie być to, co chcesz, ponieważ zaciemnia sam algorytm. Kod z pewnością nie jest tak czysty, jak mógłby / powinien być.Parametr kształtu dla rozkładu Gamma jest równy 1 / parametr dyspersji dla rodziny Gamma. Współczynniki i inne rzeczy, do których możesz chcieć uzyskać dostęp programowy, są dostępne w poszczególnych elementach listy wartości zwracanych:
Prognozowania można dokonać przy użyciu wyniku procedury. Oto trochę więcej kodu R, który pokazuje, jak wygenerować oczekiwane wartości i kilka innych informacji:
I przykładowy przebieg:
Teraz dla wyodrębnienia współczynnika i kontrastów:
źródło
foo.pred$fit
wartość szacuje punkt E (y), ale składnikfoo.pred$pred.ygt0$pred
da E (y | y> 0). Dodałem w standardowym obliczeniu błędu dla y, BTW, zwrócone jako se.fit. Współczynniki można uzyskać ze składników za pomocą współczynników (foo.pred$pred.ygt0
) i współczynników (foo.pred$pred.p.ygt0
); Niedługo napiszę procedurę ekstrakcji i procedurę kontrastu.