Store r-class values and tabulate later

[Anselm Schulz]

Hello all,

I am new to stata and my work is about medical data.
My data are about diagnostic performance of different imaging modalities (I am using random data in the example). I have installed a stata extension called diagnostic test (the command is “diagt”, first variable is reference standard second the test, both are binary). In my analysis diagt is executed many times and I want to save/store the r-class data/values provided by the test e.g. sensitivity and specificity in order to summarize them later (at the end of the “do” file) in a table. I guess that diagt is not important and that it would be the same for all commands providing r-class data.

At the moment I store the values in macros and put the macros at the end in a matrix which I display with “matlist”.

My questions:
Is there an easier/more elegant way to store the r-class values and show them in a table afterwards?
It seems not optimal that I must specify each value I want to use later and put it into a macro. Is it possible to store all the r-class data at once and display specific values later?

// run diagt
diagt fasitb MRtb

// save some r-class data
global MRt_sens = r(sens)
global MRt_spec = r(spec)

// run diagt for another modality/test
diagt fasitb MRnb

// save the same r-class data into different macros
global MRn_sens = r(sens)
global MRn_spec = r(spec)

// now I generate my matrix and the summary table

matrix A = $MRn_sens, $MRn_spec\ $MRt_sens, $MRt_spec \ $MCC_overall, . \ $MCC_sens, $MCC_spec
matrix rownames A = "MR nativ" "MR toilax" "McNemar Overall, Nativ vs Toilax" "McNemar sens spec"
matrix colnames A = sensitivity specificity
matlist A, format(%15.1f) twidth(40) title(Summary of Results)

My table looks like this:

Here a list over all r-class data provided by diagt and I will probaly use several of them.

return list

scalars:
               r(sens) =  50.12003200853561
            r(sens_lb) =  48.50661567935583
            r(sens_ub) =  51.73326103478671
               r(spec) =  49.96000639897616
            r(spec_lb) =  48.71277950932287
            r(spec_ub) =  51.20727063140053
               r(prev) =  37.49
            r(prev_lb) =  36.53994481639038
            r(prev_ub) =  38.44734619397738
                r(ppv) =  37.52746155382465
             r(ppv_lb) =  36.18391913702044
             r(ppv_ub) =  38.8855612589619
                r(npv) =  62.54756659323052
             r(npv_lb) =  61.1880343889589
             r(npv_ub) =  63.89241263219665
                r(lrp) =  1.001599488763926
             r(lrp_lb) =  .9619277032697625
             r(lrp_ub) =  1.042907416515918
                r(lrn) =  .9983979504151261
             r(lrn_lb) =  .9587131869858447
             r(lrn_ub) =  1.039725416239468
                 r(or) =  1.003206675602117
              r(or_lb) =  .9251808504343817
              r(or_ub) =  1.087812866382799
                r(roc) =  .500400185585022
             r(roc_lb) =  .490277306910667
             r(roc_ub) =  .510523064259377
              r(oprev) =  .
           r(oprev_lb) =  .
           r(oprev_ub) =  .
               r(opos) =  .
            r(opos_lb) =  .
            r(opos_ub) =  .
               r(oneg) =  .
            r(oneg_lb) =  .
            r(oneg_ub) =  .

Any help would be appreciated very much!

[Mike Lacy]

I'd say you're doing pretty good as a new Stata user to do what you have done so far:

Anyway, because -diagt- returns so much material in r(), I'd suspect that any way of handling these results will be fairly tedious. I hope someone else has an easier solution, but the best Stata tool I can think is the built-in suite of -postfile- commands, which allows you to store a set of values into a Stata data set. Below, I've provided a sketch of a way to save all of the r-class results returned by a series of executions of -diagt-. There are likely some small mistakes in what I show, as I don't have any good material on which to test my syntax, but it should be pretty close to being right:

