Ubah kode SAS NLMIXED untuk regresi gamma nol-meningkat menjadi R

Saya mencoba menjalankan regresi nol-naik untuk variabel respon kontinu dalam R. Saya menyadari implementasi gamls, tetapi saya benar-benar ingin mencoba algoritma ini oleh Dale McLerran yang secara konsep sedikit lebih mudah. Sayangnya, kodenya ada di SAS dan saya tidak yakin bagaimana menulisnya kembali untuk sesuatu seperti nlme.

Kode tersebut adalah sebagai berikut:

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;

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

MENAMBAHKAN:

Catatan: Tidak ada efek campuran hadir di sini - hanya diperbaiki.

Keuntungan dari pemasangan ini adalah (walaupun koefisiennya sama seperti jika Anda secara terpisah memasukkan regresi logistik ke P (y = 0) dan regresi kesalahan gamma dengan tautan log ke E (y | y> 0)) Anda dapat memperkirakan fungsi gabungan E (y) yang mencakup nol. Seseorang dapat memprediksi nilai ini dalam SAS (dengan CI) menggunakan garis predict (1 - p_yEQ0)*mu.

Lebih lanjut, seseorang dapat menulis pernyataan kontras khusus untuk menguji signifikansi variabel prediktor pada E (y). Sebagai contoh, ini adalah versi lain dari kode SAS yang telah saya gunakan:

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; 

Kemudian untuk memperkirakan "gift1" versus "gift2" (b1 versus b2) kita dapat menulis pernyataan estimasi ini:

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

Bisakah R melakukan ini?

Setelah menghabiskan beberapa waktu pada kode ini, tampaknya bagi saya seolah-olah pada dasarnya:

1) Melakukan regresi logistik dengan sisi kanan b0_f + b1_f*x1dan y > 0sebagai variabel target,

2) Untuk pengamatan yang y> 0, lakukan regresi dengan sisi kanan b0_h + b1_h*x1, kemungkinan Gamma dan link=log,

3) Juga memperkirakan parameter bentuk distribusi Gamma.

Ini memaksimalkan kemungkinan secara bersama, yang bagus, karena Anda hanya perlu melakukan panggilan fungsi satu. Namun, kemungkinannya terpisah, jadi Anda tidak mendapatkan perkiraan parameter yang ditingkatkan sebagai hasilnya.

Berikut adalah beberapa kode R yang memanfaatkan glmfungsi untuk menghemat upaya pemrograman. Ini mungkin bukan yang Anda inginkan, karena mengaburkan algoritma itu sendiri. Kode itu juga tidak sebersih yang seharusnya / seharusnya.

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

Parameter bentuk untuk distribusi Gamma sama dengan 1 / parameter dispersi untuk keluarga Gamma. Koefisien dan hal-hal lain yang mungkin ingin Anda akses secara programatik dapat diakses pada elemen individual dari daftar nilai pengembalian:

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

Prediksi dapat dilakukan dengan menggunakan output dari rutin. Berikut ini beberapa kode R yang menunjukkan cara menghasilkan nilai yang diharapkan dan beberapa informasi lainnya:

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

Dan contoh dijalankan:

> 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 

Sekarang untuk ekstraksi koefisien dan kontrasnya:

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