Programming for Applications

Chapter 16: Analyzing Data

Yu-You Liou (NTU)

Shih Chien University

2026-08-02

First Analyses

Summary Statistics: mean, min, max

The building blocks for everything that follows. With the dow30 quote data from Chapter 12 (packaged in nutshell):

temp <- tempfile(fileext = ".rda")
download.file("https://raw.githubusercontent.com/cran/nutshell/master/data/dow30.rda",
              temp, mode = "wb")
load(temp)
mean(dow30$Open); min(dow30$Open); max(dow30$Open)
[1] 36.24574
[1] 0.99
[1] 122.45

Two arguments worth memorizing:

mean(c(1, 2, 3, 4, 5, NA))                  # any NA poisons the result...
[1] NA
mean(c(1, 2, 3, 4, 5, NA), na.rm=TRUE)      # ...unless you remove it
[1] 3
mean(c(-1, 0:100, 2000))                    # outliers drag the mean...
[1] 68.43689
mean(c(-1, 0:100, 2000), trim=0.1)          # ...trim filters that fraction
[1] 50

range, quantile, fivenum, IQR

range(dow30$Open)
[1]   0.99 122.45
quantile(dow30$Open, probs=c(0, 0.25, 0.5, 0.75, 1.0))
     0%     25%     50%     75%    100% 
  0.990  19.655  30.155  51.680 122.450 
fivenum(dow30$Open)
[1]   0.990  19.650  30.155  51.680 122.450
IQR(dow30$Open)
[1] 32.025

All of these combine naturally with apply/tapply for per-group statistics (Chapter 12).

summary and str

The most convenient overview is generic summary — quartiles + mean for numerics, top counts for factors (nothing meaningful for plain character vectors):

summary(dow30[, c("Open", "Volume")])
      Open            Volume         
 Min.   :  0.99   Min.   :1.336e+06  
 1st Qu.: 19.66   1st Qu.:1.111e+07  
 Median : 30.16   Median :1.822e+07  
 Mean   : 36.25   Mean   :5.226e+07  
 3rd Qu.: 51.68   3rd Qu.:4.255e+07  
 Max.   :122.45   Max.   :2.672e+09

The popular alternative str displays an object’s structure:

str(dow30)
'data.frame':   7482 obs. of  8 variables:
 $ symbol   : Factor w/ 30 levels "MMM","AA","AXP",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Date     : Factor w/ 252 levels "2008-09-22","2008-09-23",..: 252 251 250 249 248 247 246 245 244 243 ...
 $ Open     : num  73.9 75.1 75.3 74.8 74.6 ...
 $ High     : num  74.7 75.2 75.5 75.5 74.9 ...
 $ Low      : num  73.9 74.5 74.5 74.5 74 ...
 $ Close    : num  74.5 74.6 74.9 75.4 74.7 ...
 $ Volume   : num  2560400 4387900 3371500 2722500 3566900 ...
 $ Adj.Close: num  74.5 74.6 74.9 75.4 74.7 ...

stem: a Text Histogram

stem(x, scale=1, width=80, atom=1e-08) draws a stem-and-leaf plot right in the console. Distances of missed NFL field goals, 2005:

temp <- tempfile(fileext = ".rda")
download.file("https://raw.githubusercontent.com/cran/nutshell/master/data/field.goals.rda",
              temp, mode = "wb")
load(temp)
stem(subset(field.goals, play.type=="FG no")$yards)

  The decimal point is at the |

  20 | 0
  22 | 
  24 | 
  26 | 00
  28 | 0000000
  30 | 0000000
  32 | 00000000
  34 | 000
  36 | 0000
  38 | 00000000000000
  40 | 0000000000
  42 | 0000000000000000
  44 | 000000000000
  46 | 000000000000000000
  48 | 000000000000000000
  50 | 000000000000
  52 | 0000000000000000000
  54 | 0000
  56 | 000
  58 | 00
  60 | 00
  62 | 0

