### Generalized Latent Variable Modeling by Skrondal and Rabe-Hesketh Section 9.5: Nicotine gum and smoking cessation Meta-analysis

The data set used in this section is gum.dat from Silagy (2003).

The do-file is gum.do.

The programs used are gllamm and gllapred. Use the command ssc describe gllamm and follow instructions to download these GLLAMM commands. For more information on GLLAMM see http://www.gllamm.org.

```
insheet using gum.dat, clear
rename study studynam
gen study=_n
list, clean

studynam    d1    n1    d0    n0   study
1.      Blondal 1989    37    92    24    90       1
2.     Campbell 1991    21   107    21   105       2
3.   Fagerstrom 1982    30    50    23    50       3
4.          Fee 1982    23   180    15   172       4
5.       Garcia 1989    21    68     5    38       5
6.       Garvey 2000    75   405    17   203       6
7.        Gross 1995    37   131     6    46       7
8.         Hall 1985    18    41    10    36       8
9.         Hall 1987    30    71    14    68       9
10.         Hall 1996    24    98    28   103      10
11.   Hjalmarson 1984    31   106    16   100      11
12.        Huber 1988    31    54    11    60      12
13.       Jarvis 1982    22    58     9    58      13
14.       Jensen 1991    90   211    28    82      14
15.       Killen 1984    16    44     6    20      15
16.       Killen 1990   129   600   112   617      16
17.      Malcolm 1980     6    73     3   121      17
18.     McGovern 1992    51   146    40   127      18
19.     Nakamura 1990    13    30     5    30      19
20.       Niaura 1994     5    84     4    89      20
21.       Niaura 1999     1    31     2    31      21
22.        Pirie 1992    75   206    50   211      22
23.        Puska 1979    29   116    21   113      23
24.    Schneider 1985     9    30     6    30      24
25.     Tonnesen 1988    23    60    12    53      25
26.        Villa 1999    11    21    10    26      26
27.       Zelman 1992    23    58    18    58      27
```

Prepare data for analysis, creating variables
gum=1 if in treatment arm, 0 otherwise
quit=1 if quit smoking, 0 otherwise
wt1: frequency weight
```
reshape long d n, i(study) j(gum)

(note: j = 0 1)

Data                               wide   ->   long
-----------------------------------------------------------------------------
Number of obs.                       27   ->      54
Number of variables                   6   ->       5
j variable (2 values)                     ->   gum
xij variables:
d0 d1   ->   d
n0 n1   ->   n
-----------------------------------------------------------------------------

gen y0 = n-d
rename d y1
gen i=_n
reshape long y, i(i) j(quit)

(note: j = 0 1)

Data                               wide   ->   long
-----------------------------------------------------------------------------
Number of obs.                       54   ->     108
Number of variables                   7   ->       7
j variable (2 values)                     ->   quit
xij variables:
y0 y1   ->   y
-----------------------------------------------------------------------------

drop i
rename y wt1

sort study gum quit
list in 1/12, clean

quit   study   gum          studynam   wt1     n
1.      0       1     0      Blondal 1989    66    90
2.      1       1     0      Blondal 1989    24    90
3.      0       1     1      Blondal 1989    55    92
4.      1       1     1      Blondal 1989    37    92
5.      0       2     0     Campbell 1991    84   105
6.      1       2     0     Campbell 1991    21   105
7.      0       2     1     Campbell 1991    86   107
8.      1       2     1     Campbell 1991    21   107
9.      0       3     0   Fagerstrom 1982    27    50
10.      1       3     0   Fagerstrom 1982    23    50
11.      0       3     1   Fagerstrom 1982    20    50
12.      1       3     1   Fagerstrom 1982    30    50
```

