Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Numerical integration with Mata gives different answer to R

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

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