topic Meta analysis using Glimmix in Statistical Procedures

Meta analysis using Glimmix

G_Le_Teuff — Thu, 12 Jan 2012 17:04:45 GMT

Dear SAS users,

I would like to have some information about Proc glimmix in order to perform a meta-analysis using count data (incidence rate).

There are two main approach in meta analysis : fixed and random effect. The first one is based on the inverse-variance weighted.

Some examples using proc mixed for normal data exist (See Normand, Stat Med). Based on these exemples I will try to apply on cound data using proc glimmix. Two programs are possible. One based on weighted option and the second consists to fix the estimated variance of incidence rate.

proc glimmix data=count_data /*empirical*/;

class study;

model tox = / dist=nb link=log offset=logPA2 solution; * tox = number of adverse event - PA2 = sum of person-years;

random intercept / subject=study ; * to test heterogeneity between study;

weight inv_varln; * weight = inverse variance;

run;

proc glimmix data=count_data;

class study;

model tox = / dist=nb link=log offset=logPA2 solution;

random intercept / subject=study;

random intercept / subject=study group=study; * I use it but not sure, the objective is to sepcify the variance of incidence rate of each study;

parms (0)

(0.125) (0.333) (0.031) (0.2) (0.0625) * I fixe the variance of ln(incidence rate) for each study;

(0.0169) (0.005) (0.007) (0.055) (0.0476)

(0.083) (0.033) (0.1) (0.0625) (0.0212)

(1) / hold=(2 to 15); * the last value is the scale parameter;

run;

These two programs did not give similar results. Could someone explain me why ? While using proc mixed with weight and repeated give similar result in norma data. Do I use the variance of incidence and not the variance of ln incidence ?

I am not sure to the second random statement using group=study statement. I insert it in order to have a variance -covariance matrix with size equal to the number of studies.

Best regards

Gwénaël

Meta analysis using Glimmix

SteveDenham — Fri, 13 Jan 2012 12:42:03 GMT

I've never attempted anything like this, but I did notice one curiosity in your code for the second model. You have two random statements, but they are much the same, so that the first statement results in an odd occurrence. What about a statement that looks like:

random intercept/subject=study type=un(1):

I think this would give a diagonal matrix, with the diagonal having the variances for each study. Perhaps lvm will comment, and help out on this one.

Steve Denham

Meta analysis using Glimmix

G_Le_Teuff — Fri, 13 Jan 2012 14:27:42 GMT

Dear Steve , I used the code you indicated (See below) but it doesn't work.

proc glimmix data=count_data;

class study;

model tox = / dist=nb link=log offset=logPA2 solution;

random intercept / subject=study;

random intercept / subject=study type=un(1);

parms (0)

(0.125) (0.333) (0.031) (0.2) (0.0625) * I fixe the variance of ln(incidence rate) for each study;

(0.0169) (0.005) (0.007) (0.055) (0.0476)

(0.083) (0.033) (0.1) (0.0625) (0.0212)

(1) / hold=(2 to 15); * the last value is the scale parameter;

run;

ERROR: 3 PARMS must be given instead of 20

I tried

random _residual_ / group=study; but did not converge

random

_residual_ / subject study; ERROR: 3 PARMS must be given instead

random

intercept / subject=study type=un(1 ) group=study; it works but Matrix G not defnite positive

Cov Parm Subject Group Estimate Error

Intercept study 0.000499 .

UN(1,1) study study 1 0.1250 .

UN(1,1) study study 2 0.3330 .

UN(1,1) study study 3 0.03100 .

UN(1,1) study study 4 0.2000 .

UN(1,1) study study 5 0.06250 .

UN(1,1) study study 6 0.01690 .

UN(1,1) study study 7 0.005000 .

UN(1,1) study study 8 0.007000 .

UN(1,1) study study 9 0.05500 .

UN(1,1) study study 10 0.04760 .

UN(1,1) study study 11 0.08300 .

UN(1,1) study study 12 0.03300 .

UN(1,1) study study 13 0.1000 .

UN(1,1) study study 14 0.06250 .

