# A tibble: 5 × 4
id date loc temp
<dbl> <dttm> <chr> <dbl>
1 1 2010-01-03 00:00:00 New York 42
2 2 2010-02-03 00:00:00 New York 41.4
3 3 2010-03-03 00:00:00 New York 38.5
4 4 2010-04-03 00:00:00 New York 41.1
5 5 2010-05-03 00:00:00 New York 39.8ETC5521: Diving Deeply into Data Exploration
Initial data analysis and model diagnostics
IDA in government statistics
The purpose of data cleaning is to bring data up to a level of quality such that it can reliably be used for the production of statistical models or statements.
A statistical value chain is constructed by defining a number of meaningful intermediate data products, for which a chosen set of quality attributes are well described.
van der Loo & de Jonge (2018) Statistical Data Cleaning with Applications in R
IDA in health and medical data
Huebner et al (2018)’s six steps of IDA: (1) Metadata setup, (2) Data cleaning, (3) Data screening, (4) Initial reporting, (5) Refining and updating the analysis plan, (6) Reporting IDA in documentation.
Example: Checking the data type (1/2)
lecture3-example.xlsx
- What problems are there with the computer’s interpretation of data type?
- What context specific issues indicate incorrect computer interpretation?
Example: Checking the data type visually
Case study: World development indicators (4/7)
Show code
wdi <- raw_dat |>
select(`Country Code`, `Series Code`, `1969 [YR1969]`:`2018 [YR2018]`) |>
rename_all(janitor::make_clean_names) |>
pivot_longer(x1969_yr1969:x2018_yr2018,
names_to = "year",
values_to = "value") |>
mutate(year = as.numeric(str_sub(year, 2, 5)) ) |>
pivot_wider(names_from = series_code,
values_from = value)
wdi2017 <- wdi |> filter(year == 2017)- Organise data into tidy form
- Check missing value distribution
Case study: World development indicators (5/7)
Check missings by
- variable
- country
Show code
wdi_cnt_miss <- wdi |>
filter(!is.na(country_code)) |>
bind_shadow() |>
select(country_code, year,
SP.ADO.TFRT_NA:SP.URB.GROW_NA) |>
pivot_longer(SP.ADO.TFRT_NA:SP.URB.GROW_NA,
names_to="var",
values_to="value") |>
group_by(country_code) |>
count(value) |>
mutate(value = fct_recode(value,
miss="NA",
not="!NA")) |>
pivot_wider(names_from = value, values_from = n) |>
mutate(p_miss = miss/(miss+not)) |>
select(country_code, p_miss)
wdi_cnt_p <- wdi_cnt_miss |>
ggplot(aes(x=1, y=p_miss,
label=country_code)) +
geom_quasirandom() +
ylim(c(0,1)) + ylab("Prop miss")
ggplotly(wdi_cnt_p)Case study: World development indicators (6/7)
Look at Costa Rica (CRI), most complete country
Case study: World development indicators (6/7)
Missings imputed using imputeTS using the moving average method.
- Don’t have to impute before scrutinizing data
- What are these numbers supposed to be?
SE.PRM.CMPT.ZS is “Primary completion rate, total (% of relevant age group)”
Do we have any problems?
Yes. The explanation of the variable suggests the numbers should range between 0-100.
Your turn: check another country
Your turn
- Pick a different
country_codefromwdi(not"CRI"). - Produce a
vis_miss()plot for that country, sorted by missingness. - Is it more or less complete than Costa Rica? Would you be comfortable imputing missing values for it?
Case study: Employment Data in Australia (2/3)
What do you notice?
Show code
employed_2000 <- employed |>
filter(between(year, 2005, 2014)) |>
mutate(year = factor(year)) |>
ungroup()
labels_2000 <- employed_2000 |>
group_by(year) |>
slice_tail(n = 1) |>
ungroup()
ggplot() +
geom_line(data = employed_2000,
aes(x = month, y = value/1000, color = year)) +
geom_label_repel(
data = labels_2000,
aes(x = month, y = value/1000,
label = year, color = year),
fill = alpha("white", 0.6),
direction = "y",
hjust = 0,
nudge_x = 0.5,
segment.color = "grey70"
) +
ylab("employed (mil)") +
scale_x_continuous("month", breaks=seq(1, 12, 1)) +
scale_color_viridis_d("", end = 0.8) +
theme(aspect.ratio = 1.5,
legend.position = "none")
Why do you think the number of people employed is going up each year?
- Australian population is 25.39 million in 2019
- 1.5% annual increase in population
- Vic population is 6.681 million (Sep 2020) - 26%
- NSW population is 8.166 (Sep 2020) - 32%
Case study: Employment Data in Australia (3/3)
- Also notice a change in August numbers from 2014.
- A potential explanation for this is that there was a change in the survey from 2014.
See discussion on this at Hyndsight blog (10 October 2014).
Your turn: check other month pairs
How do we provide evidence that the August-to-September is a genuine problem, that is worth investigating to understand why.
Your turn
- Repeat the same check for each consecutive different pair of adjacent months (e.g. June to July).
- Make the boxplots for each set of differences. You can plot them against month, or use facetting.
- Does the same kind of step-change show up?
- What does this tell you about whether the August/September pattern is a specific, real event, versus something you’d see between any two months?
Show code
employed_diffs <- employed_2000 |>
arrange(date) |>
mutate(diff = value - dplyr::lag(value)) |>
mutate(month_c = factor(month, levels = c(1:12),
labels = c("jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec"))) |>
group_by(month_c) |>
mutate(diff_m = diff -
median(diff, na.rm = TRUE)) |>
ungroup()
Case study: 2014 Data Mining Cup winners
Ugly plot of all observations provided in training sample, with response variable in colour, and test sample to predict.
What does this tell you about the test sample?
Case study: Warranty claims
Show code
Rows: 4,561
Columns: 14
$ Region <chr> "East", "West", "North …
$ State <chr> "Delhi", "Gujarat", "We…
$ Area <chr> "Urban", "Rural", "Urba…
$ City <chr> "New Delhi", "Ahmedabad…
$ Consumer_profile <chr> "Personal", "Personal",…
$ TV_2001_Issue <dbl> 1, 1, 0, 0, 0, 0, 1, 1,…
$ TV_2002_Issue <dbl> 1, 1, 1, 0, 0, 0, 1, 1,…
$ TV_2003_Issue <dbl> 1, 0, 1, 0, 0, 0, 1, 0,…
$ Claim_Value <dbl> 25000, 4216, 4000, 5000…
$ Service_Centre <dbl> 13, 10, 10, 12, 10, 10,…
$ Product_Age <dbl> 60, 672, 275, 10, 4, 34…
$ Purchased_from <chr> "Dealer", "Dealer", "De…
$ Call_details <dbl> 1.3, 25.0, 11.0, 1.6, 0…
$ Purpose <chr> "Complaint", "Other", "…TV_2001_Issue: failure of power supplyTV_2002_Issue: failure of inverterTV_2003_Issue: failure of motherboard
- What is the population that this data is measuring?
- What is not measured?
Example 2: El Nino
Explore missings
- plotting on edge of plots, or
- using simple imputation like mean
Impute and check
- Impute using regression or simulation
- Check distribution relative to complete cases
Show code
library(simputation)
ocean_imp_yr <- oceanbuoys %>%
bind_shadow() %>%
impute_lm(air_temp_c ~ wind_ew + wind_ns + year + longitude + latitude) %>%
impute_lm(humidity ~ wind_ew + wind_ns + year + longitude + latitude) %>%
impute_lm(sea_temp_c ~ wind_ew + wind_ns + year + longitude + latitude) %>%
add_label_shadow()
ggplot(ocean_imp_yr,
aes(x = air_temp_c,
y = humidity,
color = any_missing)) +
geom_point() +
theme(legend.title = element_blank())
Example 1: Checking variance and distribution (2/3)
Cushny, A. R. and Peebles, A. R. (1905) The action of optical isomers: II hyoscines. The Journal of Physiology 32, 501–510.
Few observations. Nothing strongly suggests violation of normality and spread of points is similar for each group.
Paired t-test
data: extra[group == 1] and extra[group == 2]
t = -4, df = 9, p-value = 0.003
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
-2.5 -0.7
sample estimates:
mean difference
-1.6 Example 2: Checking distribution and variance (3/3)
Is it plausible that the samples are from a normal population? Do they have equal variance?
Df Sum Sq Mean Sq F value Pr(>F)
spray 5 2669 534 34.7 <2e-16 ***
Residuals 66 1015 15
---
Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1What hypothesis being tested? What would the decision be?
Why does equal variance matter in this test?