* tuning_frames.do
*
* Table 5 of "Which H2O frame should you use to explain your model?"
*
* Three models fit to the same 8 data draws: a deliberately overfit one, a
* hand-regularized one, and one tuned by random grid search. Reports the
* generalization gap and x5's share of within-frame permutation importance
* on the training and testing frames.
*
* Importance is expressed as a share of the within-frame total, because with
* reps() > 1 the raw values are on the loss scale and the frames do not have
* the same column total.

clear all
set more off
local R 8

capture program drop makedata
program define makedata
    args n seed
    clear
    set seed `seed'
    set obs `n'
    gen double x1 = runiform(-3,3)
    gen double x2 = rnormal()
    gen double x3 = 0.9*x2 + sqrt(1-0.9^2)*rnormal()
    gen double x4 = runiform(-2,2)
    gen byte   g  = runiform()<0.5
    gen double x5 = runiform()
    gen byte   x6 = runiform()<0.5
    gen double y  = 4*sin(1.5*x1) + 2*x2 + 3*x4*(2*g-1) + rnormal(0,1.5)
end

* grab MSE / R2 on a named frame for the model currently in e()
capture program drop framemse
program define framemse, rclass
    args frm
    quietly h2omlpostestframe `frm'
    quietly h2omlestat metrics
    matrix M = r(metrics)
    return scalar mse = M[2,1]
    return scalar r2  = M[6,1]
end

h2o init
tempfile res
postfile pf int rep str5 model str2 frame str2 vname double imp ///
    double mse_tr double mse_te using `res'

forvalues r = 1/`R' {
    makedata 4000 `=20261001 + `r''
    quietly _h2oframe put, into(d`r')
    quietly _h2oframe split d`r', into(dt`r' dv`r' ds`r') split(0.6 0.2 0.2) rseed(19)
    quietly _h2oframe change dt`r'

    forvalues m = 1/3 {
        if `m'==1 {                                   // A: deliberately overfit
            local mname "over"
            quietly h2oml gbregress y x1 x2 x3 x4 g x5 x6, validframe(dv`r') ///
                ntrees(500) maxdepth(12) minobsleaf(1) lrate(0.1) h2orseed(19)
        }
        else if `m'==2 {                              // B: hand-regularized
            local mname "reg"
            quietly h2oml gbregress y x1 x2 x3 x4 g x5 x6, validframe(dv`r') ///
                ntrees(500) maxdepth(4) minobsleaf(20) lrate(0.05) ///
                samprate(0.8) predsamprate(0.8) ///
                stop(5, metric(mse)) scoreevery(1) h2orseed(19)
        }
        else {                                        // C: tuned by grid search
            local mname "tuned"
            quietly h2oml gbregress y x1 x2 x3 x4 g x5 x6, validframe(dv`r') ///
                ntrees(300) maxdepth(2 4 6 8) lrate(0.02 0.05 0.1) ///
                minobsleaf(5 20) stop(5, metric(mse)) scoreevery(1) ///
                h2orseed(19) ///
                tune(metric(mse) grid(random, h2orseed(20)) ///
                     maxmodels(10) nooutput)
        }

        framemse _train
        local mtr = r(mse)
        framemse ds`r'
        local mte = r(mse)

        * equal samples() on both frames -> isolates bias, not sample size
        foreach fr in dt ds {
            quietly h2omlgraph permimp, frame(`fr'`r') reps(10) samples(700) ///
                h2orseed(19) savedata(_pi, replace) nodraw
            preserve
                quietly use _pi, clear
                quietly collapse (mean) importance, by(variable)
                local tag = cond("`fr'"=="dt","tr","te")
                forvalues i = 1/`=_N' {
                    post pf (`r') ("`mname'") ("`tag'") (variable[`i']) ///
                        (importance[`i']) (`mtr') (`mte')
                }
            restore
        }
        di "rep `r' model `mname' : trainMSE=`mtr' testMSE=`mte'"
    }
}
postclose pf

use `res', clear
save tunedres, replace

di _n "===== Generalization gap by model (mean over `R' reps) ====="
preserve
    quietly keep if frame=="tr" & vname=="x1"
    table model, statistic(mean mse_tr) statistic(mean mse_te) nformat(%9.4f)
restore

di _n "===== permimp: share of within-frame total importance (%) ====="
bysort rep model frame: egen double tot = total(imp)
gen double share = 100*imp/tot
table vname model frame, statistic(mean share) nformat(%7.2f)

di _n "===== NOISE ATTRIBUTION: x5 train/test share ratio by model ====="
preserve
    keep if vname=="x5"
    keep rep model frame share
    quietly reshape wide share, i(rep model) j(frame) string
    gen double ratio = sharetr/sharete
    di "  model   x5 share TRAIN   x5 share TEST   ratio   reps tr>te"
    levelsof model, local(ms)
    foreach mm in `ms' {
        quietly summarize sharetr if model=="`mm'", meanonly
        local a = r(mean)
        quietly summarize sharete if model=="`mm'", meanonly
        local b = r(mean)
        quietly summarize ratio if model=="`mm'", meanonly
        local c = r(mean)
        quietly count if model=="`mm'" & sharetr>sharete
        local w = r(N)
        quietly count if model=="`mm'"
        di "   `mm'  " %12.3f `a' "  " %14.3f `b' "  " %8.2f `c' "     `w'/" r(N)
    }
restore

di _n "===== same for x3 (correlated proxy) and x6 (binary noise) ====="
foreach v in x3 x6 {
    preserve
        keep if vname=="`v'"
        keep rep model frame share
        quietly reshape wide share, i(rep model) j(frame) string
        gen double ratio = sharetr/sharete
        di "  `v':"
        levelsof model, local(ms)
        foreach mm in `ms' {
            quietly summarize sharetr if model=="`mm'", meanonly
            local a = r(mean)
            quietly summarize sharete if model=="`mm'", meanonly
            local b = r(mean)
            quietly summarize ratio if model=="`mm'", meanonly
            local c = r(mean)
            di "     `mm'   train=" %7.3f `a' "  test=" %7.3f `b' "  ratio=" %6.2f `c'
        }
    restore
}