UN(1,1) study study 15 0.02120 .

UN(1,1) study study 16 0.02500 .

UN(1,1) study study 17 0.06250 .

UN(1,1) study study 21 6E-19 .

Scale 0.1579 0.08449

Solutions for Fixed Effects

Standard

Effect Estimate Error DF t Value Pr > |t|

Intercept -2.6138 0.1188 0 -22.00 .

So I am wondering how fix the covariance matrix using random statement

Gwénaël

Meta analysis using Glimmix

SteveDenham — Fri, 13 Jan 2012 15:19:31 GMT

I worry about the three parms instead of 20 error. I think it is coming from the following: there are three parms: first random statement = variance between studies, second random statement = also variance between studies, and residual/scale. This would also explain why you get G matrix not positive definite, as the matrix contains a redundant column when both random statements are included.

What about:

proc glimmix data=count_data;

class study;

model tox = / dist=nb link=log offset=logPA2 solution;

random intercept/subject=study group=study type=chol(1) solution;

parms

(0.125) (0.333) (0.031) (0.2) (0.0625)

(0.0169) (0.005) (0.007) (0.055) (0.0476)

(0.083) (0.033) (0.1) (0.0625) (0.0212)

(1) / hold=(1to 15);

run;

Going to a single random statement and using the Cholesky root. I offer type=chol, as it ensures that the varainace-covariance matrix is at least positive semidefinite, so that should help, as will removing the redundant random intercept/subject=study line.

Good luck.

Steve Denham

Meta analysis using Glimmix

SteveDenham — Fri, 13 Jan 2012 15:24:53 GMT

Or, since it looks like there is some variability accounted for by intercept study (0.000429), you might try:

proc glimmix data=count_data;

class study;

model tox = / dist=nb link=log offset=logPA2 solution;

random intercept/subject=study:

random intercept/subject=study group=study type=chol(1) solution;

parms (0)

(0.125) (0.333) (0.031) (0.2) (0.0625)

(0.0169) (0.005) (0.007) (0.055) (0.0476)

(0.083) (0.033) (0.1) (0.0625) (0.0212)

(1) / hold=(2 to 16);

run;

This will estimate the first and last parameter, and hold the parameters for the fifteen studies constant. I would still stay with the Cholesky parameterization of the G matrix.

Steve Denham

Meta analysis using Glimmix

G_Le_Teuff — Sat, 14 Jan 2012 11:24:26 GMT

Steve,

I ran your 2 proposed models and compare with a Negative Binomial random model using weight statement(=inverse variance) and a NB random model with type=simple.

Note: I updated the number of fixed variance parameters.

1 - Model with weight (inverse variance)

proc glimmix data=count_data /*empirical*/;

class study;

model tox = / dist=nb link=log offset=logPA2 solution;

random intercept / subject=study ;

weight inv_varln;

where groupe=1;

run;

-2 Res Log Pseudo-Likelihood 29.74

Generalized Chi-Square 16.87

Gener. Chi-Square / DF 0.99

Covariance Parameter Estimates

Standard

Cov Parm Subject Estimate Error

Intercept study 0.2245 0.1309

Scale 0.8371 1.8751

Solutions for Fixed Effects

Standard

Effect Estimate Error DF t Value Pr > |t|

Intercept -2.6600 0.1244 17 -21.38 <.0001

2 - random Model fixing variance covariance (R)

proc glimmix data=count_data;

class study;

model tox = / dist=nb link=log offset=logPA2 solution;

random intercept / subject=study;

random _residual_ / subject=study group=stud y type=simple;

parms (0)

(0.125) (0.333) (0.031) (0.2) (0.0625)

(0.0169) (0.005) (0.007) (0.055) (0.0476)

(0.083) (0.033) (0.1) (0.0625) (0.0212)

(0.025) (0.0625) (0.0169) (0.0060) / hold=(2 to 19);

where groupe=1;

run;

-2 Res Log Pseudo-Likelihood 29.74

Generalized Chi-Square 16.87

Gener. Chi-Square / DF 0.99

