TECH3100 · Data Visualisation in R

Missing values

data frames, imputation, interpolation, and data cleaning
Lesson 8
Kaplan Business School Australia
All figures and numbers computed in base R 4.3.3 from airquality
0.1

0.1 Copyright notice

Commonwealth of Australia · Copyright Regulations 1969

This material has been reproduced and communicated to you by or on behalf of Kaplan Business School pursuant to Part VB of the Copyright Act 1968 (the Act).

The material in this communication may be subject to copyright under the Act. Any further reproduction or communication of this material by you may be the subject of copyright protection under the Act.

Do not remove this notice.

0.2

0.2 Where this week sits

WeekTopic
1–4Variable types. Descriptive statistics. Histograms, boxplots, scatter plots, correlation.
5Group presentations. Linear trend line.
6Visualisation for time series.
7Visualisation for multivariate data.
8Data frames. Missing values: imputation and interpolation. Data cleaning.
9Aggregates and joining tables.
10–11Data visualisation with ggplot.
12Individual presentations.
Picking up from last week Week 7 ended with an exercise on airquality and the question of what cor(airquality) returns before you do anything about the missing values. This lesson answers it, and then works out what to do next.
0.3

0.3 What you should be able to do by the end

  1. Build, import, inspect and subset a data frame, and perform the five standard operations on one in base R.
  2. Explain what NA means in R and predict how it propagates through a calculation.
  3. Produce three visual displays of missingness and read what each one shows.
  4. Distinguish missing completely at random, missing at random, and missing not at random, and say which evidence would support each.
  5. Apply five imputation methods and state what each does to the mean, the spread and the correlations.
  6. Apply three interpolation methods to a gapped time series, and test them against known values before trusting one.
The one that matters most Outcome 5 and outcome 6 both end in the same place. Filling a gap adds no information. Every method in this lesson makes an assumption, and your job is to name it.
0.4

0.4 One dataset for the whole lesson

airquality ships with R. Daily air-quality measurements from New York, 1 May to 30 September 1973.

ColumnMeaningMissing
OzoneMean ozone in parts per billion, 13:00 to 15:00 37 of 153
Solar.RSolar radiation in langleys, 08:00 to 12:00 7 of 153
WindAverage wind speed in miles per hour0
TempMaximum daily temperature in degrees Fahrenheit0
MonthMonth, 5 to 90
DayDay of month, 1 to 310

Why this one

  • It is genuinely incomplete. The gaps were not manufactured for teaching.
  • It is a time series, so the rows have an order that carries information. That is what makes interpolation available.
  • It is small enough to print and large enough that the choices matter.
Load it
data(airquality)
aq <- airquality
dim(aq)
#> [1] 153   6
1
Section 1

Working with data frames

A data frame is the object every dataset in this unit lives in. Five operations cover almost everything you will do to one, and all five exist in base R.

1.1

1.1 Anatomy of a data frame

Rows are observations. Columns are variables. Each column holds one type, and that is what separates a data frame from a spreadsheet.

Ozone<int>Solar.R<int>Wind<num>Temp<int>Month<int>row 1411907.4675row 2361188.0725row 31214912.6745row 41831311.5625row 5NANA14.3565aq$Ozoneone whole column,returned as a vectoraq[2, "Temp"]one cellaq[5, ]one rowOne column = one variable, and every value in it shares one type.One row = one observation. Here, one day of measurements at one station.NA marks a value that was not recorded. It is not zero, and it is not an empty string.
1.2

1.2 Building one and reading one in

Build from vectors

students <- data.frame(
  name    = c("John Smith", "Jane Doe", "Joe Schmo"),
  age     = c(34, 28, 51),
  address = c("123 Main St.", "456 Maple Ave.",
              "789 Broadway"),
  stringsAsFactors = FALSE
)

Read and write a CSV

df <- read.csv("my_data.csv")
write.csv(df, "new_data.csv", row.names = FALSE)

# read.csv arguments worth knowing
read.csv("f.csv",
         stringsAsFactors = FALSE,
         na.strings = c("", "NA", "-999"),
         strip.white = TRUE)
na.strings earns its place Datasets in the wild mark absence with an empty cell, the text "N/A", a full stop, or a sentinel number such as −999. If you do not declare these, R reads them as text or as a real value, and −999 will quietly drag your mean down.

The packages the source slides use

install.packages(c("readr", "dplyr", "tidyr"))
library(readr)
library(dplyr)

df <- read_csv("my_data.csv")
write_csv(df, "new_data.csv")

What differs

  • read_csv returns a tibble, prints ten rows, and never converts text to a factor.
  • read.csv returns a plain data frame and is available with no installation.
  • Both accept na or na.strings for custom missing markers.
For this unit Every example in these slides runs in base R, so nothing has to be installed and nothing can break in a fresh Colab session. The dplyr version is here because you will meet it constantly outside this room.
1.3

1.3 Inspect before you touch anything

Six functions, in this order

FunctionWhat it answers
dim(aq)How many rows and columns?
str(aq)What type is each column, and what do the first few values look like?
head(aq, 10)Does the top of the file look like data or like a stray header?
summary(aq)Range, quartiles, and the count of NAs per column.
colnames(aq)Are the names usable, or full of spaces?
colSums(is.na(aq))Exactly how many values are missing where?

What summary() gives you free

summary(aq$Ozone)
#>  Min. 1st Qu. Median  Mean 3rd Qu.  Max.  NA's
#>  1.00   18.00  31.50 42.13   63.25 168.00   37

The final column appears only when the variable contains missing values. It is the fastest missingness check in R, and it comes with the summary you wanted anyway.

Do this every time Running str() and summary() before any analysis takes ten seconds and catches the numeric column that arrived as text, the date that arrived as a factor, and the column that is entirely NA.
1.4

1.4 The five operations

select()keeps columnsOzoneTempMonth4167536725127451862523656filter()keeps rowsOzoneTempMonth4167536725127451862523656arrange()reorders rowsOzoneTempMonth1274518625236563672541675sorted on the first columnmutate()adds a columnOzoneTempMonthlevel41675hi36725hi12745lo18625lo23656hisummarise()collapses to one rowmeanmeann26.068.055 rows in, 1 row outEach verb takes a data frameand returns a data frame.That is what makes themchainable.
Figure 1.4. Each operation takes a data frame and returns a data frame. Highlighted cells show what each one keeps.
OperationBase R
selectaq[, c("Ozone", "Temp")]
filteraq[aq$Month == 5, ]
subset(aq, Month == 5)
arrangeaq[order(-aq$Temp), ]
mutateaq$TempC <- (aq$Temp - 32) * 5/9
transform(aq, TempC = ...)
summariseaggregate(Ozone ~ Month, aq, mean)
The one trap aq[aq$Month == 5, ] keeps rows where the test is NA and fills them with NA. subset() drops them. The two give different row counts on incomplete data, so choose one deliberately.
# select: keep columns sel <- aq[, c("Ozone", "Temp", "Month")] # filter: keep rows fil <- subset(aq, Month == 5 & !is.na(Ozone)) # arrange: reorder rows arr <- aq[order(-aq$Temp), ] # mutate: add a column mut <- transform(aq, TempC = round((Temp - 32) * 5/9, 1)) # summarise: collapse to one row per group aggregate(Ozone ~ Month, data = aq, FUN = mean) #> Month Ozone #> 1 5 23.61538 #> 2 6 29.44444 #> 3 7 59.11538 #> 4 8 59.96154 #> 5 9 31.44828 # aggregate() drops NA rows by default, which is # why the group sizes below are not 31, 30, 31...
1.5

