Hi everybody:
I?m trying to rewrite some routines originally written for SAS?s PROC
NLMIXED into LME4's glmer.
These examples came from a paper by Nelson et al. (Use of the
Probability Integral Transformation to Fit Nonlinear Mixed-Models
with Nonnormal Random Effects - 2006). Firstly the authors fit a
Poisson model with canonical link and a single normal random effect
bi ~ N(0;Sigma^2).The SAS routine was:
log_s =log(SURVT)
cens=1
proc nlmixed data=liver;
parms logsig2 = 0 b0 = -2.8 btrt = -0.54 bhrt =0.18;
xb= log_s + b0 + btrt * treat + bhrt * heart+ bi;
lambda=exp(xb);
model cens ~ poisson(lambda);
random bi ~ normal(0,exp(logsig2)) subject=INST;
run;
I obtained the same results for parameters estimates and
standarderrors, by using:
glmer(cens ~ treat + heart + offset(log(SURVT) + (1 | INST),
data=liver, family=poisson)
After that, the authors present the same model, but now bi = log(gi)
where gi has the following gamma distribution:
f(gi | theta1) = gi ^ (1/theta1 - 1) exp(-gi/theta1)/
GAMMA(1/theta1)(theta1^(1/theta1)
They used a set of transformations proposed in the paper, and with the
following rotine fitted the model achieving
estimates for the same fixed parameters and for theta1:
proc nlmixed data=liver;
parms theta1 = 1 b0 = -2:8 btrt = -0.54 bhrt =0:18;
pi = CDF('NORMAL',ai);
if(pi > 0:999999) then pi =0:999999;
gi2 = quantile('GAMMA',pi, 1/theta1);
gi = theta1 * gi2;
xb= log_s + b0 + btrt * treat + bhrt * heart+ log(gi);
lambda=exp(xb);
model cens ~ poisson(lambda);
random ai ~ normal(0,1) subject=INST;
run;
This time I' m lost. Could anyone please give me a hint?
The data set is:
liver <- as.data.frame(matrix(c(1,1,1,1,2,2,2,3,3,3,3,4,4,4,4,
+ 5,5,6,6,7,7,23.286,6.429,26.857,11.143,3.857,9.000,8.714,
+ 1.143,23.143,2.571,2.857,76.429,35.857,25.857,52.286,25.143,
+ 25.143,29.286,28.857,1.857,14.286,1,0,1,0,0,1,0,0,1,1,0,1,
+ 1,0,0,1,0,1,0,0,1,1,0,0,1,1,0,0,1,0,0,1,0,1,0,0,1,0,1,0,1,0),
+
ncol=4,dimnames=list(NULL,c("INST","SURVT","treat","heart"))))
Thanks
Cl?dio Almeida
Federal University of Minas Gerais - Brazil (55 31 33727884)
Clodio Almeida <clodioalm <at> gmail.com> writes:> > Hi everybody: >[snip]> > I obtained the same results for parameters estimates and > standarderrors, by using: > > glmer(cens ~ treat + heart + offset(log(SURVT) + (1 | INST), > data=liver, family=poisson)It's nice that this check works. I would like to see a lot more NLMIXED/glmer comparison ...> After that, the authors present the same model, but now bi = log(gi) > where gi has the following gamma distribution: > > f(gi | theta1) = gi ^ (1/theta1 - 1) exp(-gi/theta1)/ > GAMMA(1/theta1)(theta1^(1/theta1) > > They used a set of transformations proposed in the paper, and with the > following rotine fitted the model achieving > estimates for the same fixed parameters and for theta1: > > proc nlmixed data=liver; > parms theta1 = 1 b0 = -2:8 btrt = -0.54 bhrt =0:18; > pi = CDF('NORMAL',ai); > if(pi > 0:999999) then pi =0:999999; > gi2 = quantile('GAMMA',pi, 1/theta1); > gi = theta1 * gi2; > xb= log_s + b0 + btrt * treat + bhrt * heart+ log(gi); > lambda=exp(xb); > model cens ~ poisson(lambda); > random ai ~ normal(0,1) subject=INST; > run; > > This time I' m lost. Could anyone please give me a hint? > > The data set is: > liver <- as.data.frame(matrix(c(1,1,1,1,2,2,2,3,3,3,3,4,4,4,4, > + 5,5,6,6,7,7,23.286,6.429,26.857,11.143,3.857,9.000,8.714, > + 1.143,23.143,2.571,2.857,76.429,35.857,25.857,52.286,25.143, > + 25.143,29.286,28.857,1.857,14.286,1,0,1,0,0,1,0,0,1,1,0,1, > + 1,0,0,1,0,1,0,0,1,1,0,0,1,1,0,0,1,0,0,1,0,1,0,0,1,0,1,0,1,0), > + ncol=4,dimnames=list(NULL,c("INST","SURVT","treat","heart"))))It's not immediately obvious that you can do this in glmer; in fact, I suspect you can't. I don't know of an R package that can fit non-normal random effects 'out of the box'; the only solutions I know of other than SAS PROC NLMIXED are alternative programs like WinBUGS and AD Model Builder (which do have R interfaces, but have their own fairly significant learning curves). You might try asking this question on r-sig-mixed-models instead.