Covariance Parameter Estimates

Standard

Cov Parm Subject Group Estimate Error

Intercept study 0.2244 0.1308

Residual (VC) study study 1 0.1250 .

Residual (VC) study study 2 0.3330 .

Residual (VC) study study 3 0.03100 .

Residual (VC) study study 4 0.2000 .

Residual (VC) study study 5 0.06250 .

Residual (VC) study study 6 0.01690 .

Residual (VC) study study 7 0.005000 .

Residual (VC) study study 8 0.007000 .

Residual (VC) study study 9 0.05500 .

Residual (VC) study study 10 0.04760 .

Residual (VC) study study 11 0.08300 .

Residual (VC) study study 12 0.03300 .

Residual (VC) study study 13 0.1000 .

Residual (VC) study study 14 0.06250 .

Residual (VC) study study 15 0.02120 .

Residual (VC) study study 16 0.02500 .

Residual (VC) study study 17 0.06250 .

Residual (VC) study study 21 0.01690 .

Scale 0.8399 1.8777

Solutions for Fixed Effects

Standard

Effect Estimate Error DF t Value Pr > |t|

Intercept -2.6600 0.1244 17 -21.38 <.0001

3 - Model with no random effect and type=chol(1) for matrix R

proc glimmix data=count_data;

class study;

model tox = / dist=nb link=log offset=logPA2 solution;

random intercept/subject=study group=study type=chol(1);

parms

(0.125) (0.333) (0.031) (0.2) (0.0625)

(0.0169) (0.005) (0.007) (0.055) (0.0476)

(0.083) (0.033) (0.1) (0.0625) (0.0212)

(0.025) (0.0625) (0.0169) (0.0060) / hold=(1 to 19);

where groupe=1;

run;

-2 Res Log Pseudo-Likelihood 99.51

Generalized Chi-Square 113.25

Gener. Chi-Square / DF 6.66

Covariance Parameter Estimates

Standard

Cov Parm Subject Group Estimate Error

CHOL(1,1) study study 1 0.1250 .

CHOL(1,1) study study 2 0.3330 .

CHOL(1,1) study study 3 0.03100 .

CHOL(1,1) study study 4 0.2000 .

CHOL(1,1) study study 5 0.06250 .

CHOL(1,1) study study 6 0.01690 .

CHOL(1,1) study study 7 0.005000 .

CHOL(1,1) study study 8 0.007000 .

CHOL(1,1) study study 9 0.05500 .

CHOL(1,1) study study 10 0.04760 .

CHOL(1,1) study study 11 0.08300 .

CHOL(1,1) study study 12 0.03300 .

CHOL(1,1) study study 13 0.1000 .

CHOL(1,1) study study 14 0.06250 .

CHOL(1,1) study study 15 0.02120 .

CHOL(1,1) study study 16 0.02500 .

CHOL(1,1) study study 17 0.06250 .

CHOL(1,1) study study 21 0.01690 .

Scale 0.006000 .

Solutions for Fixed Effects

Standard

Effect Estimate Error DF t Value Pr > |t|

Intercept -2.7003 0.04557 17 -59.26 <.0001

4 - random Model and type=chol(1) for R matrix

proc glimmix data=count_data;

class study;

model tox = / dist=nb link=log offset=logPA2 solution;

random intercept/subject=study;

random intercept/subject=study group=study type=chol(1);

parms (0)

(0.125) (0.333) (0.031) (0.2) (0.0625)

(0.0169) (0.005) (0.007) (0.055) (0.0476)

(0.083) (0.033) (0.1) (0.0625) (0.0212)

(0.025) (0.0625) (0.0169) (0.0060) / hold=(2 to 19);

where groupe=1;

run;

-2 Res Log Pseudo-Likelihood 29.73

Generalized Chi-Square 16.89

Gener. Chi-Square / DF 0.99

Covariance Parameter Estimates

Standard

Cov Parm Subject Group Estimate Error

Intercept study 0.1043 .

CHOL(1,1) study study 1 0.1250 .

CHOL(1,1) study study 2 0.3330 .

