Stata output language_2

library(tidyverse)
library(htmltools)
library(evaluate)
library("Statamarkdown")
library(githubinstall)
suppressMessages(gh_install_packages("konfound", ref = "newitcv_output", force = TRUE))

## Warning in fread(download_url, sep = "\t", header = FALSE, stringsAsFactors =
## FALSE, : Found and resolved improper quoting out-of-sample. First healed line
## 4848: <<Puriney honfleuR "Evening, honfleuR" by Seurat>>. If the fields are not
## quoted (e.g. field separator does not appear within any field), try quote="" to
## avoid this warning.

## Do you want to install the package (Y/n)?  
## 
## ── R CMD build ─────────────────────────────────────────────────────────────────
##   
   checking for file ‘/private/var/folders/52/my2sfsnn1qs78nfvl14zyynh0000gn/T/RtmpL4XKEG/remotesfe4a49d24b50/konfound-project-konfound-0aa960d/DESCRIPTION’ ...
  
✔  checking for file ‘/private/var/folders/52/my2sfsnn1qs78nfvl14zyynh0000gn/T/RtmpL4XKEG/remotesfe4a49d24b50/konfound-project-konfound-0aa960d/DESCRIPTION’
## 
  
─  preparing ‘konfound’:
## 
  
   checking DESCRIPTION meta-information ...
  
✔  checking DESCRIPTION meta-information
## 
  
─  checking for LF line-endings in source and make files and shell scripts
## 
  
─  checking for empty or unneeded directories
## 
  
─  building ‘konfound_0.4.0.tar.gz’
## 
  
   
## 

library(konfound)

0) Summary

• The updated R konfound command output language has been integrated into STATA.
• The output language task is finished, but there are a few minor tasks to complete within the next two weeks

1. Add signsuppresion & eff_thr options in STATA version (which requires integration with Xuesen’s recent update in pkonfound_v2 model 1)
2. Develop an updated loop to address discrepancies in output values (e.g., RIR, p-value) between R and STATA

1) STATA Output Language for RIR in Linear Model (model 0 in STATA)

• It appears that the STATA version output simultaneously displays both RIR and ITCV.
• User-entered effect threshold (eff_thr != NA) in Stata should be updated soon to align with Xuesen’s updated version.

1.1) Invalidate

• There is a little difference in RIR value (= 6) in Stata and R

// Stata, eff_thr = NA
pkonfound_v2 -0.3 0.01 5000 10, nu(0.00) sig(0.05) onetail(0) model_type(0) switch_trm(1) 

pkonfound is working
------------------
Impact Threshold for Omitted Variable

An omitted variable would have to be correlated at 0.611 with the outcome and at -.611 with the predictor
of interest (conditioning on observed covariates. Signs are interchangeable) to invalidate an inference.
Correspondingly the impact of an omitted variable (as defined in Frank 2000) must be
0.611 x -.611=-0.3736 to invalidate an inference.
------------------
The Threshold for % Bias to Invalidate/Sustain the Inference

You entered an estimated effect of -0.3. To invalidate the inference of an effect using the threshold
of -.02 for statistical significance with alpha = .05, 93.33% of the (-0.3) estimate would have
to be due to bias. This implies that to invalidate the inference one would expect to have to
replace 4667 (93.33%) observations with cases for which the treatment effect is 0 (RIR = 4667)

Note:
For non-linear models, the impact threshold should not be used.
The % bias calculation is based on the original coefficient,
compare with the use of average partial effects as in the [konfound] command.

See Frank et al. (2013) for a description of the method.

Citation: Frank, K.A., Maroulis, S., Duong, M., and Kelcey, B. (2013).
What would it take to change an inference?
Using Rubin's causal model to interpret the robustness of causal inferences.
Education, Evaluation and Policy Analysis, 35 437-460.

### R, invalidate
## eff_thr = NA
pkonfound(est_eff = -0.3, std_err = 0.01, n_obs = 5000, n_covariates = 10, alpha = .05, tails = 2, index = "RIR", nu = 0, model_type = "ols", eff_thr = NA, to_return = "print")

Robustness of Inference to Replacement (RIR):
TO INVALIDATE:

