Convertir le code SAS NLMIXED pour une régression gamma gonflée à zéro en R

11

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?

a11msp
la source
2
L'utilisateur 779747 a noté dans sa publication croisée à Rhelp que cela avait été publié ici en premier. Je n'ai pas vu de demande spécifique d'afficher un tel avis dans SO, mais certains (la plupart?) D'entre nous nous y attendons, car c'est l'attente indiquée dans les listes de diffusion R.
DWin

Réponses:

9

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*x1et y > 0comme 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 et link=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 glmfonction 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.

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

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:

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

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:

# 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)
}

Et un exemple d'exécution:

> 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 

Maintenant pour l'extraction des coefficients et les contrastes:

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
la source
2
Vous avez raison en ce qui concerne les "parties" (c'est-à-dire la régression logit pour PR (y> 0) et la régression gamma pour E (y | y> 0) mais il s'agit de l'estimation combinée (et des erreurs types, CI) qui sont d'intérêt principal - c'est-à-dire E (y). Les prévisions de cette quantité sont faites dans le code SAS par (1 - p_yEQ0) * mu. Cette formulation vous permet d'effectuer des contrastes sur les coefficients de cette valeur combinée.
B_Miner
@B_Miner - J'ai ajouté du code + des exemples qui résolvent partiellement le problème de prédiction, merci de l'avoir signalé.
jbowman le
N'est-ce pas seulement des estimations distinctes? En SAS, NLMIXED donnera la possibilité d'estimer l'estimation ponctuelle de E (y) ainsi que d'un CI (en utilisant la méthode delta je crois). En outre, vous pouvez écrire des contrastes définis par l'utilisateur des paramètres comme je l'ai montré ci-dessus pour tester l'hypothèse linéaire. Il doit y avoir une alternative R?
B_Miner
Eh bien, oui et non. Pour utiliser l'exemple, le résultat foo.pred$fitdonne l'estimation ponctuelle de E (y), mais le composant foo.pred$pred.ygt0$predvous 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.
jbowman
Pouvez-vous décrire d'où cela vient: se.fit <- sqrt (((1-p0) * se.ev) ^ 2 + (ev * se.p0) ^ 2 + (se.p0 * se.ev) ^ 2)
B_Miner