*Modeling how PAC beta changes without versus with controlling for MAC
*Assuming no MAC effect
{
clear
program drop changePAC_w_MAC_noMAC // this line causes problems ion first run and needs to be omitted 
program define changePAC_w_MAC_noMAC, rclass
	version 14.1
	drop _all
	matrix input m = (1, 0.101473869, 1, 0.076398164, 0.752885097, 1) // making matrix with set correlation structure
	drawnorm TL F M, n(3241) corr(m) cstorage(lower) clear 
	*drawnorm above draws a normal distrubtion of variables with above m correlation structure and n sample size
	*corr * // to double check correlation structure entered correctly--increase n to gain more precision - silenced to not slowdown simulation 
	reg TL F
	gen F1=_b[F]
	reg TL F M
	gen F2_minus_F1=_b[F]-F1
	sum F2_minus_F1 
	return scalar mean = r(mean)
	end

timer clear 1 // timer lines are just to see how long simulations take to run
timer on 1
simulate mean=r(mean), reps(10): changePAC_w_MAC_noMAC //Set for 10 for diagnostics and then pushed up for 10,000 for authoritative simulation analysis
timer off 1
timer list 1
*233 s for 10k reps on laptop, 267s for sim3, 2390.62 s for ts2

sum mean
sum mean if mean>0 // number of times greater than zero
display r(N)/10000
// above is percent of times PAC beta controling for MAC is greater than PAC without MAC control
}



*Modeling how PAC beta changes without versus with controling for MAC
*Assuming no PAC effect
{
clear
program drop changePAC_w_MAC_noPAC
program define changePAC_w_MAC_noPAC, rclass
	version 14.2
	drop _all
	matrix input m = (1, 0.064650816, 1, 0.085870761, 0.752885097, 1)
	drawnorm TL F M, n(3241) corr(m) cstorage(lower) clear
	*corr // to double check correlation
	reg TL F
	gen F1=_b[F]
	reg TL F M
	gen F2_minus_F1=_b[F]-F1
	sum F2_minus_F1 
	return scalar mean = r(mean)
	end

timer clear 1
timer on 1
*simulate _b _se, reps(10): myprog2
simulate mean=r(mean), reps(10): changePAC_w_MAC_noPAC
timer off 1
timer list 1
*257 s on sim3

sum mean
sum mean if mean>0 // number of times 
display r(N)/10000
// above is percent of times PAC beta controling for MAC is greater than PAC without MAC control
}


*Modeling how PAC beta changes without versus with controling for MAC
*Assuming negative MAC effect half the magnitude of PAC effect
{
clear
program drop progtemp
program define progtemp, rclass
	version 14.2
	drop _all
	matrix input m = (1, 0.101473869, 1, 0.056866588, 0.752885097, 1)
	drawnorm TL F M, n(3241) corr(m) cstorage(lower) clear
	*corr // to double check correlation
	reg TL F
	gen F1=_b[F]
	reg TL F M
	gen F2_minus_F1=_b[F]-F1
	sum F2_minus_F1 
	return scalar mean = r(mean)
	end

timer clear 1
timer on 1
simulate mean=r(mean), reps(10): progtemp
timer off 1
timer list 1
*233 s for 10k reps on laptop, 267s for sim3, 2390.62 s for ts2

sum mean
sum mean if mean>0 // number of times 
display r(N)/10000*100
// above is percent of times PAC beta controling for MAC is greater than PAC without MAC control
}




*Modeling PAC beta versus MAC beta when together in same model
*Assuming no MAC effect
{
clear
program drop PACvMAC_beta_noMAC
program define PACvMAC_beta_noMAC
	drop _all
	matrix input m = (1, 0.101473869, 1, 0.076398164, 0.752885097, 1)
	drawnorm TL F M, n(3241) corr(m) cstorage(lower) clear
	reg TL F M
	end

timer clear 1
timer on 1
simulate _b _se, reps(10): PACvMAC_beta_noMAC
timer off 1
timer list 1
*120 s on sim3, 316s on ts2, 169 on laptop

sum _b_F _b_M if _b_F>_b_M

foreach v in F M {
 gen t_`v' = _b_`v'/_se_`v'
 gen p_`v' = 2*(1-normal(abs(t_`v')))
}
hist _b_F
hist _b_M
hist p_F
sum _b_F _b_M p_F p_M 
sum 


sum p_F if p_F<0.05
sum p_M if p_M<0.05
hist p_M
}






