ETC5521 Worksheet Week 8

Going beyond two variables, exploring high dimensions

Author

Prof. Di Cook

Exercise 1: Risk taking on vacation

The book “Market Segmentation Analysis - Understanding It, Doing It, and Making It Useful” by Sara Dolnicar, Bettina Grün and Friedrich Leisch has a selection of useful data sets. These can be used by installing the package:

install.packages("https://statistik.boku.ac.at/nachlass_leisch/MSA/packages/MSA_0.3-1.tar.gz", repos = NULL, type = "source")

If you can’t install the package, you can download the data from the book website.

Load the risk data. This data has 563 observations and 6 variables.

The data was collected by academic researchers using a permission based online panel.

The sample was taken from adult Australian residents who have undertaken at least one holiday in the last year which involved staying away from home for at least four nights.

The respondents were asked: “Which risks have you taken in the past?” and answered on a 5-point scale with options:

Never (1)
Rarely (2)
Quite often (3)
Often (4)
Very often (5)

The six types of risk, which form the variables were:

Recreational risks (e.g., rock-climbing, scuba diving)
Health risks (e.g., smoking, poor diet, high alcohol consumption)
Career risks (e.g., quitting a job without another to go to)
Financial risks (e.g., gambling, risky investments)
Safety risks (e.g., speeding)
Social risks (e.g., standing for election, publicly challenging a rule or decision)

Use AI to find an explanation of market segmentation. Write a few sentences to explain it.

Solution

Market segmentation divides customers into groups with similar characteristics so a business can serve each group more effectively. It is sometimes achieved by using cluster analysis to discover groups in the data.

Then so what is cluster analysis?

What type of plots should you make for this data?

Solution

Because this is categorical (ordinal) data a good choice would be to make Likert plots.

Code

data(risk)
risk_f <- as_tibble(risk) |>
  mutate_if(is.numeric, as.factor)
risk_likert <- gglikert(
  risk_f,
  add_labels = FALSE,
  add_totals = FALSE,
  sort = "ascending"
) +
  scale_fill_brewer(palette = "RdBu", direction = -1) 
risk_likert

What this doesn’t consider is whether some customers respond similarly to other customers. How would you check this?

Solution

This data can be considered to be numeric, and so we would compute the correlation.

Code

cor(risk)

             Recreational Health Career Financial Safety Social
Recreational         1.00   0.27   0.40      0.29   0.37   0.44
Health               0.27   1.00   0.28      0.37   0.40   0.33
Career               0.40   0.28   1.00      0.40   0.32   0.44
Financial            0.29   0.37   0.40      1.00   0.42   0.35
Safety               0.37   0.40   0.32      0.42   1.00   0.32
Social               0.44   0.33   0.44      0.35   0.32   1.00

Another way to look at this is using a tour.

Code

animate_xy(risk)

There is moderate association between the responses. It is strictly linear, also.

A common approach to do market segmentation is to conduct a cluster analysis (outside the scope of this class), which will divide the observations into groups. The question becomes “how many groups would be suitable to summarise the customer behaviour?” We are going to use the guided tour to help make the decision.

This is code that will do the clustering:

Code

risk_d <- apply(risk, 2, function(x) (x - mean(x)) / sd(x))

# Clustering
nc <- 5 # Set the number of clusters
set.seed(1145)
r_km <- kmeans(risk_d, centers = nc, iter.max = 500, nstart = 5)

r_km_d <- risk_d |>
  as_tibble() |>
  mutate(cl = factor(r_km$cluster)) |>
  bind_cols(model.matrix(~ as.factor(r_km$cluster) - 1))
colnames(r_km_d)[(ncol(r_km_d) - nc + 1):ncol(r_km_d)] <- paste0(
  "cluster",
  1:nc
)
r_km_d <- r_km_d |>
  mutate_at(vars(contains("cluster")), function(x) x + 1)

Change the number of clusters nc starting with 2 up to maybe 5. Examine how it is breaking the data, and what the groups would mean using the 6 types of risk.

Code