CHOL(1,1) study study 3 0.03100 .

CHOL(1,1) study study 4 0.2000 .

CHOL(1,1) study study 5 0.06250 .

CHOL(1,1) study study 6 0.01690 .

CHOL(1,1) study study 7 0.005000 .

CHOL(1,1) study study 8 0.007000 .

CHOL(1,1) study study 9 0.05500 .

CHOL(1,1) study study 10 0.04760 .

CHOL(1,1) study study 11 0.08300 .

CHOL(1,1) study study 12 0.03300 .

CHOL(1,1) study study 13 0.1000 .

CHOL(1,1) study study 14 0.06250 .

CHOL(1,1) study study 15 0.02120 .

CHOL(1,1) study study 16 0.02500 .

CHOL(1,1) study study 17 0.06250 .

CHOL(1,1) study study 21 0.01690 .

Scale 0.1103 0.09373

Solutions for Fixed Effects

Standard

Effect Estimate Error DF t Value Pr > |t|

Intercept -2.6504 0.1239 0 -21.39 .

The random model with weight statement (model 1) and the random model with type=simple (model 2) give similar results. Note: When I use a poisson model the scale parameter is not defined. I need to insert a random _residual_ while scale parameter estimation is by default with Negative binomial. The last 3 models give a Pearson/df close to 1 (control of overdispersion).

The model 3 indicates the need to adjust for overdispersion chiSquare/df=6.66 (>>1).But I realize I make a mistake because I wrongly fixed the parm hold=(1 to 19) leading to fix the scale parameter.

The model 4 (random effect type=chol(1)) give no SE estimation for random effect and DF=0 for fixed effect. How explain it ? Very strange. Thank per advance for you help.

Gwénaël

Meta analysis using Glimmix

lvm — Mon, 16 Jan 2012 00:18:32 GMT

There are several things to keep in mind with meta-analysis. First of all, with fixed or random effects meta-analysis, you need to fix the within-study variances. With normal distributions, this is straight-forward by first fixing the residual at 1 and then using weights. (This is easily accomplished, using either mixed or glimmix (see my other posts on a similar topic)). Or, one specifies a separate within-study (i.e., "residual") variance for each study and holds all of these fixed. You are trying this approach here, but there are some complications, so I don't think you are getting what you want. With a discrete distribution (Poisson or negative binomial (NB)), the within-study ("residual") variance is defined based on the distribution (=mean for Poisson, =(mean+(scale)*mean^2 for NB). Here, "scale" means the within-study scale parameter for the NB (equivalent to 1/k in other parameterizations). This residual scale term is on the scale of the raw data, not the link (but the other random statements refer to the scale of the link). You have the pre-specified within-study variances, but I don't believe they are being substituted for the above-described within-study variances, using your coding. Rather, I think with the first two approaches (your latest post) you are getting the internally-defined residual variance multiplied by the listed variance parameters (or the equivalent with the weights). This gets more complicated because, by using the NB distribution, the program is also estimating a NB scale parameter based on data across the studies (that is the scale parameter being displayed). Thus, it appears that none of your choices is really fixing your within-study variances at the values you want.

I have done a lot with meta-analysis, but either with normal data, or with non-normal data but without fixing the within-study variances (using inidiviudal data within the studies). I have definitely not tried this using the NB distribution (with its complication of estimating a separate scale across the studies). I don't have an exact suggestion right now, but will think about it more. Here is another complication. If you substitute a variance for "residual", then you no longer have a true Poisson or NB, since these are defined, in part, by the variance:mean relation.

Here is one approach to consider. Choose Poisson, and determine the "residual" variance based on the mean. Then re-scale your listed within-study variances, so that the product gets you back to the desired within-study variances. Then use your method 2 (last post). I need to think about this some more. By the way, your approaches 1 and 2 are closest to what you want (with the qualifiers that I list here).

Meta analysis using Glimmix

G_Le_Teuff — Mon, 16 Jan 2012 15:19:21 GMT

