Is there an easy way to feed an adjusted lp problem back into the ivmte() function?

[jongohkim91]

Dear Mr. Shea,

Currently I am trying to restrict some variables in MTR functions(m0 and m1) to have an identical coefficient regardless of the treatment status. Before, you have given me this reply: Hence, I have run the ivmtefunction first without restricting the variables to obtain an initial lp problem model(lp.result$model), and I adjusted the lp model to have an additional constraint row. More specifically, I have added a row having mostly zero values but c(-1,1) values for the columns with a name associated with the variable of interest. For instance, if we use the same AE example from Angrist and Evans (1998), I will run the ivmte function with parameters specified below.

results_ate <- ivmte(data = dataset,
                      target = "ate",
                      m0 = ~ u + I(u^2) + morekids,
                      m1 = ~ u + I(u^2) + morekids,
                      propensity = morekids ~ samesex,
                      ivlike = c(worked ~ morekids | samesex,
                                 worked ~ morekids + samesex + samesex*morekids),
                      solver = 'lpsolveapi',
                      noisy = TRUE)

I would like to make sure the MTR terms for morekids in m0 & m1 will have the same coefficient and insert an additional constraint row into a new lp model. The row has values like these in the A matrix.

 slack1-         slack1+         slack2-         slack2+         slack3-         slack3+ 
              0               0               0               0               0               0 
        slack4-         slack4+         slack5-         slack5+         slack6-         slack6+ 
              0               0               0               0               0               0 
[m0](Intercept)           [m0]u      [m0]I(u^2)    [m0]morekids [m1](Intercept)           [m1]u 
              0               0               0              -1               0               0 
     [m1]I(u^2)    [m1]morekids 
              0               1 

It has = & 0 for senseand rhsrespectively.

Let's say the new lp model is stored in the variable new_lp. Since I would like to check whether the change I have made results in an outcome I desired, I must ,as far as I know, feed new_lp back into the ivmtefunction to obtain all necessary information(results_ate$gstar.coef, results_ate$bounds and etc.).

However, the ivmte() function is deeply connected with few big wrapper functions(ivmteEstiamte() & audit()) and its "daughter" functions(lpSetup(), bound(), and etc.). I, therefore, should edit various call variables throughout the code or edit the parameters of many functions to assure .

Is there any much easier way to use the new_lp as an input in the ivmte() or should I go through all the steps from ivmte() function so that all relevant parameters are set as it should be?

Thank you very much in advance for your reply.

Best regards, Jongoh

[a-torgovitsky]

@jongohkim91 I'm interested in your use case, could you explain your motivation a bit more? You are using morekids as the treatment variable, but then also including it as a covariate in m0 and m1. What is your goal from doing that?

(@jkcshea this might actually cause some problems with our weighting expressions. It might be something we should expressly disallow.)

[jongohkim91]

@a-torgovitsky Sorry for causing confusion. Yes, as you pointed it out, it does not make sense to restrict morekids in Angrist and Evans setting.

What I exactly want is to restrict some control variables which enters in the MTR function to have the same coefficient in both m0& m1. There are two major reasons why I would like to do this.

First, I think some variables do not vary according to a treatment status so I would like them to be restricted. Second, I would like to check whether I have correctly applied the MTR method by comparing a result from MTR analysis with a result from Stata using traditional MTE framework with restricted variables. I have already a result produced by Stata and I would like to expand the analysis to the MTR setting. I should get a similar result when I am using the MTR framework as in the MTE case.

After I obtain a coherent result, then I would like to take advantage of the flexibility of MTR analysis by adding more variables into the MTR function to see which instrument(or intervention) affects people the most.

This is the motivation behind this issue and I hope I gave you a clear explanation.

[a-torgovitsky]

@jongohkim91 It makes more sense. But if you restrict the coefficient to be the same in the two MTRs, then the MTE will not depend on this control variable, right? So your previous MTE-based analysis was not conditioning on these control variables? Or was it?

