SDaA_using_survey.R
In SDaA: Sampling: Design and Analysis

### R code from vignette source 'SDaA_using_survey.Rnw'

###################################################
### code chunk number 1: config
###################################################
	options(width=65)


###################################################
### code chunk number 2: loadPkg
###################################################
  library(SDaA)
  library(survey)


###################################################
### code chunk number 3: SDaA_using_survey.Rnw:45-49
###################################################
  ### Example 3.2, p. 63 
  agsrsDesign <- svydesign(ids=~1, weights = ~1, data = agsrs)
  svyratio(numerator = ~acres92, denominator = ~acres87, 
      design = agsrsDesign) # proportion B hat


###################################################
### code chunk number 4: acreagetwoyears
###################################################
  ### part of Example 3.2, p. 64
  plot(I(acres92/10^6) ~ I(acres87/10^6), 
      xlab = "Millions of Acres Devoted to Farms (1987)", 
      ylab = "Millions of Acres Devoted to Farms (1992)", data = agsrs)
  abline(lm(I(acres92/10^6) ~ 0 + I(acres87/10^6), # through the origin 
          data = agsrs), col = "red", lwd = 2)


###################################################
### code chunk number 5: seedlingData
###################################################
  ### Example 3.5, p. 72, table 3.3
  seedlings <- data.frame(tree = 1:10, 
      x = c(1, 0, 8, 2, 76, 60, 25, 2, 1, 31),
      y = c(0, 0, 1, 2, 10, 15, 3, 2, 1, 27))
  names(seedlings) <- c("tree", "x", "y")


###################################################
### code chunk number 6: seedlingsplot
###################################################
  plot(y ~ x, data = seedlings, xlab = "Seedlings Alive (March 1992)",
      ylab = "Seedlings That Survived (February 1994)")
  # abline(lm(y ~ 0 + x, data = seedlings), lwd = 2, col = "red")
  # TODO: add proper abline


###################################################
### code chunk number 7: SDaA_using_survey.Rnw:93-99
###################################################
  ### Example 3.6, p. 75
  pf <- data.frame(photo = c(10, 12, 7, 13, 13, 6, 17, 
          16, 15, 10, 14, 12, 10, 5,
          12, 10, 10, 9, 6, 11, 7, 9, 11, 10, 10),
      field = c(15, 14, 9, 14, 8, 5, 18, 15, 13, 15, 11, 15, 12,
          8, 13, 9, 11, 12, 9, 12, 13, 11, 10, 9, 8))


###################################################
### code chunk number 8: modelreg
###################################################
	### Example 3.9, p. 83
  recacr87 <- agsrs$acres87
  recacr87[recacr87 > 0] <- 1/recacr87[recacr87 > 0] # cf. p. 450
  model1 <- lm(acres92 ~ 0 + acres87, weights = recacr87, data = agsrs)
  summary(model1)


###################################################
### code chunk number 9: modelregplot
###################################################
  ### Figure 3.6, p. 85
  wtresid <- resid(model1) / sqrt(agsrs$acres87) 
  plot(wtresid ~ I(agsrs$acres87/10^6), 
      xlab = "Millions of Acres Devoted to Farms (1987)",
      ylab = "Weighted Residuals")


###################################################
### code chunk number 10: acresboxplot
###################################################
	boxplot(acres92/10^6 ~ region, xlab = "Region", 
      ylab = "Millions of Acres", data = agstrat)


###################################################
### code chunk number 11: gpaex
###################################################
  ### Example 5.2, p. 137 middle
  GPA <- cbind(expand.grid(1:4, 1:5), 
      gpa = c(3.08, 2.60, 3.44, 3.04, 2.36, 3.04, 3.28, 2.68, 2.00, 2.56, 
          2.52, 1.88, 3.00, 2.88, 3.44, 3.64, 2.68, 1.92, 3.28, 3.20))
  names(GPA)[1:2] <- c("person_num", "cluster")
  GPA$pwt <- 100/5
  
  clusterDesign <- svydesign(ids = ~ cluster, weights = ~ pwt, data = GPA)
  
  svytotal(~ gpa, design = clusterDesign)
  #      total     SE 
  # gpa 1130.4 67.167
  
  # Stata results: 1130.4   67.16666 ---> corresponds perfectly


###################################################
### code chunk number 12: SDaA_using_survey.Rnw:172-175
###################################################
  ### Figure 5.3
	plot(volume ~ clutch, xlim = c(0,200), pch=19, data = coots,
      xlab = "Clutch Number", ylab = "Egg Volume")


###################################################
### code chunk number 13: SDaA_using_survey.Rnw:178-181
###################################################
  ### Figure 5.3
  plot(volume ~ clutch, xlim = c(0,200), pch=19, data = coots,
      xlab = "Clutch Number", ylab = "Egg Volume")


###################################################
### code chunk number 14: readStatepop
###################################################
  data(statepop)
  statepop$psi <- statepop$popn / 255077536


###################################################
### code chunk number 15: SDaA_using_survey.Rnw:197-201
###################################################
	### page 191, figure 6.1
  plot(phys ~ psi, data = statepop, 
       xlab = expression(paste(Psi[i], " for County")),
       ylab = "Physicians in County (in thousands)")


###################################################
### code chunk number 16: SDaA_using_survey.Rnw:210-215
###################################################
	### Figure 7.1
  data(htpop)
  popecdf <- ecdf(htpop$height)
  plot(popecdf, do.points = FALSE, ylab = "F(y)", 
       xlab = "Height Value, y")