1.5 Chaining operations

R 4.1 added a pipe to the language itself. It needs no package.

# without a pipe: nested, read inside out
head(aq[order(-aq$Temp), c("Ozone", "Temp")], 5)

# with the native pipe: read left to right
aq |>
  subset(Month == 7) |>
  transform(hot = Temp > 85) |>
  head(3)

# x |> f(y) is exactly f(x, y)
# The value on the left becomes the FIRST argument.
What a pipe buys Nothing that intermediate variables could not do. It removes the naming step and puts the operations in the order you would say them out loud, which makes a long chain easier to check.

The magrittr pipe %>% predates the native one and behaves almost identically for these purposes.

library(dplyr)

aq %>%
  filter(Month == 7) %>%
  select(Ozone, Temp) %>%
  arrange(desc(Temp)) %>%
  mutate(hot = Temp > 85) %>%
  head(3)

# Differences worth knowing:
#   %>% allows the placeholder  . anywhere
#   |>  puts the value first, and nowhere else
#   %>% needs magrittr or dplyr loaded
#   |>  is part of the language
In an exam or a submission The native pipe works in any R 4.1 or later session with nothing installed. Use it unless you have been told otherwise.
1.6

1.6 Tidy data: long and wide

Tidy means one variable per column and one observation per row. Neither layout below is more correct than the other; each is tidy for a different question.

Wideone row per day, one column per variabledayOzoneTemp141672367231274Longone row per measurementdayvariablevalue1Ozone411Temp672Ozone362Temp723Ozone123Temp74reshapedirection ="long"Tidy means one variable per column and one observation per row. Which layout is tidydepends on what you call an observation. For a chart of Ozone against Temp, wide is tidy.For a chart that facets by variable, long is tidy. Reshaping loses nothing; it renames.
1.7

1.7 Duplicates, types and strings

Duplicates

sum(duplicated(df))     # how many repeats
df[!duplicated(df), ]   # keep the first of each
unique(df)              # the same thing

# duplicated on chosen columns only
df[!duplicated(df[, c("id", "date")]), ]

Types

str(df)                 # what R thinks each is
df$price <- as.numeric(df$price)
df$grp   <- as.factor(df$grp)
df$date  <- as.Date(df$date, "%d/%m/%Y")

Strings into numbers

p <- c("$1,200", "$980", "$1,050")
as.numeric(gsub("[$,]", "", p))
#> [1] 1200  980 1050

# splitting a fixed-width field
bd <- c("03151998", "11021985")
substr(bd, 1, 2)   #> "03" "11"
substr(bd, 3, 4)   #> "15" "02"
substr(bd, 5, 8)   #> "1998" "1985"

# splitting on a separator
strsplit("guest_AU", "_")[[1]]
#> [1] "guest" "AU"
Order matters Strip the characters first, then convert. Calling as.numeric() on "$1,200" returns NA with a warning, and you have just created a missing value that was never missing.
1.Q

1.Q Knowledge check: Section 1

Q1. aq[aq$Ozone > 100, ] and subset(aq, Ozone > 100) return different numbers of rows. Why?
A comparison against NA returns NA, and NA used as a row index produces a row of NAs. subset() treats NA as FALSE. On this dataset the difference is 37 phantom rows.
Q2. Which layout is tidy for a scatterplot of Ozone against Temperature?
The observation is a day, and the two variables the plot needs are two columns of that row. For a chart with one panel per variable, the long layout becomes the tidy one.
Q3. A price column arrives as the text "$1,200". What does as.numeric(price) return?
Nothing stops. You now hold a missing value that the source data never had. Strip the dollar sign and the comma with gsub() first, then convert.
2
Section 2

Seeing what is missing

Before you fill a single value, find out how much is missing, where it sits, and whether the pattern is related to anything you can observe.

2.1

2.1 What NA is, and how it spreads

x <- c(41, 36, NA, 18)mean(x)→NAone missing value poisons the whole resultmean(x, na.rm = TRUE)→31.667the NA is dropped, and n falls from 4 to 3sum(is.na(x))→1counts how many are missingx == NA→NA NA NA NAcomparison with NA is never TRUEis.na(x)→F F T Fthis is the test that worksNA means unknown. R refuses to guess, so any arithmetic touching an NA returns NA.
Figure 2.1. Every line was run in R 4.3.3.
  1. It is not zeroZero is a recorded measurement. NA is the absence of one. Averaging a column that treats absence as zero pulls the mean towards zero.
  2. It is not an empty string"" is a value of length zero. It is a string that R will happily sort, count and paste.
  3. It is not NULLNULL has length zero and disappears from a vector. NA occupies a position and keeps the length.
  4. It is not NaNNaN is the result of an undefined calculation such as 0/0. is.na(NaN) is TRUE, but is.nan(NA) is FALSE.
  5. It propagatesAny arithmetic touching an NA returns NA, because R will not guess. na.rm = TRUE drops it and silently reduces n.
  6. It cannot be comparedx == NA returns NA, never TRUE. The only test that works is is.na(x).
x <- c(41, 36, NA, 18)

mean(x)                 #> NA
mean(x, na.rm = TRUE)   #> 31.66667
length(x)               #> 4
sum(!is.na(x))          #> 3

is.na(x)                #> F F T F
x == NA                 #> NA NA NA NA

# Distinguishing the four
is.na(NA);   is.na(NaN)     #> TRUE TRUE
is.nan(NA);  is.nan(NaN)    #> FALSE TRUE
length(c(1, NULL, 3))       #> 2
length(c(1, NA,   3))       #> 3

# na.rm changes the denominator. Report it.
mean(aq$Ozone, na.rm = TRUE)  #> 42.12931
sum(!is.na(aq$Ozone))         #> 116
2.2

2.2 Look at the whole table at once

Ozone37 NASolar.R7 NAWind0 NATemp0 NAMonth0 NADay0 NAMay 1Jun 1Jul 1Aug 1Sep 1One thin band = one day. 153 days, 6 columns, 918 cells.missingrecordedThe June block in the Ozone column is the thing to notice. Missingness hereis not scattered at random through the year.
Figure 2.2. Every cell of airquality, coloured by whether it was recorded. 153 rows by 6 columns.
  1. UnitOne thin horizontal band is one day. The six vertical strips are the six columns.
  2. EncodingRed marks a missing value, pale blue marks a recorded one. Nothing else is encoded, so the whole figure is one binary variable.
  3. ScaleRows run in date order from 1 May at the top to 30 September at the bottom. That ordering is what makes the block structure visible.
  4. StructureMissingness is concentrated in two columns. Wind, Temp, Month and Day are complete.
  5. Groups and exceptionsThe Ozone column carries a solid block through the second half of June. Solar.R carries one shorter block in early August. The two blocks do not line up.
  6. Claim and limitOzone is missing in a sustained run rather than at scattered single days. The figure shows where values are missing. It says nothing about why, and nothing about whether the missing days differ from the recorded ones.
