*! version 1.2.0 MLB 12Dec2013 program define mfxrcspline, rclass if c(stata_version) >= 11 { version 11 } else { version 10 } syntax [if] [in], /// [ /// at(str asis) /// link(str) /// manually specify link if e(cmd) isn't recognized CUSTomdydxb(string) /// manually specify a custom link ciopts(str asis) /// lineopts(str asis) /// level(passthru) /// EQuation(str asis) /// SHowknots /// YTItle(passthru) /// legend(passthru) /// addplot(str asis) /// extra overlaid graph GENerate(namelist max=3) /// noci /// * /// other -twoway options ] marksample touse local splines : char _dta[rcsplines] local oldvar : char _dta[oldvar] local knots : char _dta[knots] if "`splines'" == "" | "`oldvar'" == "" | "`knots'" == "" { di as err "adjustrcspline can only be used when the spline was created with mkspine2 with the cubic option" exit 198 } capture confirm variable `oldvar' if _rc { di as err "the original variable on which the spline is based needs to be present in the dataset" exit 198 } tempname b matrix `b' = e(b) local indepvars : colnames `b' local ok : list splines in indepvars if !`ok' { di as err "all variables created in the last call to mkspline2 must be" di as err "independent variables in the last estimation command" exit 198 } capture assert `:word count `generate'' != 2 if _rc { di as err "only 1 or 3 new variables can be specified in the generate option" exit 198 } if "`generate'" != "" { confirm new variable `generate' } if `"`equation'"' != "" { local eq `"[`equation']"' } tempname b matrix `b' = e(b) local indepvars : colnames `b' // Collect and check the link function if "`link'" != "" { // allow case insensitive link names local link = lower("`link'") // allow abreviations local l = max(4,length("`link'")) if "`link'" == substr("log", 1, `l') { local link = "probit" } else if "`link'" == substr("identity", 1, `l') { local link = "identity" } else if "`link'" == substr("logit", 1, `l') { local link = "logit" } else if "`link'" == substr("probit", 1, `l') { local link = "probit" } else if "`link'" == substr("logcomplement",1,`l') { local link = "logcomplement" } else if "`link'" == substr("loglog", 1, `l') { local link = "loglog" } else if "`link'" == substr("cloglog", 1, `l') { local link = "cloglog" } else if "`link'" == substr("reciprocal", 1, `l') { local link = "reciprocal" } else if "`: word 1 of `link''" == substr("power", 1, `l') { capture confirm numeric `: word 2 of `link'' if _rc { di as error "the second element specified in the link() option should be a number when specifying the power link " exit 198 } if `: word 2 of `link'' == 0 { local link "log" } else { local link = "power `: word 2 of `link''" } } else if "`: word 1 of `link''" == substr("opower", 1, `l') { capture confirm numeric `: word 2 of `link'' if _rc { di as error "the second element specified in the link() option should be a number when specifying the odds power link " exit 198 } if `: word 2 of `link'' == 0 { local link "logit" } else { local link = "opower `: word 2 of `link''" } } else { di as err "link `link' not recognized" exit 198 } } else if "`customdydxb'" == "" { if inlist("`e(cmd)'", "regress") | /// ("`e(cmd)'" == "glm" & "`e(linkt)'" == "Identity") { local link "identity" } if inlist("`e(cmd)'", "logit", /// "logistic", /// "betafit", /// "seqlogit") | /// ("`e(cmd)'" == "glm" & "`e(linkt)'" == "Logit") { local link "logit" } if inlist("`e(cmd)'", "probit") | /// ("`e(cmd)'" == "glm" & "`e(linkt)'" == "Probit") { local link = "probit" } if inlist("`e(cmd)'", "poisson", "nbreg") | /// ("`e(cmd)'" == "glm" & "`e(linkt)'" == "Log") { local link = "log" } if "`e(cmd)'" == "glm" & "`e(linkt)'" == "Log complement" { local link = "logcomplement" } if "`e(cmd)'" == "glm" & "`e(linkt)'" == "Log-log" { local link = "loglog" } if ("`e(cmd)'" == "glm" & "`e(linkt)'" == "Complementary log-log") | /// "`e(cmd)'" == "cloglog" { local link = "cloglog" } if "`e(cmd)'" == "glm" & "`e(linkt)'" == "Reciprocal" { local link = "reciprocal" } if "`e(cmd)'" == "glm" & strpos("`e(linkt)'", "Power") == 1 { local link = "power `e(power)'" } if "`e(cmd)'" == "glm" & "`e(linkt)'" == "Odds power" { local link = "opower `e(power)'" } } if "`link'`customdydxb'" == "" & "`e(cmd)'" != "glm" { di as err "last estimation command (`e(cmd)') was not recognized" di as err "specify the link() or customdydxb() option to manually specify the link function" exit 198 } if "`link'`customdydxb'" == "" & "`e(cmd)'" == "glm" { di as err "the `e(linkt)' link function used in the last glm command was not recognized" di as err "specify the link() or customdydxb() option to manually specify the link function" exit 198 } // dxb/dx local i = 1 local n : word count `knots' foreach k of local knots { local k`i++' = `k' } gettoken sp1 rest: splines local dxbdx "`eq'_b[`sp1']" local a1 = 1/((`k`n''-`k1')^2) local n1 = `n' - 1 local i = 1 foreach var of varlist `rest' { local a2 = (`k`n''-`k`i'')/((`k`n''-`k`n1'')*(`k`n''-`k1')^2) local a3 = (`k`n1''-`k`i'')/((`k`n''-`k`n1'')*(`k`n''-`k1')^2) local va "cond(`oldvar'-`k`i'' >0, 3*`eq'_b[`var']*`a1'*(`oldvar'-`k`i'' )^2, 0)" local vb "cond(`oldvar'-`k`n1''>0, 3*`eq'_b[`var']*`a2'*(`oldvar'-`k`n1'')^2, 0)" local vc "cond(`oldvar'-`k`n'' >0, 3*`eq'_b[`var']*`a3'*(`oldvar'-`k`n'' )^2, 0)" local dxbdx "`dxbdx' + `va' - `vb' + `vc'" local `i++' } // dy/dxb // set other covariates at the mean or as specified in at() if "`link'" != "identity" { if c(stata_version) >= 11 { foreach coef of local indepvars { _ms_parse_parts `coef' if inlist("`r(type)'", "variable", "factor") & !r(omit) local other "`other' `r(name)'" if "`r(type)'" == "factor" & !r(omit) local factor "`factor' `r(name)'" if "`r(type)'" == "interaction" { local intvar "" forvalues i = 1/`r(k_names)' { local intvar "`r(name`i')'" if `: list intvar in splines' { di as err "the spline variables may not be part of an interaction" exit 198 } } } } local cons "_cons" local other : list other - cons local factor : list uniq factor } else { local cons "_cons" local other : list indepvars - cons } local other : list uniq other local other : list other - splines // local notis created because "=" sign might be attached to variable name without space local notis : subinstr local at "=" " ", all local other : list other - notis quietly { // gen id var to merge the new vars into the original dataset tempvar id qui gen long `id' = _n preserve tempname atmat if `"`e(offset)'"' != "" { local xb `"(xb() + `e(offset)')"' tempvar offvar qui gen double `offvar' = `e(offset)' local expvar `e(offset)' // remove "ln(" and ")" in case -exposure()- was specified instead of offset local expvar1 : subinstr local expvar "ln(" "", all if "`expvar1'" != "`expvar'" { local exposure "exposure" } local expvar : subinstr local expvar1 ")" "", all if !`: list offvar in notis' { sum `offvar' if `touse', meanonly if "`exposure'" != "" { qui replace `expvar' = exp(r(mean)) matrix `atmat' = exp(r(mean)) } else { qui replace `expvar' = r(mean) matrix `atmat' = r(mean) } local atmatrow "`expvar'" } } else { local xb "xb()" } if "`other'" != "" { foreach var of local other { if `: list var in factor' { Leftmode `var' if `touse' tempvar t qui gen double `t' = r(lmode) qui drop `var' rename `t' `var' matrix `atmat' = nullmat(`atmat') \ `=r(lmode)' local atmatrow "`atmatrow' `var'" } else { sum `var' if `touse', meanonly tempvar t qui gen double `t' = r(mean) qui drop `var' rename `t' `var' matrix `atmat' = nullmat(`atmat') \ `=r(mean)' local atmatrow "`atmatrow' `var'" } } } local at : subinstr local at " =" "=", all local at : subinstr local at "= " "=", all local k_at : word count `at' tokenize `at' forvalues i = 1/`k_at' { gettoken var val : `i', parse("=") qui drop `var' qui gen double ``i'' matrix `atmat' = nullmat(`atmat') \ ``val'' local atmatrow "`atmatrow' `var'" } } } // link specific derivatives if "`link'" == "identity" { local dydxb "1" } if "`link'" == "logit" { local dydxb "invlogit(`xb')*invlogit(-`xb')" } if "`link'" == "probit" { local dydxb "normalden(`xb')" } if "`link'" == "log"{ local dydxb "exp(`xb')" } if "`link'" == "logcomplement" { local dydxb "-1*exp(`xb')" } if "`link'" == "loglog"{ local mu "exp(-1*exp(-1*(`xb')))" local dydxb "-1*(`mu')*ln(`mu')" } if "`link'" == "cloglog" { local mu "invcloglog(`xb')" local dydxb "(`mu'-1)*ln(1-`mu')" } if "`link'" == "reciprocal"{ local dydxb "-1*`xb'^-2" } if "`: word 1 of `link''" == "power" { local power : word 2 of `link' local mu "`xb'^(1/`power')" local dydxb "1/`power'*(`mu')^(1-`power')" } if "`: word 1 of `link''" == "opower" { local power : word 2 of `link' local mu "(1 + `power'*`xb')^(1/`power') / (1 + (1 + `power'*`xb')^(1/`power'))" local dydxb "((`mu')*(1-(`mu')))/(1+`power'*xb())" } if "`customdydxb'" != "" { local dydxb "`customdydxb'" } // dy/dx tempvar eff lb ub mis mark qui predictnl double `eff' = (`dydxb')*(`dxbdx') if `touse' , ci(`lb' `ub') `level' qui gen byte `mis' = missing(`oldvar') qui bys `oldvar' `mis' : gen byte `mark' = _n == 1 if `mis' == 0 if "`showknots'" != "" { local shknots "xline(`knots')" } if `"`ytitle'"' == "" { if c(stata_version) >= 11 { local ytitle `"ytitle("({&part} `e(depvar)')/({&part} `oldvar')")"' } else { local ytitle `"ytitle("(d `e(depvar)')/(d `oldvar')")"' } } if `"`legend'"' == "" { local legend "legend(nodraw)" } if "`ci'" == "" { twoway rarea `lb' `ub' `oldvar' if `mark', /// sort /// `ciopts' || /// line `eff' `oldvar' if `mark', /// sort /// lpattern(solid) /// `ytitle' /// `shknots' /// `lineopts' /// `legend' /// `options' /// || `addplot' } else { twoway line `eff' `oldvar' if `mark', /// sort /// lpattern(solid) /// `ytitle' /// `shknots' /// `lineopts' /// `legend' /// `options' /// || `addplot' } if "`link'" != "identity" & "`generate'" == ""{ qui restore } if "`link'" != "identity" & "`generate'" != ""{ qui keep `id' `eff' `lb' `ub' sort `id' tempfile newvars qui save `newvars' qui restore sort `id' tempvar merge qui merge `id' using `newvars', _merge(`merge') } if `: word count `generate'' == 1 { gen double `generate' = `eff' } if `: word count `generate'' == 3 { tokenize `generate' gen double `1' = `eff' gen double `2' = `lb' gen double `3' = `ub' } capture confirm matrix `atmat' if !_rc { matrix rownames `atmat' = `atmatrow' matrix colnames `atmat' = "value" return matrix atmat = `atmat' } end program define Leftmode, sort rclass syntax varname [if] marksample touse tempvar freq minus quietly { bys `touse' `varlist' : gen `freq' = sum(`touse')*`touse' by `touse' `varlist' : replace `freq' = (_n == _N)*`freq'[_N] gen `minus' = -`varlist' sort `touse' `freq' `minus' } return scalar lmode = `varlist'[_N] end