@jkcshea I wonder how hard it would be to put in this functionality. The first thing that comes to mind is letting the user pass a list of control variables through a new option. Then all coefficients in that list are restricted to be identical across the two MTRs. This requires parsing variable names, which we know is a lot of fun, but maybe all of the structure is in place already to be able to do that fairly easily?

[jkcshea]

@jongohkim91 It will take some time to update the package to allow for the constraints you describe. So attached is an example of how you can perform the exercise you describe using the current output from ivmte. There, I estimate the bounds and the MTR coefficients when I include your new equality constraint. Although you will be unable to perform an audit as in the package, you can set up the problem as though an audit was performed (more detailed description in the R file). The other objects returned by the package (e.g. propensity score model, S-weights) will be as before. Let me know if you have any questions.

example.zip

As a reminder, the MTR coefficient estimates returned by ivmte are not estimates of the true coefficient values. I mention this simply because your post suggested you were interested in the MTR coefficients.

@a-torgovitsky Ah, yes, parsing variable names---good times. Yes, I think we've paid the fixed cost, so (hopefully!) it will not be very difficult to implement. At the moment, I think the user must pass the control variables as a formula. The reason is that there may be interactions and factor variables that need to be expanded. So the arguments will look something like:

ivmte(m0 = ~ x1 + I(x2^2)+ factor(x3):x4 + u,
      m1 = ~ x1 + I(x2^2)+ factor(x3):x4 + u
      controls = ~ I(x^2) + factor(x3):x4,
      ...)

I'll let you know if that changes, or if I come up with something better.

[a-torgovitsky]

@jkcshea Makes sense and that seems totally fine to me. We probably want to call the option equalcoef or equal_coef or equal.coef depending on our naming conventions. Open to other names as well, that's just the first thing that popped in my head.

[jongohkim91]

@a-torgovitsky Yes, we were conditioning the restricted control variables in the MTE analysis. It is true that theses variables do not influence the MTE directly but we included them since we are more interested in looking into the effects of the unrestricted control variables on various treatment effects. To achieve more precise estimates of these unrestricted variables, we have included multiple restricted control variables in the MTE function and also planning to include them in MTR functions as well.

@jkcshea Thank you very much for the example code! It helped me a lot. Could you please elaborate a bit more of the MTR coefficients? What do you mean that the MTR coefficient estimates returned by ivmte are not estimates of the true coefficient values? I and my fellow researchers are not specifically interested in the MTR coefficients but we need to compare the MTR coefficients with the MTE coefficients from from STATA to assure both results are consistent.

[jkcshea]

What do you mean that the MTR coefficient estimates returned by ivmte are not estimates of the true coefficient values?

If you are partially identified, then there are many possible values of the MTR coefficients that are observationally equivalent. The set of coefficient estimates returned by the package is just one of those possible values. Even if the true value of the target parameter was equal to, say, its upper bound, it is possible for the coefficients in results$gstar.coef$max.g0 and result$gstar.coef$max.g1 to be very different from the true MTR coefficients, e.g. because the solution to the LP problem is on an edge rather than a vertex.

[a-torgovitsky]

@a-torgovitsky Yes, we were conditioning the restricted control variables in the MTE analysis. It is true that theses variables do not influence the MTE directly but we included them since we are more interested in looking into the effects of the unrestricted control variables on various treatment effects. To achieve more precise estimates of these unrestricted variables, we have included multiple restricted control variables in the MTE function and also planning to include them in MTR functions as well.

Maybe I don't understand, but I think what you want to do here is the opposite: do not constrain the MTR functions to have the same coefficients.

For example, consider MTR(d, u | x) = a{d} + b{d}u + c_{d}x

If you constraint c{0} = c{1}, then you will have MTE(u | x) = (a{1} - a{0}) + (b{1} - b{0})u which does not depend on x.

So if you previously estimated an MTE specification that depends on x, and you want to reproduce that specification with MTRs, then you want to let your MTRs have different coefficients on x, i.e. allow for c{0} != c{1}.

[jongohkim91]

@jkcshea Ah, I see. Thank you! I will keep that in mind.

@a-torgovitsky Please allow me to explain more in detail. In our setting, we have multiple control variables. Most of them will be restricted and a few of them will be allowed to have different coefficient according to one's treatment status.

