CSU Hayward

Statistics Department

Appendix D: Main Parts of S-plus
Simulation Programs

Introduction

The main parts of the S-plus programs we used to perform the simulations described in Sections 4 and 8 of this paper are shown below. Those who are familiar with reading computer programs may find it easier to understand the simulation procedure by reading the code than by reading the descriptions in the main article.

Here we have omitted code for initializations, declarations, and printing of graphics. The full versions of these programs are also available on this site. See the links at the end of each section below.

Section 4

The fragment of S-plus code below shows the main loop for generating values of the D-process. Variable names are the same as those in the main part of the article, except that D(m) is called Disease[m] here, and T(m) is called Test[m]. The function rbinom(1,1,eta) randomly generates 1 observation from a binomial distribution with 1 trial and success probability h. The vector pi.prop[m] contains the running proportions used to make Figure 1. The quantity pi.post is the prevalence estimate of Table 2. 

#Generate Values
Sum <- 0
for (m in 2:M2)
{      
      if (Disease[m-1] == 1)
            Test[m-1] <- rbinom(1,1,eta)
      else
            Test[m-1] <- rbinom(1,1,1-theta)
                                                         
      if (Test[m-1] == 1) 
	    Disease[m] <- rbinom(1,1,gamma)
      else 
	    Disease[m] <- rbinom(1,1,1-delta)
                  
      Sum <- Sum + Disease[m]
      Pi.prop[m] <- Sum / m
}

The full program for the simulation of section 4 is available on this site.

Section 8

The fragment of code below shows the main loop for simulating values of prevalence P. Again here, the variable names in the code follow closely the notation used in the article. The last half of the values in the vector pi[m] are used to compute the estimates shown in Table 4 and to make the histogram in Figure 3. All of its values are used in Figures 4 and 5. The function rbeta samples from the beta distribution to obtain values of pi[m] . 

#Generate the chain of prevalence values 
Sum <- 0
for (m in 2:M2)                         
{
        X[m]    <-      rbinom(1, A, (pi[m-1]*eta)/
                                ((pi[m-1]*eta)+(1-pi[m-1])*(1-theta)))

        Y[m]    <-      rbinom(1, B, (pi[m-1]*(1-eta))/
                                ((pi[m-1]*(1-eta))+(1-pi[m-1])*theta))

        pi[m]   <-      rbeta(1, (X[m] + Y[m] + alpha), 
                                (A + B - X[m] - Y[m] + beta))

        Sum     <-      Sum + pi[m]

        Average[m] <-   Sum / m
}

The full program for the simulation of section 8 is available on this site.