arcsin transformation of P is normal as follow:Y = arcsin(P^(1/2)) ~ N(arcsin(x^(1/2)),1/4n), x = true population proportion, n = number of individuals in a sample.

Stata code:

/**

Program to produce simulated "prevalence samples"

Input parameters (arguments):

nsam number of samples or studies to simulate

n sample size of each "study"

pi population prevalence

cl confidence level (2-sided) (default value 0.95)

Result:

leaves dataset of simulated prevalence values with 100*`cl'% CIs in memory

(as current data)

**/

program prevsim, rclass

version 12

syntax [, nsam(integer 1) n(integer 100) pi(real 0.5) cl(real 0.95)]

drop _all

set obs `nsam'

tempvar obst

local zn=-invnorm((1-`cl')/2)

gen `obst' = rbinomial(`n',`pi')

gen prev = `obst'/`n'

gen ci_upp = prev+`zn'*sqrt(prev*(1-prev)/`n')

gen ci_low = prev-`zn'*sqrt(prev*(1-prev)/`n')

end

a. Modify the code to calculate confidence interval based on the arcsin transformation (calculate confidence interval for w = arcsin(x^(1/2)), and use this to calculate confidence interval for x).

add code before end:

gen arcsinp = asin(sqrt(prev))

gen ci_u_as = arcsinp+`zn'*sqrt(0.25/`n') gen ci_l_as = arcsinp-`zn'*sqrt(0.25/`n') gen ci_pa_upp = sin(ci_u_as)^2

gen ci_pa_low = sin(ci_l_as)^2

Assume x = 0.03, n = 90. We have 10 hypothetical repeated studies. Calculate sample proportion, and 90% confidence interval by arcsin transformation and standard way. How many confidence intervals in each case include x? How many should expect?

1. b. repeat part (a) by 10,000 simulations. Calculate and compare the coverage probabilities of the two confidence intervals, by modifying the Stata code:

/*

do-file showing use of 'prevsim' command defined in prevsim.ado

*/

* simulate 1000 samples of size 10 with true prevalence 0.2

version 12

prevsim, nsam(1000) n(10) pi(0.2) cl(0.95)

* look at resulting sampling distribution of (estimated) prevalences

hist prev, discrete

* calculate coverage

count if ci_low<=0.2 & ci_upp>=0.2

* return list/* command shows the saved results produced by 'count' command */

display r(N)/1000

* display the CIs to visualise coverage

gen x=_n

twoway pccapsym ci_low x ci_upp x if x<=100, horiz xline(0.2) xline(0,lpatt(dash)) msym(oh) msiz(small)

* repeat the above with larger sample size (n) to see better properties (normal sampling distribution)

Explain how the results are relevant to2 cases in a sample size of 90?

c. Repeat part(b) for 900 sample size and compare with part(b).