I will simply put x1 as a control variable which we would restrict its MTR coefficient to have the same values in any treatment status, and x2 is a control variable which could have different coefficient values.

Let MTR(d, u | x1, x2) = a{d} + b{d}u + c{d}x1 + e{d}x2.

We are interested in looking into the effects of x2 and believe the effect of x1 wouldn't depend on the treatment status, d. Hence I restrict x1to have identical coefficient for MTR(0, u | x) and MTR(1, u | x).

Then, MTE(u | x1, x2) = (a{1} - a{0}) + (b{1} - b{0})u + (e_{1} - e{0})x2

It is true that the MTE does not depend on x1. However more accurate estimation of e{1} - e{0} is possible when the term c{d}x1 is included in the MTR function.

I hope this makes sense.

[a-torgovitsky]

It is true that the MTE does not depend on x1. However more accurate estimation of e{1} - e{0} is possible when the term c{d}x1 is included in the MTR function.

That motivation certainly makes sense.

But what specification for the MTE did you estimate in Stata? Did it include both x1 and x2? If so, then restricting the coefficient on x1 to be the same in the MTRs will imply a more restrictive MTE specification than what you had previously. Maybe that's a good thing for your research goals -- it can certainly help increase precision. I was just reacting to your goal of trying to get the specification of MTRs to be the same as the previous MTE specification.

[jongohkim91]

@a-torgovitsky Thank you for your comments. I truly appreciate your feedback.

But what specification for the MTE did you estimate in Stata? Did it include both x1 and x2?

Yes, we included both x1 and x2 in STATA. Let x1 be a vector of restricted variables and x2 be a vector of unrestricted variables. Our MTE specification is as follows. MTE(u | x1, x2) = (a{1} - a{0}) + (b{1} - b{0})u + (e{1} - e{0})x2 + (f{1} - f{0})x2*u

If so, then restricting the coefficient on x1 to be the same in the MTRs will imply a more restrictive MTE specification than what you had previously. Maybe that's a good thing for your research goals -- it can certainly help increase precision.

Since we are aiming to estimate an effect of a policy, I would say it is better for us to have more precise estimation. However, I am not sure whether I fully understand why such phenomenon occurs.

I guess it is because MTE(u | x1, x2) is simply MTE(u | x2), and thus a corresponding MTR specification is
MTR(d, u | x2) = a{d} + b{d}u + e{d}x2 + f{d}x2*u.

Hence MTR(d, u | x1, x2) = a{d} + b{d}u + c{d}x1 + e{d}x2 + f_{d}x2*u would be more restrictive, and it will result in more precise estimation of MTE.

In order to obtain a similar result from the STATA analysis using MTE(u | x1, x2) = MTE(u | x2),

I should specify MTR functions as MTR(d, u | x2) = a{d} + b{d}u + e{d}x2 + f{d}x2*u

rather than MTR(d, u | x1, x2) = a{d} + b{d}u + c{d}x1 + e{d}x2 + f_{d}x2*u.

Did I understand your point correctly?

[a-torgovitsky]

I didn't follow the premise:

Yes, we included both x1 and x2 in STATA. Let x1 be a vector of restricted variables and x2 be a vector of unrestricted variables. Our MTE specification is as follows. MTE(u | x1, x2) = (a{1} - a{0}) + (b{1} - b{0})u + (e{1} - e{0})x2 + (f{1} - f{0})x2*u

This looks like you are only including x2, not x1. Do you mean that you are also including x1 in propensity score estimation?

[jongohkim91]

@a-torgovitsky I am sorry for the late reply and confusing you.

MTE(u | x1, x2) = (a{1} - a{0}) + (b{1} - b{0})u + (e{1} - e{0})x2 + (f{1} - f{0})x2u since (c{1} - c{0})x1= 0 x1 = 0.

Yes we include both x1 and x2 in propensity score estimation.

[a-torgovitsky]

That sounds like you are using x1 as an instrument. i.e. it is in the propensity score, but excluded from the outcome equation (MTE specification). That doesn't seem like what you want.

