Meta analysis using Glimmix

G_Le_Teuff · Posted 01-12-2012 12:04 PM

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

SteveDenham · Posted 01-13-2012 07:42 AM

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

G_Le_Teuff · Posted 01-13-2012 09:27 AM

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

SteveDenham · Posted 01-13-2012 10:19 AM

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

SteveDenham · Posted 01-13-2012 10:24 AM

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

G_Le_Teuff · Posted 01-14-2012 06:24 AM

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

lvm · Posted 01-15-2012 07:18 PM

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).

G_Le_Teuff · Posted 01-16-2012 10:19 AM

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

G_Le_Teuff · Posted 01-16-2012 11:15 AM

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)

lvm · Posted 01-16-2012 02:24 PM

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.

G_Le_Teuff · Posted 01-17-2012 04:47 AM

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

G_Le_Teuff · Posted 01-17-2012 05:47 AM

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

lvm · Posted 01-17-2012 10:43 AM

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.

G_Le_Teuff · Posted 01-17-2012 10:58 AM

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 ?

G_Le_Teuff · Posted 01-17-2012 11:26 AM

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.

SAS Innovate 2025: Save the Date