Konvertieren Sie den SAS NLMIXED-Code für die Gamma-Regression ohne Inflation in R.

11

Ich versuche, eine Null-Inflations-Regression für eine kontinuierliche Antwortvariable in R auszuführen. Mir ist eine Gamlss-Implementierung bekannt, aber ich möchte diesen Algorithmus von Dale McLerran wirklich ausprobieren, der konzeptionell etwas einfacher ist. Leider ist der Code in SAS und ich bin nicht sicher, wie ich ihn für so etwas wie nlme neu schreiben soll.

Der Code lautet wie folgt:

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;

Von: http://listserv.uga.edu/cgi-bin/wa?A2=ind0805A&L=sas-l&P=R20779

HINZUFÜGEN:

Hinweis: Hier sind keine gemischten Effekte vorhanden - nur behoben.

Der Vorteil dieser Anpassung besteht darin, dass Sie (obwohl die Koeffizienten dieselben sind, als ob Sie eine logistische Regression an P (y = 0) und eine Gammafehlerregression mit logarithmischer Verknüpfung an E (y | y> 0) separat anpassen) dies können Schätzen Sie die kombinierte Funktion E (y), die die Nullen enthält. Diesen Wert kann man in SAS (mit einem CI) anhand der Linie vorhersagenpredict (1 - p_yEQ0)*mu .

Ferner kann man benutzerdefinierte Kontrastanweisungen schreiben, um die Signifikanz von Prädiktorvariablen auf E (y) zu testen. Hier ist zum Beispiel eine andere Version des SAS-Codes, den ich verwendet habe:

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; 

Um dann "Geschenk1" gegen "Geschenk2" (b1 gegen b2) zu schätzen, können wir diese Schätzungserklärung schreiben:

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)) ; 

Kann R das tun?

a11msp
quelle
2
user779747 hat in seinem Crossposting zu Rhelp festgestellt, dass dies zuerst hier gepostet wurde. Ich habe keine spezielle Anfrage gesehen, eine solche Mitteilung in SO zu veröffentlichen, aber einige (die meisten?) Von uns Cross-HelpeRs erwarten dies, da dies die in den R-Mailinglisten angegebene Erwartung ist.
DWin

Antworten:

9

Nachdem ich einige Zeit mit diesem Code verbracht habe, scheint es mir, als ob es im Grunde genommen:

1) Führt eine logistische Regression mit der rechten Seite b0_f + b1_f*x1und y > 0als Zielvariable durch,

2) Für die Beobachtungen für die y> 0 ist , eine Regression mit dem rechten Seite durchführt b0_h + b1_h*x1, eine Gamma - Wahrscheinlichkeit und link=log,

3) Schätzt auch den Formparameter der Gammaverteilung.

Es maximiert die Wahrscheinlichkeit gemeinsam, was sehr schön ist, da Sie nur einen Funktionsaufruf ausführen müssen. Die Wahrscheinlichkeit trennt sich jedoch trotzdem, sodass Sie keine verbesserten Parameterschätzungen erhalten.

Hier ist ein R-Code, der die glmFunktion nutzt , um Programmieraufwand zu sparen. Dies ist möglicherweise nicht das, was Sie möchten, da es den Algorithmus selbst verdeckt. Der Code ist sicherlich auch nicht so sauber, wie er sein könnte / sollte.

McLerran <- function(y, x)
{
  z <- y > 0
  y.gt.0 <- y[y>0]
  x.gt.0 <- x[y>0]

  m1 <- glm(z~x, family=binomial)
  m2 <- glm(y.gt.0~x.gt.0, family=Gamma(link=log))

  list("p.ygt0"=m1,"ygt0"=m2)
}

# Sample data
x <- runif(100)
y <- rgamma(100, 3, 1)      # Not a function of x (coef. of x = 0)
b <- rbinom(100, 1, 0.5*x)  # p(y==0) is a function of x
y[b==1] <- 0

foo <- McLerran(y,x)
summary(foo$ygt0)

Call:
glm(formula = y.gt.0 ~ x.gt.0, family = Gamma(link = log))

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.08888  -0.44446  -0.06589   0.28111   1.31066  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.2033     0.1377   8.737 1.44e-12 ***
x.gt.0       -0.2440     0.2352  -1.037    0.303    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

(Dispersion parameter for Gamma family taken to be 0.3448334)

    Null deviance: 26.675  on 66  degrees of freedom
Residual deviance: 26.280  on 65  degrees of freedom
AIC: 256.42

Number of Fisher Scoring iterations: 6

Der Formparameter für die Gammaverteilung ist gleich 1 / der Dispersionsparameter für die Gammafamilie. Auf Koeffizienten und andere Dinge, auf die Sie möglicherweise programmgesteuert zugreifen möchten, kann über die einzelnen Elemente der Rückgabewertliste zugegriffen werden:

