I know my question might not be exactly on the usual SO topics, but I'm really desperate. I've been trying to find what's wrong with my code for months.
I need the SAS equivalent with PROC NLMIXED to get the complete variance-covariance matrix of this model, written on R with lme4::glmer :
(after "poissonizing" my survival data with the library surrosurv, it looks like the snippet below. var1, var2 and cov are factor variables)
data_pois = surrosurv::poissonize(data, interval.width = interval.width, factors = c ('var1', 'var2', 'cov'), compress = FALSE)
head(data_pois)
id event var1 var2 cov time interval
1225 1 0 -0.5 13 1 0.0005669899 0
822 2 1 -0.5 9 1 0.0947011826 0
611 3 1 -0.5 7 1 0.0951096974 0
235 4 1 -0.5 3 0 0.1028146920 0
1847 5 0 -0.5 19 1 0.1122793044 0
1321 6 1 -0.5 14 1 0.1298256694 0
mod = lme4::glmer(
formula = event ~ -1 + interval : factor(var1) + cov + offset(log(time))+ (0 + var1 | var2 ) + (1 | var2 ),
data = data_pois,
family = poisson(link = "log"),
control=glmerControl(optCtrl=list(maxfun=2e5), optimizer="bobyqa"))
The details of the covariate are:
I'm assuming you want an overall independent random intercept and another set of correlated ones for each non-reference level of var1.
From the information that you provided and with my admittedly pretty superficial knowledge of both languages and the underlying statistics, I think the matching NLMIXED call could look something like this:
proc nlmixed data=data_pois;
/* Create dummies as needed/desired - more will be required if >2 levels exist. */
dummy_var1 = (var1 eq -0.5);
dummy_cov = (cov eq 0);
/* Random effects, see caveats below. */
random r1 r2 ~ normal([0, 0], [rsig1, 0, rsig2]) subject=var2;
/* Linear predictor. */
eta = b1*interval + b2*dummy_var1*interval + b3*dummy_cov1 + r1*dummy_var1 + r2 + log(time);
mu = exp(eta);
/* Option 1: use built-in distribution. */
model event ~ poisson(mu);
/* Option 2: write out log-likelihood by hand. */
ll = event * log(mu) - mu - lgamma(event + 1);
model event ~ general(ll); /* The DV can be anything here, its value doesn't matter. */
run;
This is quite bare-bones, I would almost certainly add a parms statement to initialize model parameters to more sensible values such as those precomputed from a fixed effects-only model. It's usually also a good idea to express your variances in the log scale (which will double as a bounds statement) and add extra boundary checks before making any transformations. You'll probably want to include additional tuning and/or convergence options in the NLMIXED call as well.
Finally, returning to the random effects, it really depends on what you want to do here exactly. A limitation of PROC NLMIXED is that random effects have to be nested if you're adding more than one - this is not an issue here if you're using a single blocking level. It's still up to you to define how the additional levels of r1 (to be added) covary if there's more than 2 levels of var1, I believe glmer makes this unstructured by default.