--- title: "Discriminant Analysis" editor: markdown: wrap: 72 --- ## Discriminant analysis - ANOVA and MANOVA: predict a (counted/measured) response from group membership. - Discriminant analysis: predict group membership based on counted/measured variables. - Covers same ground as logistic regression (and its variations), but emphasis on classifying observed data into correct groups. ## ... continued - Does so by searching for linear combination of original variables that best separates data into groups (canonical variables). - Assumption here that groups are known (for data we have). If trying to "best separate" data into unknown groups, see *cluster analysis*. ## Packages ```{r bDiscrim-1, message=FALSE, warning=FALSE} library(MASS, exclude = "select") library(tidyverse) library(ggrepel) library(ggord) # installation instructions below library(MVTests) # for Box M test ``` - `ggrepel` allows labelling points on a plot so they don't overwrite each other. ## About `select` - Both `dplyr` (in `tidyverse`) and `MASS` have a function called `select`, and *they do different things*. - When you load `MASS`, make sure to load it *without* its `select`, so that when you use `select`, you get the one you're used to. - If you forget, and you intend to use the `tidyverse` `select`, you will get a problem that is almost impossible to debug unless you have seen it before. ## Installing `ggord` - `ggord` (the package) contains a function, also called `ggord`, that makes a nice picture of a discriminant analysis. - It lives on `r-universe`, so to install: ```{r} #| eval: false install.packages("ggord", repos = "https://fawda123.r-universe.dev") ``` ## Example 1: seed yields and weights ```{r bDiscrim-2, message=FALSE} my_url <- "http://ritsokiguess.site/datafiles/manova1.txt" hilo <- read_delim(my_url, " ") g <- ggplot(hilo, aes(x = yield, y = weight, colour = fertilizer)) + geom_point(size = 4) ``` ::: columns ::: {.column width="40%"} Recall data from MANOVA: needed a multivariate analysis to find difference in seed yield and weight based on whether they were high or low fertilizer. ::: ::: {.column width="60%"} ```{r} #| echo = FALSE g ``` ::: ::: ## Basic discriminant analysis ```{r bDiscrim-4 } hilo.1 <- lda(fertilizer ~ yield + weight, data = hilo) ``` - Uses `lda` from package MASS. - "Predicting" group membership from measured variables. ## Output (in `hilo.1`) \small ```{r bDiscrim-5} #| echo: false hilo.1 ``` \normalsize ## Things to take from output 1/2 - Group means: high-fertilizer plants have (slightly) higher mean yield and weight than low-fertilizer plants. - "Coefficients of linear discriminants": \texttt{LD1, LD2,}\ldots are scores constructed from observed variables that best separate the groups. - For any plant, get LD1 score by taking $-0.76$ times yield plus $-1.25$ times weight, add up, standardize. ## Things to take from output 1/2 - the LD1 coefficients are like slopes: - if yield higher, LD1 score for a plant lower - if weight higher, LD1 score for a plant lower - High-fertilizer plants have higher yield and weight, thus low (negative) LD1 score. Low-fertilizer plants have low yield and weight, thus high (positive) LD1 score. - One LD1 score for each observation. Plot with actual groups, because the LD1 scale should do the best job of distinguishing the fertilizer groups. ## How many linear discriminants? - Smaller of these: - Number of variables - Number of groups *minus 1* - Seed yield and weight: 2 variables, 2 groups, $\min(2,2-1)=1$. ## Getting LD scores Feed output from LDA into `predict`: ```{r bDiscrim-6 } p <- predict(hilo.1) hilo.2 <- cbind(hilo, p) ``` ## the LD scores \scriptsize ```{r} hilo.2 ``` \normalsize ## LD1 scores in order ```{r bDiscrim-7} hilo.2 %>% select(fertilizer, yield, weight, LD1) %>% arrange(desc(LD1)) ``` ## LD1 scores and fertilizer Most positive LD1 score is most obviously low fertilizer, most negative is most obviously high. High fertilizer have yield and weight high, negative LD1 scores. ## Plotting LD1 scores With one LD score, plot against (true) groups, eg. boxplot: ```{r bDiscrim-8, fig.height=3.4} ggplot(hilo.2, aes(x = fertilizer, y = LD1)) + geom_boxplot() ``` ## What else is in `hilo.2`? - `class`: predicted fertilizer level (based on values of `yield` and `weight`). - `posterior`: predicted probability of being low or high fertilizer given `yield` and `weight`. - `LD1`: scores for (each) linear discriminant (here is only LD1) on each observation. ## Predictions and predicted groups \ldots based on `yield` and `weight`: \footnotesize ```{r bDiscrim-12} hilo.2 %>% select(yield, weight, fertilizer, class) ``` \normalsize ## Count up correct and incorrect classification ```{r} with(hilo.2, table(obs = fertilizer, pred = class)) ``` - Each predicted fertilizer level is exactly same as observed one (perfect prediction). - Table shows no errors: all values on top-left to bottom-right diagonal. ## Posterior probabilities \footnotesize show how clear-cut the classification decisions were: ```{r bDiscrim-13} hilo.2 %>% mutate(across(starts_with("posterior"), \(p) round(p, 4))) %>% select(-LD1) ``` ## Comments Only obs. 7 has any doubt: `yield` low for a high-fertilizer, but high `weight` makes up for it. ## Example 2: jobs and survey scores 244 people who do one of three different jobs also took a survey that gave them scores on three different traits called Outdoor, Social, and Conservative. Can we use these survey scores to distinguish the people who do the three different jobs? Data in . ## Read in ```{r} my_url <- "https://datafiles.ritsokiguess.site/jobs.txt" jobs0 <- read_table(my_url) jobs0 %>% slice_sample(n = 10) ``` ## Problem - The jobs are numbered 1, 2, and 3, but we actually know that 1 is customer service, 2 is mechanic, 3 is dispatcher. It would be better to have those names in the dataset. - `recode_values`: ```{r} jobs0 %>% mutate(job = recode_values( job, 1 ~ "cs", 2 ~ "mech", 3 ~ "disp" )) -> jobs ``` ## Check ```{r} jobs %>% slice_sample(n = 10) ``` ## Discriminant analysis ```{r} jobs.1 <- lda(job ~ outdoor + social + conservative, data = jobs) jobs.1 ``` ## Comments - Now have 3 variables and 3 groups, so $\min(3, 3-1) = 2$ LDs. - "Proportion of trace" says LD1 does best at separating jobs, but LD2 might do something as well. - High score on LD1: high on `social`, low on other two. - High score on LD2: *low* on `outdoor`. ## Discriminant scores - via `predict`: ```{r} p <- predict(jobs.1) jobs.2 <- cbind(jobs, p) glimpse(jobs.2) ``` ## Plotting discriminant scores - now have two of them (quantitative), plus true job (categorical): ```{r} ggplot(jobs.2, aes(x = x.LD1, y = x.LD2, colour = job)) + geom_point() ``` ## The bi-plot Individuals *and* original variables on same plot: ```{r} #| warning: false ggord(jobs.1, jobs$job, ellipse = FALSE, vec_ext = 5) ``` ## Comments - on first plot: - customer service mostly to right (high LD1) - dispatcher mostly to left (low LD1) - mechanic mostly at bottom (low LD2) - but some mixing - on bi-plot: - similar story, but also: - high LD1 goes with high `social`, low `conservative` - low LD1 goes with high `outdoor` ## Inference from bi-plot - customer service high on `social`, low on `conservative` - dispatchers opposite of customer service - mechanics high on `outdoor` ## How good was the classification? ```{r} with(jobs.2, table(obs = job, pred = class)) ``` - Most of the jobs were gotten correct - but mechanics in particular were often confused with people from other jobs ## Some posterior probabilities for people gotten wrong ```{r} jobs.2 %>% filter_out(job == class) %>% select(job, class, starts_with("posterior")) %>% slice_sample(n = 10) ``` ## Example 3: professions and leisure activities - 15 individuals from three different professions (politicians, administrators and belly dancers) each participate in four different leisure activities: reading, dancing, TV watching and skiing. After each activity they rate it on a 0--10 scale. - How can we best use the scores on the activities to predict a person's profession? - Or, what combination(s) of scores best separate data into profession groups? ## The data \scriptsize ```{r} #| message = FALSE my_url <- "http://ritsokiguess.site/datafiles/profile.txt" active <- read_delim(my_url, " ") active ``` \normalsize ## Discriminant analysis \tiny ```{r bDiscrim-30, message=F} active.1 <- lda(job ~ reading + dance + tv + ski, data = active) active.1 ``` \normalsize ## Comments - Two discriminants, first fair bit more important than second. - `LD1` depends (negatively) most on `dance`, a bit on `tv`. - `LD2` depends mostly (positively) on `tv`. ## Misclassification ```{r bDiscrim-31 } p <- predict(active.1) active.2 <- cbind(active, p) active.2 with(active.2, table(obs = job, pred = class)) ``` Everyone correctly classified. ## Plotting LDs \small ```{r bDiscrim-32, fig.height=4} g <- ggplot(active.2, aes(x = x.LD1, y = x.LD2, colour = job, label = job)) + geom_point() + geom_text_repel() + guides(colour = "none") g ``` \normalsize ## Biplot ```{r bDiscrim-33, fig.height=4.3} # ggbiplot(active.1, groups = active$job) ggord(active.1, active$job, ellipse = FALSE, vec_ext = 2) ``` ## Comments on plot - Groups well separated: bellydancers bottom left, administrators bottom right, politicians upper middle. - Bellydancers most negative on `LD1`: like dancing most. - Administrators most positive on `LD1`: like dancing least. - Politicians most negative on `LD2`: like TV-watching most. ## Posterior probabilities \scriptsize ```{r bDiscrim-35} active.2 %>% mutate(across(starts_with("posterior"), \(p) round(p, 3))) %>% rename_with(\(p) str_remove(p, "osterior"), starts_with("posterior")) %>% select(job, class, starts_with("p.")) ``` \normalsize Not much doubt. ## Example 4: remote-sensing data - View 25 crops from air, measure 4 variables `x1-x4`. - Go back and record what each crop was. - Can we use the 4 variables to distinguish crops? ## The data \footnotesize ```{r bDiscrim-39 } #| message = FALSE my_url <- "http://ritsokiguess.site/datafiles/remote-sensing.txt" crops <- read_table(my_url) crops ``` \normalsize ## Discriminant analysis \tiny ```{r bDiscrim-40 } crops.1 <- lda(crop ~ x1 + x2 + x3 + x4, data = crops) crops.1 ``` \normalsize ## Assessing - 3 LDs (four variables, four groups). - 1st two important. - `LD1` mostly `x1` (plus) - `LD2` `x3` (minus) ## Predictions - Thus: \footnotesize ```{r bDiscrim-43} p <- predict(crops.1) crops.2 <- cbind(crops, p) with(crops.2, table(obs = crop, pred = class)) ``` \normalsize - Not very good, eg. only half the Soybeans and Sugarbeets classified correctly. ## Plotting the LDs ```{r piacentini,fig.height=3.4} ggplot(crops.2, aes(x = x.LD1, y = x.LD2, colour = crop)) + geom_point() ``` Corn (red) mostly left, cotton (green) sort of right, soybeans and sugarbeets (blue and purple) mixed up. ## Biplot ```{r bDiscrim-45, fig.height=6} ggord(crops.1, crops$crop, ellipse = FALSE, vec_ext = 10) # ggbiplot(crops.1, groups = crops$crop) ``` ## Comments - Corn low on LD1 (left), hence low on `x1` - Cotton tends to be high on LD1 (high `x1`) - one cotton very low on LD2 (high `x3`?) - Rather mixed up. ## Posterior probs (some) \scriptsize ```{r bDiscrim-52 } crops.2 %>% mutate(across(starts_with("posterior"), \(p) round(p, 3))) %>% rename_with(\(p) str_remove(p, "osterior"), starts_with("posterior")) %>% filter(crop != class) %>% select(crop, class, starts_with("p.")) ``` \normalsize ## Comments - These were the misclassified ones, but the posterior probability of being correct was not usually too low. - The correctly-classified ones are not very clear-cut either.