J'essaie d'exécuter une régression zéro gonflée pour une variable de réponse continue dans R. Je connais une implémentation gamlss, mais j'aimerais vraiment essayer cet algorithme de Dale McLerran qui est conceptuellement un peu plus simple. Malheureusement, le code est en SAS et je ne sais pas comment le réécrire pour quelque chose comme nlme.
Le code est comme suit:
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;
De: http://listserv.uga.edu/cgi-bin/wa?A2=ind0805A&L=sas-l&P=R20779
AJOUTER:
Remarque: Il n'y a pas d'effets mixtes ici - seulement fixes.
L'avantage de cet ajustement est que (même si les coefficients sont les mêmes que si vous ajustez séparément une régression logistique à P (y = 0) et une régression d'erreur gamma avec un lien logarithmique à E (y | y> 0)), vous pouvez estimer la fonction combinée E (y) qui comprend les zéros. On peut prédire cette valeur en SAS (avec un CI) en utilisant la ligne predict (1 - p_yEQ0)*mu
.
De plus, on est capable d'écrire des déclarations de contraste personnalisées pour tester la signification des variables prédictives sur E (y). Par exemple, voici une autre version du code SAS que j'ai utilisé:
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;
Ensuite, pour estimer "cadeau1" par rapport à "cadeau2" (b1 par rapport à b2), nous pouvons écrire cette déclaration d'estimation:
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)) ;
Est-ce que R peut faire ça?
Réponses:
Ayant passé un peu de temps sur ce code, il me semble que c'est essentiellement:
1) Effectue une régression logistique avec le côté droit
b0_f + b1_f*x1
ety > 0
comme variable cible,2) Pour les observations pour lesquelles y> 0, effectue une régression avec le côté droit
b0_h + b1_h*x1
, une probabilité gamma etlink=log
,3) Estime également le paramètre de forme de la distribution gamma.
Cela maximise la probabilité conjointement, ce qui est bien, car vous n'avez qu'à effectuer l'appel d'une fonction. Cependant, la probabilité se sépare de toute façon, vous n'obtenez donc pas d'estimations de paramètres améliorées.
Voici un code R qui utilise la
glm
fonction pour économiser l'effort de programmation. Ce n'est peut-être pas ce que vous aimeriez, car cela obscurcit l'algorithme lui-même. Le code n'est certainement pas aussi propre qu'il pourrait / devrait l'être non plus.Le paramètre de forme pour la distribution Gamma est égal à 1 / le paramètre de dispersion pour la famille Gamma. Les coefficients et autres éléments auxquels vous pourriez souhaiter accéder par programmation sont accessibles sur les éléments individuels de la liste de valeurs de retour:
La prédiction peut être effectuée en utilisant la sortie de la routine. Voici un peu plus de code R qui montre comment générer des valeurs attendues et d'autres informations:
Et un exemple d'exécution:
Maintenant pour l'extraction des coefficients et les contrastes:
la source
foo.pred$fit
donne l'estimation ponctuelle de E (y), mais le composantfoo.pred$pred.ygt0$pred
vous donnera E (y | y> 0). J'ai ajouté dans le calcul de l'erreur standard pour y, BTW, retourné comme se.fit. Les coefficients peuvent être obtenus à partir des composants par les coefficients (foo.pred$pred.ygt0
) et les coefficients (foo.pred$pred.p.ygt0
); J'écrirai une routine d'extraction et une routine de contraste dans peu de temps.