# A missingness matrix in one call.
# t() and the reversed index put day 1 at the top.
image(1:ncol(aq), 1:nrow(aq),
      t(is.na(aq))[, nrow(aq):1],
      col  = c("#dfeaf2", "#CC0000"),
      axes = FALSE, xlab = "", ylab = "Day")
axis(1, 1:ncol(aq), colnames(aq), las = 2)
box()

# The same information as text
colSums(is.na(aq))
#>   Ozone Solar.R    Wind    Temp   Month     Day
#>      37       7       0       0       0       0

# Read this chart BEFORE you decide anything.
# A block and a scatter of single days call for
# completely different treatments.
2.3

2.3 Count it two ways

Missing values per columnOzone37 (24.2%)Solar.R7 (4.6%)Wind0 (0.0%)Temp0 (0.0%)Month0 (0.0%)Day0 (0.0%)What listwise deletion costsCells missing4.79%44 of 918 valuesRows lost27.45%42 of 153 daysUnder 5 per cent of the values are missing. Dropping every incomplete rowremoves more than a quarter of the days.
Figure 2.3. Missing values per column, and what dropping every incomplete row would cost.
  1. UnitThe upper panel counts cells. The lower panel counts rows. They answer different questions.
  2. EncodingBar length gives the count in the upper panel and the percentage in the lower one.
  3. ScaleThe lower bars share a common 0 to 100 per cent scale, so the two are directly comparable.
  4. StructureOnly two of the six columns contain any missing values at all. Ozone accounts for 37 of the 44 missing cells.
  5. Groups and exceptions44 missing cells out of 918 is 4.79 per cent. Those 44 cells are spread across 42 different days, which is 27.45 per cent of the rows.
  6. Claim and limitListwise deletion costs far more than the raw missing rate suggests. The multiplier depends entirely on how the missing cells are spread across rows, so this figure does not generalise to another dataset.
# cells
sum(is.na(aq))        #> 44
prod(dim(aq))         #> 918
100 * 44 / 918        #> 4.793028

# rows
sum(!complete.cases(aq))          #> 42
100 * mean(!complete.cases(aq))   #> 27.45098

# per column
colSums(is.na(aq))
round(100 * colMeans(is.na(aq)), 1)
#>   Ozone Solar.R    Wind    Temp   Month     Day
#>    24.2     4.6     0.0     0.0     0.0     0.0

barplot(colSums(is.na(aq)), las = 2,
        col = "#CC0000", ylab = "Missing values")

# The gap between 4.8% and 27.5% is the whole
# argument against reaching for complete.cases()
# as a first move.
2.4

2.4 How the columns fail together

Every combination of missingness across Ozone and Solar.ROzoneSolar.RDaysokok11172.5%NAok3522.9%okNA53.3%NANA21.3%Only 2 days lose both. The two columns fail almost independently, soan imputation that borrows Solar.R to predict Ozone will work on 35 of the 37 gaps.The other 2 stay missing whatever you do, unless you impute Solar.R first.
Figure 2.4. Every combination of missing and recorded across the two incomplete columns.
  1. UnitOne row of the figure is one pattern. The count is the number of days showing that pattern.
  2. EncodingRed marks a column that is missing for that pattern. Bar length repeats the count.
  3. ScaleCounts are out of 153 days. The four patterns are exhaustive and mutually exclusive.
  4. Structure111 days are complete. 35 lose Ozone alone, 5 lose Solar.R alone, and 2 lose both.
  5. Groups and exceptionsIf the two columns failed independently you would expect about 1.7 days losing both. The observed figure is 2, so the two failures are close to independent.
  6. Claim and limitSolar.R is available on 35 of the 37 days when Ozone is missing. This is what makes model-based imputation of Ozone practical. The remaining 2 days stay missing under any method that needs Solar.R.
table(Ozone   = is.na(aq$Ozone),
      Solar.R = is.na(aq$Solar.R))
#>        Solar.R
#> Ozone   FALSE TRUE
#>   FALSE   111    5
#>   TRUE     35    2

# Expected count under independence
37 * 7 / 153      #> 1.692810

# Which days lose both?
which(is.na(aq$Ozone) & is.na(aq$Solar.R))
#> [1]  5 27

# For more than two columns, count the distinct
# patterns directly.
pat <- apply(is.na(aq), 1, paste, collapse = "")
sort(table(pat), decreasing = TRUE)
2.5

2.5 Three mechanisms, three consequences

NameWhat it means Example in this datasetConsequence
MCAR Missing completely at random. Whether a value is missing has nothing to do with anything, observed or not. An instrument that fails on days chosen by a coin toss. Deleting incomplete rows is unbiased. You lose precision only.
MAR Missing at random. Whether a value is missing depends on other variables you did record. Ozone missing far more often in June than in any other month. Deletion is biased. Imputation that conditions on the observed variable can recover the structure.
MNAR Missing not at random. Whether a value is missing depends on the value itself. A monitor that saturates and records nothing on the highest-ozone days. No method in this lesson fixes it. You need the reason for the absence.
What the data can and cannot tell you You can gather evidence against MCAR from the data alone, by finding an observed variable that predicts missingness. You cannot distinguish MAR from MNAR from the data alone, because the values that would settle it are the ones you do not have. That distinction has to come from knowing how the measurements were taken.

Naming the mechanism is a claim about the world, not a calculation. Say which one you are assuming and why.

2.6

2.6 Testing for MCAR in this dataset

Days with a missing Ozone reading, by month5Mayof 31 days21Junof 30 days5Julof 31 days5Augof 31 days1Sepof 30 daysJune is missing 21 of its 30 Ozone readings. The other four months losebetween 1 and 5 days each.Whether a value is missing depends on Month, which you can observe. Thatrules out missing completely at random.
Figure 2.6. Days with a missing Ozone reading, counted by month.
  1. UnitOne bar is one month. The height is the number of days in that month with no Ozone reading.
  2. EncodingBar height gives the count. The June bar is coloured differently because it is the finding, not because of anything in the data.
  3. ScaleCounts run from 0 to 24. Month lengths differ by one day, which does not affect this comparison.
  4. StructureMay, July and August each lose 5 days. September loses 1.
  5. Groups and exceptionsJune loses 21 of its 30 days. That is four times the next highest month.
  6. Claim and limitWhether an Ozone value is missing depends on Month, which is fully observed, so the data are not missing completely at random. This rules out MCAR. It cannot separate MAR from MNAR, because a June instrument fault and a June ozone level that the instrument could not record leave identical traces.
Procedure For each observed variable, compare its distribution on the rows where your target is missing against the rows where it is present. A variable that shifts is evidence against MCAR.
# Does missingness depend on an observed variable?
tapply(is.na(aq$Ozone), aq$Month, sum)
#>  5  6  7  8  9
#>  5 21  5  5  1

table(aq$Month)
#>  5  6  7  8  9
#> 31 30 31 31 30

