Want to convert SAS code using proc nlmixed to R code using nlme

Question

The question is similar to the one in the following post:

" troubles converting proc nlmixed (SAS) to nlme (R) "

Answer 1

This doesn't look like a mixed model at all -- what is the random effect? It looks like a slightly unusual generalized linear model: I would think this would do it.

glm((1-p)~d1+d2+d3+d4-1,family=binomial(link="log"))

(you'll have to flip the sign of the parameters yourself).

Tiwari et al. 2006 fit a model of this type.

R equivalent:

dd <- read.table("sas_vals.txt",header=TRUE,na.string=".")
g1 <- glm((1-Y)~d1+d2+d3+d4-1,family=binomial(link="log"),data=dd,
    start=rep(-0.1,4))

This is a bit difficult, probably in part because of the tiny data set (I assume this is a subset of the real data -- fitting 4 parameters to 21 binary observations would be a very bad idea ...)

Or:

library("bbmle")
g2 <- mle2(Y~dbinom(prob=1-exp(-p),size=1),
           parameters=list(p~d1+d2+d3+d4-1),
           start=list(p=0.1),
           data=dd)

The log-likelihoods suggest that the mle2 fit is better, and refitting the glm fit with those starting values works better:

g3 <- update(g1,start=-coef(g2))
g4 <- mle2(Y~dbinom(prob=1-(exp(-p)^exp(b*Treat)),size=1),
           parameters=list(p~d1+d2+d3+d4-1),
           start=list(p=0.1,b=0),
           data=dd)
summary(g4)
##      Estimate Std. Error z value  Pr(z)
## p.d1 0.106163   0.088805  1.1955 0.2319
## p.d2 0.241029   0.178263  1.3521 0.1763
## p.d3 0.105970   0.113455  0.9340 0.3503
## p.d4 0.179232   0.189951  0.9436 0.3454
## b    0.559624   0.740500  0.7557 0.4498

This seems to match the SAS results pretty well.