RIR = 4673

You entered an estimated effect of -0.3. To invalidate
the inference of an effect using the threshold of -0.02 for
statistical significance with alpha = 0.05, 93.465% of the
(-0.3) estimate would have to be due to bias. This implies
that to invalidate the inference one would expect to have to
replace 4673 (93.465%) observations with cases for which the
treatment effect is 0 (RIR = 4673).

See Frank et al. (2013) for a description of the method.

Citation: Frank, K.A., Maroulis, S., Duong, M., and Kelcey, B. (2013).
What would it take to change an inference?
Using Rubin's causal model to interpret the robustness of causal inferences.
Education, Evaluation and Policy Analysis, 35 437-460.

For other forms of output, run ?pkonfound and inspect the to_return argument

For models fit in R, consider use of konfound().

1.2) Sustain

• There is no significant difference between Stata and R

// Stata, eff_thr = NA
pkonfound_v2 0.003 0.1 5000 10, nu(0.00) sig(0.05) onetail(0) model_type(0) switch_trm(1) 

pkonfound is working
------------------
Your estimate is not statistically significant.
------------------
------------------
Impact Threshold for Omitted Variable

An omitted variable would have to be correlated at 0.163 with the outcome and at -.163 with the predictor
of interest (conditioning on observed covariates. Signs are interchangeable) to sustain an inference.
Correspondingly the impact of an omitted variable (as defined in Frank 2000) must be
0.163 x -.163=-0.0266 to sustain an inference.
------------------
The Threshold for % Bias to Invalidate/Sustain the Inference

You entered an estimated effect of 0.003. To invalidate the inference of an effect using the threshold
of .196 for statistical significance with alpha = .05, -6433.33% of the (0.003) estimate would have
to be due to bias. This implies that to invalidate the inference one would expect to have to
replace 4924 (98.47%) observations with cases for which the treatment effect is 0 (RIR = 4924)

Note:
For non-linear models, the impact threshold should not be used.
The % bias calculation is based on the original coefficient,
compare with the use of average partial effects as in the [konfound] command.

See Frank et al. (2013) for a description of the method.

Citation: Frank, K.A., Maroulis, S., Duong, M., and Kelcey, B. (2013).
What would it take to change an inference?
Using Rubin's causal model to interpret the robustness of causal inferences.
Education, Evaluation and Policy Analysis, 35 437-460.

### R, sustain
## eff_thr = NA
pkonfound(est_eff = 0.003, std_err = 0.1, n_obs = 5000, n_covariates = 10, alpha = .05, tails = 2, index = "RIR", nu = 0, model_type = "ols", eff_thr = NA, to_return = "print")

Robustness of Inference to Replacement (RIR):
TO SUSTAIN:

RIR = 4923

You entered an estimated effect of 0.003. The threshold value for
statistical significance is 0.196 (alpha = 0.05). To reach that
threshold, 98.47% of the (0.003) estimate would have to be due to
bias. This implies that to sustain an inference one would expect to
have to replace 4923 (98.47%) observations with effect of 0 with
cases with effect of 0.196 (RIR = 4923).

See Frank et al. (2013) for a description of the method.

Citation: Frank, K.A., Maroulis, S., Duong, M., and Kelcey, B. (2013).
What would it take to change an inference?
Using Rubin's causal model to interpret the robustness of causal inferences.
Education, Evaluation and Policy Analysis, 35 437-460.

For other forms of output, run ?pkonfound and inspect the to_return argument

For models fit in R, consider use of konfound().

2) STATA & R Output Language for RIR in Non-linear Model (model 1 in STATA)

2.1) Invalidate, changeSE = TRUE

// Case 2-1: Stata, changeSE = TRUE, Invalidate
pkonfound_v2 -0.3 0.01 5000 0 2500, nu(0.00) sig(0.05) onetail(0) model_type(1) switch_trm(1) replace(1)

RIR = 221

The table implied by the parameter estimates and sample sizes you entered:

             |      Fail    Success 
-------------+----------------------
     Control |      1156            
   Treatment |      1344       1156 

