*! version 18 26Nov2025

help gpsmdpolest

Title

gpsmdpolest -- The command performs the estimation of the dose-response function

Syntax

gpsmdpolest outcome treatment_dimensions , gpsmd(string) model(string) exogenous(varlist) file_pred(string) numboot(numlist integer max=1)[dividingint(numlist integer max=1) matrtreat(string) level(numlist max=1) cutpoints(numlist integer max=1) index(string) ln(varlist) matrixwithresults(string)]

outcome: the outcome variable of interest.

treatment_dimensions: the variables storing the dimension of the treatment.

Description

gpsmdpolest The command performs the estimation of the dose-response function according to Egger and von Ehrlich (2013).
gpsmdpolest supports linear models where terms can be at a user-defined power. It further supports the logarithmic transformation of propensity score as well as the interactions.

Since the model includes the propensity score, a generated regressor (Wooldridge 2010), the t-approximation is not reliable for estimating the confidence intervals. Therefore, gpsmdpolest estimates confidence intervals by bootstrap (BC) (Carpenter and Bithell 2000).

The program does not produce any graphs. Drawing graphs in dimensions higher than two in Stata is not easy, and any graph would require some adjustments.
Therefore, the main program's output consists of a dataset that the user can further process with programs like graph3d (Rostam-Afschar and Jessen 2014) to obtain a graphical representation of the results – see also graph twoway contour.

The command should be used together with gpsmd, gpsmdcomsup, and gpsmdbal to estimate the dose-response function.

Note: gpsmd must be invoked before invoking gpsmdpolest.

Options

model(string): a string with the right side of the model.
The right side of the model must be explicitly written due to how the program parses inputs (e.g. "T1 + T2 + gps + T1*gps + T2*gps + T1^(2) + T2^2 + (gps^2) + ((T1*gps)^(2)) + (T2*gps)^2 + ln(gps) + (ln(gps))^2 + (T2*ln(gps))^2 + T2*ln(gps) "). The program does not support interactions between the dimensions of the treatment.

