#Written by Erik Gregory# #egregory2007@yahoo.com# #California State University Sacramento# #19 July 2010# #Dr. Taylor wants to estimate the proportion of teenage patients in her practice# #currently experiencing flu symptoms. On one day she collects all of the data for# #her study to determine the mean length of time with flu symptoms.# #Dr. Taylor observes that patients seem to be more likely to have flu symptoms for# #a longer period of time than the CDC predicts. She wants to know if this measurement# #is statistically significant or not.# #Suppose the true number of individuals in Dr. Taylor's practice with flu symptoms# #starting within the last week is 50. The simulation proceeds as follows.# #Random number generator seed...(anyone get the reference?)# set.seed(48151623) #variable definitions# categories <- 0*(1:8) l <- 0*(1:50) m <- 0*(1:50) #the 8 possible durations# averages <- 1:8 #weights on our categories# weights <- 8/averages #sample size# N = 50 #start dates between 0 and 7 days# start.date <- round(runif(N, 0, 7)-0.5) #durations of symptoms between 1 and 7 days# duration <- round(runif(N, 1, 8)) patients <- data.frame(start.date = start.date, duration = duration, end.date = start.date + duration) #sort out our sampled individuals# patients <- patients[order(patients\$end.date),] row.names(patients) <- 1:N j <- 1 while(patients\$end.date[j] < 7) { j <- j+1 } w <- j sampled <- patients[w:N,] row.names(sampled) <- 1:(N-w+1) #sorts the indicies of our sampled/unsampled individuals (to be used later)# patients <- patients[order(patients\$start.date),] row.names(patients) <- 1:N j <- 1 while(j<=N) { if(patients\$end.date[j] >= 7) { l[j] <- j } else { m[j] <- j } j <- j+1 } l <- l[l>0] m <- m[m>0] #Sorts our sampled individuals into categories by duration# sampled <- sampled[order(sampled\$duration),] row.names(sampled) <- 1:(N-w+1) j <- 1 i <- 1 while(j<=N-w+1) { k <- 0 while(sampled\$duration[j] == i && j<=N-w+1) { k <- k+1 categories[i] <- k j <- j+1 } i <- i+1 } #The uncorrected mean duration# uncorrected <- mean(sampled\$duration) #Actual mean duration# true <- mean(patients\$duration) #Corrected mean duration# corrected <- sum(categories)/sum(categories/averages) #Plot of the longitudinal lines# plot(patients\$start.date, 1:N, main = "Patients", xlim = c(0, 16), xlab = "Time", cex = .01, ylab = "Patient") points(patients\$end.date, 1:N, cex = 0.01) segments(patients\$start.date[l], l , patients\$end.date[l], l, col = "blue") segments(patients\$start.date[m], m, patients\$end.date[m], m) abline(v = 7, col = "red")