The reported effect size =-0.300, SE = 0.010, p-value = 0.000.
The SE has been adjusted to 0.057 to generate real numbers in the implied table.
Numbers in the table cells have been rounded to integers, which may slightly alter
the estimated effect from the value originally entered.

To invalidate the inference that the effect is different from 0 (alpha = .05)
one would need to replace 221 (16.44%) treatment failure cases
with cases for which the probability of failure in the control group (46.24%) applies (RIR = 221).
This is equivalent to transferring 119 cases from treatment failure to treatment success
(Fragility = 119).

Note that RIR = Fragility/[1-P(failure in the control group)]

The transfer of 119 cases yields the following table:

             |      Fail    Success 
-------------+----------------------
     Control |      1156       1344 
   Treatment |      1225       1275 
Effect size =-0.111, SE = 0.057, p-value = 0.051.
This is based on t = estimated effect/standard error

## Case 2-1: R, changeSE = TRUE, Invalidate
pkonfound(-0.3, 0.01, 5000, n_covariates = 0, alpha = .05, tails = 2, nu = 0, n_treat = 2500, switch_trm = TRUE, model_type = "logistic", to_return = "print")

RIR = 220

The table implied by the parameter estimates and sample sizes you entered:

          Fail Success
Control   1157    1343
Treatment 1344    1156

The reported effect size = -0.300, SE = 0.010, p-value = 0.000. 
The SE has been adjusted to 0.057 to generate real numbers in the 
implied table. Numbers in the table cells have been rounded 
to integers, which may slightly alter the estimated effect from 
the value originally entered.

To invalidate the inference that the effect is different from 0 
(alpha = 0.050) one would need to replace 220 (16.37%) treatment failure 
cases with cases for which the probability of failure in the control 
group (46.28%) applies (RIR = 220). This is equivalent to transferring 
118 cases from treatment failure to treatment success (Fragility = 118).

Note that RIR = Fragility/[1-P(failure in the control group)]

The transfer of 118 cases yields the following table:

          Fail Success
Control   1157    1343
Treatment 1226    1274
Effect size = -0.111, SE = 0.057, p-value = 0.051. 
This is based on t = estimated effect/standard error

For other forms of output, run ?pkonfound and inspect the to_return argument

For models fit in R, consider use of konfound().

2.2) Invalidate, changeSE = FALSE

// Case 2-2: Stata, changeSE = FALSE, Invalidate
pkonfound_v2 -0.3 0.1 5000 10 2500, nu(0.00) sig(0.05) onetail(0) model_type(1) switch_trm(1) replace(1)

RIR = 207

The table implied by the parameter estimates and sample sizes you entered:

             |      Fail    Success 
-------------+----------------------
     Control |      2246        254 
   Treatment |      2307        193 

The reported effect size =-0.300, SE = 0.100, p-value = 0.003.
Values have been rounded to the nearest integer. This may cause little change
to the estimated effect for the table.

To invalidate the inference that the effect is different from 0 (alpha = .05)
one would need to replace 207 (8.970000000000001%) treatment failure cases
with cases for which the probability of failure in the control group (89.84%) applies (RIR = 207).
This is equivalent to transferring 21 cases from treatment failure to treatment success
(Fragility = 21).

Note that RIR = Fragility/[1-P(failure in the control group)]

The transfer of 21 cases yields the following table:

             |      Fail    Success 
-------------+----------------------
     Control |      2246        254 
   Treatment |      2286        214 
Effect size =-0.189, SE = 0.097, p-value = 0.052.
This is based on t = estimated effect/standard error

## Case 2-2: R, changeSE = FALSE, Invalidate
pkonfound(-0.3, 0.1, 5000, n_covariates = 10, alpha = .05, tails = 2, n_treat = 2500, switch_trm = TRUE, model_type = "logistic", to_return = "print")

RIR = 207

The table implied by the parameter estimates and sample sizes you entered:

          Fail Success
Control   2246     254
Treatment 2307     193

The reported effect size = -0.300, and SE = 0.100, p-value = 0.003. 
Values have been rounded to the nearest integer. This may cause 
a little change to the estimated effect for the table.