Thank. but I am not sure how to determine the "residual" variance based on the mean. Do I have to take the total number of toxicities (tox) or the incidence rate i.e total number of toxicities divided by person-years. Here A=mu for Poisson model.

This question because the weigth is based on the inverse-variance and an estimation of the variance is c/person-years**2 for incidence rate and 1/c for log(IR).

In the case of binomial data, I have to use the same approach ? with A=pi(1-pi).

Thank par advance

Meta analysis using Glimmix

G_Le_Teuff — Mon, 16 Jan 2012 16:15:44 GMT

To be sure to well understand your proposal and following my last post:

In a data set, I add 2 supplementary lines

mu=tox;

varlnbis=varln/mu;w_varlnbis=w_varln/mu;

When I used mu=log(tox) I obtain the following results

* weighted approach -2.6739(0.1286) random=0.2824(0.1030)

* fixed the new variance -2.6778(0.1302) random=0.2980(0.1050)

and with mu=tox the following results

* weighted approach -2.6765(0.1275) random=0.2352(0.099)

* fixed the new variance -2.6778(0.1307) random=0.3047(0.1056)

Meta analysis using Glimmix

lvm — Mon, 16 Jan 2012 19:24:39 GMT

You are on the right track based on my previous comments. The difference is still probably due to the estimated scale for the chosen distribution. But the more I think about it, I don't think that Poisson or NB are appropriate for your analysis. Even if the original data in each study are discrete, either Poisson or NB, you are doing the meta-analysis on an estimated parameter from each study. This parameter estimate would be continuous, and probably approximated described by a normal distribution. Your within-study variance (or its inverse, the weight) is based on large-sample (normal) theory. Thus, I think you could avoid a lot of confusion by directly analyzing the parameter estimate for the studies using normal-based approaches.

You can continue with use of GLIMMIX, but base the analysis on LOGtox. The weight should be the inverse of the within-study variance for the logtox for each study. You still need a parms statement, so that the residual is held fixed at 1. The combination of the weight and the parms statements gives you fixed within-study variances and an estimated between-study variance. The output says that the residual is 1, but because of the weights, the residual is actually different for each study (=1/weight).

proc glimmix data=count_data;

class study;

model LOGtox = / offset=logPA2 solution;

random intercept / subject=study ; * to test heterogeneity between study;

weight inv_varln; * weight = inverse variance;

parms (1) (1) / hold=2;

run;

You can read more at:

http://www.lexjansen.com/pharmasug/2000/stats/st09.pdf

(this was written before GLIMMIX was available, but it is very useful). There are also other ways of using a long parms statement and no weight.

Meta analysis using Glimmix

G_Le_Teuff — Tue, 17 Jan 2012 09:47:14 GMT

Thank you for this information but I have some difficulties to choose what variance used to compute the weight. I am only able to calculate the variance of incidence rate or log(incidence rate) but I do not know which used.

For a poisson regression, the within-study variance is defined as mean but on which scale : incidence rate or log(incidence rate) ? Could you indicate me on which scale the variance of a poisson model is calculated. The model is log(tox/PA^2) = intercept + random effect. Variance=mean but mean is what ? log(tox/PA^2) or tox/PA^2 or tox or log(tox). Sorry if my questions are stupid.

For a normal regression (y=log(tox)), do I calculate the variance of log(tox/PA2) ?

After that, I will be able to compare the 2 approches (poisson and normal).

Thank you per advance

Meta analysis using Glimmix

G_Le_Teuff — Tue, 17 Jan 2012 10:47:29 GMT

Please below the results of different approaches

Fixed Random -2ResLL

Poisson using weight = inverse-variance of incidence rate -2.6779 (0.1309) 0.3078(0.1058) 31.13

Poisson using weight = inverse-variance of log incidence rate -2.6739(0.1286) 0.2824(0.1030) 30.78

Normal using weight=inverse-variance of log incidence rate* -2.7431(0.1395) 0.3131(0.1197) 36.05

* because the response is y=log(tox)

The -2ResLL is different between the first 2 models and normal model. Poisson model seem better fit than Normal model. Based on these results, could we select the poisson model ?

