Konwertuj kod SAS NLMIXED dla zerowej regresji gamma na R

11

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ć?

a11msp
źródło
2
użytkownik779747 zauważył w swoim krzyżowym poście do Rhelpa, że ​​zostało to tutaj zamieszczone jako pierwsze. Nie widziałem konkretnej prośby o opublikowanie takiego powiadomienia w SO, ale niektórzy (większość?) Z naszych cross-helpeR'ów tego oczekują, ponieważ takie jest oczekiwane na listach dyskusyjnych R.
DWin

Odpowiedzi:

9

Po spędzeniu trochę czasu na tym kodzie wydaje mi się, że zasadniczo:

1) Dokonuje regresji logistycznej z prawą stroną b0_f + b1_f*x1i y > 0jako zmienną docelową,

2) W przypadku obserwacji, dla których y> 0, wykonuje regresję z prawą stroną b0_h + b1_h*x1, prawdopodobieństwo gamma i link=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ę glmfunkcję, 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ć.

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

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:

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

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:

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

I przykładowy przebieg:

> 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 

Teraz dla wyodrębnienia współczynnika i kontrastów:

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 
łucznik
źródło
2
Masz rację co do tego, co dzieje się z „częściami” (tj. Regresja logit dla PR (y> 0) i regresja gamma dla E (y | y> 0), ale jest to łączna ocena (i błędy standardowe, CI) które są głównym przedmiotem zainteresowania - tj. E (y). Prognoz tej ilości dokonuje się w kodzie SAS przez (1 - p_yEQ0) * mu. Ten preparat pozwala przeprowadzać kontrasty na współczynniki tej łącznej wartości.
B_Miner
@B_Miner - Dodałem trochę kodu + przykładów, które częściowo rozwiązują problem z prognozowaniem, dziękuję za zwrócenie na to uwagi.
łucznik
Czy to jednak nie tylko osobne szacunki? W SAS NLMIXED da możliwość oszacowania punktowego E (y), a także CI (przy użyciu metody delta, jak sądzę). Możesz również napisać zdefiniowane przez użytkownika kontrasty parametrów, jak pokazano powyżej, aby przetestować hipotezę liniową. Musi istnieć alternatywa R.
B_Miner
Cóż, tak i nie. Aby skorzystać z tego przykładu, zwrócona foo.pred$fitwartość szacuje punkt E (y), ale składnik foo.pred$pred.ygt0$predda 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.
jbowman
Czy możesz opisać, skąd pochodzi: se.fit <- sqrt (((1-p0) * se.ev) ^ 2 + (ev * se.p0) ^ 2 + (se.p0 * se.ev) ^ 2)
B_Miner