*Modeling PAC beta versus MAC beta when together in same model
*Assuming no PAC effect
{
clear
program drop progtemp
program define progtemp
	drop _all
	matrix input m = (1, 0.064650816, 1, 0.085870761, 0.752885097, 1)
	drawnorm TL F M, n(3241) corr(m) cstorage(lower) clear
	reg TL F M
	end

timer clear 1
timer on 1
simulate _b _se, reps(10): progtemp
timer off 1
timer list 1
*

sum _b_F _b_M if _b_F>_b_M
sum _b_F if _b_F>0
sum _b_M if _b_M>0
}




*Modeling PAC beta versus MAC beta when together in same model
*Assuming negative MAC effect half the magnitude of PAC effect
{
clear
program drop progtemp
program define progtemp
	drop _all
	matrix input m = (1, 0.101473869, 1, 0.056866588, 0.752885097, 1)
	drawnorm TL F M, n(3241) corr(m) cstorage(lower) clear
	reg TL F M
	end

timer clear 1
timer on 1
simulate _b _se, reps(10): progtemp
timer off 1
timer list 1

sum _b_F _b_M if _b_F>_b_M
sum _b_F if _b_F>0
sum _b_M if _b_M>0
}






*modeling where PAC controlled for when including GPAC to see what happens to coefficients
*true PAC should be 0.101473869 and GPAC effect should be 0.507

*first controlling for PAC and assuming GPAC 0.5 of PAC
{
clear
program drop prog2
program define prog2
	drop _all
	matrix input m = (1, 0.101473869, 1, 0.05091268, 0.041652468, 1)
	drawnorm TL F FF, n(2913) corr(m) cstorage(lower) clear
	reg TL F FF
	end

timer clear 1
timer on 1
simulate _b _se, reps(10): prog2
timer off 1
timer list 1
*

sum _b_F _b_FF
centile _b_FF, centile(1 2.5 5 25 50 75 95 97.5 99)

/* results
                                                      -- Binom. Interp. --
    Variable |       Obs  Percentile    Centile        [95% Conf. Interval]
-------------+-------------------------------------------------------------
       _b_FF |    10,000          1     .004692        .0030358    .0062236
             |                  2.5     .011402        .0104258    .0122466
             |                    5    .0164529        .0159307    .0170899
             |                   25    .0342207        .0337241    .0346961
             |                   50    .0465902          .04617    .0470011
             |                   75    .0592479        .0587162    .0597335
             |                   95    .0759821        .0753539    .0768838
             |                 97.5    .0813579         .080743    .0821146
             |                   99    .0882748        .0866443    .0891206


*/

}



*now controlling for PAC 
*and assuming no GPAC--correlation only evident due to correlation with PAC
*True GPAC of 0.004226637
{
clear
program drop prog2
program define prog2
	drop _all
	matrix input m = (1, 0.101473869, 1, 0.004226637, 0.041652468, 1)
	drawnorm TL F FF, n(2913) corr(m) cstorage(lower) clear
	reg TL F FF
	end

timer clear 1
timer on 1
simulate _b _se, reps(10): prog2
timer off 1
timer list 1
*

sum _b_F _b_FF
centile _b_FF, centile(1 2.5 5 25 50 75 95 97.5 99)

/* results
                                                       -- Binom. Interp. --
    Variable |       Obs  Percentile    Centile        [95% Conf. Interval]
-------------+-------------------------------------------------------------
       _b_FF |    10,000          1   -.0427735       -.0435801   -.0418779
             |                  2.5   -.0360036       -.0371715    -.034839
             |                    5   -.0306063       -.0313663   -.0297657
             |                   25    -.012794       -.0132639   -.0123001
             |                   50   -.0000496       -.0005197     .000457
             |                   75    .0121378        .0116895    .0126861
             |                   95     .030653        .0299391     .031367
             |                 97.5     .036122        .0352277     .037488
             |                   99    .0441277        .0421036    .0454494


*/

}


