\documentclass[11pt]{article} \usepackage[margin=1in,bottom=.5in,includehead,includefoot]{geometry} \usepackage{hyperref} \usepackage{language} \usepackage{alltt} \usepackage{fancyhdr} \pagestyle{fancy} \fancyhf{} %% Now begin customising things. See the fancyhdr docs for more info. \chead{} \lhead[\sf \thepage]{\sf \leftmark} \rhead[\sf \leftmark]{\sf \thepage} \lfoot{} \cfoot{Statistical Sleuth in R: Chapter 2} \rfoot{} \newcounter{myenumi} \newcommand{\saveenumi}{\setcounter{myenumi}{\value{enumi}}} \newcommand{\reuseenumi}{\setcounter{enumi}{\value{myenumi}}} \pagestyle{fancy} \def\R{{\sf R}} \def\Rstudio{{\sf RStudio}} \def\RStudio{{\sf RStudio}} \def\term#1{\textbf{#1}} \def\tab#1{{\sf #1}} \usepackage{relsize} \newlength{\tempfmlength} \newsavebox{\fmbox} \newenvironment{fmpage}[1] { \medskip \setlength{\tempfmlength}{#1} \begin{lrbox}{\fmbox} \begin{minipage}{#1} \vspace*{.02\tempfmlength} \hfill \begin{minipage}{.95 \tempfmlength}} {\end{minipage}\hfill \vspace*{.015\tempfmlength} \end{minipage}\end{lrbox}\fbox{\usebox{\fmbox}} \medskip } \newenvironment{boxedText}[1][.98\textwidth]% {% \begin{center} \begin{fmpage}{#1} }% {% \end{fmpage} \end{center} } \newenvironment{boxedTable}[2][tbp]% {% \begin{table}[#1] \refstepcounter{table} \begin{center} \begin{fmpage}{.98\textwidth} \begin{center} \sf \large Box~\expandafter\thetable. #2 \end{center} \medskip }% {% \end{fmpage} \end{center} \end{table} % need to do something about exercises that follow boxedTable } \newcommand{\cran}{\href{http://www.R-project.org/}{CRAN}} \title{The Statistical Sleuth in R: \\ Chapter 2} \author{ Ruobing Zhang\and Kate Aloisio \and Nicholas J. Horton\thanks{Department of Mathematics, Amherst College, nhorton@amherst.edu} } \date{\today} \begin{document} \maketitle \tableofcontents %\parindent=0pt <>= opts_chunk$set( dev="pdf", fig.path="figures/", fig.height=3, fig.width=4, out.width=".47\\textwidth", fig.keep="high", fig.show="hold", fig.align="center", prompt=TRUE, # show the prompts; but perhaps we should not do this comment=NA # turn off commenting of ouput (but perhaps we should not do this either ) @ <>= print.pval = function(pval) { threshold = 0.0001 return(ifelse(pval < threshold, paste("p<", sprintf("%.4f", threshold), sep=""), ifelse(pval > 0.1, paste("p=",round(pval, 2), sep=""), paste("p=", round(pval, 3), sep="")))) } @ <>= require(Sleuth2) require(mosaic) trellis.par.set(theme=col.mosaic()) # get a better color scheme set.seed(123) # this allows for code formatting inline. Use \Sexpr{'function(x,y)'}, for exmaple. knit_hooks$set(inline = function(x) { if (is.numeric(x)) return(knitr:::format_sci(x, 'latex')) x = as.character(x) h = knitr:::hilight_source(x, 'latex', list(prompt=FALSE, size='normalsize')) h = gsub("([_#$%&])", "\\\\\\1", h) h = gsub('(["\'])', '\\1{}', h) gsub('^\\\\begin\\{alltt\\}\\s*|\\\\end\\{alltt\\}\\s*$', '', h) }) showOriginal=FALSE showNew=TRUE @ \section{Introduction} This document is intended to help describe how to undertake analyses introduced as examples in the Second Edition of the \emph{Statistical Sleuth} (2002) by Fred Ramsey and Dan Schafer. More information about the book can be found at \url{http://www.proaxis.com/~panorama/home.htm}. This file as well as the associated \pkg{knitr} reproducible analysis source file can be found at \url{http://www.amherst.edu/~nhorton/sleuth}. This work leverages initiatives undertaken by Project MOSAIC (\url{http://www.mosaic-web.org}), an NSF-funded effort to improve the teaching of statistics, calculus, science and computing in the undergraduate curriculum. In particular, we utilize the \pkg{mosaic} package, which was written to simplify the use of R for introductory statistics courses. A short summary of the R needed to teach introductory statistics can be found in the mosaic package vignette (\url{http://cran.r-project.org/web/packages/mosaic/vignettes/MinimalR.pdf}). To use a package within R, it must be installed (one time), and loaded (each session). The package can be installed using the following command: <>= install.packages('mosaic') # note the quotation marks @ Once this is installed, it can be loaded by running the command: <>= require(mosaic) @ This needs to be done once per session. In addition the data files for the \emph{Sleuth} case studies can be accessed by installing the \pkg{Sleuth2} package. <>= install.packages('Sleuth2') # note the quotation marks @ <>= require(Sleuth2) @ We also set some options to improve legibility of graphs and output. <>= trellis.par.set(theme=col.mosaic()) # get a better color scheme for lattice options(digits=3, show.signif.stars=FALSE) @ The specific goal of this document is to demonstrate how to calculate the quantities described in Chapter 2: Inference Using \emph{t}-Distributions using R. \section{Bumpus's Data on Natural Selection} Is humerus length related to whether the bird would survive or perish? That's the question being addressed by Case Study 2.1 in the \emph{Sleuth}. \subsection{Statistical summary and graphical display} We begin by reading the data and summarizing the variables. <<>>= summary(case0201) fav=favstats(Humerus ~ Status, data=case0201); fav @ A total of \Sexpr{nrow(case0201)} subjects are included in the data: \Sexpr{nrow(subset(case0201, Status=="Survived"))} are adult male sparrows that survived and \Sexpr{nrow(subset(case0201, Status=="Perished"))} that perished. The following figure replicates Display 2.1 on page 29. <<>>= bwplot(Status ~ Humerus, data=case0201) @ <<>>= densityplot(~ Humerus, groups=Status, auto.key=TRUE, data=case0201) @ Both distributions are approximately normally distributed. \subsection{Inferential procedures (two-sample t-test)} First, we calculate the pooled SD and the standard error between these two different sample average (page 40, Display 2.8). <<>>= # Calculate Pooled SD n1 = fav["Perished", "n"]; n1 n2 = fav["Survived", "n"]; n2 s1 = fav["Perished", "sd"]; s1 s2 = fav["Survived", "sd"]; s2 Sp = sqrt(((n1-1)*(s1)^2+(n2-1)*(s2)^2)/(n1+n2-2)); Sp # Calculate standard error SE = Sp*sqrt(1/n1+1/n2); SE @ So the pooled SD is \Sexpr{round(Sp, 2)} and the standard error is \Sexpr{round(SE, 1)}. Based on this information, we can construct a 95\% confidence interval (page 41, Display 2.9). <<>>= Y1 = fav["Perished", "mean"]; Y1 Y2 = fav["Survived", "mean"]; Y2 Yd = Y2-Y1; Yd df = n1+n2-2; df qt = qt(0.975, df); qt hw = qt*SE; hw lower = Yd-hw; lower upper = Yd+hw; upper @ So the 95\% confidence interval of the difference between means is (\Sexpr{round(lower, 1)}, \Sexpr{round(upper, 1)}) Now we want to calculate the $t$-statistic and $p$-value (as shown on page 44, Display 2.10). <<>>= tstats = (Yd-0)/SE; tstats # The hypothesis difference=0 onepval = 1-pt(tstats, df); onepval twopval = 2*onepval; twopval @ The one-sided $p$-value is \Sexpr{round(onepval, 2)} and the two-sided $p$-value is \Sexpr{round(twopval, 2)}. We can get the results of ``Summary of Statistical Findings" (page 29) by using the following code: <<>>= t.test(Humerus ~ Status, var.equal=TRUE, data=case0201) confint(lm(Humerus ~ Status, data=case0201)) @ \section{Anatomical Abnormalities Associated with Schizophrenia} Is the area of brain related to the development of schizophrenia? That's the question being addressed by case study 2.2 in the \emph{Sleuth}. \subsection{Statistical summary and graphical display} We begin by reading the data and summarizing the variables. <<>>= summary(case0202) @ A total of \Sexpr{nrow(case0202)} subjects are included in the data. There are \Sexpr{nrow(case0202[ "Affected"])} pairs of twins; one of the twins has schizophrenia, and the other does not. So there are \Sexpr{nrow(case0202["Affected"])} affected subjects and \Sexpr{nrow(case0202["Unaffect"])} unaffected subjects. The difference in area of left hippocampus of these pairs of twins is: <<>>= DIFF = case0202[, "Unaffect"]-case0202[, "Affected"] favstats(DIFF) @ This matches the results on page 30, Display 2.2. <<>>= densityplot(DIFF) @ \subsection{Inferential procedures (two-sample t-test)} We want to calculate the paired t-test and 95\% confidence interval. <<>>= # Calculate t-statistics difmean = favstats(DIFF)[, "mean"]; difmean difsd = favstats(DIFF)[, "sd"]; difsd difSE = difsd/sqrt(15); difSE tscore = (difmean-0)/difSE; tscore # hypothesis difference=0 twopvalue = 2*(1-pt(tscore, 15-1)); twopvalue # Construct confidence interval q = qt(0.975, 15-1); q schizolower = difmean-q*difSE; schizolower schizoupper = difmean+q*difSE; schizoupper @ So the two-sided $p$-value is \Sexpr{round(twopvalue, 3)} and the 95\% confidence interval is (\Sexpr{round(schizolower, 2)}, \Sexpr{round(schizoupper, 2)}). Or we can get the results displayed on page 31 by conducting a paired $t$-test. <<>>= t.test(case0202[, "Unaffect"], case0202[, "Affected"], paired=TRUE) @ \end{document}