> coefficients(foo$p.ygt0)
(Intercept)           x 
   2.140239   -2.393388 

Die Vorhersage kann über die Ausgabe der Routine erfolgen. Hier ist ein weiterer R-Code, der zeigt, wie erwartete Werte generiert werden, sowie einige andere Informationen:

# Predict expected value
predict.McLerren <- function(model, x.new)
{
  x <- as.data.frame(x.new)
  colnames(x) <- "x"
  x$x.gt.0 <- x$x

  pred.p.ygt0 <- predict(model$p.ygt0, newdata=x, type="response", se.fit=TRUE)
  pred.ygt0 <- predict(model$ygt0, newdata=x, type="response", se.fit=TRUE)  

  p0 <- 1 - pred.p.ygt0$fit
  ev <- (1-p0) * pred.ygt0$fit

  se.p0 <- pred.p.ygt0$se.fit
  se.ev <- pred.ygt0$se.fit

  se.fit <- sqrt(((1-p0)*se.ev)^2 + (ev*se.p0)^2 + (se.p0*se.ev)^2)

  list("fit"=ev, "p0"=p0, "se.fit" = se.fit,
       "pred.p.ygt0"=pred.p.ygt0, "pred.ygt0"=pred.ygt0)
}

Und ein Probelauf:

> x.new <- seq(0.05,0.95,length=5)
> 
> foo.pred <- predict.McLerren(foo, x.new)
> foo.pred$fit
       1        2        3        4        5 
2.408946 2.333231 2.201889 2.009979 1.763201 
> foo.pred$se.fit
        1         2         3         4         5 
0.3409576 0.2378386 0.1753987 0.2022401 0.2785045 
> foo.pred$p0
        1         2         3         4         5 
0.1205351 0.1733806 0.2429933 0.3294175 0.4291541 

Nun zur Koeffizientenextraktion und den Kontrasten:

coef.McLerren <- function(model)
{
  temp1 <- coefficients(model$p.ygt0)
  temp2 <- coefficients(model$ygt0)
  names(temp1) <- NULL
  names(temp2) <- NULL
  retval <- c(temp1, temp2)
  names(retval) <- c("b0.f","b1.f","b0.h","b1.h")
  retval
}

contrast.McLerren <- function(b0_f, b1_f, b2_f, b0_h, b1_h, b2_h)
{
  (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))
}


> coef.McLerren(foo)
      b0.f       b1.f       b0.h       b1.h 
 2.0819321 -1.8911883  1.0009568  0.1334845 
jbowman
quelle
2
Sie haben Recht, was mit den "Teilen" passiert (dh Logit-Regression für PR (y> 0) und Gamma-Regression für E (y | y> 0), aber es ist die kombinierte Schätzung (und Standardfehler, CI). das ist von Hauptinteresse - dh E (y). Vorhersagen dieser Menge werden im SAS-Code durch (1 - p_yEQ0) * mu gemacht. Diese Formulierung ermöglicht es Ihnen, Kontraste auf den Koeffizienten für diesen kombinierten Wert
durchzuführen
@B_Miner - Ich habe einige Code + Beispiele hinzugefügt, die das Vorhersageproblem teilweise beheben, danke, dass Sie darauf hingewiesen haben.
Jbowman
Ist dies nicht nur eine separate Schätzung? In SAS gibt NLMIXED die Möglichkeit, die Punktschätzung von E (y) sowie einen CI zu schätzen (unter Verwendung der Delta-Methode, die ich glaube). Sie können auch benutzerdefinierte Kontraste der Parameter schreiben, wie oben gezeigt, um die lineare Hypothese zu testen. Es muss eine R-Alternative geben?
B_Miner
Ja und nein. Um das Beispiel zu verwenden, foo.pred$fitgibt die Rückgabe die Punktschätzung von E (y) an, aber die Komponente foo.pred$pred.ygt0$predwürde Ihnen E (y | y> 0) geben. Ich habe in der Standardfehlerberechnung für y, BTW, als se.fit zurückgegeben. Die Koeffizienten können aus den Komponenten durch Koeffizienten ( foo.pred$pred.ygt0) und Koeffizienten ( foo.pred$pred.p.ygt0) erhalten werden; Ich werde in Kürze eine Extraktionsroutine und eine Kontrastroutine schreiben.
Jbowman
Können Sie bitte beschreiben, woher dies kommt: se.fit <- sqrt (((1-p0) * se.ev) ^ 2 + (ev * se.p0) ^ 2 + (se.p0 * se.ev) ^ 2)
B_Miner