[jkcshea]

The package has been updated to handle equality constraints on the MTR coefficients. The user passes a one-sided formula to the argument equal.coef indicating which terms in m0 and m1 should have equal coefficients. There are checks in place to ensure that equal.coef is passed correctly, e.g. the function checks whether all the terms in equal.coef are contained in m0 and m1.

Here is an example where I impose the equality constraint on the term x1,

> set.seed(10L)
> dtcf <- ivmte:::gendistCovariates()$data.full
> args <- list(ivlike = ey ~ 1 + d + x1,
+              data = dtcf,
+              propensity = d ~ x1 + z1,
+              m0 = ~ x1 + u,
+              m1 = ~ x1 + u,
+              uname = 'u',
+              target = "ate",
+              criterion.tol = 0.01,
+              initgrid.nu = 5,
+              initgrid.nx = 5,
+              audit.nx = 5,
+              audit.nu = 5,
+              noisy = FALSE)

> ## Estimate without equality constraints
> result <- do.call(ivmte, args)
> result$gstar.coef
$min.g0
[m0](Intercept)          [m0]x1           [m0]u 
      0.2306349       0.4413234      -0.4149248 

$max.g0
[m0](Intercept)          [m0]x1           [m0]u 
     -0.4120800       0.2502855       0.6530651 

$min.g1
[m1](Intercept)          [m1]x1           [m1]u 
      0.6719583       0.0000000      -0.6287006 

$max.g1
[m1](Intercept)          [m1]x1           [m1]u 
      0.3618246       0.3101337       0.0000000 

> ## Impose equality constraint on x1 coefficient
> args$equal.coef <- ~ 0 + x1
> result <- do.call(ivmte, args)
> result$gstar.coef
$min.g0
[m0](Intercept)          [m0]x1           [m0]u 
      0.4260714       0.2458869      -0.5238869 

$max.g0
[m0](Intercept)          [m0]x1           [m0]u 
     -0.4120800       0.2777644       0.6296761 

$min.g1
[m1](Intercept)          [m1]x1           [m1]u 
      0.4260714       0.2458869      -0.1302407 

$max.g1
[m1](Intercept)          [m1]x1           [m1]u 
     0.37305113      0.27776443      0.02114278 

Note that if I instead declared equal.coef = ~ x1, I would have constrained both the intercept term and the coefficient on x1 to be equal for treated and control units.

These equality constraints may be imposed on any term common to both MTRs. If the constraint is imposed on a term that breaks out into many terms, e.g. a spline or a factor variable, then all the sub-terms will be constrained to have equal coefficients across the MTRs. For example, equal.coef = ~ factor(x2) will constrain the coefficients of each of the dummy variables generated to have common coefficients for treatment and control.

The equality constraints may also be imposed when the direct regression approach is used.

[a-torgovitsky]

Looks great, thanks! Could you make an issue reminding us to update the documentation to include both this feature and the direct regression approach? I can do that at some point.

[jkcshea]

Ah yes, I had completely forgotten about that. Sure, I'll post the issue. Also, I realized I forgot to update the GMM procedure to allow for these equality constraints, so I'll update that as well.

[jongohkim91]

That sounds like you are using x1 as an instrument. i.e. it is in the propensity score, but excluded from the outcome equation (MTE specification). That doesn't seem like what you want.

No we have a separate vector Z for instruments.

Let X= [x1 x2] and Z is a vector of instruments. (x1 is a vector of restricted variables and x2 is a vector of the unrestricted)

In the STATA MTE analysis, we use local IV approach and MTE(X=x, U=u) = x(\beta{1} - \beta{0}) + k(p) where k(p) = E[U{1} - U{0} | U=u] & p is a propensity score. (p ∈ P(X,Z)).

Hence x1, x2 & Z enters in the propensity score estimation process but Z is not included in the MTE function.

We have set a constraint for x1 to have beta values, \beta{1} = \beta{0}, and set k(p) to be a 2nd order polynomial. Thus, MTE = x2(\beta{1} - \beta{0}) + k(p)