Continuous random intercept and slope model for Table 9.9 using gllamm (iteration log not shown):
```
gen treat = gum -0.5
gen cons=1
eq int: cons
eq slop: treat

gllamm quit treat, i(study) nrf(2) eqs(int slop) l(logit) /*
*/ f(binom) weight(wt) adapt nocor trace allc

number of level 1 units = 5908
number of level 2 units = 27

Condition Number = 1.7685027

gllamm model

log likelihood = -3074.1459

------------------------------------------------------------------------------
quit |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
treat |   .5670943   .0881729     6.43   0.000     .3942786      .73991
_cons |    -1.1658   .1401855    -8.32   0.000    -1.440559   -.8910418
------------------------------------------------------------------------------

Variances and covariances of random effects
------------------------------------------------------------------------------

***level 2 (study)

var(1): .48050074 (.15333759)
cov(2,1): fixed at 0

var(2): .05190697 (.04639595)
------------------------------------------------------------------------------

------------------------------------------------------------------------------
quit |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
quit         |
treat |   .5670943   .0881729     6.43   0.000     .3942786      .73991
_cons |    -1.1658   .1401855    -8.32   0.000    -1.440559   -.8910418
-------------+----------------------------------------------------------------
stu1_1       |
cons |   .6931816   .1106042     6.27   0.000     .4764014    .9099619
-------------+----------------------------------------------------------------
stu1_2       |
treat |    .227831    .101821     2.24   0.025     .0282656    .4273964
------------------------------------------------------------------------------
```

Slight discrepancy for the treatment effect standard deviation disappears when nip(20) is used.

Empirical Bayes predictions of study-specific effect sizes using gllapred:
```
gllapred eff, u

(means and standard deviations will be stored in effm1 effs1 effm2 effs2)
-3101.5793 -3092.8436 -3088.3938 -3084.9338 -3081.8962 -3079.2417
-3077.5891 -3076.3852 -3074.9084 -3074.1634 -3074.1459 -3074.1459
log-likelihood:-3074.1459
```

95% intervals
```
gen eff = effm2 + _b[treat]
gen lower = effm2 + _b[treat] - 1.96*effs2
gen upper = effm2 + _b[treat] + 1.96*effs2

```

List results:
```
sort study
qui by study: gen last=_n==_N
list studynam lower eff upper if last==1, clean

studynam       lower        eff      upper
4.      Blondal 1989     .222363   .5866519   .9509408
8.     Campbell 1991    .0146321   .3882984   .7619646
12.   Fagerstrom 1982    .1783771   .5678892   .9574013
16.          Fee 1982    .1462535   .5193934   .8925332
20.       Garcia 1989    .2424062   .6489648   1.055523
24.       Garvey 2000    .3575259   .6972256   1.036925
28.        Gross 1995    .2477593   .6437063   1.039653
32.         Hall 1985      .19287   .5980018   1.003134
36.         Hall 1987    .3139485   .6961169   1.078285
40.         Hall 1996   -.0332232   .3325011   .6982254
44.   Hjalmarson 1984     .259717   .6314499   1.003183
48.        Huber 1988    .4692236   .8577137   1.246204
52.       Jarvis 1982    .3085245   .7037886   1.099053
56.       Jensen 1991    .1502759   .4940711   .8378663
60.       Killen 1984    .1223252   .5393925   .9564598
64.       Killen 1990    .0739083   .3133026   .5526969
68.      Malcolm 1980    .2148292   .6320365   1.049244
72.     McGovern 1992    .0535688   .3903244   .7270799
76.     Nakamura 1990    .2542217   .6702684   1.086315
80.       Niaura 1994    .1086343   .5270553   .9454764
84.       Niaura 1999    .0755295   .5087444   .9419592
88.        Pirie 1992    .2841303   .5926718   .9012133
92.        Puska 1979     .139473   .5054876    .871502
96.    Schneider 1985    .1465653   .5650316    .983498
100.     Tonnesen 1988    .2204977   .6133725   1.006247
104.        Villa 1999    .1476256   .5656016   .9835777
108.       Zelman 1992    .1351462    .522485   .9098238
```

Plot results using forest plot (similar to Figure 9.6):
```
encode studynam, gen(s)
twoway (rcap lower upper s, horizontal blpattern(dash)) (scatter s eff),   /*
*/ ylabel(1(1)27,valuelabels angle(horizontal) labs(small)) legend(off)  /*
*/ xscale(range(-.5 1.5)) ysize(7) xlab(-.5(.5)1.5) xline(0) saving(gumf)

```