*now not controlling for PAC and assuming GPAC 0.5 of PAC
{
clear
program drop prog2
program define prog2
	drop _all
	matrix input m = (1, 0.101473869, 1, 0.05091268, 0.041652468, 1)
	drawnorm TL F FF, n(2913) corr(m) cstorage(lower) clear
	reg TL FF
	end

timer clear 1
timer on 1
simulate _b _se, reps(10): prog2
timer off 1
timer list 1
*120 s on sim3, 316s on ts2, 169 on laptop

sum _b_F _b_FF
centile _b_FF, centile(1 2.5 5 25 50 75 95 97.5 99)

/* results


    Variable |       Obs  Percentile    Centile        [95% Conf. Interval]
-------------+-------------------------------------------------------------
       _b_FF |    10,000          1    .0084801         .007033    .0097083
             |                  2.5    .0151792        .0142968    .0164998
             |                    5    .0207839        .0201747    .0215638
             |                   25    .0385307         .038023    .0389614
             |                   50    .0507274        .0502685    .0512164
             |                   75    .0633202        .0629105    .0638656
             |                   95     .080692        .0799301     .081528
             |                 97.5    .0864983        .0855864    .0872704
             |                   99    .0937368        .0923194    .0950752


*/
}




*now controlling for PAC 
*and assuming no GPAC--correlation only evident due to correlation with PAC
*True GPAC of 0.004226637
{
clear
program drop prog2
program define prog2
	drop _all
	matrix input m = (1, 0.101473869, 1, 0.004226637, 0.041652468, 1)
	drawnorm TL F FF, n(2913) corr(m) cstorage(lower) clear
	reg TL FF
	end

timer clear 1
timer on 1
simulate _b _se, reps(10): prog2
timer off 1
timer list 1
*

sum _b_F _b_FF
centile _b_FF, centile(1 2.5 5 25 50 75 95 97.5 99)

/* results
                                                       -- Binom. Interp. --
    Variable |       Obs  Percentile    Centile        [95% Conf. Interval]
-------------+-------------------------------------------------------------
       _b_FF |    10,000          1   -.0383353       -.0397184   -.0372395
             |                  2.5    -.032073       -.0329736   -.0312471
             |                    5   -.0261623       -.0269513   -.0254788
             |                   25   -.0079633       -.0084519   -.0074662
             |                   50    .0045785        .0040651    .0050416
             |                   75    .0166616        .0161219    .0170751
             |                   95    .0347956         .033981    .0354604
             |                 97.5    .0402862         .039537    .0415178
             |                   99    .0473858        .0462049    .0490042

*/

}








*Modeling PAC beta versus PAYC beta when together in same model
*Assuming no PAYC effect
{
clear
program drop progtemp
program define progtemp
	drop _all
	matrix input m = (1, 0.080404349, 1, 0.052202135, 0.649245169, 1)
	drawnorm TL PAC PAYC, n(2606) corr(m) cstorage(lower) clear
	reg TL PAC PAYC
	end

timer clear 1
timer on 1
simulate _b _se, reps(10): progtemp
timer off 1
timer list 1
*120 s on sim3, 316s on ts2, 169 on laptop

sum _b_PAC if _b_PAC>_b_PAYC
display r(N)/10000*100 // 95.96
sum _b_PAC if _b_PAC>0
sum _b_PAYC if _b_PAYC>0
}