###################################################
### code chunk number 17: SDaA_using_survey.Rnw:218-223
###################################################
  ### Figure 7.2
  minht <- min(htpop$height)
  breaks <- c(minht-1, seq(from = minht, to = max(htpop$height), by = 1))
  hist(htpop$height, ylab = "f(y)", breaks = breaks, 
       xlab = "Height Value, y", freq = FALSE)


###################################################
### code chunk number 18: SDaA_using_survey.Rnw:226-230
###################################################
  ### Figure 7.3
  data(htsrs)
  hist(htsrs$height, ylab = "Relative Frequency", 
       xlab = "Height (cm)", freq = FALSE)


###################################################
### code chunk number 19: SDaA_using_survey.Rnw:235-239
###################################################
  ### Figure 7.4
  data(htstrat)
  hist(htstrat$height, ylab = "Relative Frequency", 
      xlab = "Height (cm)", freq = FALSE)


###################################################
### code chunk number 20: SDaA_using_survey.Rnw:242-247
###################################################
	### Figure 7.5 (a)
  minht <- min(htstrat$height)
  breaks <- c(minht-1, seq(from = minht, to = max(htstrat$height), by = 1))
  hist(htstrat$height, ylab = expression(hat(f)(y)), breaks = breaks, 
       xlab = "Height Value, y", freq = FALSE)


###################################################
### code chunk number 21: SDaA_using_survey.Rnw:250-254
###################################################
  ### Figure 7.5 (b)
  stratecdf <- ecdf(htstrat$height)
  plot(stratecdf, do.points = FALSE, ylab = expression(hat(F)(y)), 
       xlab = "Height Value, y")


###################################################
### code chunk number 22: SDaA_using_survey.Rnw:259-262
###################################################
	### Figure 7.6
  data(syc)
  hist(syc$age, freq = FALSE, xlab = "Age")


###################################################
### code chunk number 23: SDaA_using_survey.Rnw:273-286
###################################################
  ### Figure 7.8
  sycdesign <- svydesign(ids= ~ psu, strata = ~ stratum,
     data = syc, weights=~finalwt)
  # p. 235: "Each of the 11 facilities with 360 or more youth
  # formed its own stratum (strata 6-16)", so in order
  # to avoid a lonely.psu error message
  #  Error in switch(lonely.psu, certainty = scale * crossprod(x), remove = scale *  : 
  #          Stratum (6) has only one PSU at stage 1
  # we set the option to "certainty" for this example
  # to see the problem, use: by(syc$psu, syc$stratum, unique) 
  oo <- options(survey.lonely.psu = "certainty")
  svyboxplot(age ~ factor(stratum), design = sycdesign) # mind the factor
  options(oo)


###################################################
### code chunk number 24: SDaA_using_survey.Rnw:293-298
###################################################
  ### Figure 7.9
  library(ggplot2)
  p <- ggplot(syc, aes(x = factor(stratum), y = factor(age)))
  g <- p + stat_sum(aes(group=1, weight = finalwt, size = ..n..)) 
  print(g)


###################################################
### code chunk number 25: SDaA_using_survey.Rnw:310-315
###################################################
	### Figure 7.10
  oo <- options(survey.lonely.psu = "certainty")
  sycstrat5 <- subset(sycdesign, stratum == 5)
  svyboxplot(age ~ factor(psu), design = sycstrat5)
  options(oo)


###################################################
### code chunk number 26: SDaA_using_survey.Rnw:318-323
###################################################
  ### Figure 7.11
  sycstrat5df <- subset(syc, stratum == 5)
  p <- ggplot(sycstrat5df, aes(x = factor(psu), y = factor(age)))
  g <- p + stat_sum(aes(group=1, weight = finalwt, size = ..n..)) 
  print(g)


###################################################
### code chunk number 27: SDaA_using_survey.Rnw:341-348
###################################################
  ### Example 10.1
  hh <- rbind(c(119, 188),
              c(88, 105))
  rownames(hh) <- c("cableYes", "cableNo")
  colnames(hh) <- c("computerYes", "computerNo")
  addmargins(hh)
  chisq.test(hh, correct = FALSE)  # OK      


###################################################
### code chunk number 28: SDaA_using_survey.Rnw:354-364
###################################################
	### Example 10.2 (nursing students and tutors)
  nst <- rbind(c(46, 222),
               c(41, 109),
               c(17, 40),
               c(8, 26))
  colnames(nst) <- c("NR", "R")
  rownames(nst) <- c("generalStudent", "generalTutor", "psychiatricStudent", 
      "psychiatricTutor")
  addmargins(nst)
  chisq.test(nst, correct = FALSE) # OK        


###################################################
### code chunk number 29: SDaA_using_survey.Rnw:367-376
###################################################
	### Example 10.3 (Air Force Pilots)
  afp <- data.frame(nAccidents = 0:7, 
                    nPilots = c(12475, 4117, 1016, 269, 53, 14, 6, 2))
  # estimate lambda
  lambdahat <- sum(afp$nAccidents * afp$nPilots  / sum(afp$nPilots))
  # expected counts
  observed <- afp$nPilots
  expected <- dpois(0:7, lambda = lambdahat) * sum(afp$nPilots)
  sum((observed - expected)^2 / expected) # NOT OK


###################################################
### code chunk number 30: SDaA_using_survey.Rnw:383-390
###################################################
	### Example 10.4
  hh2 <- rbind(c(238, 376),
               c(176, 210))
  rownames(hh2) <- c("cableYes", "cableNo")
  colnames(hh2) <- c("computerYes", "computerNo")
  addmargins(hh2)
  chisq.test(hh2, correct = FALSE)  # OK


###################################################
### code chunk number 31: SDaA_using_survey.Rnw:399-400
###################################################
	### example 10.5