// Assemble long list of r-class result names into a macro
local rlist "sens  sens_lb  sens_ub  spec  spec_lb  spec_ub prev"
local rlist "`rlist' prev_lb  prev_ub  ppv  ppv_lb  ppv_ub npv"
local rlist "`rlist' npv_lb  npv_ub  lrp  lrp_lb  lrp_ub  lrn"
local rlist "`rlist' lrn_ub  or  or_lb  or_ub  roc  roc_lb roc_ub"
local rlist "`rlist' oprev  oprev_lb  oprev_ub  opos  opos_lb"
local rlist "`rlist' opos_ub  oneg  oneg_lb  oneg_ub"
// Create parenthesized list of items to be stored in the format that -post- requires
local postlist ""
foreach item of local rlist {
  local postlist "`postlist' (r(`item'))"
}
// Set up for using  -postfile-
postutil clear
tempname memhold
tempfile results
postfile `memhold' `rlist' using "`results'"
// Loop over all your calls to -diagt-  and store each set of r() results
// into a temporary data file called `results'
For all of your calls to AllYourCallsTo_diagt {
   diagt .....
   post `memhold' `postlist'
}   
postclose `memhold'
clear
use `results'

I have always found the documentation on -postfile- a bit thin, so searching online for material on it is helpful. Here's something I just found:
https://datatoday.blogspot.com/2011/...-in-stata.html

[daniel klein]

Concerning the approach in #1, I would create the matrix first, then fill in the values. You do not need to use globals (which you should always avoid). Something along the lines

// create the results matrix
matrix A = J(4, 3, .z)