dividingint(numlist integer max=1): if matrtreat(string) is not specified the program generates a matrix by dividing the dimensions in dividingint number of intervals.
The Cartesian product of the extremes of the intervals in the different dimensions constitutes the set of treatment points for which the program estimates the response.
The set of treatment points will be stored in r(matrtreat) as a matrix with (`dividingint'+1)^|dimensions| (where |dimensions| is the number of dimensions) rows and columns equal |dimensions|.

matrtreat(string): The user can specify the treatment points for which she is interested in estimating the response.
Treatment points must be stored in a Stata matrix named as specified in matrtreat(string). The matrix must have one column for each treatment dimension.
A row of the matrix identifies a single point. The user can specify only one option between dividingint and matrtreat.

exogenous(varlist): the exogenous variables the user wants to use in the reduced equations.

file_pred(string): As explained above, the program does not generate any graphs. It generates a dataset with the necessary information for the user to generate the desired graphs (see below for a more detailed description of the file generated).
In file_pred(string), the user must specify the first characters of the name for the files .dta storing the results.

level(numlist max=1): the confidence level for the confidence intervals (default 0.05).

numboot(numlist integer max=1): the number of bootstrap samples. Since bootstrapping is the only way to obtain the confidence intervals, this is not an optional argument.

cutpoints(numlist integer max=1): the chosen number of discrete intervals of the dimensions of the treatment for the common support estimation.
It is worth noticing that when common support is required, the program estimates the dose-response function by using only those observations that lie on the common support.
It is suggested to use the same number used to calculate the common support.

index(string): the point t_d where the user wants to calculate the GPS. It can be "mean" or "p50": "mean" for the mean, and "p50" for the median.

ln(varlist): the treatment dimensions that have to be log-transformed.

matrixwithresults(string): if the argument is “T”, the program returns a matrix called r(returnresults)that includes all the results as well as the chosen doses.
The default is “T”. If “F”, the matrix is not generated. This option can be helpful when the number of treatment points exceeds Stata matrix limits.

Examples

Setup

Setting the dataset:

. clear all

. set obs 1200

Generating independent variables for the propensity score estimation:

. seed 13131

. gen X1 = 1* rnormal(0,1)

. gen X2 = 2* rnormal(0,1)

. gen X3 = 3* rnormal(0,1)

. gen X4 = 4* rnormal(0,1)

. gen X5 = 5* rnormal(0,1)

. gen X6 = 6* rnormal(0,1)

. gen X7 = 7* rnormal(0,1)

Generating the treatment dimensions:

. matrix R = (25, 2 \2, 25)

. drawnorm V1 V2, cov(R)

. gen T1= 1*X1 + .5*X2 + 1*X3 + .5*X4 + 1*X5 + .5*X6 + 1*X7 + V1

. gen T2= .5*X1 + 1*X2 + .5*X3 + 1*X4 + .5*X5 + 1*X6 + .5*X7 + V2

Generating the outcomes:

. gen res= rnormal(0, 25)

. gen Y = 1+ 2*T1 + 1.5*T2 + 1*X1 + 1.5*X2 + 2*X3 + 1*X4 + 1.5*X5 + 2*X6 + 1*X7 + res

The estimation:

. gpsmd T1 T2, exogenous(X1 X2 X3 X4 X5 X6 X7) gpsmd(GPS)

. gpsmdcomsup T1 T2, exogenous(X1 X2 X3 X4 X5 X6 X7) index("p50") cutpoints(2) obs_notsup(Commonsupport)

. gpsmdbal X1 X2 X3 X4 X5 X6 X7, index("p50") cutpoints(2) nq_gpsmd(4) discrtreat(Discretetreat) obs_notsup(Commonsupport)

. gpsmdpolest Y T1 T2, gpsmd(GPS) exogenous(X1 X2 X3 X4 X5 X6 X7) model("T1 + T2 + GPS + T1*GPS + T2*GPS") file_pred(ExampleStata) numboot(1000) dividingint(3) index("p50") cutpoints(2)

Stored results

Variables

gpsmdpolest generates the variables specified in model but the treatment dimensions and the GPS.
All the variables are named starting with __I_, __P_, or __LN_. I_ represent interaction variables and P_ variables with power (variables like ((T1*GPS)^(2)) are named with both, e.g., __P_2_I_T1_GPS).
LN_ represents logarithmic transformation. The user should check whether in her dataset there are variables whose names begin with these characters.
If it is the case, it is preferable to change their names before running the program.

Datasets

gpsmdpolest generates the following dataset:

- One dataset named as specified in file_pred(). It includes one row for each treatment point.
The columns store the response, the partial derivatives, and the upper and lower bound of the confidence intervals.

gpsmdbal stores the following in r():

Macros

r(gpsmd) the name of the variable with the GPS estimates

r(exogenous) the exogenous variables for the reduced equation estimation

r(Dimensions) the dimensions of the treatment(gpsmdbal)

r(listgenvar)) the program generates the variables as specified in model(string). This macro reports the list of the variables generated

r(regmodel) the command for the regression for the polynomial estimation

r(cmd) macro with the command

r(cmdline) macro with the cdmline. This macro reports the command just invoked, including options and specifications

r(Outcome) macro containing the name of the outcome variable

Matrices

r(matrtreat) a matrix with the treatment points for which the dose-response has been estimated

r(returnresults) if matrixwithresults(string) is equal to "T", the program returns a matrix with the same information included in `file_pred’.dta.

Bibliography and Sources

Carpenter, James, and John Bithell. 2000. ‘Bootstrap Confidence Intervals: When, Which, What? A Practical Guide for Medical Statisticians’. Statistics in Medicine 19 (9). Wiley Online Library: 1141–64.

Egger, Peter H., and Maximilian von Ehrlich. 2013. ‘Generalized Propensity Scores for Multiple Continuous Treatment Variables’. Economics Letters 119 (1): 32–34.

Egger, Peter H., and Peter Egger. 2016. ‘Heterogeneous Effects of Tariff and Nontariff Policy Barriers in General Equilibrium’. Beiträge zur Jahrestagung des Vereins für Socialpolitik 2016: Demographischer Wandel - Session: Trade Barriers, No. G18-V3, ZBW - Deutsche Zentralbibliothek für Wirtschaftswissenschaften, Leibniz-Informationszentrum Wirtschaft, Kiel und Hamburg.

Egger, Peter H., Maximilian v. Ehrlich, and Douglas R. Nelson. 2020. ‘The Trade Effects of Skilled versus Unskilled Migration’. Journal of Comparative Economics 48 (2): 448–64.

Leuven, Edwin, and Barbara Sianesi. 2003. ‘PSMATCH2: Stata Module to Perform Full Mahalanobis and Propensity Score Matching, Common Support Graphing, and Covariate Imbalance Testing’.

Präg, Patrick. 2019. ‘Visualizing Individual Outcomes of Social Mobility Using Heatmaps’. Socius 5 (January). SAGE Publications: 1-2.

Rostam-Afschar, Davud, and Robin Jessen. 2014. ‘GRAPH3D: Stata Module to Draw Colored, Scalable, Rotatable 3D Plots’. Boston College Department of Economics.

Wooldridge, Jeffrey M. 2010. Econometric Analysis of Cross Section and Panel Data. MIT press.