Additionnaly, can you confirm that the random effect as estimated in log scale. THANK

Meta analysis using Glimmix

lvm — Tue, 17 Jan 2012 15:43:33 GMT

It looks like you are on the right track. The results are not expected to be the same. You can't compare log-likelihoods because pseudo-likelihood is used for the Poisson distributions but true restricted likelihoods are used for the normal distribution.

Meta analysis using Glimmix

G_Le_Teuff — Tue, 17 Jan 2012 15:58:18 GMT

But the results can be more different : I applied these 3 approaches on other data (specific toxicities and not total toxicity). So, no simple solution to this problem. Yes, it was what I am thinking about the comparison of Poisson and Normal distribution.

I am wondering whether you will be interest to collaborate on this topic ?

Meta analysis using Glimmix

G_Le_Teuff — Tue, 17 Jan 2012 16:26:38 GMT

For information, please find enclosed a reference

Mixed-effects Poisson regression models for meta-analysis of follow-up studies with constant or varying durations. The international journal of Biostatistics. Vol 5, Issue 1, 2009. P.G .Bagos, G.K. Nikolopoulos. Examples are given but using Strata.

Meta analysis using Glimmix

G_Le_Teuff — Wed, 18 Jan 2012 12:44:45 GMT

Dear Lvm,

I read a second time the paper you cited above (Meta-analysis with linear and Nonlinear multilevel models using proc mixed and proc nlmixed). The authors used a proc mixed based on log(OR) with weight statement and compare the results using a proc nlmixed with a random effect. In the SAS code of proc nlmixed, I did not seen how they fixed the within-study variance. Very strange ! because the variance function, here, is mu(1-mu)/n for a proportion. Have you some idea about it !

Meta analysis using Glimmix

lvm — Wed, 18 Jan 2012 14:20:32 GMT

I think your last posting gets at the critical point. Using the NLMIXED code (or equivalent now in GLIMMIX for this example), the within-study variance is automatically determined based on the specified distribution. One does not fix the within-study variance. If the data have a Poisson distribution, the variance is determined by the mean (by definition); i.e., the variance is completely known, and there is no reason to fix it in each study. The software takes care of this. This is why the approach you are taking is difficult -- you are specifying a discrete distribution (NB or Poisson), but you want to fix the within-study variances at something other than what is defined for those discrete distributions. You are essentially getting a overdispersed distribution scenario, where the degree of overdispersion (or possibly underdispersion) varies with the study. I attemped one workaround, where you take into account the internal variance structure (determined by the software), and then rescaled it to achieve YOUR fixed variances. However, I think you would be better off using the normal approximations, because this is easier to control the within-study variances.

Meta analysis using Glimmix

G_Le_Teuff — Wed, 18 Jan 2012 18:04:43 GMT

Thank. It's clearer. But based on the paper for log ddd ratios (normal data) and binomial data, is Poisson with internal residual variance determined by software (wihout weight statement) a correct and reasonable approach for incidence rate and gives results similar to the normal distribution with variance of log incidence.

Meta analysis using Glimmix

G_Le_Teuff — Wed, 18 Jan 2012 18:25:17 GMT

I am thinking that my confusion come from the 2 following responses : y= (#event / Person-year) and y=#event. In your last post, you propose the model

proc glimmix data=count_data;

class study;

model LOGtox = / offset=logPA2 solution;

random intercept / subject=study ; * to test heterogeneity between study;

weight inv_varln; * weight = inverse variance;

parms (1) (1) / hold=2;

run;

But I have an approximation of the variance of tox/person-year or log(tox/person-year). So I have to peform a model where the response is log(tox/person-year) and delete the offset.

proc glimmix data=count_data;

class study;

model LOG(tox/person-years) = / solution;

random intercept / subject=study ; * to test heterogeneity between study;

weight inv_varln; * weight = inverse variance;

parms (1) (1) / hold=2;

run;

In my last posts I was wondering whether I have to consider the denominator of the response (offset) in the function variance determined by the software or whether the software takes into account only the numerator of incidence rate i.e #events.