// now fill the matrix
local row = 1
foreach y in MRtb MRn {
    diagt fastib `y'
    matrix A[`row', 1] = r(sens)
    matrix A[`row', 2] = r(spec)
    local ++row
}

Alternatively, you can use postfile, as Mike has pointed out.

Concerning the more specific question

Is it possible to store all the r-class data at once and display specific values later?

Yes. There are extended functions that return the names of r-class results. Here is an example

. sysuse auto , clear
(1978 Automobile Data)

. summarize price

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       price |        74    6165.257    2949.496       3291      15906

. return list

scalars:
                  r(N) =  74
              r(sum_w) =  74
               r(mean) =  6165.256756756757
                r(Var) =  8699525.974268789
                 r(sd) =  2949.495884768919
                r(min) =  3291
                r(max) =  15906
                r(sum) =  456229

. local r_scalars : r(scalars)

. display "`r_scalars'"
N sum_w mean Var sd min max sum

I trust you can figure out how to sue this approach in loop.

[Mike Lacy]

r(scalars), I was looking for something like that! <grin>

[Anselm Schulz]

Great! Thank you so much for your advices and fast response!

However, I am struggling using the r_scalar macro. It seems that it nicely creates a "sub macro/variable" (sorry I do not know how to call them) for each r-class - but how do I use/display the specific values e.g. for "min" only?

Thanks in advance for your support!

PS why is using global macros problematic - I naively assumed it was a good idea since they are everywhere available?

[daniel klein]

However, I am struggling using the r_scalar macro. It seems that it nicely creates a "sub macro/variable" (sorry I do not know how to call them) for each r-class - but how do I use/display the specific values e.g. for "min" only?

That part of my answer was referring more to a general concept, not to your specific problem. The rest of this post will elaborate a bit on these basic concepts and might not be relevant to solving your specific problem.

The local macro r_scalars merely holds the names of all scalars that are currently stored in r(). There is nothing special about the name r_scalars; I choose that name because I feel it describes the locals' contents. I should have named the local names_of_r_scalars. To store the contents of those scalars for later use, you will have to loop through the names. Here is an illustration:

// get the names of all scalars in r()
local names_of_r_scalars : r(scalars)

// now store the contents
foreach name of local names_of_r_scalars {
    local `name' = r(`name')
}

After the loop has finished, all scalars that have been returned in r() are now stored in local macros. Using my example with summarize, the local macro min would not hold the value that has been returned in r(min). Here is how this looks like:

. sysuse auto , clear
(1978 Automobile Data)

. summarize price

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       price |        74    6165.257    2949.496       3291      15906

. return list

scalars:
                  r(N) =  74
              r(sum_w) =  74
               r(mean) =  6165.256756756757
                r(Var) =  8699525.974268789
                 r(sd) =  2949.495884768919
                r(min) =  3291
                r(max) =  15906
                r(sum) =  456229

. local names_of_r_scalars : r(scalars)

. foreach name of local names_of_r_scalars {
  2.     local `name' = r(`name')
  3. }

. display "`min'"
3291

If you wanted this for more than one variable, you could do something like this:

sysuse auto , clear

foreach var in price mpg {
    
    summarize `var'
    
    local names_of_r_scalars : r(scalars)
    foreach name of local names_of_r_scalars {
        local `var'_`name' = r(`name')
    }
}

display "`price_min'"
display "`mpg_min'"

The relevant output is

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       price |        74    6165.257    2949.496       3291      15906

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
         mpg |        74     21.2973    5.785503         12         41

. display "`price_min'"
3291

. display "`mpg_min'"
12

This approach will fail if the created names for local macros are longer than 32 (actually, the limit for local macros is 31) characters. For example, you could not store r(minimum) of variable abcdefghijklmnopqrstuvwxyz. I will not show how to do that.

PS why is using global macros problematic - I naively assumed it was a good idea since they are everywhere available?

That is precisely the reason why you should not use global macros: it is too easy to lose track of the stuff you create in one place that will then mess up stuff in another. Even if you manage to keep track of the global macros that you define, you usually do not know about the global macros that others, including StataCorp define in their code. Here is one example where defining a global goes wrong

sysuse auto
mi set wide
mi register imputed rep78
global MI_IMPUTE_obj "this is my global"
display "$MI_IMPUTE_obj"
mi impute mvn rep78 = price , add(1)
display "$MI_IMPUTE_obj"

The relevant output is this

. display "$MI_IMPUTE_obj"
this is my global

. mi impute mvn rep78 = price , add(1)

[...]

. display "$MI_IMPUTE_obj"

(note blank line at the end).

[Anselm Schulz]

Wow and thank you very much for your very extensive and detailed response That really helps a lot!!!

[Anselm Schulz]

Dear Daniel,

I have now intensely tested your first approach - which I believe would be best suited for my task (so r(scalar) will also work and perhaps I have to use this approach anyway).
The problem that occurs is probably related to the command "diagt" which is no original stata command. It seems that "diagt" deletes somehow all matrices and some macros.
Your approach works nicely with e.g. "sum" command.

I found also this old post by Nick Cox using a similar approach, sorry that I did not see it earlier, but I did not know what to look for...
https://www.stata.com/statalist/arch.../msg01216.html

Here is my code

sysuse auto
// here I generate binary varibles for later use with "diagt", "price_b_RF" will be the reference standard. The output is of course nonsense...
recode price (0/5006.5=0) (5006.5/16000=1), gen(price_b_RF)
recode mpg (0/20=0) (20/16000=1), gen(mpg_b)
recode headroom (0/3=0) (3/5=1), gen(headroom_b)

matrix mean_sd = J(4, 2, .)
local test 1

local row = 1
foreach var in price mpg {
    sum `var'
    matrix mean_sd[`row', 1] = r(mean)
    matrix mean_sd[`row', 2] = r(sd)
    //matrix rownames sens_spec[`row'] = "`var'"
    local ++row
}

matrix list mean_sd
// the matrix is displayed nicely :)

// Now I am trying "diagt"
matrix sens_spec = J(4, 2, .)
matrix dir
macro dir
// two matrices are available

local row = 1
foreach var in mpg_b headroom_b {
    diagt price_b_RF `var'
    // it seems as if the matrix rownames sens_spec was deleted/dropped
    // the do file stops with this error message:
    //matrix sens_spec not found
    //r(111);
    matrix sens_spec[`row', 1] = r(sens)
    matrix sens_spec[`row', 2] = r(spec)
    local ++row
}

matrix list sens_spec
matrix dir
macro dir
// it seems as if all matrices were deleted/dropped :(
// the local macro "_test" is still there, however "_row" has disappeared
// Why does this happen?

The do file stopps with these lines (I have described it in the do file as well):
matrix sens_spec not found
r(111);
Do you have any idea why this is happening and how to solve it?

Thanks in advance for your help and effords!!

[daniel klein]

Thank you for providing example code. Your code runs fine for me. I am using

which diagt
[...]
*! diagt 2.032, 30 June 2003
*! by PT Seed ([email protected])
*! based on diagtest.ado (Aurelio Tobias, STB-56: sbe36)
*! and further suggestions from Tom Steichen

i.e., diagt from SSC. As you have correctly stated, diagt is not an official Stata command. We call such commands "community-contributed commands" these days. You will also come across the term "user-written programs". As explained in the FAQs, you are supposed to state where such commands come from (e.g., SSC, SJ, STB). It is usually best to also include the results of

which command

where you substitute command with the community-contributed command that you are using; I have done this for diagt above. The reason for this specific FAQ is that there are often multiple versions of community-contributed commands out there and we need to make sure that we are talking about the same thing. Often, but not always, the latest version of community-contributed commands is hosted on the SSC.

Let's move on. In diagt.ado, I find the comment

* bugfix to preserve matrices.

which makes me believe that your guess is correct: older versions of diagt have indeed deleted matrices. This bug seems to have been fixed. You probably have an out-of-date version of diagt installed.

I should mention that I have used Stata 14.1 to run your code.

[Anselm Schulz]

Thank you again for your excellent advice and I will try to improve my posts in the future. Sorry that I did not sees the obvious – I always thought that the error had to be in my code…

However I can not find your version of diagt 2.032 when I search in stata.
My version of stata is 16.1 and diagt is

which diagt
[...]
*! diagt 2.0.5, 30 June 2003
*! by PT Seed ([email protected])
*! based on diagtest.ado (Aurelio Tobias, STB-56: sbe36)
*! and further suggestions from Tom Steichen

I have tried an earlier version of diagt

*! diagt 1.1.2, 5 Oct 2000                       (STB-59: sbe36.1)
*! by PT Seed ([email protected])
*! based on diagtest.ado (Aurelio Tobias, STB-56: sbe36)

with this version everything works just fine, however the older version provides not all diagnostic test values included in the 2.0.5 version.
Is it perhaps possible that you could attach your diagt files to your post if they are difficult to download/find elsewhere or you post the content of diagt.ado and I replace/addust my version?

[Paul T Seed]

Thank you for the information that user-written programs. are now known as "community-contributed commands".
The latest version of diagti is below. I am glad it continues to be useful, but I do not have the time to maintain it properly.
If users want more modern-looking graphs, they can use the other "community-contributed commands" batplot, with the option notrend.

*! diagti 2.053, 7 April 2004
*! by PT Seed ([email protected])
*! Immediate version of diagt,
*! based on diagtest.ado (Aurelio Tobias, STB-56: sbe36)
*
* Likelihood ratios added
* Odds ratios and ROC areas added
* CI based on csi, or

* add 0.5 when zeros appear of with -sf- option when calculating LRs
* bugfix to preserve matrices.

program define diagti , rclass
version 8.0

set trace off

tokenize "`*'", parse(" ,")
local i= 1
while `i' <= 4 {
confirm integer number ``i''
local i = `i' + 1
}
local a = `1'
local b = `2'
local c = `3'
local d = `4'

local n1_ = `a' + `b'
local n2_ = `c' + `d'
local n_1 = `a' + `c'
local n_2 = `b' + `d'

local n = `n_1' + `n_2'
mac shift 4
local options "Prev(string) Level(real $S_level) noTable sf sf0 tb Woolf Bamber Hanley Odds Star dp(integer 1)"
parse "`*'"

if "`prev'" ~= "" {
local _prev_ "(lr)"
}

if "`star'" == "star" {
local star "di ""
}
else {
local star "*"
}

local dp2 = `dp' + 1
local dp3 = `dp' + 2

`star' tb affects the quantities based on risk ratios:

if index("`prev'", "%") {
local prev = real(subinstr("`prev'","%","",1))
cap assert `prev' >=0 & `prev' < 100
if _rc {
di in red "Prevalence must be between 0 and 100%
exit _rc
}
}
else if "`prev'" ~= ""{
local prev = `prev'*100
cap assert `prev' >=0 & `prev' < 100
if _rc {
di in red "Prevalence must be between 0 and 100%
di in red "Percentages must be indicated by ''%'' sign
exit _rc
}
}


if ((`a'*`b'* `c'*`d' == 0) & ("`sf0'" ~= "")) | ("`sf'" ~= "") {
local sf "(sf)"
local sk_sf = 1
local tb
local woolf
* NOTE: sf overrides tb and woolf
}
else {
local sk_sf = 0
if "`tb'" ~= "" {
local _tb_ "(tb)"
local _woolf_ "(tb)"
}

if "`woolf'" ~= "" {
local _woolf_ "(Woolf)"
}
}

* set trace on

`star' Prevalence
if "`_prev_'" == ""{
qui cii `n' `n1_', level(`level')
local prev = 100*r(mean)
local prev_lb = 100*r(lb)
local prev_ub = 100*r(ub)
}

`star' Sensitivity
qui cii `n1_' `a' , level(`level')
local sens = 100*r(mean)
local sens_lb = 100*r(lb)
local sens_ub = 100*r(ub)

`star' Specificity
qui cii `n2_' `d', level(`level')
local spec = 100*r(mean)
local spec_lb = 100*r(lb)
local spec_ub = 100*r(ub)

`star' Without substitution formula
if "`sf'" == "" {
qui csi `a' `c' `b' `d' , or `tb' `woolf' level(`level')
`star' Likelihood ratios
local lrp = r(rr)
local lrp_ub = r(ub_rr)
local lrp_lb = r(lb_rr)

`star' Odds ratio
local or = r(or)
local or_lb = r(lb_or)
local or_ub = r(ub_or)

qui csi `b' `d' `a' `c' , `tb' `woolf' level(`level')
local lrn = r(rr)
local lrn_ub = r(ub_rr)
local lrn_lb = r(lb_rr)

* If prevalence not given, PVs and CIs worked out directly.
if "`_prev_'" == "" {
`star' PPV
qui cii `n_1' `a', level(`level')
local ppv = 100*r(mean)
local ppv_lb = 100*r(lb)
local ppv_ub = 100*r(ub)
`star' NPV
qui cii `n_2' `d', level(`level')
local npv = 100*r(mean)
local npv_lb = 100*r(lb)
local npv_ub = 100*r(ub)
}

`star' If prevalence given, PVs and CIs worked out from LRs.
if "`_prev_'" ~= "" {
`star' Give CI for PVs
local odds = `prev'/(100-`prev')
local ppv = (`sens'*`prev') / ((`sens'*`prev')+((100-`spec')*(100-`prev')))*100
local ppv_lb = `lrp_lb'*`odds'
local ppv_ub = `lrp_ub'*`odds'

local npv = ((`spec'*(100-`prev')) / (((100-`sens')*`prev')+(`spec'*(100-`prev'))))*100
local npv_lb = `lrn_ub'*`odds'
local npv_ub = `lrn_lb'*`odds'
foreach v in ppv_lb ppv_ub {
local `v' = 100*``v''/(1+``v'')
}
foreach v in npv_lb npv_ub {
local `v' = 100/(1+``v'')
}
}
}