# Compare the other columns on missing vs
# observed days
obs <- !is.na(aq$Ozone)
round(rbind(
  observed = colMeans(aq[ obs, 2:4], na.rm = TRUE),
  missing  = colMeans(aq[!obs, 2:4], na.rm = TRUE)), 2)
#>          Solar.R  Wind  Temp
#> observed  184.80  9.86 77.87
#> missing   189.51 10.26 77.92

# Solar.R, Wind and Temp barely move.
# Month moves a great deal. That is the signal.

boxplot(Temp ~ is.na(Ozone), data = aq,
        names = c("recorded", "missing"),
        ylab = "Temperature (F)")
2.7

2.7 What cor() does with missing values

This is the question Week 7 finished on. The default is to refuse.

cor(aq)["Ozone", "Temp"]
#> [1] NA

cor(aq, use = "complete.obs")["Ozone", "Temp"]
#> [1] 0.6985

cor(aq, use = "pairwise.complete.obs")["Ozone", "Temp"]
#> [1] 0.6984

sum(complete.cases(aq))                      #> 111
sum(!is.na(aq$Ozone) & !is.na(aq$Temp))      #> 116

The three settings

use =Behaviour
"everything"The default. Any NA in a pair returns NA for that pair.
"complete.obs"Drops any row with an NA anywhere, then computes every correlation on the same 111 rows.
"pairwise.complete.obs"Uses whatever rows each pair has. Ozone against Temp gets 116 rows; Wind against Temp gets all 153.
The catch with pairwise Each cell of the matrix is computed on a different sample, so the matrix as a whole can fail to be a valid correlation matrix. For a single coefficient it is fine. For anything that consumes the whole matrix at once, use complete.obs and report the n.

Here the two settings agree to three decimal places, which happens when the extra 5 rows resemble the other 111. That agreement is a result on this dataset rather than a rule.

2.Q

2.Q Knowledge check: Section 2

Q1. In airquality, 4.79 per cent of the cells are missing but 27.45 per cent of the rows are incomplete. What explains the gap?
44 missing cells land on 42 distinct days. A row is lost as soon as one cell in it is missing, so a thin scatter of missing cells destroys a large share of rows.
Q2. Ozone is missing on 21 of 30 June days and on 1 to 5 days in every other month. What does this establish?
Missingness depends on Month, which you observed, so MCAR is ruled out. Separating MAR from MNAR needs knowledge of why the June readings were not taken, and the data cannot supply it.
Q3. x <- c(4, NA, 6). Which expression returns c(FALSE, TRUE, FALSE)?
Comparing anything with NA yields NA, so x == NA gives three NAs. is.na() is the only test that answers the question.
3
Section 3

Imputation

Filling a gap adds no information. Every method below trades one distortion for another, and the trade is your decision to defend.

3.1

3.1 Four questions before you fill anything

QUESTION 01

How much

What fraction of cells, and what fraction of rows? A column that is 60 per cent missing is not a candidate for imputation; it is a candidate for deletion.

Ozone: 24.2 per cent of cells
QUESTION 02

Where

Scattered single values, or sustained blocks? A block cannot be filled from its neighbours, because it has none nearby.

Ozone: one run of 10 days
QUESTION 03

Why

Does missingness depend on an observed variable, on the unobserved value itself, or on neither? This is the mechanism from Section 2.

Ozone: depends on Month
QUESTION 04

What for

A mean, a correlation, a chart, a model? Mean imputation preserves the mean and damages every correlation. The right method depends on which quantity you are about to report.

Choose the method last
AND THEN

Do nothing, sometimes

Reporting a statistic on the complete cases, with the n stated, is a valid answer. It is the baseline every other method has to beat.

n = 111 of 153
ALWAYS

Keep the flag

Add a logical column recording which values you filled. It costs one column and it makes every downstream result auditable.

aq$Ozone_imputed
3.2

3.2 Listwise deletion

What it does

Drops every row containing a missing value anywhere, then proceeds as though the data were complete.

cc <- aq[complete.cases(aq), ]
nrow(cc)          #> 111
nrow(aq)          #> 153

# na.omit() is the same thing
nrow(na.omit(aq)) #> 111

# Restricted to the columns you actually need
keep <- complete.cases(aq[, c("Ozone", "Temp")])
sum(keep)         #> 116
The cheapest improvement available Apply complete.cases() to the columns your analysis uses, not to the whole frame. On this dataset that recovers 5 days for an Ozone against Temp analysis, at no cost.

When it is the right choice

  • The data are missing completely at random, so the surviving rows are a fair sample.
  • The loss is small. A few per cent of rows costs precision and nothing else.
  • You want a defensible baseline that invents nothing.

What it costs here

Rows kept111 of 153
Rows lost42, or 27.45 per cent
June days kept9 of 30
The bias, stated plainly Section 2 showed that missingness depends on Month. Deleting incomplete rows therefore removes most of June. Any seasonal statistic computed on the surviving 111 days describes a year with almost no June in it.
3.3

3.3 Mean and median imputation

Observed onlyn = 116, mean 42.13, sd 32.99050100150Ozone (ppb)After mean imputationn = 153, mean 42.13, sd 28.69050100150Ozone (ppb)37 identical values land in one binThe mean is preserved exactly. The spread is not.It falls from 32.99 to 28.69, a drop of 13 per cent, and no new information entered the data.
Figure 3.3. The distribution of Ozone before and after the 37 gaps are filled with the observed mean of 42.13.
  1. UnitOne bar is a 10 ppb bin. Bar height is the number of days falling in it.
  2. EncodingHeight gives the count. The gold bar marks the bin that receives all 37 filled values.
  3. ScaleBoth panels share the same bins and the same vertical scale, so the change in shape is the change in the data.
  4. StructureThe observed distribution is right-skewed, with a long tail towards 168 ppb.
  5. Groups and exceptionsAfter imputation, 37 identical values sit in a single bin. The mean is unchanged at 42.13. The standard deviation falls from 32.99 to 28.69.
  6. Claim and limitMean imputation preserves the mean exactly and shrinks every other moment. The shrinkage follows from the arithmetic: adding values at the mean adds zero to the sum of squared deviations while increasing n. Any standard error computed afterwards is too small.
Why the shrinkage is guaranteed $$s^2 = \frac{1}{n-1}\sum_{i=1}^{n}\left(x_i - \bar{x}\right)^2$$

A value filled at \(\bar{x}\) contributes \((\bar{x}-\bar{x})^2 = 0\) to the sum and 1 to \(n\). Filling \(m\) gaps therefore multiplies \(s^2\) by \(\frac{n-1}{n+m-1}\), which is always below 1. Here \(n = 116\) and \(m = 37\), giving a factor of 0.757 on the variance and 0.870 on the standard deviation.

oz <- aq$Ozone

# mean imputation
mi <- oz
mi[is.na(mi)] <- mean(oz, na.rm = TRUE)

# median imputation, better for a skewed column
mdi <- oz
mdi[is.na(mdi)] <- median(oz, na.rm = TRUE)

c(observed = sd(oz, na.rm = TRUE),
  mean_imp = sd(mi),
  med_imp  = sd(mdi))