To invalidate the inference that the effect is different from 0 
(alpha = 0.050) one would need to replace 207 (8.97%) treatment failure 
cases with cases for which the probability of failure in the control 
group (89.84%) applies (RIR = 207). This is equivalent to transferring 
21 cases from treatment failure to treatment success (Fragility = 21).

Note that RIR = Fragility/[1-P(failure in the control group)]

The transfer of 21 cases yields the following table:

          Fail Success
Control   2246     254
Treatment 2286     214
Effect size = -0.189, SE = 0.097, p-value = 0.052. 
This is based on t = estimated effect/standard error

For other forms of output, run ?pkonfound and inspect the to_return argument

For models fit in R, consider use of konfound().

2.3) Sustain, changeSE = TRUE, RIR exceeds 100%

• There is a slight difference in p-value in implied table between R (=.779) and Stata (=.930)
• This difference could be significant.
• The calculation of the p-value is identical in both R and Stata. However, it appears that the t_start values in R and Stata are slightly different, which might explain the observed discrepancy.
• Currently woring on debug

// Case 2-3: Stata, changeSE = TRUE, Sustain
pkonfound_v2 0.05 0.3 50 5 25, nu(0.00) sig(0.01) onetail(0) model_type(1) switch_trm(1) replace(1)

RIR = 20

The table implied by the parameter estimates and sample sizes you entered:

             |      Fail    Success 
-------------+----------------------
     Control |        13            
   Treatment |        12         13 

The reported effect size = 0.050, SE = 0.300, p-value = 0.930.
The SE has been adjusted to 0.566 to generate real numbers in the implied table.
Numbers in the table cells have been rounded to integers, which may slightly alter
the estimated effect from the value originally entered.

To reach the threshold that would sustain an inference that the effect is different
from 0 (alpha = .01) one would need to replace 20 (166.67%) treatment failure cases
with cases for which the probability of failure in the control group (52%) applies (RIR = 20).
This is equivalent to transferring 9 cases from treatment failure to treatment success
(Fragility = 9).

Note that RIR = Fragility/[1-P(failure in the control group)]

Note the RIR exceeds 100%. Generating the transfer of 9 cases would
require replacing more cases than are in the treatment failure condition.

The transfer of 9 cases yields the following table:

             |      Fail    Success 
-------------+----------------------
     Control |        13         12 
   Treatment |         3         22 
Effect size = 2.072, SE = 0.734, p-value = 0.007.
This is based on t = estimated effect/standard error

## Case 2-3: R, changeSE = TRUE, Sustain
pkonfound(0.05, 0.3, 50, n_covariates = 5, alpha = .01, tails = 2, nu = 0, n_treat = 25, switch_trm = TRUE, model_type = "logistic", to_return = "print")

RIR = 19

The table implied by the parameter estimates and sample sizes you entered:

          Fail Success
Control     13      12
Treatment   12      13

The reported effect size = 0.050, SE = 0.300, p-value = 0.779. 
The SE has been adjusted to 0.566 to generate real numbers in the 
implied table. Numbers in the table cells have been rounded 
to integers, which may slightly alter the estimated effect from 
the value originally entered.

To reach the threshold that would sustain an inference that the 
effect is different from 0 (alpha = 0.010) one would need to replace 19 
(158.33%) treatment failure cases with cases for which the probability of 
failure in the control group (52.00%) applies (RIR = 19). This is equivalent 
to transferring 9 cases from treatment failure to treatment success
(Fragility = 9).

Note that RIR = Fragility/[1-P(failure in the control group)]

Note the RIR exceeds 100%. Generating the transfer of 9 cases would
require replacing more cases than are in the treatment failure condition.

The transfer of 9 cases yields the following table:

          Fail Success
Control     13      12
Treatment    3      22
Effect size = 2.072, SE = 0.734, p-value = 0.007. 
This is based on t = estimated effect/standard error

For other forms of output, run ?pkonfound and inspect the to_return argument

For models fit in R, consider use of konfound().

2.4) Sustain, changeSE = FALSE

• There is a similar slight difference in p-value in implied table between R (=.841) and Stata (=.881)
• It might or might not be significant. But debug should be conducted.