*Modeling PAC beta alone versus PAYC beta when together in same model
*Assuming no PAYC effect  
{
clear
program drop progtemp
program define progtemp, rclass
	version 14.2
	drop _all
	matrix input m = (1, 0.080404349, 1, 0.052202135, 0.649245169, 1)
	drawnorm TL PAC PAYC, n(2606) corr(m) cstorage(lower) clear
	corr // to double check correlation
	reg TL PAC
	gen F1=_b[PAC]
	reg TL PAC PAYC
	gen F2_minus_F1=_b[PAC]-F1
	sum F2_minus_F1 
	return scalar mean = r(mean)
	end

timer clear 1
timer on 1
simulate mean=r(mean), reps(10): progtemp
timer off 1
timer list 1


sum mean
sum mean if mean>0 // number of times 
display r(N)/10000*100 // 50.09
// above is percent of times PAC beta controling for PAYC is greater than PAC without PAYC control
}



*Modeling PAC beta alone versus PAYC beta when together in same model
*Assuming no PAC effect  
{
clear
program drop progtemp
program define progtemp, rclass
	version 14.2
	drop _all
	matrix input m = (1, 0.019025664, 1, 0.029304282, 0.649245169, 1)
	drawnorm TL PAC PAYC, n(2606) corr(m) cstorage(lower) clear
	*corr // to double check correlation
	reg TL PAC
	gen F1=_b[PAC]
	reg TL PAC PAYC
	gen F2_minus_F1=_b[PAC]-F1
	sum F2_minus_F1 
	return scalar mean = r(mean)
	end

timer clear 1
timer on 1
simulate mean=r(mean), reps(10): progtemp
timer off 1
timer list 1


sum mean
sum mean if mean>0 // number of times 
display r(N)/10000*100 // 12.52%
// above is percent of times PAC beta controling for PAYC is greater than PAC without PAYC control
}



*Modeling PAC beta versus PAYC beta when together in same model
*Assuming no PAC effect  
{
clear
program drop progtemp
program define progtemp
	drop _all
	matrix input m = (1, 0.019025664, 1, 0.029304282, 0.649245169, 1)
	drawnorm TL PAC PAYC, n(2606) corr(m) cstorage(lower) clear
	reg TL PAC PAYC
	end

timer clear 1
timer on 1
simulate _b _se, reps(10): progtemp
timer off 1
timer list 1
*120 s on sim3, 316s on ts2, 169 on laptop

sum _b_PAC if _b_PAC>_b_PAYC
display r(N)/10000*100 // 26.48
sum _b_PAC if _b_PAC>0
sum _b_PAYC if _b_PAYC>0
}












*Modeling PAC beta versus PAOC beta when together in same model
*Assuming no PAOC effect
{
clear
program drop progtemp
program define progtemp
	drop _all
	matrix input m = (1, 0.065765208, 1, 0.045030346, 0.684713805, 1)
	drawnorm TL PAC PAOC, n(2378) corr(m) cstorage(lower) clear
	reg TL PAC PAOC
	end

timer clear 1
timer on 1
simulate _b _se, reps(10): progtemp
timer off 1
timer list 1
*120 s on sim3, 316s on ts2, 169 on laptop

sum _b_PAC if _b_PAC>_b_PAOC
display r(N)/10000*100 // 
sum _b_PAC if _b_PAC>0
sum _b_PAOC if _b_PAOC>0
}

*Modeling PAC beta alone versus PAOC beta when together in same model
*Assuming no PAOC effect  
{
clear
program drop progtemp
program define progtemp, rclass
	version 14.2
	drop _all
	matrix input m = (1, 0.065765208, 1, 0.045030346, 0.684713805, 1)
	drawnorm TL PAC PAOC, n(2378) corr(m) cstorage(lower) clear
	*corr // to double check correlation
	reg TL PAC
	gen F1=_b[PAC]
	reg TL PAC PAOC
	gen F2_minus_F1=_b[PAC]-F1
	sum F2_minus_F1 
	return scalar mean = r(mean)
	end

timer clear 1
timer on 1
simulate mean=r(mean), reps(10): progtemp
timer off 1
timer list 1


sum mean
sum mean if mean>0 // number of times 
display r(N)/10000*100 // 
// above is percent of times PAC beta controling for PAOC is greater than PAC without PAOC control
}



