{smcl}
{txt}{sf}{ul off}{.-}
      name:  {res}<unnamed>
       {txt}log:  {res}/Users/eric/Work/MB/ATCR/200/Labs/lab 2.smcl
  {txt}log type:  {res}smcl
 {txt}opened on:  {res}18 Sep 2018, 06:55:03
{txt}
{com}. 
. ********************************************************************************
. /* Binomial distribution.
> 
> The probability of being selected by random digit dialing for a political survey 
> in your area code is 0.025. The surveys are conducted randomly each month for 
> one year.  What is the probability that you will not be selected in the entire 
> year (i.e. the probability you will not be selected in any of 12 months)?  Use 
> the binomial formula.
> 
> We will do this two ways using built-in three built-in Stata functions.  
> Briefly check the help files for these functions.  The help files should appear 
> in a viewer window. */
. 
. help comb
{txt}
{com}. help binomialtail
{txt}
{com}. help binomialp
{txt}
{com}. 
. /* Now we implement the binomial formula directly using comb(N,k), which 
> evaluates the combinatoric (N choose k), and the exponentiation operator ^. */
. 
. display as result 1-comb(12,0)*0.025^0*(1-.025)^12
{res}.26200165
{txt}
{com}. 
. /* The "as result" option in the display command bolds the output. An easier 
> way to do this is to use the binomialtail function. */
. 
. di as result binomialtail(12,1,.025)
{res}.26200165
{txt}
{com}. 
. /* Now we calculate the probability of being chosen once or twice, first using 
> the binomial formula. */
. 
. di as result comb(12,1)*0.025^1*(1-.025)^11+comb(12,2)*0.025^2*(1-.025)^10
{res}.25910001
{txt}
{com}. 
. * and now using the binomialp function
. di as result binomialp(12,1,.025)+binomialp(12,2,.025)
{res}.25910001
{txt}
{com}. 
. /* Now we use binomialtail to get the probability of being selected 2 times 
> or fewer: */
. 
. di as result 1-binomialtail(12,3,.025)
{res}.99709836
{txt}
{com}. 
. /* Now we will get summary statistics and plot the binomial distribution with 
> first 10, then 100, trials and varying probabilities of success.  Using loops, 
> we'll first make a dataset with 11 or 101 observations, one for each possible 
> number of successes, then calculate the corresponding probability. */
. 
. foreach n in 10 100 {c -(}
{txt}  2{com}.         clear
{txt}  3{com}.         set obs `=`n'+1' // note syntax to plug in the sample size + 1
{txt}  4{com}.         gen x = _n-1 // _n is the observation number
{txt}  5{com}.         foreach p in .5 .18 .82 {c -(}
{txt}  6{com}.         
.                 * calculate mean and SD
.                 local mean = round(`n'*`p',.01)
{txt}  7{com}.                 local sd = round(sqrt(`n'*`p'*(1-`p')),.01)
{txt}  8{com}.                 di as result "N = `n', P = `p', mean = `mean', SD = `sd'"
{txt}  9{com}.                 
.                 * calculate probability of each number of successes
.                 local pct = `p'*100 // variable and plot names can't include "."
{txt} 10{com}.                 gen p`pct' = binomialp(`n',x,`p')
{txt} 11{com}.                 * plot the distribution
.                 twoway (bar p`pct' x, fcolor(blue)), ///
>                         title("N = `n', P = `pct'%") caption("mean = `mean', SD = `sd'") ///
>                         name(b`n'_`pct', replace)
{txt} 12{com}.                 {c )-}
{txt} 13{com}.         {c )-}
{txt}{p}
number of observations (_N)  was 0,
now 11
{p_end}
{res}N = 10, P = .5, mean = 5, SD = 1.58
N = 10, P = .18, mean = 1.8, SD = 1.21
N = 10, P = .82, mean = 8.199999999999999, SD = 1.21
{txt}{p}
number of observations (_N)  was 0,
now 101
{p_end}
{res}N = 100, P = .5, mean = 50, SD = 5
N = 100, P = .18, mean = 18, SD = 3.84
N = 100, P = .82, mean = 82, SD = 3.84
{txt}
{com}. 
. /* Note that Stata has trouble rounding one of the means; this comes from the 
> conversion back from binary to decimal notation.  Now combine the plots to make 
> them easier to compare. */
. 
. graph combine b10_50 b10_18 b10_82 b100_50 b100_18 b100_82, rows(2)
{res}{txt}
{com}.         
. * Are the plots consistent with the calculated means and SDs?  
. 
. /* According to the CDC, 20% of all adult men in the U.S. smoke cigarettes.  If 
> you selected repeated samples of size 12 from all adult males in the U.S., what
> would be the mean number of men who do not smoke? What would be the mean number 
> of men who smoke? What would be the standard deviation? */
. 
. di as result "Mean # non-smokers = " round(12*.8,.01)
{res}Mean # non-smokers = 9.6
{txt}
{com}. di as result "Mean # smokers = " round(12*.2,.01)
{res}Mean # smokers = 2.4
{txt}
{com}. di as result "SD of # of smokers = " round(sqrt(12*.8*.2),.01)
{res}SD of # of smokers = 1.39
{txt}
{com}. 
. /* Suppose you select a sample of 12 men and find that 6 of them smoke.  
> Assuming .20 is the true probability of smoking among men, what is the 
> probability that you would have a sample this extreme? */
. 
. di as result ///
>         "Probability of >=6 smokers in sample of 12, assuming prevalence of 20%: " ///
>         round(binomialtail(12,6,.2),.001)
{res}Probability of >=6 smokers in sample of 12, assuming prevalence of 20%: .019
{txt}
{com}. 
. ********************************************************************************
. /* Hypergeometric distribution
> 
> A binomial outcome k can be thought of as the number of successes in a sample 
> of size n, sampled with replacement from a population of size N, in which the 
> proportion of successes K/N=p. The corresponding hypergeometric outcome k is 
> the number of successes in a sample of size n, in this case sampled without 
> replacement from the same population.  Commonly we assume that the population 
> is so large that sampling with or without replacement is of no practical 
> importance.  But what if N is smaller?
> 
> Dr. Allen gives this New York City-based example:  Having bought a bag of 30 
> roasted chestnuts, you walk home eating them with gusto, and arrive home having
> eaten 20.  In opening the remaining 10, you find that 7 contain worms.  What 
> is the probability that none if the first 20 contains worms? 
> 
> First check the hypergeometric probability functions. */
. 
. help hypergeometricp
{txt}
{com}. help hypergeometric
{txt}
{com}. 
. /* We can rephrase the problem as follows:  In a sample of 20 sampled without 
> replacement from a population of size 30 including 7 "successes", what is the 
> probability of drawing 0 successes? */
. 
. di as result "probability you recently ate no wormy chestnuts = " ///
>         round(hypergeometric(30,7,20,0),.000001)
{res}probability you recently ate no wormy chestnuts = .000059
{txt}
{com}. 
. * We can also calculate this from the formula  
. di as result "probability you recently ate no wormy chestnuts = " ///
>         round(comb(7,0)*comb(23,20)/comb(30,20),.000001)
{res}probability you recently ate no wormy chestnuts = .000059
{txt}
{com}. 
. * Note that we could have used either hypergeometric function in this case.
. 
. ********************************************************************************
. /* Normal distribution
> 
> First, familiarize yourself with the Stata normal cumulative probability  and 
> inverse cumulative probability functions. */
. 
. help normal
{txt}
{com}. help invnormal
{txt}
{com}.         
. /* The lengths of adult men's feet are approximately normal, with mean 11 
> inches and standard deviation 1.5 inches. 
> 
> What is the probability of a foot length of more than 13 inches? Because the 
> argument to the normal function is Z-score (i.e., the observed value minus the 
> mean, divided by the SD), we calculate that first.  Using a local variable 
> avoids typing in a calculated number, which is prone to error. Before running 
> this next bit of code, try to put in words what it is doing. */
. 
. local z = (13-11)/1.5 
{txt}
{com}. di as result "probability of foot length >=13 inches: " ///
>         round(1-normal(`z'),.001)
{res}probability of foot length >=13 inches: .091
{txt}
{com}. 
. /* We could also do this in a single step, because the normal function can 
> calculate its argument on the fly.  (Stata is inconsistent about allowing 
> this.) */
. 
. di as result "probability of foot length >=13 inches: " ///
>         round(1-normal((13-11)/1.5),.001)
{res}probability of foot length >=13 inches: .091
{txt}
{com}. 
. * What is the probability of a foot length between 10 and 12 inches?
. 
. local z10 = (10-11)/1.5
{txt}
{com}. local z12 = (12-11)/1.5
{txt}
{com}. local p1012 = round(normal(`z12')-normal(`z10'),.01)
{txt}
{com}. di as result "probability of foot length 10-12 inches: `p1012'"
{res}probability of foot length 10-12 inches: .5
{txt}
{com}. 
. /* This could also have been done in a single step; try that if you like.
> 
> What are the foot lengths at the 25th percentile and 75th percentile? The 
> argument of the invnormal function is a probability and the output is a 
> Z-score, from which we calculate the observed value using the mean and SD. */
. 
. di as result "25th pctile = " round(11+invnormal(.25)*1.5,.01) ///
>         "  inches"  _n "75th pctile = " round(11+invnormal(.75)*1.5,.01) " inches"
{res}25th pctile = 9.99  inches
75th pctile = 12.01 inches
{txt}
{com}. 
. /* The probability of 0.04 that a randomly chosen adult male foot length will 
> be less than how many inches? */
. 
. di as result "4th pctile = " round(11+invnormal(.04)*1.5,.01) "inches"
{res}4th pctile = 8.37inches
{txt}
{com}. 
. /* Stata had trouble with the rounding again.  
> 
> We could also have done this using simulation, which is commonly used when a 
> distribution isn't well understood (in contrast to the Normal distribution).  
> To get reasonably reliable estimates, we'll use a large simulated sample of 
> 100,000, and also set the seed for the random number generator, to make the 
> results reproducible. */
. 
. clear
{txt}
{com}. set obs 100000
{txt}{p}
number of observations (_N)  was 0,
now 100,000
{p_end}

{com}. set seed 200
{txt}
{com}. gen fl = rnormal(11,1.5)
{txt}
{com}. label var fl "foot length"
{txt}
{com}. 
. /* Following lab 1, we'll use the ordering of the elements of the recode 
> command to ensure that any values of 13 are counted in the >= 13 category. */
. 
. recode fl (13/max=1 ">=13 inches") (min/13=0 "<13 inches"), gen(fl_13up)
{txt}(100000 differences between fl and fl_13up)

{com}. label var fl_13up "foot length"
{txt}
{com}. tab fl_13up

{txt}foot length {c |}      Freq.     Percent        Cum.
{hline 12}{c +}{hline 35}
 <13 inches {c |}{res}     90,812       90.81       90.81
{txt}>=13 inches {c |}{res}      9,188        9.19      100.00
{txt}{hline 12}{c +}{hline 35}
      Total {c |}{res}    100,000      100.00
{txt}
{com}. 
. recode fl (10/12=1 "10-12 inches") (min/10 12/max=0 "<10, >12 inches"), gen(fl_1012)
{txt}(100000 differences between fl and fl_1012)

{com}. label var fl_1012 "foot length"
{txt}
{com}. tab fl_1012

    {txt}foot length {c |}      Freq.     Percent        Cum.
{hline 16}{c +}{hline 35}
<10, >12 inches {c |}{res}     50,779       50.78       50.78
{txt}   10-12 inches {c |}{res}     49,221       49.22      100.00
{txt}{hline 16}{c +}{hline 35}
          Total {c |}{res}    100,000      100.00
{txt}
{com}. 
. /* The results are close to the theoretical values of 9.1% and 50%.  We can 
> also get the quantiles of the distribution as local variables, then use the 
> command return list to see the results. */
. 
. _pctile fl, percentiles(4 25 75)
{txt}
{com}. return list

{txt}scalars:
                 r(r1) =  {res}8.369423389434814
                 {txt}r(r2) =  {res}9.983631134033203
                 {txt}r(r3) =  {res}12.02052927017212
{txt}
{com}. 
. /* Results are again close to the theoretical values. To get a sense of 
> simulation error, try redoing this with a smaller simulated sample. */
. 
. ********************************************************************************
. /* Poisson and negative binomial distributions
> 
> Count outcomes are commonly encountered in epidemiology.  Examples include the 
> number of cases of incident disease in a given number of person-years of 
> follow-up, number of binge drinking episodes, and so on.  
> 
> Check out the four Poisson probability functions: */
. 
. help poissonp 
{txt}
{com}. help poisson 
{txt}
{com}. help poissontail
{txt}
{com}. help invpoisson
{txt}
{com}. help invpoissontail
{txt}
{com}. 
. /* Suppose you accrue 10,545 person-years of follow-up in your cohort, and the 
> population incidence of myocardial infarction is 2.5% (i.e., 2.5 cases per 100 
> person-years).  Assuming that the distribution is Poisson, which is reasonable 
> in this case, what is the expected number of cases?  We'll save the expected 
> number of cases as a global variable, which will be remembered even if you 
> run the commands below one at at time; note that we call global variables 
> using a different syntax. */
. 
. global mean = 10545*.025
{txt}
{com}. di as result "Expected number of cases = $mean"
{res}Expected number of cases = 263.625
{txt}
{com}. 
. * What is the probability of observing at least 300 MI cases?
. 
. di as result "probability of observing >=300 MI cases: " ///
>         round(poissontail($mean,300),.001)
{res}probability of observing >=300 MI cases: .015
{txt}
{com}. 
. * What is the chance of observing between 250 and 275 cases?
. 
. di as result "probability of observing 250-275 MI cases: " ///
>         round(poisson($mean,275)-poisson($mean,249),.001)
{res}probability of observing 250-275 MI cases: .576
{txt}
{com}. 
. * What is the IQR of the distribution?
. 
. local q25 = round(invpoisson($mean, .25),.1)
{txt}
{com}. local q75 = round(invpoisson($mean, .75),.1)
{txt}
{com}. di as result ///
>         "IQR of Poisson distribution with mean `=round($mean,.1)': `q25'-`q75'"
{res}IQR of Poisson distribution with mean 263.6: 274.8-252.9
{txt}
{com}.         
. /* Although the Poisson model might be adequate for the number of MIs, the 
> number of binge drinking episodes might be "over-dispersed", with variance 
> greater than its mean.  
> 
> The negative binomial distribution with parameters n and p gives the 
> probabililty of observing k failures before n successes are observed, where 
> p is the probability of success. The extra parameter n effectively allows the 
> variance to be larger than the mean.
> 
> That being said, the parameters n and p are a counter-intuitive way to think 
> of the distribution. Consider that the mean of the distribution is n(1-p)/p 
> and its variance is n(1-p)/p^2.  A little algebra shows that p is the ratio 
> of the mean to the variance, and n is given by mean^2/(variance-mean). 
> Consider a negative binomial distribution with mean 2 and variance 5 
> (so it is over-dispersed, relative to the Poisson model).  We can check 
> to see that these formulas are correct. */
. 
. local m = 2
{txt}
{com}. local v = 5
{txt}
{com}. local p = `m'/`v'
{txt}
{com}. local n = `m'^2/(`v'-`m')
{txt}
{com}. clear
{txt}
{com}. set obs 100000
{txt}{p}
number of observations (_N)  was 0,
now 100,000
{p_end}