#> observed mean_imp  med_imp
#> 32.98788 28.69214 29.05375

mean(mi)   #> 42.12931  (unchanged)
mean(mdi)  #> 39.55556  (pulled towards the median)

# Why the sd must fall:
#   var = sum((x - xbar)^2) / (n - 1)
#   each filled value adds 0 to the numerator
#   and 1 to n, so the ratio can only shrink.
3.4

3.4 Group-wise imputation

If missingness depends on a variable you observed, condition on it. Fill each gap with the mean of its own group rather than the mean of everything.

# ave() computes a group statistic and returns it
# aligned to the original rows
grp <- ave(oz, aq$Month,
           FUN = function(v) mean(v, na.rm = TRUE))

gi <- oz
gi[is.na(gi)] <- grp[is.na(gi)]

round(tapply(oz, aq$Month, mean, na.rm = TRUE), 2)
#>     5     6     7     8     9
#> 23.62 29.44 59.12 59.96 31.45

c(mean = mean(gi), sd = sd(gi))
#>     mean       sd
#> 40.85098 29.59461

Why it beats the overall mean here

  • The monthly means run from 23.6 to 60.0. Filling a June gap with the annual mean of 42.13 overstates it by about 13 ppb.
  • The correlation with Temperature falls to 0.646 rather than 0.609, because the filled values now vary with season.
  • The spread is slightly better preserved: 29.59 against 28.69.
The limit of the idea Every filled value within a group is still identical, so within-group variance is still destroyed. Group-wise imputation moves the values to the right neighbourhood; it does not give them any spread.
Choosing the grouping Use the variable you found evidence for in Section 2. Here that was Month. Grouping on a variable unrelated to missingness adds complexity and buys nothing.
3.5

3.5 Regression imputation

Mean imputationr falls to 0.60960708090100050100150Temperature (°F)Regression imputationr rises to 0.717607080901000501001502 fills below zero(−21.6 and −2.3 ppb)Temperature (°F)Ozone (ppb)observed (116 days)filled by the meanfilled by the fitted modelMean imputation puts every filled value on one horizontal line.Regression puts every filled value on the fitted surface, and two of them below zero.
Figure 3.5. Ozone against Temperature. Blue points were recorded; coloured points are the 35 values filled by each method.
  1. UnitOne point is one day. 116 days were recorded, and up to 37 are filled.
  2. EncodingPosition gives Temperature and Ozone. Colour separates recorded values from filled ones.
  3. ScaleBoth panels use identical axes. The vertical axis extends below zero so that impossible fills stay visible.
  4. StructureThe recorded points show a positive, widening relationship. Ozone rises with Temperature and its spread rises with it.
  5. Groups and exceptionsMean imputation places all 37 fills on one horizontal line. Regression places 35 fills exactly on the fitted surface, with no scatter, and two of them below zero at −21.6 and −2.3 ppb. Ozone was never observed below 1.
  6. Claim and limitRegression imputation raises the reported correlation from 0.698 to 0.717, above the complete-case value. That rise is an artefact: filled points carry no residual, so they sit more tightly on the line than real data ever does. The model also has no knowledge that a concentration cannot be negative.
fit <- lm(Ozone ~ Solar.R + Wind + Temp, data = aq)
summary(fit)$r.squared   #> 0.6059

pr <- predict(fit, newdata = aq)
ri <- oz
k  <- is.na(ri) & !is.na(pr)
ri[k] <- pr[k]

sum(k)            #> 35   filled
sum(is.na(ri))    #> 2    still missing
sum(ri < 0, na.rm = TRUE)   #> 2  impossible

# Two problems, two fixes:

# 1. The model needs Solar.R, which is itself
#    missing on 2 of those days.
# 2. Nothing constrains the prediction to the
#    physically possible range.
ri <- pmax(ri, min(oz, na.rm = TRUE))

# Stochastic regression restores the scatter by
# adding a residual drawn from the model.
set.seed(3100)
sri <- oz
sri[k] <- pr[k] + rnorm(sum(k), 0, sigma(fit))
3.6

3.6 Hot deck imputation

Fill each gap with a value drawn from the observed values of the same column. The filled column then has the right marginal distribution.

set.seed(3100)
hi <- oz
hi[is.na(hi)] <- sample(oz[!is.na(oz)],
                        size = sum(is.na(oz)),
                        replace = TRUE)

c(observed = sd(oz, na.rm = TRUE), hotdeck = sd(hi))
#> observed  hotdeck
#> 32.98788 32.22754

cor(hi, aq$Temp)   #> 0.5683
cor(oz, aq$Temp, use = "complete.obs")   #> 0.6984

# Nearest-neighbour hot deck: draw from donors
# that resemble the recipient
donors <- aq$Temp > 80 & !is.na(aq$Ozone)
sample(aq$Ozone[donors], 1)

The trade it makes

SpreadCorrelation
Observed32.990.698
Mean28.690.609
Hot deck32.230.568

Hot deck preserves the spread better than any other method here, and damages the correlation with Temperature more than any other method here. The two facts have the same cause: a randomly drawn value has a realistic magnitude and no relationship to the day it lands on.

The repair Draw from donors that resemble the recipient on the variables you did observe. Restricting the donor pool to days with similar temperature keeps the spread and recovers much of the relationship. That is the idea behind every nearest-neighbour imputation method.

Hot deck is random, so it needs set.seed() to be reproducible. Report the seed.

3.7

3.7 What each method does to the numbers

Standard deviation of Ozoneobserved value 32.99, marked by the rule32.99Listwisen = 11128.69Meann = 15329.05Mediann = 15329.60Month meann = 15330.57Regressionn = 15132.23Hot deckn = 153Correlation between Ozone and Temperaturecomplete-case value 0.698, marked by the rule0.6980.6090.6010.6460.7170.568
Figure 3.7. Ozone standard deviation and its correlation with Temperature, under six treatments of the same 37 gaps.
  1. UnitOne bar is one treatment applied to the same column. Nothing else changes between bars.
  2. EncodingBar height gives the statistic. The dashed rule marks the complete-case value, which is the honest baseline.
  3. ScaleThe upper panel is truncated below 26 and the lower panel below 0.54, so differences are exaggerated. Read the printed numbers, not the bar heights.
  4. StructureEvery single-value method shrinks the spread. Every method except regression weakens the correlation.
  5. Groups and exceptionsHot deck keeps the spread closest to the truth (32.23 against 32.99) and damages the correlation most (0.568). Regression does the reverse and pushes the correlation above the baseline to 0.717.
  6. Claim and limitNo method here preserves both the spread and the relationship. The right choice depends on which of the two your analysis reports. Multiple imputation, which is beyond this unit, exists precisely to escape this trade.
comp <- data.frame(
  method = c("listwise", "mean", "median",
             "month", "regression", "hot deck"),
  n  = c(116, 153, 153, 153, 151, 153),
  sd = round(c(sd(oz, na.rm = TRUE), sd(mi), sd(mdi),
               sd(gi), sd(ri, na.rm = TRUE), sd(hi)), 2),
  r  = round(c(
    cor(oz,  aq$Temp, use = "complete.obs"),
    cor(mi,  aq$Temp), cor(mdi, aq$Temp),
    cor(gi,  aq$Temp),
    cor(ri,  aq$Temp, use = "complete.obs"),
    cor(hi,  aq$Temp)), 3))
