{smcl}
{com}{sf}{ul off}{txt}{.-}
      name:  {res}<unnamed>
       {txt}log:  {res}/Users/eric/Work/MB/ATCR/200/Labs/lab 3.smcl
  {txt}log type:  {res}smcl
 {txt}opened on:  {res} 9 Jul 2018, 04:54:55
{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
. 
. 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 forumla
. 
. 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 distribution = " round(sqrt(12*.8*.2),.01)
{res}SD of distribution = 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 binomialtail(12,6,.2)
{res}.01940528
{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 smalle?
. 
. * 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 help for 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 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 result is close to the theoretical value of 9.1%.  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, the number of binge drinking episodes reported between study visits, 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 = $En"
{res}Expected number of cases = 
{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, k, and p, gives the probabililty of observing k failures before n successes are observed, where p is the probability of success.  
. 
. * That of course is a counter-intuitive way to think of the distribution. Consider that the mean of the distribution of k 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}. 
. * 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),.01)
{res}Probability of 10 failures before 5th success = .06
{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),.01)
{res}Probability of 5-10 failures before 7th success = .43
{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 stat 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 time units?
. 
. forvalues t = 1/10 {c -(}
{txt}  2{com}.         local p_`t' = round(exponential(10,`t'),.01)
{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: .1
Probability of failure within 2 years: .18
Probability of failure within 3 years: .26
Probability of failure within 4 years: .33
Probability of failure within 5 years: .39
Probability of failure within 6 years: .45
Probability of failure within 7 years: .5
Probability of failure within 8 years: .55
Probability of failure within 9 years: .59
Probability of failure within 10 years: .63
{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),.01)
{res}Probability of surviving >8 years: .45
{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 3.smcl
  {txt}log type:  {res}smcl
 {txt}closed on:  {res} 9 Jul 2018, 04:55:05
{txt}{.-}
{smcl}
{txt}{sf}{ul off}