In the MTR frame work, to estimate the same MTE, I initially thought I should set the MTR function as follows. MTR = a{d} + b{d}u + c{d}u^2 + e{d}x1 + f_{d}x2. Then, the derived MTE from MTR is MTEmtr = (a{1} - a{0}) + (b{1} - b{0})u + (c{1} - c{0})u^2 + (e{1} - e{0})x1 + (f{1} - f{0})x2 but with additional constraint e_{1} - e{0} = 0, making MTEmtr = (a{1} - a{0}) + (b{1} - b{0})u + (c{1} - c{0})u^2 + (f_{1} - f{0})x2

The reason why for such complicated MTR setting is that we are interested in seeking which variable has statistically significant effect on our outcome variable Y in the outcome equation: E[Y | X=x, U=u] = x\beta{0} + x(\beta{1} - \beta_{0})p + K(p) where K(p) is 3rd order polynomial.

I hope that made sense and my explanation is clear.

[jongohkim91]

@jkcshea Thank you very much for the update. I truly appreciate it!

[a-torgovitsky]

Ok, let us know if the functionality works the way you expected.

[jongohkim91]

Sorry for giving you such a late feedback. Since I have received great help from you all, I owe you a detailed explanation of how our MTR analysis went.

Before giving you details, let me explain how my research environment is. I am currently using a remote environment from the data center which we acquired our data from. Since there is no internet connection in the environment, I had to source all the original code files. I tried to install the github version locally but it was unsuccessful. I am waiting for the data center to install the github version on next Tuesday. Also lpSolveAPI is being used to solve LP problems because the data center is not eligible to apply for the free academic license from Gurobi.

Here is my experience so far:

1. The new update seems to work fine. I have not encountered any errors yet related to restricting the variables.

2. I have encountered a rather rare case. I estimated the bounds for ATE, ATT, & ATU in the setting of our research. In case of ATE & ATT, the package worked perfectly. It is not surprising because we have a large share of compliers(more than 70%) in our dataset. Bounds were produced and I was able to get all MTR coefficients. However, in the ATU case, I was only able to get the bounds but the upper bound was negative infinity(-1.0e30). On the other hand the lower bound was -0.22. By checking the messages produced by the ivmte function, I was able to realize that the first audit was unsuccessful because the LP model was unbounded. The package restarted the audit with new settings and then I was able to get the weird bounds. Maybe due to the fact the upper bound is negative infinity, many values in the returned list from ivmte() have NULL values. bounds, call.options, gstar, gstar.weights, moments, propensity, s.set are the only variables which I could get real values from. It did not return any specifics regarding the LP model. The lp.result list has simply a single NULL value.

3. This is more of a question. Is the package designed to stop the analysis if one bootstrap attempt fails to solve the LP even after maximum times of resampling? Since I was able to get the bounds, I tried to include bootstrapping procedure for 200 times with same specification but with smaller initial & audit grid for faster computation time. During the bootstrap procedure, the analysis got interrupted because the LP generated from one bootstrap attempt was unbounded. I am guessing the code stopped running because the LP was still unbounded after 3 times of resampling.

If you have any questions or need further clarification, please let me know.

[a-torgovitsky]

For problem 2) you could try setting the initial grid to be much larger using the initgrid.nx and initgrid.nu parameters.

[jkcshea]

Thanks for posting this @jongohkim91. Issues 2 and 3 should not be happening... Is it possible to post an example code? It will be much easier to diagnose if I can replicate the error, but I understand if the work isn't ready to be shared.

For 2:

• Can you show the full set of messages?
• Are any warnings returned?

For 3:

• Bootstraps should keep going even when the estimation procedure fails.
• Again, are there no error messages or warnings?

[jongohkim91]

In our research, we aim to measure the impact of different citizenship policies on immigrant children's educational outcome by taking advantage of the introduction of birthright citizenship in Germany in 2000. Generally immigrant children could be naturalized if one of their parents have lived in Germany for at least 8 years(parental eligibility) or children themselves stayed in Germany for at least years 8 years(individual eligibility) before the year 2000 reform. After the reform, immigrant children who were born in Germany and by one foreign-born parent who has legally resided in Germany for more than 8 years could become a German citizen when parents declare their child's birth.