{com}. gen y = rnbinomial(`n',`p')
{txt}
{com}. tabstat y, s(mean var)

{txt}{ralign 12:variable} {...}
{c |}      mean  variance
{hline 13}{c +}{hline 20}
{ralign 12:y} {...}
{c |}{...}
 {res}  1.99457   4.95259
{txt}{hline 13}{c BT}{hline 20}

{com}. clear
{txt}
{com}. 
. * Now familiarize with the negative binomial probability functions
. 
. help nbinomialp
{txt}
{com}. help nbinomial
{txt}
{com}. help nbinomialtail
{txt}
{com}. help invbinomial
{txt}
{com}. help invbinomialtail
{txt}
{com}. 
. /* Now use these functions to calculate the probability of observing exactly 
> 10 failures before the 5th success, if the probability of success is .4. 
> Also calculate the mean and variance of the distribution. */
. 
. di as result "Probability of 10 failures before 5th success = " ///
>         round(nbinomialp(5,10,.4),.001)
{res}Probability of 10 failures before 5th success = .062
{txt}
{com}. di as result "Mean = " 5*.6/.4 "  Variance = " 5*.6/.4^2
{res}Mean = 7.5  Variance = 18.75
{txt}
{com}. 
. /* Now calculate the probability of observing 5-10 failures before the 7th 
> success, with success probabiity .6. */
. 
. di as result  "Probability of 5-10 failures before 7th success = " ///
>         round(nbinomial(7,10,.6)-nbinomial(7,4,.6),.001)
{res}Probability of 5-10 failures before 7th success = .432
{txt}
{com}. 
. ********************************************************************************
. /* Exponential failure time distribution
> 
> The exponential distribution describes time to an endpoint if the failure rate 
> is a constant, and is closely related to the Poisson distribution, as explained 
> in class. Stata parameterizes the distribution in terms of the its mean or 
> scale, which is the inverse of the failure rate per unit time.  Have a look 
> at the functions for this distribution: */
. 
. help exponential()
{txt}
{com}. help invexponential()
{txt}
{com}. 
. /* Suppose the failure rate of 10% per year, so the mean failure time is 10 
> years.  What is the probability of failing within 2, 5, and 10 years? */
. 
. forvalues t = 1/10 {c -(}
{txt}  2{com}.         local p_`t' = round(exponential(10,`t'),.001)
{txt}  3{com}.         di as result "Probability of failure within `t' years: `p_`t''"
{txt}  4{com}.         {c )-}
{res}Probability of failure within 1 years: .095
Probability of failure within 2 years: .181
Probability of failure within 3 years: .259
Probability of failure within 4 years: .33
Probability of failure within 5 years: .393
Probability of failure within 6 years: .451
Probability of failure within 7 years: .503
Probability of failure within 8 years: .551
Probability of failure within 9 years: .593
Probability of failure within 10 years: .632
{txt}
{com}.         
. /* Note that the cumulative risk is not a linear function of time.  Why?
> 
> What is the probability of surviving more than 8 years?  Note that the exponentialtail function is not what we want for this purpose. */
. 
. di as result "Probability of surviving >8 years: " ///
>         1-round(exponential(10,8),.001)
{res}Probability of surviving >8 years: .449
{txt}
{com}. 
. * Finally, get the IQR of the distribution
. local q25 = round(invexponential(10,.25),.01)
{txt}
{com}. local q75 = round(invexponential(10,.75),.01)
{txt}
{com}. di as result ///
>         "IQR of exponential distribution with mean 10: `q25'-`q75' years"
{res}IQR of exponential distribution with mean 10: 2.88-13.86 years
{txt}
{com}. 
. log close
      {txt}name:  {res}<unnamed>
       {txt}log:  {res}/Users/eric/Work/MB/ATCR/200/Labs/lab 2.smcl
  {txt}log type:  {res}smcl
 {txt}closed on:  {res}18 Sep 2018, 06:55:11
{txt}{.-}
{smcl}
{txt}{sf}{ul off}