Correlation and Covariance

Three Correlation Measures

Correlation answers “when x rises, does y rise — and how tightly?”, measuring linear dependence on a −1…1 scale. Three statistics:

  • Pearson — Excel’s CORREL; rooted in normal-distribution theory, best for normal-ish data;
  • Spearman — nonparametric, no distributional assumptions (correlates the ranks);
  • Kendall’s tau — compares concordant vs. discordant pairs of rankings.
cor(x, y = NULL, use = "everything",
    method = c("pearson", "kendall", "spearman"))

x,y may be two vectors, or one data frame/matrix (→ a correlation matrix). use handles NAs: "all.obs" (error), "everything" (NA result), "complete.obs" / "na.or.complete" (drop incomplete rows), "pairwise.complete.obs" (drop per pair).

Example: Weight Gain and Birth Weight

Does the mother’s weight gain correlate with the baby’s weight? Clean first (valid values, single births, gestation > 35 weeks), then look — 100k+ points need smoothScatter, not plot:

temp <- tempfile(fileext = ".rda")
download.file("https://raw.githubusercontent.com/cran/nutshell/master/data/births2006.smpl.rda",
              temp, mode = "wb")
load(temp)
births2006.cln <- births2006.smpl[
  !is.na(births2006.smpl$WTGAIN) & !is.na(births2006.smpl$DBWT) &
  births2006.smpl$DPLURAL == "1 Single" & births2006.smpl$ESTGEST>35,]
smoothScatter(births2006.cln$WTGAIN, births2006.cln$DBWT)
cor(births2006.cln$WTGAIN, births2006.cln$DBWT)
[1] 0.1751866
cor(births2006.cln$WTGAIN, births2006.cln$DBWT, method="spearman")
[1] 0.1776192

A modest correlation by either measure — the blob is angled, just barely. (The book used the full 3.2M-record file; our 10% sample gives near-identical values.)

Covariance

Covariance is the Pearson numerator — same arguments as cor:

cov(births2006.cln$WTGAIN, births2006.cln$DBWT)
[1] 1176.469

Relatives: cov2cor converts a covariance matrix to a correlation matrix; cov.wt(x, wt=, cor=, center=, method=c("unbiased","ML")) computes weighted covariance.

Principal Components Analysis

prcomp and princomp

PCA rewrites possibly-correlated variables as uncorrelated components. Two implementations: prcomp (preferred; formula method prcomp(~..., data=, subset=, na.action=), default method with retx, center=TRUE, scale., tol) and the S-PLUS-compatible princomp (eigendecomposition of the correlation/covariance matrix) — the book demonstrates the latter.

Data: team batting statistics 2000–2008, extracted from the Baseball Databank by SQL (renamed to readable play names) and packaged as team.batting.00to08:

temp <- tempfile(fileext = ".rda")
download.file("https://raw.githubusercontent.com/cran/nutshell/master/data/team.batting.00to08.rda",
              temp, mode = "wb")
load(temp)
batting.pca <- princomp(~singles+doubles+triples+homeruns+
                          walks+hitbypitch+sacrificeflies+
                          stolenbases+caughtstealing,
                        data=team.batting.00to08)
summary(batting.pca)
Importance of components:
                           Comp.1     Comp.2     Comp.3      Comp.4      Comp.5
Standard deviation     74.9009809 61.8710858 31.8113983 27.98819003 23.78885885
Proportion of Variance  0.4610727  0.3146081  0.0831687  0.06437897  0.04650949
Cumulative Proportion   0.4610727  0.7756807  0.8588494  0.92322841  0.96973790
                            Comp.6      Comp.7      Comp.8      Comp.9
Standard deviation     12.88429066 9.150840397 8.283972499 7.060503344
Proportion of Variance  0.01364317 0.006882026 0.005639904 0.004096998
Cumulative Proportion   0.98338107 0.990263099 0.995903002 1.000000000