animate_xy(r_km_d[, 1:6], guided_tour(lda_pp(r_km_d$cl)), col = r_km_d$cl)

Solution

Two and three clusters breaks the data along the main direction of association, into low and high risk taker groups. Probably we’d choose the three groups here to have low, medium and high risk.

Four clusters breaks the moderate group into two, which balances recreational and social (possibly also career) against safety. This could be a useful grouping.

Five clusters breaks out this group further into health against social and recerational. Possibly this is also useful.

I’d tend to stick wih the three groups, as the clearest distinction between customers, but I’m not a tourism operator.

Exercise 2: Parkinsons

This dataset is composed of a range of biomedical voice measurements from 31 people, 23 with Parkinson’s disease (PD). Each column in the table is a particular voice measure, and each row corresponds one of 195 voice recording from these individuals (“name” column). The main aim of the data is to discriminate healthy people from those with PD, according to “status” column which is set to 0 for healthy and 1 for PD.

The data is available at The UCI Machine Learning Repository in ASCII CSV format. The rows of the CSV file contain an instance corresponding to one voice recording. There are around six recordings per patient, the name of the patient is identified in the first column. There are 24 variables in the file, including the persons name in column 1.

The data are originally analysed in: Max A. Little, Patrick E. McSharry, Eric J. Hunter, Lorraine O. Ramig (2008), ‘Suitability of dysphonia measurements for telemonitoring of Parkinson’s disease’, IEEE Transactions on Biomedical Engineering (to appear).

Code

library(cassowaryr)
# Load the data
data(pk)

How many pairwise plots would you need to look at, to look at all of them?

Solution

There are 23 numeric variables in the data set, which would require 253 pairwise plots to be made.

Compute several of the scagnostics (monotonic, outlying, clumpy2) for the first five variables of variables, except for name. (Note: We are using just five for computing speed, but the scagnostics could be calculated on all variables.)

Code

# Compute the scagnostics on the relevant variables
s <- calc_scags_wide(pk[,2:5],
                scags=c("outlying","monotonic",
                        "clumpy2"))
s

Solution

# A tibble: 6 × 5
  Var1           Var2         outlying clumpy2 monotonic
  <fct>          <fct>           <dbl>   <dbl>     <dbl>
1 MDVP:Fhi(Hz)   MDVP:Fo(Hz)     0.412   0.756    0.796 
2 MDVP:Flo(Hz)   MDVP:Fo(Hz)     0.172   0.541    0.324 
3 MDVP:Flo(Hz)   MDVP:Fhi(Hz)    0.336   0.671    0.0956
4 MDVP:Jitter(%) MDVP:Fo(Hz)     0.289   0        0.270 
5 MDVP:Jitter(%) MDVP:Fhi(Hz)    0.541   0.764    0.0978
6 MDVP:Jitter(%) MDVP:Flo(Hz)    0.430   0        0.407

Sort the scagnostics, separately by the values on (i) monotonic (ii) outlying (iii) clumpy2, and plot the pair of variables with the highest values on each.

Solution

Code

# Check the results for monotonic
s |> 
  select(Var1, Var2, monotonic) |> 
  arrange(desc(monotonic))

# A tibble: 6 × 3
  Var1           Var2         monotonic
  <fct>          <fct>            <dbl>
1 MDVP:Fhi(Hz)   MDVP:Fo(Hz)     0.796 
2 MDVP:Jitter(%) MDVP:Flo(Hz)    0.407 
3 MDVP:Flo(Hz)   MDVP:Fo(Hz)     0.324 
4 MDVP:Jitter(%) MDVP:Fo(Hz)     0.270 
5 MDVP:Jitter(%) MDVP:Fhi(Hz)    0.0978
6 MDVP:Flo(Hz)   MDVP:Fhi(Hz)    0.0956

Code

ggplot(data=pk, 
       aes(x=`MDVP:Fhi(Hz)`, y=`MDVP:Fo(Hz)`)) + 
  geom_point() + theme(aspect.ratio = 1)