// Case 2-4: Stata, changeSE = FALSE, Sustain
pkonfound_v2 -0.03 0.2 500 10 250, nu(0.00) sig(0.05) onetail(0) model_type(1) switch_trm(1) replace(1)

RIR = 25

The table implied by the parameter estimates and sample sizes you entered:

             |      Fail    Success 
-------------+----------------------
     Control |       180         70 
   Treatment |       182         68 

The reported effect size =-0.030, SE = 0.200, p-value = 0.881.
Values have been rounded to the nearest integer. This may cause little change
to the estimated effect for the table.

To reach the threshold that would sustain an inference that the effect is different
from 0 (alpha = .05) one would need to replace 25 (36.76%) treatment success cases
with cases for which the probability of failure in the control group (72%) applies (RIR = 25).
This is equivalent to transferring 17 cases from treatment success to treatment failure
(Fragility = 17).

Note that RIR = Fragility/[1-P(success in the control group)]

The transfer of 17 cases yields the following table:

             |      Fail    Success 
-------------+----------------------
     Control |       180         70 
   Treatment |       199         51 
Effect size =-0.417, SE = 0.211, p-value = 0.049.
This is based on t = estimated effect/standard error

## Case 2-4: R, changeSE = FALSE, Sustain
pkonfound(-0.03, 0.2, 500, n_covariates = 10, alpha = .05, tails = 2, nu = 0, n_treat = 250, switch_trm = TRUE, model_type = "logistic", to_return = "print")

RIR = 24

The table implied by the parameter estimates and sample sizes you entered:

          Fail Success
Control    180      70
Treatment  182      68

The reported effect size = -0.030, and SE = 0.200, p-value = 0.841. 
Values have been rounded to the nearest integer. This may cause 
a little change to the estimated effect for the table.

To reach the threshold that would sustain an inference that the 
effect is different from 0 (alpha = 0.050) one would need to replace 24 
(35.29%) treatment success cases with cases for which the probability of 
failure in the control group (72.00%) applies (RIR = 24). This is equivalent 
to transferring 17 cases from treatment success to treatment failure
(Fragility = 17).

Note that RIR = Fragility/[1-P(success in the control group)]

The transfer of 17 cases yields the following table:

          Fail Success
Control    180      70
Treatment  199      51
Effect size = -0.417, SE = 0.211, p-value = 0.049. 
This is based on t = estimated effect/standard error

For other forms of output, run ?pkonfound and inspect the to_return argument

For models fit in R, consider use of konfound().

2.5) needtworow = TRUE, Invalidate, changeSE = FALSE

// Case 2-5: Stata, changeSE = FALSE, needtworow = T, Invalidate, right
pkonfound_v2 0.25 0.5 70 20 20, nu(0.00) sig(0.001) onetail(0) model_type(1) switch_trm(1) replace(1)

RIR = 36

The table implied by the parameter estimates and sample sizes you entered:

             |      Fail    Success 
-------------+----------------------
     Control |        27         23 
   Treatment |        10         10 

The reported effect size = 0.250, SE = 0.500, p-value = 0.639.
The SE has been adjusted to 0.530 to generate real numbers in the implied table.
Numbers in the table cells have been rounded to integers, which may slightly alter
the estimated effect from the value originally entered.

The inference cannot be sustained merely by switching cases in only one treatment
condition. Therefore, cases have been switched from treatment failure to treatment success
and from control success to control failure. The final Fragility (= 16) and RIR (= 36)
reflect both sets of changes.

Please compare the after transfer table with the implied table.

             |      Fail    Success 
-------------+----------------------
     Control |        35         15 
   Treatment |         1         19 
Effect size = 3.813, SE = 1.072, p-value = 0.001.
This is based on t = estimated effect/standard error

## Case 3: R, needtworows = T, Sustain, changeSE = T, 
pkonfound(0.25, 0.5, 70, n_covariates = 20, alpha = .001, tails = 2, nu = 0, n_treat = 20, switch_trm = TRUE, model_type = "logistic", to_return = "print")

RIR = 37

The table implied by the parameter estimates and sample sizes you entered:

          Fail Success
Control     27      23
Treatment   10      10