Nonoarametric maximum likelihood (NPML) estimates using gateaux derivative:

```
gllamm quit treat, i(study) nrf(2) eqs(int slop) l(logit) /*
*/ f(binom) weight(wt) ip(f) nip(2)

number of level 1 units = 5908
number of level 2 units = 27

Condition Number = 12.81975

gllamm model

log likelihood = -3106.6101

------------------------------------------------------------------------------
quit |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
treat |   .5317816   .0666912     7.97   0.000     .4010692     .662494
_cons |  -1.086669   .0919203   -11.82   0.000    -1.266829   -.9065085
------------------------------------------------------------------------------

Probabilities and locations of random effects
------------------------------------------------------------------------------

***level 2 (study)

loc1: -.53793, .35823
var(1): .19270432

loc2: -.09286, .06184
cov(2,1): .03326538
var(2): .0057424
prob: 0.3997, 0.6003
------------------------------------------------------------------------------
```

Add one mass at a time using gateaux derivative method (output not shown):
```
local n=2
while \$HG_error==0&`n'<12{
matrix a = e(b)
local ll = e(ll)
local k = e(k)
local n = `n' + 1
gllamm quit treat , i(study) nrf(2) eqs(int slop) l(logit) /*
*/ f(binom) weight(wt) ip(f) nip(`n') lf0(`k' `ll') /*
*/ gateaux(-5 5 40) from(a) trace
}

```

Final estimates for 10 masses:
```number of level 1 units = 5908
number of level 2 units = 27

Condition Number = 111.7378

gllamm model                                      Number of obs   =       5908
LR chi2(3)      =       0.36
Log likelihood = -3061.5044                       Prob > chi2     =     0.9488

------------------------------------------------------------------------------
quit |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
treat |   .5909019   .1020571     5.79   0.000     .3908736    .7909301
_cons |  -1.197405   .1530979    -7.82   0.000    -1.497471   -.8973385
------------------------------------------------------------------------------

Probabilities and locations of random effects
------------------------------------------------------------------------------

***level 2 (study)

loc1: -.17192, -1.7707, -.74877, 1.2864, .58648, -.05122, .11831, -.95241, .57538, .39284
var(1): .57783952

loc2: -.37184, -.0051, .28878, -.02044, 1.1465, .16963, -.63535, -.16783, -.13978, .05351
cov(2,1): .01816208
var(2): .09297283
prob: 0.1216, 0.106, 0.0442, 0.0461, 0.0382, 0.1482, 0.0325, 0.0355, 0.1617, 0.2662
------------------------------------------------------------------------------
```

Store locations and log odds as variables:
```
matrix a=e(zlc2)'
svmat a
matrix b=e(zlc3)'
svmat b
matrix prob = e(zps2)'
svmat prob

```

Convert log odds to probabilities and calculate corresponding rounded percentages:
```
replace prob1=exp(prob1)
gen percent=round(prob1*100,1)

```

Produce a plot similar to the upper panel of Figure 9.7:
```
twoway (scatter b1 a1 [w=percent], msymbol(Oh))   /*
*/ (scatter b1 a1, mlabel(percent) mlabpos(0) msymbol(none)), /*
*/ xtitle(Intercept) ytitle(Slope) legend(off) saving(blobs, replace)

```

References

Silagy, C. (2003). Nicotine replacement therapy for smoking cessation (Cochrane review). The Cochrane Library, Issue 4. Chichester: Wiley.

Rabe-Hesketh, S., Skrondal, A. and Pickles, A. (2002). Reliable estimation of generalised linear mixed models using adaptive quadrature. The Stata Journal 2, 1-21.

Skrondal, A. and Rabe-Hesketh, S. (2004). Generalized Latent Variable Modeling: Multilevel, Longitudinal and Structural Equation Models. Boca Raton, FL: Chapman & Hall/ CRC Press.

Outline
Datasets and do-files