################################################### ### chunk number 1: ################################################### options(SweaveHooks=list(fig=function(){par(cex.main=1.1, mar=c(4.1,3.6,1.6,1.1), pty="s", mgp=c(2.25,0.5,0), tck=-0.02)})) ################################################### ### chunk number 2: norm-samp ################################################### y <- rnorm(100) plot(density(y), type="l") ################################################### ### chunk number 3: sampdist-means ################################################### av <- numeric(100) # will take 100 means for (i in 1:100){ av[i] <- mean(rnorm(4)) } lines(density(av), col="red") ################################################### ### chunk number 4: survey-sample ################################################### library(MASS) y <- na.omit(survey[survey$Sex=="Female", "Height"]) av <- numeric(100) for(i in 1:100)av[i] <- mean(sample(y, 4, replace=TRUE)) ################################################### ### chunk number 5: miss-vs-nomiss ################################################### library(MASS) library(lattice) complete <- complete.cases(Pima.tr2) completeF <- factor(c("oneORmore", "none")[as.numeric(complete)+1]) Pima.tr2$completeF <- completeF densityplot(~bmi, groups=completeF, data=Pima.tr2, auto.key=list(columns=2)) # Use completeF in place of complete, for better labeling ################################################### ### chunk number 6: boot-miss ################################################### rownum <- seq(along=complete) # generate row numbers allpresSample <- sample(rownum[complete], replace=TRUE) # By default, sample size is the number of vales of the first argument NApresSample <- sample(rownum[!complete], replace=TRUE) densityplot(~bmi, groups=completeF, data=Pima.tr2, auto.key=list(columns=2), subset=c(allpresSample, NApresSample)) ################################################### ### chunk number 7: boot-diff-of-means ################################################### twot <- function(n=99){ complete <- complete.cases(Pima.tr2) rownum <- seq(along=complete) # generate row numbers d2 <- numeric(n+1) d2[1] <- with(Pima.tr2, mean(bmi[complete], na.rm=TRUE)- mean(bmi[!complete], na.rm=TRUE)) for(i in 1:n){ allpresSample <- sample(rownum[complete], replace=TRUE) NApresSample <- sample(rownum[!complete], replace=TRUE) d2[i+1] <- with(Pima.tr2, mean(bmi[allpresSample], na.rm=TRUE)- mean(bmi[NApresSample], na.rm=TRUE)) } d2 } # NB: There are n+1 values in all; the mean difference for the # actual sample values, plus differences for n bootstrap samples ## Now estimate and plot the density function d2 <- twot(n=999) dens <- density(d2) plot(dens) ## Check the proportion of values of d2 that are less than or equal to 0 sum(d2<0)/length(d2) ################################################### ### chunk number 8: 95% CI ################################################### round(sort(d2)[c(26, 975)], 2) ################################################### ### chunk number 9: ex-5a ################################################### library(DAAGxtras) sampvalues <- simulateSampDist(numINsamp=c(4,9)) plotSampDist(sampvalues=sampvalues, graph="density", titletext=NULL) ################################################### ### chunk number 10: sim-density-plot ################################################### getOption("SweaveHooks")[["fig"]]() sampvalues <- simulateSampDist(numINsamp=c(4,9)) plotSampDist(sampvalues=sampvalues, graph="density", titletext=NULL) ################################################### ### chunk number 11: density-qq eval=FALSE ################################################### ## sampvalues <- simulateSampDist() ## plotSampDist(sampvalues=sampvalues, graph=c("density", "qq")) ################################################### ### chunk number 12: density-qq-plot ################################################### getOption("SweaveHooks")[["fig"]]() sampvalues <- simulateSampDist() plotSampDist(sampvalues=sampvalues, graph=c("density", "qq")) ################################################### ### chunk number 13: Adelaide ################################################### library(MASS) y <- na.omit(survey[survey$Sex=="Female", "Height"]) sampvalues <- simulateSampDist(y) plotSampDist(sampvalues=sampvalues)