The top pair of variables on monotonic has a strong positive association with the majority of points, and a few outliers. The top pair of variables on outlying, is also the top pair on clumpy2, and has outliers with some clumpiness in the mass of points with low values.

Make an interactive scatterplot matrix. Browse over it to choose other interesting pairs of variables and make the plots.

Solution

Code

# Create an interactive splom
s <- s |>
  mutate(vars = paste(Var1, Var2))
highlight_key(s) |>
  GGally::ggpairs(columns = 3:5, mapping = aes(label = vars)) |>
  ggplotly(tooltip = "all", 
           width=500, 
           height=500) |>
  highlight("plotly_selected")

The scagnostics help us to find interesting associations between pairs of variables. However, the problem here is to detect differences between Parkinsons’ patients and normal patients. How would you go about that? Think about some ideas long the line of scagnostics but look for differences between the two groups.

Solution

Code

# One way to examine difference between Parkinsons and healthy
pk_med <- pk |> 
  select(-name) |> 
  group_by(status) |>
  summarise_all(list(median, sd)) |>
  pivot_longer(
    cols=`MDVP:Fo(Hz)_fn1`:`PPE_fn2`,
    names_to="var", 
    values_to="value") |>
  separate(var, c("var","stat"), "_") |>
  mutate(stat = fct_recode(stat,
                           "m"="fn1",
                           "s"="fn2")) |>
  pivot_wider(names_from=stat,
              values_from=value) |>
  group_by(var) |>
  summarise(
    d = (m[status==0]-m[status==1])/sqrt(s[status==0]^2+s[status==1]^2))
pk_med |> arrange(desc(d)) |> head()

# A tibble: 6 × 2
  var               d
  <chr>         <dbl>
1 MDVP:Fo(Hz)   0.870
2 HNR           0.647
3 MDVP:Fhi(Hz)  0.518
4 MDVP:Flo(Hz)  0.211
5 NHR          -0.243
6 MDVP:RAP     -0.356

Code

ggplot(pk, aes(x=factor(status), y=`MDVP:Fo(Hz)`)) + 
  geom_boxplot()

Generally we are looking for variables where the differences between the Parkinsons and normal patients are big. You need to measure big, relative to the variance of each group. Doing a two sample t-test for each variable is one approach. Here, I’ve computed the median for each group of patients and compared the difference in medians relative to the pooled standard deviation in each group.

Now try to do this using the scagnostics.

Solution

Code

# Check the pair of variables already examined
ggplot(data=pk, 
       aes(x=`MDVP:Fhi(Hz)`, 
           y=`MDVP:Fo(Hz)`, 
           colour=factor(status))) + 
  geom_point(alpha=0.5) + 
  scale_colour_discrete_divergingx(palette="Zissou 1") +
  theme(aspect.ratio = 1,
        legend.title = element_blank())

Code

# Compute scagnostics separately for each group and compare differences
pk_0 <- pk |>
  filter(status == 0)
pk_1 <- pk |>
  filter(status == 1)
s_0 <- calc_scags_wide(pk_0[,2:5],
                scags=c("outlying","monotonic",
                        "clumpy2"))
s_1 <- calc_scags_wide(pk_1[,2:5],
                scags=c("outlying","monotonic",
                        "clumpy2"))
s_0 <- s_0 |>
  mutate(status = 0)
s_1 <- s_1 |>
  mutate(status = 1)
s <- bind_rows(s_0, s_1)  
s_wide <- s |>
  pivot_longer(outlying:monotonic, names_to="scag", values_to="value") |>
  pivot_wider(names_from = status, values_from = "value") |>
  mutate(d = abs(`1`-`0`)) |>
  arrange(desc(d))

Code

# Plot top pair
ggplot(data=pk, 
       aes(x=`MDVP:Jitter(%)`, 
           y=`MDVP:Fo(Hz)`, 
           colour=factor(status))) + 
  geom_point(alpha=0.5) + 
  scale_colour_discrete_divergingx(palette="Zissou 1") +
  theme(aspect.ratio = 1,
        legend.title = element_blank())