Generalized ordered probit: "Looping over program" problem

[Inna Semel]

Dear Statalist User,

I have a piece of code for estimation generalized ordered probit model from the book A.M.Jones (2007) Applied Econometrics for Health Economists (Ch.4 Bias in sell-reported data).
I use Stata 14.2. This code works fine:

reshape long vig_mobility_, i(id) j(vignum)
tabulate vignum, generate(vigdum)
drop vigdum1

/* Program for GOP  - cutpoints are functions of individual characteristics, health index function of vignette dummies */

cap program drop gop
program define gop

args lnf b m1 m2 m3 m4

    tempvar p1 p2 p3 p4 p5
quietly {
    gen double `p1' = 0
    gen double `p2' = 0
    gen double `p3' = 0
    gen double `p4' = 0
    gen double `p5' = 0

    replace `p1' = normal(`m1'-`b')
    replace `p2' = normal(`m2'-`b') - normal(`m1'-`b')
    replace `p3' = normal(`m3'-`b') - normal(`m2'-`b')
    replace `p4' = normal(`m4'-`b') - normal(`m3'-`b')
    replace `p5' = 1                   - normal(`m4'-`b')

    replace `lnf' = (vig_mobility_==1)*ln(`p1')+(vig_mobility_==2)*ln(`p2')+(vig_mobility_==3)*ln(`p3')+(vig_mobility_==4)*ln(`p4')+(vig_mobility_==5)*ln(`p5')  
}
end

*** Estimate GOP for the vignettes                                                       ***
ml model lf gop (xb: vigdum*, nocons) (mu1: $cvars) (mu2: $cvars) (mu3: $cvars) (mu4: $cvars) if vig_mobility_!=.
ml search
ml maximize

 drop vigdum*
 reshape wide vig_mobility_, i(id) j(vignum)

I need to use this code several times. So I try to use the loop:

cap program drop gop
program define gop

args lnf b m1 m2 m3 m4

    tempvar p1 p2 p3 p4 p5
quietly {
    gen double `p1' = 0
    gen double `p2' = 0
    gen double `p3' = 0
    gen double `p4' = 0
    gen double `p5' = 0

    replace `p1' = normal(`m1'-`b')
    replace `p2' = normal(`m2'-`b') - normal(`m1'-`b')
    replace `p3' = normal(`m3'-`b') - normal(`m2'-`b')
    replace `p4' = normal(`m4'-`b') - normal(`m3'-`b')
    replace `p5' = 1                   - normal(`m4'-`b')

    replace `lnf' = (vig_`v'_==1)*ln(`p1')+(vig_`v'_==2)*ln(`p2')+(vig_`v'_==3)*ln(`p3')+(vig_`v'_==4)*ln(`p4')+(vig_`v'_==5)*ln(`p5')  
}
end


local varlist = "mobility vigact depressed worry pain discomfort relationships conflicts visionm visionf  sleepn sleepm memory learning care appear"

foreach v of local varlist {
    reshape long vig_`v'_, i(id) j(vignum)
    tabulate vignum, generate(vigdum)
    drop vigdum1        

ml model lf gop (xb: vigdum*, nocons) (mu1: $cvars) (mu2: $cvars) (mu3: $cvars) (mu4: $cvars) if vig_`v'_!=.
ml search
ml maximize

 drop vigdum*
 reshape wide vig_`v'_, i(id) j(vignum)
 
}

I see the result of the execution

reshape long vig_`v'_, i(id) j(vignum)
tabulate vignum, generate(vigdum)

in the loop. And an error "vig__ not found".

I can't handle it myself. Any help would be very highly appreciated.

[Phil Bromiley]

In debugging this kind of problem, set trace on is extremely helpful - it will tell you how Stata is interpreting your macros. You can also use display [ di " `v' " ] to see what value is in the macro. It looks like you have a value of `v' that is missing. I wonder if the double space in the local varlist= statement is causing the problem.

[Richard Williams]

I haven't read the article so I don't know exactly how the authors define a goprobit model. But the gologit2 program, available from SSC, has an option for probit link. You might see if it does what you want rather than try to program it yourself.

[Inna Semel]

Dear Phill,

Thank you for your advice!
Unfortunately, "double space in the local varlist = statement" isn't a the source of the problem, I checked.
I'm not really strong at debugging programs. I will try to figure it out. Thanks!

[Inna Semel]

Dear Richard,

Thank you for your reply!

It seems to me that I have a problem not in estimation the model as such (it could be any model in the loop , not necessarily goprobit), but in transferring control from loop to program (calling program inside the loop ) — somewhere in this place I make a mistake.

Thanks again!

[Jorrit Gosens]

This is fairly difficult to help with without a data example. Can you share some example data with dataex or otherwise a csv or txt file?

If I just do:

clear
set trace on 
local varlist = "  mobility vigact depressed worry pain discomfort relationships conflicts visionm visionf  sleepn sleepm memory learning care appear"
foreach v of local varlist {
gen `v'=""
}