Two components already carry ~78% of the variance; four carry ~92%.

loadings

loadings shows each variable’s contribution to each component:

loadings(batting.pca)

Loadings:
               Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9
singles         0.313  0.929  0.136         0.136                            
doubles                       0.437  0.121 -0.877                0.100       
triples                                                  -0.424 -0.775  0.449
homeruns       -0.235         0.383  0.825  0.324                            
walks          -0.914  0.328 -0.150 -0.182                                   
hitbypitch                                         0.989                     
sacrificeflies                                           -0.321 -0.330 -0.882
stolenbases            0.131 -0.758  0.502 -0.307         0.232              
caughtstealing               -0.208  0.104               -0.813  0.521  0.105

               Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9
SS loadings     1.000  1.000  1.000  1.000  1.000  1.000  1.000  1.000  1.000
Proportion Var  0.111  0.111  0.111  0.111  0.111  0.111  0.111  0.111  0.111
Cumulative Var  0.111  0.222  0.333  0.444  0.556  0.667  0.778  0.889  1.000

Scree Plot and Biplot

plot(batting.pca)                                      # "scree" plot

biplot(batting.pca, cex=0.5, col=c("gray50", "black")) # variables + observations

The biplot draws variable contributions to components 1–2 and the observations on one scale: singles and walks dominate the first two components. (This data returns for linear modeling in Chapter 20.)

Factor Analysis

factanal

Some quantities are observable (test scores); some are not (intelligence). Factor analysis infers hidden factors from observed variables — factanal in stats:

factanal(x, factors, data = NULL, covmat = NULL, n.obs = NA,
         subset, na.action, start = NULL,
         scores = c("none", "regression", "Bartlett"),
         rotation = "varimax", control = NULL, ...)
Argument Description
x / data formula or numeric matrix (+ data frame for the formula)
factors how many factors to fit
covmat / n.obs alternatively: a covariance matrix (or cov.wt output) + observation count
subset / na.action / start the usual suspects
scores "none", "regression" (Thompson), "Bartlett" (weighted least squares)
rotation factor rotation function, default "varimax"

Bootstrap Resampling

The Idea, and boot

Is your statistic sensitive to a few outliers? What’s its sampling range? Bootstrapping answers for any statistic: repeatedly resample observations with replacement, recompute the statistic, study the spread — formally, an estimate of an estimator’s bias.

library(boot)
boot(data, statistic, R, sim="ordinary", stype="i",
     strata=rep(1,n), L=NULL, m=0, weights=NULL,
     ran.gen=function(d, p) d, mle=NULL, simple=FALSE, ...)

Key arguments: data; statistic — a function of (data, selection) where stype defines the second argument ("i" indices, "f" frequencies, "w" weights); R — number of replicates; sim — "ordinary", "parametric", "balanced", "permutation", "antithetic"; strata for multisample problems; ran.gen/mle for parametric simulation; simple=TRUE trades time for memory.

Example: How Biased Is the Median?

Media report median home prices. Bootstrapping June 2008 San Francisco sale prices:

temp <- tempfile(fileext = ".rda")
download.file("https://raw.githubusercontent.com/cran/nutshell/master/data/sanfrancisco.home.sales.rda",
              temp, mode = "wb")
load(temp)
home.sale.prices.june2008 <- subset(sanfrancisco.home.sales,
  format(date, "%Y-%m")=="2008-06")$price
library(boot)
b <- boot(data=home.sale.prices.june2008,
          statistic=function(d, i) { median(d[i]) },
          R=1000)
b

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = home.sale.prices.june2008, statistic = function(d, 
    i) {
    median(d[i])
}, R = 1000)


Bootstrap Statistics :
    original  bias    std. error
t1*   845000 -3990.5     22636.2

Verdict: the median is a very slightly biased estimator here — and now we also know its standard error, free of any normality assumption.