comp
#>       method   n    sd     r
#> 1   listwise 116 32.99 0.698
#> 2       mean 153 28.69 0.609
#> 3     median 153 29.05 0.601
#> 4      month 153 29.60 0.646
#> 5 regression 151 30.57 0.717
#> 6   hot deck 153 32.23 0.568

# Build this table for your own data before you
# choose. It takes ten lines and it is the whole
# justification for whichever row you pick.
3.8

3.8 Flag what you filled

One extra column turns an unverifiable dataset into an auditable one.

aq2 <- aq
aq2$Ozone_imputed <- is.na(aq2$Ozone)
aq2$Ozone <- mi

sum(aq2$Ozone_imputed)   #> 37

# Now every downstream result can be checked
# with and without the filled values
cor(aq2$Ozone[!aq2$Ozone_imputed],
    aq2$Temp[!aq2$Ozone_imputed])   #> 0.6984
cor(aq2$Ozone, aq2$Temp)            #> 0.6087

# And every chart can show them
plot(aq2$Temp, aq2$Ozone,
     pch = ifelse(aq2$Ozone_imputed, 1, 19),
     col = ifelse(aq2$Ozone_imputed,
                  "#CC0000", "#0072B2"))
legend("topleft", bty = "n",
       pch = c(19, 1), col = c("#0072B2", "#CC0000"),
       legend = c("recorded", "imputed"))

What to write in your report

How many37 of 153 Ozone values, 24.2 per cent
Which methodMonthly mean
Why that oneMissingness depends on Month, and the monthly means differ by a factor of 2.5
What it costStandard deviation falls from 32.99 to 29.60
What is still missingNone for Ozone. 7 Solar.R values were left as NA.
Seed, if randomset.seed(3100)
The test of a good report A reader who disagrees with your choice should be able to redo the analysis their way from what you wrote. If your imputation is invisible in the output, they cannot.
3.Q

3.Q Knowledge check: Section 3

Q1. After mean imputation the mean of Ozone is unchanged and the standard deviation falls from 32.99 to 28.69. Why must the standard deviation fall?
Variance is the sum of squared deviations divided by n minus 1. A value equal to the mean contributes zero to the numerator and one to the denominator, so the ratio can only shrink. The effect is arithmetic, not empirical.
Q2. Regression imputation raises the Ozone-Temperature correlation from 0.698 to 0.717, above the complete-case value. What does that tell you?
Every predicted value sits exactly on the fitted surface. Adding 35 points with zero scatter to 116 points with real scatter must strengthen the apparent relationship. Adding a random residual to each prediction removes the artefact.
Q3. Which method here preserves the standard deviation best, and what does it cost?
Hot deck draws real values, so the marginal distribution survives at sd 32.23. Each drawn value has no connection to the day it lands on, so the relationship with Temperature is damaged more than by any other method here.
4
Section 4

Interpolation for time series

When rows have an order, the neighbours of a gap are informative. That is the one situation where filling a value can be defended on evidence rather than on convenience.

4.1

4.1 Why order changes the problem

MayJunJulAugSep050100150Ozone (ppb)37 missing days in 17 separate runs. Shaded bands mark the gaps.Run lengths: 1 1 3 6 1 2 2 10 1 1 1 2 2 1 1 1 1The longest run is 10 consecutive days starting 21 June. A method that interpolatesa single missing day is doing something very different from one that spans ten.
Figure 4.1. Ozone across the 153 days of the record. Shaded bands mark the 37 missing days.
  1. UnitOne point is one day. The line joins consecutive recorded days and breaks across every gap.
  2. EncodingHorizontal position gives the day index, vertical position gives Ozone. Shading marks days with no reading.
  3. ScaleThe horizontal axis is calendar time from 1 May to 30 September. Ozone runs from 1 to 168 ppb.
  4. StructureOzone is low in May, rises through July and August, and falls again in September. Day-to-day variation is large relative to that seasonal signal.
  5. Groups and exceptionsThe 37 missing days fall in 17 separate runs, of lengths 1 1 3 6 1 2 2 10 1 1 1 2 2 1 1 1 1. The longest is 10 consecutive days from 21 June.
  6. Claim and limitFilling a single missing day between two recorded neighbours is a very different act from filling ten. The same function call does both, and it reports no difference in confidence between them.
plot(aq$Ozone, type = "o", pch = 19, cex = 0.5,
     col = "#0072B2", xlab = "Day index",
     ylab = "Ozone (ppb)")

# where the gaps are
gap <- which(is.na(aq$Ozone))
abline(v = gap, col = "#CC000022", lwd = 2)

# how long each run of missing days is
r <- rle(is.na(aq$Ozone))
r$lengths[r$values]
#>  [1]  1  1  3  6  1  2  2 10  1  1  1  2  2  1  1  1  1
max(r$lengths[r$values])   #> 10

# where the longest run starts
start <- cumsum(c(1, head(r$lengths, -1)))
start[r$values & r$lengths == 10]   #> 52
aq[52, c("Month", "Day")]           #> 6, 21

# Interpolation is only available because the rows
# have an order. Shuffle them and every method in
# this section becomes meaningless.
4.2

4.2 Carrying the last value forward

LOCF holds the most recent recorded value until a new one arrives. NOCB does the same backwards.

locf <- function(v) {
  for (i in seq_along(v)[-1])
    if (is.na(v[i])) v[i] <- v[i - 1]
  v
}
nocb <- function(v) {
  for (i in rev(seq_along(v))[-1])
    if (is.na(v[i])) v[i] <- v[i + 1]
  v
}

x <- c(NA, NA, 30, 34, NA, NA, 41, 38, NA)
locf(x)        #> NA NA 30 34 34 34 41 38 38
nocb(x)        #> 30 30 30 34 41 41 41 38 NA
nocb(locf(x))  #> 30 30 30 34 34 34 41 38 38

Where it belongs

  • Step-like series that genuinely hold a value between changes: a price, a stock level, a policy rate, a sensor status.
  • Series where the last known value is the operationally correct answer, because that is what a system downstream would have used.

Where it fails

  • Smoothly varying quantities. Holding a temperature flat for ten days produces a plateau that no thermometer ever recorded.
  • A leading NA. LOCF has nothing to carry forward, so the first values stay missing until you also apply NOCB.
  • Long runs. The assumption gets weaker with every day it is extended, and the output gives no sign of it.
Reporting LOCF creates repeated identical values. Any subsequent measure of volatility, day-to-day change or autocorrelation will be biased towards stability.
oz <- aq$Ozone
f  <- locf(oz)

sum(is.na(f))   #> 0
# Ozone starts and ends with a recorded value,
# so LOCF alone leaves nothing missing here.

aq$Ozone[1]    #> 41
aq$Ozone[153]  #> 20

# The ten-day run from index 52
round(f[50:63])
# ten identical values in the middle of a series
# whose day-to-day variation is large