`star' With substitution formula
else {
local zval = invnorm(1-(1-`level'/100)/2)
local a0 = `a' + .5
local b0 = `b' + .5
local c0 = `c' + .5
local d0 = `d' + .5

`star' NOTE: Likelihood ratios depend on sensitivity & specificity
`star' They are not affected by prevalences
local lrp = (`a0'/(`a0'+`b0'))/(`c0'/(`c0'+`d0'))
local lrp_v = `b0'/(`a0'*(`a0'+`b0')) + `d0'/(`c0'*(`c0'+`d0'))
local lrp_lb = exp(ln(`lrp')-`zval'*sqrt(`lrp_v'))
local lrp_ub = exp(ln(`lrp')+`zval'*sqrt(`lrp_v'))

local lrn = (`b0'/(`a0'+`b0'))/(`d0'/(`c0'+`d0'))
local lrn_v = `a0'/(`b0'*(`a0'+`b0')) + `c0'/(`d0'*(`c0'+`d0'))
local lrn_lb = exp(ln(`lrn')-`zval'*sqrt(`lrn_v'))
local lrn_ub = exp(ln(`lrn')+`zval'*sqrt(`lrn_v'))

local or = (`a0'/(`b0'))/(`c0'/(`d0'))
local or_v = 1/`a0' + 1/`b0' + 1/`c0' + 1/`d0'
local or_lb = exp(ln(`or')-`zval'*sqrt(`or_v'))
local or_ub = exp(ln(`or')+`zval'*sqrt(`or_v'))

`star' If prevalence not given, PVs and CIs worked out directly.
if "`_prev_'" == "" {
`star' PPV
qui cii `n_1' `a', level(`level')
local ppv = 100*r(mean)
local ppv_lb = 100*r(lb)
local ppv_ub = 100*r(ub)
`star' NPV
qui cii `n_2' `d', level(`level')
local npv = 100*r(mean)
local npv_lb = 100*r(lb)
local npv_ub = 100*r(ub)
}

`star' If prevalence given, PVs and CIs worked out from LRs.
if "`_prev_'" ~= "" {
`star' Give CI for PVs
local odds = `prev'/(100-`prev')
local ppv_lb = `lrp_lb'*`odds'
local ppv = ((`sens'*`prev') / ((`sens'*`prev')+((100-`spec')*(100-`prev'))))*100
local npv = ((`spec'*(100-`prev')) / (((100-`sens')*`prev')+(`spec'*(100-`prev'))))*100
local ppv_ub = `lrp_ub'*`odds'
local npv_lb = `lrn_ub'*`odds'
local npv_ub = `lrn_lb'*`odds'
foreach v in ppv_lb ppv_ub {
local `v' = 100*``v''/(1+``v'')
}
foreach v in npv_lb npv_ub {
local `v' = 100/(1+``v'')
}
}
}