the loop finishes without an error about some empty variable name, so difficult to understand why it would be different when the command is reshape.

[Inna Semel]

Dear Jorrit

Thanks for your reply.

I added small data sample and short do-file to illustrate my question.

[Jorrit Gosens]

I'm fairly certain that the problem is that when you define your program your varlist has not been defined yet and `v' evaluates as "" (i.e., nothing).
I'm trying to put teh program into the loop but struggling a bit.

Edit: sorry, that doesn't seem to be it. Program should take in the value for `v' when running the loop.

Edit2: or it is the way the program is defined after all? SOrry, I'm a bit lost here. Hope someone else can help.

version  14
clear all
set more off
import delimited Inna_data_example1.csv, clear
drop vig_vigact_2-vig_learning_5
set trace on


global cvars c.age i.female i.married i.educ i.urban i.quintile_c i.relig i.work i.alcoh i.smoke

* I try looop, becouse I have 16 sets of vignettes for 16 health domians. It does't work - I think I do a mistake when I call program
 
cap program drop gop
program define gop
display vig_`v'_
end


*local varlist mobility vigact depressed worry pain discomfort relationships conflicts visionm visionf sleepn sleepm memory learning care appear

foreach v in mobility{
    reshape long vig_`v'_, i(id) j(vignum)
    tabulate vignum, generate(vigdum)
    drop vigdum1        

gop
*** Estimate GOP for the vignettes                                                       ***
*ml model lf gop (xb: vigdum*, nocons) (mu1: $cvars) (mu2: $cvars) (mu3: $cvars) (mu4: $cvars) if vig_mobility_!=. 
*ml search
*ml maximize

 *drop vigdum*
 *reshape wide vig_mobility_, i(id) j(vignum)
 
}

the above, just used to simplify and identify where the error occurs, yields:

- gop
  -------------------------------------------------------------------------------------------------------------- begin gop ---
  - display vig_`v'_
  = display vig__
vig__ not found
  ---------------------------------------------------------------------------------------------------------------- end gop ---
  }

[Inna Semel]

Thanks Jorrit! I really appreciate your time and effort.

[Inna Semel]

I understood the problem.
Here

  
 foreach v in mobility{

v - local macro thus not visible outside the loop.
Therefore, the program gop doesn' t see `v'.
But I still do not know how to solve it.

[Jorrit Gosens]

So it has something to do with how the local macro is used in the program. I have tried some more with different quotes around the local but not manage to get the program to use the string that the local (intuitively) should hold.
For some reason, it does work if you use a global macro. If anyone here is more educated in the dark art of how Stata uses global/local macros, then some further explanation would be helpful.
In any case, the below does what youre after. If you keep 'set trace on' you can see that the correct strings from your local varlist are used in the program now
Also note the space between the global macro name and the ==, otherwise Stata tries to find a global called blah==1)

version  14
clear all
set more off
import delimited Inna_data_example1.csv, clear

set trace on

global cvars c.age i.female i.married i.educ i.urban i.quintile_c i.relig i.work i.alcoh i.smoke

cap program drop gop
program define gop

args lnf b m1 m2 m3 m4

    tempvar p1 p2 p3 p4 p5
quietly {
    gen double `p1' = 0
    gen double `p2' = 0
    gen double `p3' = 0
    gen double `p4' = 0
    gen double `p5' = 0

    replace `p1' = normal(`m1'-`b')
    replace `p2' = normal(`m2'-`b') - normal(`m1'-`b')
    replace `p3' = normal(`m3'-`b') - normal(`m2'-`b')
    replace `p4' = normal(`m4'-`b') - normal(`m3'-`b')
    replace `p5' = 1                   - normal(`m4'-`b')

    replace `lnf' = (vig_$blah ==1)*ln(`p1')+(vig_$blah ==2)*ln(`p2')+(vig_$blah ==3)*ln(`p3')+(vig_$blah ==4)*ln(`p4')+(vig_$blah ==5)*ln(`p5')  
}
end

 
local varlist mobility vigact depressed worry pain discomfort relationships conflicts visionm visionf  sleepn sleepm memory learning care appear


foreach v of local varlist {
    reshape long vig_`v'_, i(id) j(vignum)
    tabulate vignum, generate(vigdum)
    drop vigdum1        
    global blah `v'_
    
*** Estimate GOP for the vignettes                                                       ***
ml model lf gop (xb: vigdum*, nocons) (mu1: $cvars) (mu2: $cvars) (mu3: $cvars) (mu4: $cvars) if vig_`v'_!=.
ml search
ml maximize

 drop vigdum*
 reshape wide vig_`v'_, i(id) j(vignum)
}

[Inna Semel]

Thanks a lot Jorrit!

I tried to redefine the global macro `v' in a loop - but I didn't succeed. Your method works great!!! (I obviously need to learn more about local and global macros in Stata).
Many thanks again!