The reported effect size = 0.250, SE = 0.500, p-value = 0.763. 
The SE has been adjusted to 0.530 to generate real numbers in the 
implied table. Numbers in the table cells have been rounded 
to integers, which may slightly alter the estimated effect from 
the value originally entered.

The inference cannot be sustained merely by switching cases in
only one treatment condition. Therefore, cases have been switched from
treatment failure to treatment success and from control success to control failure.
The final Fragility(= 17) and RIR(= 37) reflect both sets of changes. 
Please compare the after transfer table with the implied table.

          Fail Success
Control     35      15
Treatment    1      19
Effect size = 3.792, SE = 1.071, p-value = 0.001. 
This is based on t = estimated effect/standard error

For other forms of output, run ?pkonfound and inspect the to_return argument

For models fit in R, consider use of konfound().

3) tkonfound (model 2 in STATA)

3.1) chisq p

• Due to an error in R markdown, please refer to the figure that directly obtain the output from STATA for the chi-squared p-value option.

// chisq p
pkonfound_v2 35 17 17 38, sig(0.05) model_type(2) switch_trm(1) replace(0) test1(1)

Figure

3.2) fisher p

// fisher p
pkonfound_v2 35 17 17 38, sig(0.05) model_type(2) switch_trm(1) replace(0) test1(0)

tkonfound is working

Background Information:
This function calculates the number of cases that would have to be replaced 
with zero effect cases (RIR) to invalidate an inference made about the association 
between the rows and columns in a 2x2 table. 
One can also interpret this as switches from one cell to another, such as from 
the treatment success cell to the treatment failure cell.

Conclusion:
To invalidate the inference, one would need to replace 14 treatment success
cases for which the probability of failure in the control group applies (RIR = 14).
This is equivalent to transferring 9 cases from treatment success to treatment failure.

For the User-entered Table, the estimated odds ratio is 4.602, with p-value of 0

             |      Fail    Success 
-------------+----------------------
     Control |        35            
   Treatment |        17         38 

For the Transfer Table, the estimated odds ratio is 2.296 with p-value of .051

             |      Fail    Success 
-------------+----------------------
     Control |        35         17 
   Treatment |        26         29 

RIR = 14

4) Calculation of p-value in Implied/Transferred table

• Regardless of the research context (whether it involves upper or lower p-value tests), the ttail() function in STATA directly calculates the area under the curve.
• In R, the default option for the pt() command is a left-tailed test, which might be problematic for situations requiring a right-tailed test.
• Therefore, using pt(q, df, lower.tail = FALSE) in R eliminates concerns regarding the research context, similar to the ttail() function in STATA.
• This approach ensures that the calculation of p-values aligns well between STATA and R.

## previous p-calculation
  if (tails == 2) {
    p_start <- (2 * (1 - pt(abs(final_solution$t_start), n_obs - n_covariates - 2)))
    p_final <- (2 * (1 - pt(abs(final_solution$t_final), n_obs - n_covariates - 2)))
  } 
  # For a one-tailed test (assuming upper tail)
  else if (tails == 1 & t_obs > 0) {
    p_start <- (1 - pt(final_solution$t_start, n_obs - n_covariates - 2))
    p_final <- (1 - pt(final_solution$t_final, n_obs - n_covariates - 2))
  } 
  # For a one-tailed test (assuming lower tail)
  else if (tails == 1 & t_obs < 0) {
    p_start <- (pt(final_solution$t_start, n_obs - n_covariates - 2))
    p_final <- (pt(final_solution$t_final, n_obs - n_covariates - 2))
  }       

## updated p-calculation
    if (tails == 2) {
      p_start <- 2 * pt(abs(final_solution$t_start), n_obs - n_covariates - 2, lower.tail = FALSE)
      p_final <- 2 * pt(abs(final_solution$t_final), n_obs - n_covariates - 2, lower.tail = FALSE)
        }
        else if (tails == 1) {
            p_start = pt(abs(final_solution$t_start), n_obs - n_covariates - 2, lower.tail = FALSE)
            p_final = pt(abs(final_solution$t_final), n_obs - n_covariates - 2, lower.tail = FALSE)
            }