What that plateau costs

The controlled test on slide 4.6 measures it. On a smooth seasonal series LOCF reaches a mean absolute error of 3.89 degrees, against 3.43 for linear interpolation. It is the second best of the four methods tested, and it is worse than linear on exactly the long runs where the difference matters most.

A useful hybrid Use LOCF for runs of one or two, and leave longer runs as NA. A short hold is a mild assumption; a ten-day hold is a fabrication. Two lines of code separate them:
r <- rle(is.na(oz))
ok <- rep(r$lengths <= 2, r$lengths)
f2 <- oz
f2[ok] <- locf(oz)[ok]
4.3

4.3 Linear interpolation

46505458626670010020010 missing dayshighest recorded value in this window: 135Day indexOzone (ppb)LOCFlinearsplineobservedLOCF holds the last observed value flat across the whole run. Linear draws one straightsegment between the endpoints. The spline overshoots.The spline reaches 261 ppb inside the gap.Linear interpolation cannot leave the observed range, because every filled value is a weighted average of two real ones.
Figure 4.3. Days 46 to 70 of the Ozone record, with the three methods drawn across the same gaps.
  1. UnitOne circle is a recorded day. Each line is one method's reconstruction of the whole window.
  2. EncodingLine style separates the methods. All three are drawn on the same axes, so vertical distance between them is a real disagreement.
  3. ScaleOzone in ppb on the vertical axis; day index on the horizontal. The shaded band marks the ten-day run beginning at day 52.
  4. StructureOutside the gaps all three methods pass exactly through the recorded points. They are identical wherever data exists.
  5. Groups and exceptionsInside the ten-day band they separate widely. LOCF holds flat at 12 ppb, linear rises in one straight segment, and the spline climbs to 261 ppb. The highest value ever recorded in this window is 135.
  6. Claim and limitThe methods disagree most in the middle of the longest run, and only one of them stays inside the observed range. The disagreement measures how much of the reconstruction is invented rather than observed, and the run is exactly where no evidence exists to settle it.
idx <- seq_along(oz)
k   <- !is.na(oz)

# straight line between neighbouring points
lin <- approx(idx[k], oz[k], xout = idx, rule = 2)$y

sum(is.na(lin))          #> 0
range(lin)               #> 1 168
range(oz, na.rm = TRUE)  #> 1 168

# Linear interpolation can never leave the range
# of the observed data, because every filled value
# is a weighted average of two real ones.

# approx() also does the reverse: resampling a
# series onto a new grid
approx(idx[k], oz[k], xout = seq(1, 153, by = 7))

# method = "constant" turns approx() into LOCF
approx(idx[k], oz[k], xout = idx,
       method = "constant", rule = 2)$y
4.4

4.4 Spline interpolation and what it can do

-1000100200observed maximum 168Day indexOzone (ppb)The spline runs from −102.4 to 261.1. The observed data run from 1 to 168.Seven interpolated values are negative concentrations.A smooth curve is not the same as a plausible one.
Figure 4.4. A natural cubic spline fitted through every recorded Ozone value and evaluated on all 153 days.
  1. UnitBlue points are the 116 recorded days. The purple curve is the spline evaluated everywhere.
  2. EncodingPosition only. The shaded region below zero marks impossible concentrations.
  3. ScaleThe vertical axis has been extended to −120 and 280 ppb so the whole curve fits. The observed data occupy 1 to 168.
  4. StructureThrough dense stretches the spline tracks the data closely, because it is constrained to pass through every point.
  5. Groups and exceptionsAcross the wide gaps it overshoots badly. The fitted values run from −102.4 to 261.1, and seven of them are negative.
  6. Claim and limitA spline guarantees smoothness, not plausibility. Cubic polynomials joined across a long gap can swing far outside the data. Any interpolation of a bounded quantity needs the bound applied afterwards, and a method that needs that much repair is usually the wrong method.
spl <- spline(idx[k], oz[k], xout = idx)$y

range(spl)              #> -102.4  261.1
range(oz, na.rm = TRUE) #>    1     168
sum(spl < 0)            #> 7

# Enforcing the physical bound afterwards
spl2 <- pmax(spl, min(oz, na.rm = TRUE))

# A monotone spline cannot overshoot between
# points, which removes most of the damage
spl3 <- spline(idx[k], oz[k], xout = idx,
               method = "hyman")$y
range(spl3)

# splinefun() returns a function you can evaluate
# anywhere, rather than a fixed set of values
f <- splinefun(idx[k], oz[k])
f(52:61)

# Use a spline when the underlying quantity really
# is smooth and the gaps are short. Ozone over ten
# days is neither.
4.5

4.5 The edges of the series

rule = 1 (the default)outside the observed range, approx() gives NA2581114NANANANANANArule = 2the nearest observed value is carried outward2581114held flatInterpolation fills between observations. Filling beyond the first and last observation is extrapolation,and rule = 2 does it by assuming the series stopped changing. Choose it on purpose or leave the NAs in place.
Figure 4.5. Illustrative series with missing values before the first and after the last recorded point.
  1. UnitOne circle is a recorded observation. The dashed segments are the fitted values outside the observed range.
  2. EncodingPosition only. The vertical dashed rules mark the first and last recorded points.
  3. ScaleBoth panels use identical axes, so the only difference between them is the rule argument.
  4. StructureBetween the first and last recorded points the two panels are identical. That region is interpolation.
  5. Groups and exceptionsOutside that region, rule = 1 returns NA and rule = 2 carries the nearest recorded value outward as a flat line.
  6. Claim and limitFilling beyond the first and last observation is extrapolation, and rule = 2 does it by assuming the series stopped changing. On airquality both settings give the same answer, because Ozone begins at 41 and ends at 20, with no missing values at either end. Check that before you rely on it.
x <- c(NA, NA, 30, 34, NA, NA, 41, 38, NA)
i <- seq_along(x); k <- !is.na(x)

approx(i[k], x[k], xout = i, rule = 1)$y
#> NA NA 30.0 34.0 36.3 38.7 41.0 38.0 NA

approx(i[k], x[k], xout = i, rule = 2)$y
#> 30.0 30.0 30.0 34.0 36.3 38.7 41.0 38.0 38.0

# On airquality the distinction does not bite:
is.na(aq$Ozone[1])     #> FALSE
is.na(aq$Ozone[153])   #> FALSE
sum(is.na(approx(idx[k], oz[k],
                 xout = idx, rule = 1)$y))   #> 0

# rule can take two values, one per end:
#   rule = c(1, 2)  -> NA before, carried after

# Default to rule = 1 and decide about the ends
# on purpose. An NA you can see is safer than a
# fabricated value you cannot.
4.6

4.6 Test the method before you trust it