`star' ROC area
if _N >0 {
preserve
drop _all
label drop _all
}
qui {
tabi `a' `b' \ `c' `d'
expand pop
drop if pop == 0
drop pop

rename col test
label var test "Test result"
label define test 1 "Pos." 2 "Neg." 0 "Neg."
label values test test

rename row true
label var true "True disease status"
label define true 1 "Abnormal" 2 "Normal" 0 "Normal"
label values true true

recode true 2 = 0
recode test 2 = 0

roctab true test , `bamber' `hanley' level(`level')
local roc = r(area)
local roc_lb = r(lb)
local roc_ub = r(ub)
if "`bamber'" ~= "" {
local roc_opt = "(Bamber)"
}
if "`hanley'" ~= "" {
local roc_opt = "(Hanley)"
}
}


`star' Table
if "`table'" == "" {
tabulate true test, `options'
}



di _n in gr _col(51) "[`level'% Confidence Interval]"
di in gr "---------------------------------------------------------------------------"
if "`_prev_'" == "" {
di in gr "Prevalence Pr(A)" _col(44) in ye %6.`dp'f `prev' "%" _col(54) %6.`dp'f `prev_lb' "%" _col(65) %6.`dp'f `prev_ub' "%"
}
else {
di in gr "Prevalence Pr(A)" _col(44) in ye %6.`dp'f `prev' "%" in gr _col(54) "----- (given) -----"
}
di in gr "---------------------------------------------------------------------------"
di in gr "Sensitivity Pr(+|A)" _col(44) in ye %6.`dp'f `sens' "%" _col(54) %6.`dp'f `sens_lb' "%" _col(64) %6.`dp'f `sens_ub' "%"
di in gr "Specificity Pr(-|N)" _col(44) in ye %6.`dp'f `spec' "%" _col(54) %6.`dp'f `spec_lb' "%" _col(64) %6.`dp'f `spec_ub' "%"
di in gr "ROC area (Sens. + Spec.)/2" _col(46) in ye %5.`dp2'f `roc' _col(56) %5.`dp2'f `roc_lb' _col(66) %5.`dp2'f `roc_ub' _skip(1) in gr "`roc_opt'"

