[R] Fitting GLMM models with glmer
Clodio Almeida
clodioalm at gmail.com
Fri Sep 24 22:40:53 CEST 2010
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)
More information about the R-help
mailing list