Dear all,
I am attempting to code the random-effects meta-analysis model of Henmi and Copas (Stat Med 2010; DOI: 10.1002/sim.4029). They give code in R, which (after a couple of typo corrections, which I am aware of) has since been incorporated into the -metafor- package for R. I have tested this code and am confident that it works. In Mata, however, I have become stumped by the evaluation of a particular integral. Mata (using Adrian Mander's integrate() function) is giving me different answers to R, and with respect to a graph of the integrand function.
Here is the underlying example data, using -dataex- in Stata (taken from Lee and Done, Cochrane Library 2004; DOI: 10.1002/14651858.CD003281.pub2):
...and then some setup, in Mata:
Next, here is my integrand function, and a graph of it. It rises from near the origin to a peak at x=~2, at which y=~0.035, and then drops back to near zero:
Now to the central issue. Evaluating integrand() at particular values (say 1, 2 or 3) agrees with R. But integrating it does not, e.g.:
...gives ~0.016, whereas R gives ~0.042.
Worse still, ultimately I need to integrate from some value x up to infinity, but specifying an upper limit of infinity seems to make it completely unstable. By contrast, R very quickly gives the answer from, say, 0 to infinity as ~0.046.
What am I doing wrong here?
Many thanks,
David.
I am attempting to code the random-effects meta-analysis model of Henmi and Copas (Stat Med 2010; DOI: 10.1002/sim.4029). They give code in R, which (after a couple of typo corrections, which I am aware of) has since been incorporated into the -metafor- package for R. I have tested this code and am confident that it works. In Mata, however, I have become stumped by the evaluation of a particular integral. Mata (using Adrian Mander's integrate() function) is giving me different answers to R, and with respect to a graph of the integrand function.
Here is the underlying example data, using -dataex- in Stata (taken from Lee and Done, Cochrane Library 2004; DOI: 10.1002/14651858.CD003281.pub2):
Code:
* Example generated by -dataex-. To install: ssc install dataex clear input byte study float(eff se) 1 -.130 .361 2 -1.622 .556 3 .418 .649 4 .079 .287 5 -.179 .600 6 -.241 .696 7 -.381 .264 8 -1.912 .734 9 0.000 1.438 10 -1.093 .536 11 -1.302 .525 12 -.453 .681 13 -3.099 1.082 14 -1.660 .479 15 -1.169 .373 16 .161 .315 end label values study study label def study 1 "Agarwal (2000)", modify label def study 2 "Agarwal (2002)", modify label def study 3 "Alkaissi (1999)", modify label def study 4 "Alkaissi (2002)", modify label def study 5 "Allen (1994)", modify label def study 6 "Andrzejowski (1996)", modify label def study 7 "Duggal (1998)", modify label def study 8 "Dundee (1986)", modify label def study 9 "Ferrara-Love (1996)", modify label def study 10 "Gieron (1993)", modify label def study 11 "Harmon (1999)", modify label def study 12 "Harmon (2000)", modify label def study 13 "Ho (1996)", modify label def study 14 "Rusy (2002)", modify label def study 15 "Wang (2002)", modify label def study 16 "Zarate (2001)", modify
...and then some setup, in Mata:
Code:
mata:
st_view(yi=., ., "eff", .)
st_view(si=., ., "se", .)
vi = si:^2
k = length(yi)
wi = 1:/vi
eff = mean(yi, wi)
Q = crossdev(yi, eff, wi, yi, eff)
W1 = sum(wi)
W2 = sum(wi:^2)/W1
W3 = sum(wi:^3)/W1
W4 = sum(wi:^4)/W1
tausq = max((0, (Q - (k-1))/(W1 - W2))) // truncated D+L
sd_eff = sqrt((tausq*W2 + 1)/W1)
VR = 1 + tausq*W2
SDR = sqrt(VR)
real scalar EQ(real scalar r) {
external real scalar k, W1, W2, W3, tausq, VR
real scalar ans
ans = (k - 1) + tausq*(W1 - W2) + (tausq^2)*((1/VR^2)*(r^2) - 1/VR)*(W3 - W2^2)
return(ans)
} // conditional mean of Q given R=r
real scalar VQ(real scalar r) {
external real scalar k, W1, W2, W3, W4, tausq, VR
real scalar ans
ans = 2*(k - 1) + 4*tausq*(W1 - W2) + 2*(tausq^2)*(W1*W2 - 2*W3 + W2^2)
ans = ans + 4*(tausq^2)*((1/VR^2)*(r^2) - 1/VR)*(W3 - W2^2)
ans = ans + 4*(tausq^3)*((1/VR^2)*(r^2) - 1/VR)*(W4 - 2*W2*W3 + W2^3)
ans = ans + 2*(tausq^4)*((1/VR^2) - 2*(1/VR^3)*r^2)*(W3 - W2^2)^2
return(ans)
} // conditional variance of Q given R=r
real scalar finv(real scalar f) {
external real scalar k, W1, W2
real scalar ans
ans = (W1/W2 - 1)*((f^2)-1) + (k-1)
return(ans)
} // inverse function of f
Next, here is my integrand function, and a graph of it. It rises from near the origin to a peak at x=~2, at which y=~0.035, and then drops back to near zero:
Code:
mata:
real rowvector integrand(real rowvector x) {
real scalar ans
ans = gammap((EQ(x[1])^2)/VQ(x[1]), finv(x[1])/(VQ(x[1])/EQ(x[1]))) * normalden(x[1])
return(ans)
}
X=(0..1000)/100
Y=X
for(i=1; i<=length(X); i++) {
Y[i] = integrand(X[i])
}
plot=Y',X'
mm_plot(plot)
end
Now to the central issue. Evaluating integrand() at particular values (say 1, 2 or 3) agrees with R. But integrating it does not, e.g.:
Code:
mata: integrate(&integrand(), 1, 3, 500)
Worse still, ultimately I need to integrate from some value x up to infinity, but specifying an upper limit of infinity seems to make it completely unstable. By contrast, R very quickly gives the answer from, say, 0 to infinity as ~0.046.
What am I doing wrong here?
Many thanks,
David.
