Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Help w/ precision when using locals vs. mata?

    Hi!
    I'm having some trouble getting variables to equal each other when using locals vs. mata (I think...). I'm trying to run a multiverse g-computation analysis where I run a bunch of different models and compare each "reality" to a "hypothetical" scenario. The code I have so far loops through different exposure/covariate combinations while storing the model specifications, beta estimates, standard errors, and some calculations along the way as locals each in a new row in the dataset in a set of placeholder variables. At the end of each loop, as a sanity check I would like to include an "assert" statement just to make sure the calculations are correct. However, the assert statement seems to fail each time because the values are slightly different (at like the 15th decimal point) even though I'm (theoretically) using the same variables/values in the assert statements as I did in the original calculation.

    I was hoping generating all variables as "double" was the solution, however, that didn't seem to solve the problem...
    Are there some slight precision differences between storing things as locals vs. mata? If so, is there a way around that?

    Best,
    David

    Code:
    * Example generated by -dataex-. For more info, type help dataex
    clear
    input byte(counterfactual cohort) str6 studyid byte wave double(executive verbal_memory) float(alcohol_consumption CDC_DGA_moderate_binary CDC_heavy_binary) double age byte sex
    1 1 "000011" 1       .       .  0 0 0 72 2
    0 1 "000011" 1   -.683   .2927  . . . 72 2
    1 1 "000011" 2       .       .  0 0 0 74 2
    0 1 "000011" 2  -.5531   .0734  0 1 0 74 2
    0 1 "000011" 3  -.2413  1.0408  0 1 0 75 2
    1 1 "000011" 3       .       .  0 0 0 75 2
    1 1 "000011" 4       .       .  0 0 0 77 2
    0 1 "000011" 4  -.8488  -.3437  0 1 0 77 2
    1 1 "000031" 1       .       .  0 0 0 78 2
    0 1 "000031" 1   1.795  2.0329  0 1 0 78 2
    0 1 "000031" 2  1.4843   .2244  0 1 0 80 2
    1 1 "000031" 2       .       .  0 0 0 80 2
    0 1 "000031" 3  1.8837   .6061  0 1 0 82 2
    1 1 "000031" 3       .       .  0 0 0 82 2
    1 1 "000031" 4       .       .  0 0 0 83 2
    0 1 "000031" 4  1.3239   .2622  0 1 0 83 2
    1 1 "000041" 1       .       .  0 0 0 87 1
    0 1 "000041" 1  -.5227 -1.4793  0 1 0 87 1
    1 1 "000041" 2       .       .  0 0 0 89 1
    0 1 "000041" 2  -.3191 -1.0472  0 1 0 89 1
    0 1 "000041" 3       .       .  . . .  . 1
    1 1 "000041" 3       .       .  0 0 0  . 1
    1 1 "000041" 4       .       .  0 0 0  . 1
    0 1 "000041" 4       .       .  . . .  . 1
    0 1 "000111" 1   .6994   .8505  . . . 68 1
    1 1 "000111" 1       .       .  0 0 0 68 1
    1 1 "000111" 2       .       .  0 0 0 70 1
    0 1 "000111" 2   .1954  1.0394  0 1 0 70 1
    0 1 "000111" 3  1.0773   .5815  0 1 0 71 1
    1 1 "000111" 3       .       .  0 0 0 71 1
    0 1 "000111" 4  1.1861   .8663  0 1 0 74 1
    1 1 "000111" 4       .       .  0 0 0 74 1
    0 1 "000151" 1   1.755   .3126  3 1 0 74 2
    1 1 "000151" 1       .       .  0 0 0 74 2
    1 1 "000151" 2       .       .  0 0 0 76 2
    0 1 "000151" 2  1.9684  -.0793  3 1 0 76 2
    0 1 "000151" 3  2.0677   .8574  3 1 0 77 2
    1 1 "000151" 3       .       .  0 0 0 77 2
    0 1 "000151" 4  1.1186  -.7781  3 1 0 80 2
    1 1 "000151" 4       .       .  0 0 0 80 2
    1 1 "000201" 1       .       .  0 0 0 73 2
    0 1 "000201" 1  1.1164  2.1331  6 0 0 73 2
    0 1 "000201" 2   .7542  1.4678  3 1 0 75 2
    1 1 "000201" 2       .       .  0 0 0 75 2
    1 1 "000201" 3       .       .  0 0 0 77 2
    0 1 "000201" 3   .5116  1.0017  3 1 0 77 2
    1 1 "000201" 4       .       .  0 0 0  . 2
    0 1 "000201" 4       .       .  . . .  . 2
    0 1 "000221" 1  -.9294  -.5685 .5 1 0 67 2
    1 1 "000221" 1       .       .  0 0 0 67 2
    0 1 "000221" 2 -1.0842  -.7542 .5 1 0 69 2
    1 1 "000221" 2       .       .  0 0 0 69 2
    1 1 "000221" 3       .       .  0 0 0 70 2
    0 1 "000221" 3 -1.2583   -.006 .5 1 0 70 2
    1 1 "000221" 4       .       .  0 0 0  . 2
    0 1 "000221" 4       .       .  . . .  . 2
    0 1 "000301" 1  1.4466  -.2091  0 1 0 84 2
    1 1 "000301" 1       .       .  0 0 0 84 2
    0 1 "000301" 2   .9178  -.2847  0 1 0 86 2
    1 1 "000301" 2       .       .  0 0 0 86 2
    1 1 "000301" 3       .       .  0 0 0 88 2
    0 1 "000301" 3  1.0672   .2509  0 1 0 88 2
    1 1 "000301" 4       .       .  0 0 0 89 2
    0 1 "000301" 4   .8554 -1.2142  0 1 0 89 2
    1 1 "000311" 1       .       .  0 0 0 80 2
    0 1 "000311" 1   .4885   .4862  0 1 0 80 2
    1 1 "000311" 2       .       .  0 0 0 82 2
    0 1 "000311" 2   .7413   .0611  0 1 0 82 2
    0 1 "000311" 3   .4701    .676  0 1 0 83 2
    1 1 "000311" 3       .       .  0 0 0 83 2
    0 1 "000311" 4    .614   .4418  0 1 0 85 2
    1 1 "000311" 4       .       .  0 0 0 85 2
    0 1 "000351" 1  1.4937   .4288  6 0 0 71 2
    1 1 "000351" 1       .       .  0 0 0 71 2
    0 1 "000351" 2  1.3187   .2989  6 0 0 73 2
    1 1 "000351" 2       .       .  0 0 0 73 2
    0 1 "000351" 3   1.328  -.0727  7 1 0 75 2
    1 1 "000351" 3       .       .  0 0 0 75 2
    0 1 "000351" 4  1.1858  -.7051  7 1 0 77 2
    1 1 "000351" 4       .       .  0 0 0 77 2
    1 1 "000361" 1       .       .  0 0 0 78 1
    0 1 "000361" 1 -1.6545 -1.5781  9 0 0 78 1
    1 1 "000361" 2       .       .  0 0 0 80 1
    0 1 "000361" 2 -1.4666  -1.415 21 0 1 80 1
    1 1 "000361" 3       .       .  0 0 0 81 1
    0 1 "000361" 3 -1.8817 -1.6825 28 0 1 81 1
    0 1 "000361" 4       .       .  . . .  . 1
    1 1 "000361" 4       .       .  0 0 0  . 1
    1 1 "000381" 1       .       .  0 0 0 86 2
    0 1 "000381" 1   .5282  1.3993  0 1 0 86 2
    1 1 "000381" 2       .       .  0 0 0 88 2
    0 1 "000381" 2   .4111  1.1987  0 1 0 88 2
    0 1 "000381" 3   .3363   .0842  0 1 0 89 2
    1 1 "000381" 3       .       .  0 0 0 89 2
    0 1 "000381" 4  -.0744  1.1068  . . . 89 2
    1 1 "000381" 4       .       .  0 0 0 89 2
    1 1 "000401" 1       .       .  0 0 0 72 1
    0 1 "000401" 1   .4049 -1.2391  6 1 0 72 1
    1 1 "000401" 2       .       .  0 0 0 74 1
    0 1 "000401" 2   .3778  -.8323  6 1 0 74 1
    end
    label values counterfactual EXPAND
    label def EXPAND 0 "Reality", modify
    label def EXPAND 1 "Counterfactual", modify
    label values cohort COHORT
    label def COHORT 1 "Cohort #1", modify
    label values wave WAVE
    label def WAVE 1 "Wave 1", modify
    label def WAVE 2 "Wave 2", modify
    label def WAVE 3 "Wave 3", modify
    label def WAVE 4 "Wave 4", modify
    label values CDC_DGA_moderate_binary ALCOHOL_BINARY
    label values CDC_heavy_binary ALCOHOL_BINARY
    label def ALCOHOL_BINARY 0 "No", modify
    label def ALCOHOL_BINARY 1 "Yes", modify
    label values sex SEX
    label def SEX 1 "Male", modify
    label def SEX 2 "Female", modify
    
             *********************************************************************************
             ***** GENERATING PLACEHOLDER VARIABLES FOR MODEL SPECIFICATIONS & ESTIMATES *****
             *********************************************************************************
    
    browse 
    sort counterfactual cohort studyid wave
    
    capture drop outcome exposure covariates mean_difference standard_error_difference
    
    generate outcome = ""
    generate exposure = ""
    generate covariates = ""
    generate random_effects = ""
    generate variance_covariance_matrix = ""
    
    generate double beta = .
    generate double standard_error = .
    generate double mean_difference = .
    
    generate double contrast = .
    generate double standard_error_contrast = .
    
                                ***********************************************
                                ***** RUNNING MULTIVERSE (LEVEL ANALYSIS) *****
                                ***********************************************
    STOP                                        
    local COUNTER 0
    local LAST_OBS = _N
    
    tuples c.age i.sex 
    foreach outcome of varlist executive verbal_memory {
        foreach exposure of varlist alcohol_consumption *binary {
            forvalues i = 1/`ntuples' {                            
                // INCREMENTING COUNTER
                local `COUNTER++'
                
                // RUNNING MODEL(S)
                mixed `outcome' `exposure' `tuple`i'' || studyid : , covariance(unstructured)
                
                // STORING MATRIX OF RESULTS
                matrix Adjusted = r(table)
                matrix list Adjusted
                
                // CHECKING FORMAT OF EXPOSURE VARIABLE FOR BRANCHING
                quietly tab `exposure' if counterfactual == 0
                quietly return list
                
                // BINARY | CONTINUOUS EXPOSURE
                if (r(r) == 2 | r(r) > 4) {
                    // BETA ESTIMATE FROM MODEL
                    quietly display Adjusted[1,1]
                    local beta = Adjusted[1,1]
                    
                    // STANDARD ERROR ESTIMATE FROM MODEL            
                    quietly display Adjusted[2,1]
                    local standard_error = Adjusted[2,1]
                }
                
                // PRESERVING DATASET
                preserve
            
                    // PREDICTING OUTCOME
                    predict double predicted_outcome, xb
                    
                    // DROPPING UNNECESSARY VARIABLES
                    keep counterfactual cohort studyid wave predicted_outcome            
                    
                    // RESHAPING DATA --> WIDE
                    reshape wide predicted_outcome, i(studyid wave) j(counterfactual)
                    rename (predicted_outcome0 predicted_outcome1) (predicted_reality predicted_counterfactual)
    
                    // GENERATING SLOPES (BY PERSON-WAVE) --> DIFFERENCE BETWEEN REALITY AND COUNTERFACTUAL 
                    generate double slope = (predicted_reality - predicted_counterfactual)
    
                    // CALCULATING AVERAGE SLOPE BY PERSON
                    by studyid (wave), sort: egen double slope_difference = mean(slope)
    
                    // CALCULATING AVERAGE SLOPE DIFFERENCE
                    quietly summarize slope_difference
                    quietly return list
    
                    local mean_difference = `r(mean)'
    
                    // CALCULATING CONTRAST/STANDARD ERROR (CONTRAST)
                    generate double contrast = (`beta' * `mean_difference')
                    generate double standard_error_contrast = (`standard_error' * `mean_difference')
    
                    // COLLAPSING DATASET
                    collapse (mean) contrast standard_error_contrast
    
                    // STORING ESTIMATES AS LOCALS
                    quietly summarize contrast
                    quietly return list
    
                    local contrast = `r(mean)'
    
                    quietly summarize standard_error_contrast
                    quietly return list
    
                    local standard_error_contrast = `r(mean)'
    
                // RESTORING DATASET
                restore
                
                // CHECKING ROW AVAILABILITY & ADDING ADDITIONAL OBSERVATIONS
                capture assert !missing(beta) in `LAST_OBS'
                if (_rc == 0) {                                                                            // NOTE: IF FINAL ROW IS TAKEN --> ADDS SINGLE OBSERVATION TO DATASET
                    insobs 1, after(_N)                                                                    // NOTE: ONLY ADDS ADDITIONAL OBSERVATIONS ONCE ORIGINAL COLUMN IS FULL
                }
    
                // STORING SPECIFICATIONS & ESTIMATES        
                replace outcome = "`outcome'" in `COUNTER'
                replace exposure = "`exposure'" in `COUNTER'
                replace covariates = "`tuple`i''" in `COUNTER'
                
                replace random_effects = substr("`e(cmdline)'", strpos("`e(cmdline)'", "||"), strpos("`e(cmdline)'", ",") - strpos("`e(cmdline)'", "||")) in `COUNTER'
                replace variance_covariance_matrix = substr("`e(cmdline)'", (strpos("`e(cmdline)'", ",") + 2), .) in `COUNTER'
                
                mata: st_store(`COUNTER', "beta", (`beta'))
                mata: st_store(`COUNTER', "standard_error", (`standard_error'))
                mata: st_store(`COUNTER', "mean_difference", (`mean_difference'))
                mata: st_store(`COUNTER', "contrast", (`contrast'))
                mata: st_store(`COUNTER', "standard_error_contrast", (`standard_error_contrast'))
    
                // DROPPING PREDICTED OUTCOME
                capture drop predicted_outcome        
                
                // CHECKS
                generate double contrast_check = (beta * mean_difference)
                assert (contrast == contrast_check)
                
                generate double SE_contrast_check = (standard_error * mean_difference)
                assert (standard_error_contrast == SE_contrast_check)
                
                capture drop contrast_check SE_contrast_check
            }
        }
    }

  • #2
    Local macros are essentially holders for strings. It so happens that numbers can be held as strings, and Stata programmers make heavy use of that convenience, especially for holding not too large integers. but there can be a loss of precision. Use scalars to hold individual numbers where precision is vital.

    Comment


    • #3
      To add to Nick’s excellent answer, I recently went down a similar rabbit hole trying to validate a statistical package. The author made heavy use of storing numbers in locals which amounted to results that were closer to float precision because of the fact that locals are actually string representations of numbers rather than actually held as numbers. It makes a difference where the decimals matter or when having to divide or multiply by small quantities.

      Comment


      • #4
        A dictum I recall from the old days when "soil was still rocks" is "Never compare floating point numbers for exact equality." Something like the following would be more robust to irrelevantly small differences:
        Code:
        assert abs(contrast - contrast_check)/contrast < 1e-6  // Close enough?

        Comment

        Working...
        X