di in gr "---------------------------------------------------------------------------"
di in gr "Likelihood ratio (+) Pr(+|A)/Pr(+|N)" _col(45) in ye %6.`dp2'f `lrp' _col(55) %6.`dp2'f `lrp_lb' _col(65) %6.`dp2'f `lrp_ub' in gr _skip(`sk_sf') "`sf'" _skip(1) "`_tb_'"
di in gr "Likelihood ratio (-) Pr(-|A)/Pr(-|N)" _col(45) in ye %6.`dp2'f `lrn' _col(55) %6.`dp2'f `lrn_lb' _col(65) %6.`dp2'f `lrn_ub' in gr _skip(`sk_sf') "`sf'" _skip(1) "`_tb_'"
di in gr "Odds ratio LR(+)/LR(-)" _col(44) in ye %7.`dp2'f `or' _col(55) %6.`dp2'f `or_lb' _col(63) %8.`dp2'f `or_ub' in gr _skip(`sk_sf') "`sf'" _skip(1) "`_woolf_'"
di in gr "Positive predictive value Pr(A|+)" _col(44) in ye %6.`dp'f `ppv' "%" _col(54) %6.`dp'f `ppv_lb' "%" _col(64) %6.`dp'f `ppv_ub' "%" in gr _skip(1) "`_prev_'"
di in gr "Negative predictive value Pr(N|-)" _col(44) in ye %6.`dp'f `npv' "%" _col(54) %6.`dp'f `npv_lb' "%" _col(64) %6.`dp'f `npv_ub' "%" in gr _skip(1) "`_prev_'"
di in gr "---------------------------------------------------------------------------"

* set trace on
if "`odds'" ~= "" {
local oprev = `prev'/(100-`prev')
if "`_prev_'" == "" {
local oprev_lb = `prev_lb'/(100-`prev_lb')
local oprev_ub = `prev_ub'/(100-`prev_ub')
}
local opos = `ppv'/(100-`ppv')
local opos_lb = `ppv_lb'/(100-`ppv_lb')
local opos_ub = `ppv_ub'/(100-`ppv_ub')

local oneg = (100-`npv')/`npv'
local oneg_lb = (100-`npv_lb')/`npv_lb'
local oneg_ub = (100-`npv_ub')/`npv_ub'

if "`_prev_'" == "" {
di in gr "Pre-test odds prev/(1-prev)" _col(45) in ye %6.`dp2'f `oprev' _col(55) %6.`dp2'f `oprev_lb' _col(65) %6.`dp2'f `oprev_ub' in gr _skip(1) "`_prev_'"
}
else {
di in gr "Pre-test odds prev/(1-prev)" _col(45) in ye %6.`dp2'f `oprev' in gr _col(54) "----- (given) -----"
}
di in gr "Post-test odds (+) Pr(A|+)/(1-Pr(A|+))" _col(45) in ye %6.`dp2'f `opos' _col(55) %6.`dp2'f `opos_lb' _col(65) %6.`dp2'f `opos_ub' in gr _skip(2) "`_prev_'"
di in gr "Post-test odds (-) Pr(A|-)/(1-Pr(A|-))" _col(45) in ye %7.`dp3'f `oneg' _col(55) %7.`dp3'f `oneg_lb' _col(65) %7.`dp3'f `oneg_ub' in gr _skip(1) "`_prev_'"
di in gr "---------------------------------------------------------------------------"

local retlist "oneg opos oprev "
foreach ret in `retlist' {
cap ret scalar `ret'_ub = ``ret'_ub'
cap ret scalar `ret'_lb = ``ret'_lb'
cap ret scalar `ret' = ``ret''
}
}

if "`sf'" ~= "" {
di in gr _n "(sf) Likelihood ratios are estimated using the substitution formula.
di in gr " 0.5 is added to all cell frequencies before calculation."
}


if "`_prev_'" ~= "" {
di in gr _n "(lr) Values and confidence intervals are based on likelihood
di in gr " ratios" _skip(`sk_sf') "`sf', assuming that the prevalence is known exactly."
}
if (`a'*`b'* `c'*`d' == 0) & ("`sf'" == "") {
di _n in gr " Missing values or confidence intervals may be estimated
di in gr " using the -sf- or -sf0- options."
}

local retlist "`retlist' roc or lrn lrp npv ppv prev spec sens "


foreach ret in `retlist' {
cap ret scalar `ret'_ub = ``ret'_ub'
cap ret scalar `ret'_lb = ``ret'_lb'
cap ret scalar `ret' = ``ret''
}


end

exit

[daniel klein]

Sorry that I did not sees the obvious – I always thought that the error had to be in my code…

In your defense, this really is not obvious at all. After looking into this, I see that

which diagt
[...]
*! diagt 2.0.5, 30 June 2003
*! by PT Seed ([email protected])
*! based on diagtest.ado (Aurelio Tobias, STB-56: sbe36)
*! and further suggestions from Tom Steichen

is from SJ (more specifically, it is package sbe36_2 from SJ 4-4). Despite being the first result listed and despite the version number's indication, the version 2.032 is the more recent (in terms of bug fixes) than 2.0.5.

Anyway, type

ssc install diagt

to get version 2.032 (and version 2.053 of diagti) from the SSC. As I have mentioned earlier, the SSC is often the best source for community-contributed stuff. Unfortunately, there are also many exceptions.


EDIT:

As you see, sometimes you are lucky and authors of the commands that you are using are active on Statalist and chime in. The code Paul has shown is included in the diagt package on SSC.

[Anselm Schulz]

Thank you so much!
I am deeply impressed over the advice and support I have got from you all - you really helped me out. Nice work!

[Elisa Masiak]

Good day I am very new in this STATA, I am interested in positive predictive values, I used diagt command. I am confused on which of these 3 values (underlined). Thank you in advance



[95% Confidence Interval]

Prevalence Pr(A) 10% 8% 13.1%

Sensitivity Pr(+A) 40.3% 28.1% 53.6%
Specificity Pr(-N) 90.3% 87.5% 92.7%
ROC area (Sens. + Spec.)/2 .653 .59 .716

Likelihood ratio (+) Pr(+A)/Pr(+N) 4.16 2.8 6.2
Likelihood ratio (-) Pr(-A)/Pr(-N) .661 .537 .812
Odds ratio LR(+)/LR(-) 6.3 3.54 11.2
Positive predictive value Pr(A+) 32.5% 22.2% 44.1%
Negative predictive value Pr(N-) 92.9% 90.4% 95%

[Nick Cox]

Elisa Masiak Welcome to Statalist. Please post that as the start of a new thread with a title mentioning the diagt command, which should be explained as a community-contributed command. The subject of your question isn't related at all to the thread title.