We estimate the following MTE equation in STATA: MTE(X=x, U=u) = E[Y_{1} - Y_{0} | X=x, U=u] = x(β_{1} - β_{0}) + k(p) where Y is an educational outcome (grade repetition in this example), X is a vector of observables, p is the propensity score of a student becoming a German, and k(p) is a second-order polynomial function.

X includes gender, the country of origin, a linear and squared term of a child's age in months, a child's age when she immigrated, a linear and squared term of how long a parent has stayed in Germany in years, child cohort dummy variables, survey cohort dummies, and state dummy variables. Please note that we put additional constraint on X besides the child's gender and child's country of origin since we believe other observables do not have different values depending on the treatment status.

Therefore I have set up the ivmte function as below in order to estimate the same MTE equation from STATA in R:

result <- ivmte(data = dataset,
             m0 = ~1 + u + I(u^2) + female + factor(origin) + age_months + age_months2 +
               age_imm + ying + factor(yob) + factor(SC) + factor(state),
             m0 = ~1 + u + I(u^2) + Y1_female + factor(origin) + age_months + age_months2 + 
               age_imm + ying + factor(yob) + factor(SC) + factor(state),
             propensity = stata, #the propensity scores calculated by stata
             target = "atu", 
             ivlike = c(repetition ~ german + female  + factor(origin) + age_months + age_months2 + 
                          age_imm + ying + ying2 + factor(yob) + factor(SC) + factor(state),
                        repetition ~ german + female  + factor(origin) + age_months + age_months2 + 
                          age_imm + ying + ying2 + factor(yob) + factor(SC) + factor(state) |
                          z1 + z2 + z3 + female  + factor(origin) + age_months + age_months2 + 
                          age_imm + ying + ying2 + factor(yob) + factor(SC) + factor(state)),
             equal.coef = ~0 + age_months + age_months2 + age_imm + ying + ying2 + factor(yob) + factor(SC) + factor(state),
             treat = german,
             noisy = T,
             solver = "lpsolveapi"
             )

z1, z2, and z3 are birthright eligibility, parental eligibility, and individual eligibility respectively.

For 2: initgrid.nx and audit.nx were set to 3000. initgrid.nu and audit.nu were 200. No bootstrapping procedure was done.

• The full set of message was as follows: Performing audit procedure... Generating initial constraint grid...

  Audit count: 1 Minimum cirterion: 5.1936 Obtaining bounds... LP model was unbounded.

  Restarting audit with new settings: initgrid.nx = 3000 initgrid.nu = 200 Generating initial constraint grid...

  Audit count: 1 Minimum criterion: 51.936 Obtaining bounds... Bounds on the target parameter: [-0.2226752, -1.0000000e+30]

No error message was generated.

For 3: initgrid.nx, audit.nx, initgrid.nu and audit.nu were set to their default value. bootstrapp = 200. target was ATT.

• At the 110th bootstrap attempt,

Bootstrap iteration 110...
LP model was unbounded.

Restarting audit with new settings:
initgrid.nx = 30
initgrid.nu = 25
LP model was unbounded.

Restarting audit with new settings:
initgrid.nx = 45
initgrid.nu = 25
LP model was unbounded.

Restarting audit with new settings:
initgrid.nx = 68
initgrid.nu = 25
LP model was unbounded.
Called from: value[[3L]](cond)

• Yes, the error message printed on the console was:

error in abs(x) non-numeric argument to mathematical function It seems the error was generated during the tryCatch

If I write down all the rows in the Traceback pane,

-> eval(substitute(browser(skipCalls = pos), list(pos = (length(sys.frames()) -
stop(err)
value[[3L]](cond)
tryCatchOne(expr, names, parentenv, handlers[[1L]]
tryCatchList(expr, classes, parentenv, handlers[[1L]])
tryCatchList(expr, classes, parentenv, handlers)
tryCatch({

I hope this answers most of your questions. If you have more questions or need more details please let me know.