60708090100Day indexTemperature (°F)18 known temperatures hidden, then reconstructedMean absolute erroron the 18 hidden days, in °Fmean8.81RMSE 10.92LOCF3.89RMSE 4.91linear3.43RMSE 4.22spline5.53RMSE 7.18truthLOCFlinearsplineLinear interpolation wins on average error here, and LOCF is close behind. Meansubstitution is more than twice as bad as either.This ranking holds for a smooth seasonal series. Run the same test on your own data before you trust it.
Figure 4.6. 18 known Temperature values were hidden in four blocks, then reconstructed by each method and compared against the truth.
  1. UnitOne line is one reconstruction of the full 153-day Temperature series. The grey line is the recorded truth.
  2. EncodingLine style separates the methods. Shaded bands mark the 18 hidden days. The bar panel gives mean absolute error over those 18 days only.
  3. ScaleTemperature in degrees Fahrenheit. Errors are in the same units, so 3.43 means an average miss of about three and a half degrees.
  4. StructureAll four methods track the truth exactly outside the hidden blocks, because there they are using the real values.
  5. Groups and exceptionsLinear interpolation is best at 3.43, LOCF close behind at 3.89, spline at 5.53 and mean substitution worst at 8.81. The eight-day block accounts for most of the error in every method.
  6. Claim and limitOn this series, linear interpolation halves the error of mean substitution. Temperature is smooth and seasonal, which is the case linear interpolation is built for. Run the same test on your own series before importing this ranking.
# Temperature is complete, so we can hide values
# we already know and score the reconstruction.
tp   <- aq$Temp
gaps <- c(20:24, 55:57, 88:95, 120:121)

y <- tp
y[gaps] <- NA
k <- !is.na(y)

L  <- approx(idx[k], y[k], xout = idx, rule = 2)$y
S  <- spline(idx[k], y[k], xout = idx)$y
F  <- locf(y)
M  <- y; M[is.na(M)] <- mean(y, na.rm = TRUE)

sapply(list(mean = M, locf = F, linear = L, spline = S),
       function(z) mean(abs(z[gaps] - tp[gaps])))
#>     mean     locf   linear   spline
#> 8.807000 3.889000 3.426000 5.526000

# This is the only honest way to choose. It costs
# ten lines and it turns a preference into
# evidence you can put in a report.
4.Q

4.Q Knowledge check: Section 4

Q1. Why is interpolation available for airquality but not for a dataset of customer records?
Interpolation borrows from adjacent rows. That only means something when adjacency is meaningful. Shuffle the rows of a time series and every method in this section becomes arbitrary.
Q2. A spline fitted through the recorded Ozone values produces seven negative concentrations and a maximum of 261 ppb, against an observed maximum of 168. What does that show?
A cubic polynomial forced through points on either side of a long gap can swing far outside the data. The curve is smooth everywhere and impossible in seven places.
Q3. Linear interpolation scored a mean absolute error of 3.43 degrees on the hidden Temperature values, against 8.81 for mean substitution. What may you conclude?
The test measured one method on one series with one pattern of gaps. That is exactly the evidence you need for this dataset and no evidence at all for the next one.
5
Section 5

Choosing and reporting

The method is a decision, and a decision has to be written down. What you report matters as much as what you compute.

5.1

5.1 Choosing a treatment

SituationReach for WhyAvoid
Under a few per cent of rows incomplete, no pattern found Listwise deletion Invents nothing and costs only precision. Any imputation, which adds risk for no gain
Missingness depends on a variable you recorded Group-wise mean, or regression on that variable Conditioning on the observed cause removes the bias that deletion would leave. Overall mean, which ignores the very structure you found
You will report a spread, a standard error or a confidence interval Hot deck, or stochastic regression Single-value methods shrink the spread by construction, so every interval afterwards is too narrow. Mean and median imputation
You will report a correlation or fit a model Complete cases, with the n stated Every imputation distorts the relationship in a direction you have to argue about. Regression imputation, which strengthens the relationship it assumed
Ordered series, gaps of one or two Linear interpolation The neighbours genuinely constrain the value, and the result stays inside the observed range. Spline, whose extra flexibility buys nothing over two points
Ordered series, long runs Leave the NA in place Ten days of invented values look identical to ten days of measurement in every chart you will draw. Any method that fills silently
Step-like series such as a price or a stock level LOCF The value genuinely held until it changed. Linear interpolation, which invents a ramp
5.2

5.2 The four errors to avoid

The errorWhy it happens What to do instead
Reading a sentinel as data A file marks absence with −999, 0, or an empty cell. R reads it as a number and your mean drops. Declare them at import: na.strings = c("", "NA", "-999"), then check range() on every numeric column.
Deleting first and looking later na.omit() is the fastest thing to type, so it gets typed before anyone has counted what it removes. Count the rows it would drop and check whether they differ from the rest, before running it.
Imputing, then reporting a confidence interval Single-value imputation shrinks the spread, and the software has no way to know the values were invented. Report the interval from complete cases, or use a method that restores the variability you removed.
Filling a long run silently One function call fills a one-day gap and a ten-day gap, and the output looks the same either way. Check the run lengths with rle() first. Set a maximum gap you are willing to fill and leave the rest as NA.
The common thread Every one of these is a decision taken by default rather than on purpose. Looking at the missingness before acting on it costs a few minutes and removes all four.
5.3

5.3 Summary

Data frames

  • Rows are observations, columns are variables, and each column holds one type.
  • Five operations cover almost everything: select, filter, arrange, mutate, summarise. All five exist in base R.
  • Long and wide are both tidy. Which one depends on what you are calling an observation.

Missing values

  • NA means unknown. It propagates, and it cannot be compared with ==.
  • Count cells and rows separately. Here 4.79 per cent of cells cost 27.45 per cent of rows.
  • Look at the missingness before treating it: a matrix, a bar chart, and a pattern table.
  • Evidence against MCAR comes from an observed variable that predicts missingness. Here that variable is Month.

Imputation

  • Filling adds no information. Every method trades one distortion for another.
  • Mean imputation preserves the mean and shrinks everything else, by arithmetic necessity.
  • Regression imputation strengthens the relationship it assumed, and can produce impossible values.
  • Hot deck preserves the spread and damages the relationship.
  • Flag every value you filled.

Interpolation

  • Available only when the rows have an order.
  • Linear stays inside the observed range. Splines do not.
  • Run length decides the method. One missing day and ten missing days are different problems.
  • Hide values you already know, reconstruct them, and score the result. That turns a preference into evidence.
5.4

5.4 Before next week

Practise

  • Reproduce Figure 2.2 and Figure 2.3 from the code panels. Then run the same two on attitude and on a dataset of your own.
  • Build the comparison table from slide 3.7 for Solar.R rather than Ozone. Its 7 gaps behave differently from Ozone's 37.
  • Repeat the hide-and-score test from slide 4.6 on Wind instead of Temp. Wind is noisier, and the ranking of the four methods may not survive.
  • Take one chart from Week 7 and redraw it with imputed points marked. Decide whether the finding survives.

Week 9: Aggregates and joining tables

The aggregate() calls in this lesson were the summarise step used in passing. Next week they become the topic, along with what happens when two tables have to be matched on a key.

The link back A join is where missing values are created rather than found. A key present in one table and absent from the other produces a row of NAs that was never in either source. Everything in Section 2 applies to those NAs as well.
Bring Your comparison table for Solar.R, and one sentence naming which method you would use and what it costs.
TECH3100 · Lesson 8
← → navigate · T contents

Contents

Press T or Escape to close