*Modeling PAC beta alone versus PAOC beta when together in same model
*Assuming no PAC effect  
{
clear
program drop progtemp
program define progtemp, rclass
	version 14.2
	drop _all
	matrix input m = (1, 0.031830737, 1, 0.046487652, 0.684713805, 1)
	drawnorm TL PAC PAOC, n(2378) corr(m) cstorage(lower) clear
	*corr // to double check correlation
	reg TL PAC
	gen F1=_b[PAC]
	reg TL PAC PAOC
	gen F2_minus_F1=_b[PAC]-F1
	sum F2_minus_F1 
	return scalar mean = r(mean)
	end

timer clear 1
timer on 1
simulate mean=r(mean), reps(10): progtemp
timer off 1
timer list 1


sum mean
sum mean if mean>0 // number of times 
display r(N)/10000*100 //
// above is percent of times PAC beta controling for PAOC is greater than PAC without PAOC control
}



*Modeling PAC beta versus PAOC beta when together in same model
*Assuming no PAC effect  
{
clear
program drop progtemp
program define progtemp
	drop _all
	matrix input m = (1, 0.031830737, 1, 0.046487652, 0.684713805, 1)
	drawnorm TL PAC PAOC, n(2378) corr(m) cstorage(lower) clear
	reg TL PAC PAOC
	end

timer clear 1
timer on 1
simulate _b _se, reps(10): progtemp
timer off 1
timer list 1
*120 s on sim3, 316s on ts2, 169 on laptop

sum _b_PAC if _b_PAC>_b_PAOC
display r(N)/10000*100 // 26.48
sum _b_PAC if _b_PAC>0
sum _b_PAOC if _b_PAOC>0
}









*modeling CI for 2x reduction of PAC effect 
{
clear
program drop progtemp
program define progtemp, rclass
	drop _all
	matrix input m = (1, 0.0965, 1, 0.04825, 0.041652468, 1, 0.024125, 0.001734928, ///
	0.041652468, 1, 0.0120625,	0.0000722640, 0.001734928, 0.041652468,	1)
	drawnorm TL	PAC	GPAC GGPAC GGGPAC, n(3282) corr(m) cstorage(lower) clear
	reg TL PAC
	gen PACbeta=_b[PAC]
	reg TL PAC GPAC if _n<2914
	gen PACdividedGPAC=PACbeta/_b[GPAC]
	sum PACdividedGPAC
	return scalar mean = r(mean)
	end

timer clear 1
timer on 1
simulate mean=r(mean), reps(10000): progtemp
timer off 1
timer list 1
*
centile mean, centile(1 2.5 5 25 50 75 95 97.5 99)

/*


                                                       -- Binom. Interp. --
    Variable |       Obs  Percentile    Centile        [95% Conf. Interval]
-------------+-------------------------------------------------------------
        mean |    10,000          1    .8000177        .7435785    .8410536
             |                  2.5    .9603843        .9447164    .9889303
             |                    5    1.114969        1.097572    1.130675
             |                   25    1.639624        1.620073    1.659595
             |                   50     2.16064        2.137591     2.18507
             |                   75    2.996959        2.950066    3.046617
             |                   95    6.628198        6.325772    6.876419
             |                 97.5    9.614764        9.060184    10.34429
             |                   99     16.9589        14.75463     19.4426

median effect Narrows slightly if don't control for PAC in calculating GPAC effect:
                                                       -- Binom. Interp. --
    Variable |       Obs  Percentile    Centile        [95% Conf. Interval]
-------------+-------------------------------------------------------------
        mean |    10,000          1    .8065147        .7545073     .849601
             |                  2.5     .965015        .9379735    .9852216
             |                    5    1.088659        1.069613    1.103586
             |                   25     1.53999        1.524872    1.553152
             |                   50    1.982629        1.964518     2.00373
             |                   75     2.71653        2.686931    2.754178
             |                   95    5.221854         5.04483    5.408832
             |                 97.5    7.158693        6.777712    7.544891
             |                   99    11.43039        10.02576    13.99945


			 
*/



}