---
title: "Survey Analysis for Fake Client"
author: "Sean van der Merwe"
date: "`r Sys.Date()`"
format: 
  docx:
    reference-doc: wordreference01.docx
    toc: true
default-image-extension: svg
execute: 
  echo: false
  fig-format: svg
  fig-width: 6
  fig-height: 4
  out-width: "6in"
  out-height: "4in"
---

*Activate the navigation pane (View menu) for ease of navigation through the document.*

```{r}
#| include: false
library(knitr)
```

# Links

[Consultation Unit](https://www.ufs.ac.za/natagri/departments-and-divisions/mathematical-statistics-and-actuarial-science-home/statistical-consultation-unit/statistical-consultation-unit)

[Consultant](https://seanvdm.co.za)

[Frequently asked questions](https://seanvdm.co.za/post/faqconsult/)

# Introduction and disclaimer

Throughout our analysis, we must remember at all times that we are making judgements on the perceptions of people such as those who responded to this survey. [People who choose to respond to surveys may be systematically different from the general population.](https://www.dummies.com/education/math/statistics/statistics-conundrums-dealing-with-survey-nonresponders/) [See the Wikipedia page on Response Bias for more information.](https://en.wikipedia.org/wiki/Response_bias)

Similarly, the survey was only presented to a specific group of people, and we should not try to extend the results beyond that group without seriously considering and accounting for any systematic differences between that group and any broader group.

The data analysed here is inherently random. If the survey were to be repeated then the results will differ. We will try to estimate the extent to which this will be the case, but those estimates are themselves uncertain. 

While [the computer software used](https://en.wikipedia.org/wiki/R_(programming_language)) is tried and tested, the analysis involves multiple human elements. Both the client and statistician may have introduced human error at various stages of the research process.

**Thus, no guarantee can ultimately be given on the correctness of any findings below.**

Links are given occasionally to sources explaining the statistical concepts addressed. These are colloquial sources to aid understanding, not reference sources. However, most of these sources link to reference sources, should the researcher require academic references for specific topics.

```{r}
#| warning: false
#| message: false
#| results: hide
options(scipen = 12)
library(tidyverse)
library(openxlsx)
library(psych)
filecore <- 'FakeSurveyData'
```

```{r}
# Sean's functions that are useful for stuff
mapvalues <- function(x, from=unique(x), to=seq_along(from), warn_missing=FALSE) { 
    # Simplified replacement for mapvalues from plyr package, by Sean van der Merwe, UFS
    to <- unlist(to) # Make sure it's a vector. Also note that warn_missing is ignored.
    to[match(x, from)]
}
wrap_strings <- function(vector_of_strings, width){ as.character( sapply( vector_of_strings, \(x) { paste( strwrap(x, width=width), collapse="\n")}))}
sigfunc <- function(pvalue) { c('significant','not significant')[c(pvalue<=alpha, pvalue>alpha)] }
dropempty <- function(atbl) {
  emptyrows <- (rowSums(atbl,na.rm = TRUE)==0)
  emptycols <- (colSums(atbl,na.rm = TRUE)==0)
  if (any(emptyrows)) { atbl <- atbl[-which(emptyrows),] }
  if (any(emptycols)) { atbl <- atbl[,-which(emptycols)] }
  atbl }
UFScolours <- c(Blue="#0F204BFF", Red="#A71930FF", LightGray="#A7A8AAFF", DarkGray="#8D8D8EFF", NAS="#0039A7FF", EDU="#00675AFF", THEO="#C69317FF", HUM="#EA8400FF", SA="#A40084FF", SC="#9E83B7FF", QWA="#00B140FF", LAW="#BB133EFF", HLTH="#490E6FFF", Black="#000000FF", VeryLight="#D7D8DDFF")
nmatrixfunc <- function(X) {
  nc <- ncol(X)
  result <- matrix(1,nc,nc)
  for (i in 1:nc) {
    for (j in 1:nc) {
      result[i,j] <- sum(!(is.na(X[,i]) | is.na(X[,j])))
    }
  }
result }
corrplotfunction <- function(correldata, crosssize = 1.8, textsize = 0.9) {
  n <- nrow(correldata)
  corrmatsmall <- cor(correldata, use="pairwise.complete.obs")
  rownames(corrmatsmall) <- colnames(corrmatsmall) <- colnames(correldata)
  nmat <- nmatrixfunc(correldata)
  testmat <- abs(corrmatsmall)
  pmat <- (1-pt((testmat/sqrt(1-testmat^2))*sqrt(nmat-2),(nmat-2)))/2
  corrplot::corrplot(corrmatsmall, method = 'color', type = 'full', order = 'original', p.mat=pmat, insig = 'pch', pch.cex = crosssize, tl.cex = textsize, tl.col='black')
  corrmatsmall
}
cleantext <- function(mytext) { 
  mytext <- gsub("\\b\\d+\\b", "", mytext)
}
```


```{r}
# Read in data
d <- read.xlsx(paste0(filecore,'.xlsx'), 'Export')
info <- read.xlsx(paste0(filecore,'.xlsx'), 'Info')
```


```{r}
# Process data
nom_ques <- which(info$Type == 'Nominal')
multi_ques <- which(info$Type == 'MultiBox')
scale_ques <- which(info$Type == 'Scale')
text_ques <- which(info$Type == 'Text')
n_obs <- nrow(d)
n_nom <- length(nom_ques)
n_multi <- length(multi_ques)
n_scale <- length(scale_ques)
n_text <- length(text_ques)
groups <- info$Group[scale_ques] |> unique()
groups <- groups[!is.na(groups)]
n_grp <- length(groups)
alpha <- 0.05
```


```{r}
# Create factors and scales
scl_loc <- which(names(info) == 'ScaleLength')
scale_ques |> sapply(\(j) {
  mapvalues(d[[j]], info[j, (scl_loc + 1):(scl_loc + info$ScaleLength[j])])
}) -> num_data
for (j in which(!is.na(info$ScaleLength))) {
  d[[j]] <- factor(d[[j]], levels = info[j, (scl_loc + 1):(scl_loc + info$ScaleLength[j])])
}
```


# Statistical topics that might be relevant

## The boringness assumption

-   The classical (frequentist) framework of testing has many flaws and alternatives, but remains popular because of its simplicity
-   It starts by **assuming that everything is boring**
    -   This is called the null hypothesis
-   We assume that there is no difference between cases, treatments, situations, or groups
-   We assume that there is no change over time, nor space

## The uncertain experiment

-   We do an experiment and **see something that looks interesting**
-   No experiment is perfect, there are always errors and uncertainties
-   So we ask the question, "How interesting is this thing really?"
    -   Is this just a coincidence?
    -   If we did the experiment again, would the result be as interesting? More interesting? 
    -   Or would it be less interesting or even entirely different? 

## The p-value

-   The [p-value](https://en.wikipedia.org/wiki/P-value) is defined as the probability of seeing something at least as interesting as what we observed, under the assumption that nothing is really interesting in general.
-   If the p-value is large then we continue to assume that the null (boring) hypothesis is true. We **DO NOT** accept or conclude that the null hypothesis is true, that would be wrong. We simply say that *we have not found any evidence that contradicts it*.
-   If the p-value is small then we reject our null hypothesis and conclude that **there is at least one interesting thing** to see.
-   The probability is calculated by saying, "What might have happened had we repeated the experiment under the same settings an infinite number of times?"

## Significance level

-   So what is big and what is small for a p-value? 
-   **Before doing any experiments or collecting any data** we pick a number close to 0 and say that we are happy with that rate of false positives.
-   We call this number the significance level $(\alpha)$.
-   The significance level chosen is `r alpha`. [This means that we only look at results where the p-value is less than](https://en.wikipedia.org/wiki/P-value) `r alpha`, as other results could easily just be chance variation.
-   If you haven't done a statistical test then don't use the word *significant*.

## Multiple testing problem

-   The more things you try to test the more likely you are to stumble on an 'interesting' coincidence and end up talking nonsense
-   People who don't know what they want to know ask a ton of questions in the hope of getting 'lucky' and finding out something interesting, but chances are [they are just finding coincidences](https://en.wikipedia.org/wiki/Multiple_comparisons_problem)
-   In statistics we often say we are happy with a 5% rate of false positives, but that assumes you are doing 1 test - if you do 40 tests then you expect to get 2 false positives
-   It is possible to adjust for this in some ways, but you always lose power.
    -   The most common adjustment is the [Holm-Bonferroni method](https://en.wikipedia.org/wiki/Holm%E2%80%93Bonferroni_method)

## Tabulations

The most basic form of statistics is a simple tabulation of responses into frequencies and relative frequencies. 

Should one wish to compare two nominal variables to each other and test for independence then one might do a $\chi^2$ [(Chi Square) test](https://en.wikipedia.org/wiki/Pearson%27s_chi-squared_test).

## t-Test

-   [The t-test](https://en.wikipedia.org/wiki/Student%27s_t-test) is all about averages
-   The most popular statistical test for lots of reasons:
    -   People like to think in averages and ask whether things are the same on average
    -   It is easy to do and fairly reliable
    -   It flows from a Central Limit Theorem (and other theorems) that says:
        -   If your numbers follow any well-behaved pattern, then the average of a large enough sample will follow a *normal* distribution

## Two-sample t-test

-   The standard t-test compares the average of a single set of numbers to a predetermined number
-   If you have paired observations (e.g. matched before vs after) then you subtract first to end up with a single set of differences which you can compare to 0 (no difference)
-   But what if you have two separate set of numbers (e.g. Free State vs Gauteng)?
    -   You can compare the averages using a two-sample t-test
    -   The two-sample t-test is much less powerful and has a lot more assumptions though

## ANOVA {.smaller}

What if you have more than two sets of numbers?

-   [Analysis of Variance](https://en.wikipedia.org/wiki/Analysis_of_variance) (ANOVA) asks whether all the groups have the same average $(\mu_1=\mu_2=\mu_3=\dots)$
-   Suppose we ask each office of a business to rate their manager on a scale of 0 to 10. Each office operates independently so we can compare the manager ratings using an ANOVA to see if any managers are better or worse liked than the rest.
-   ANOVAs assume that the groups have the same variance, so if one manager plays favourites and has more variance in rating then the ANOVA will not work right.
-   If the ANOVA gives a small p-value then we reject the null hypothesis that all the averages are equal and conclude that **at least 1 average is different from the rest**.
    -   It does *not* imply that all the averages are different. To determine which averages are different people like to do a *post-hoc* test, such as the [Tukey HSD procedure](https://en.wikipedia.org/wiki/Tukey%27s_range_test)

## Rank-sum and U-test {.smaller}

If your data is badly behaved, skewed, unusual, or just doesn't really meet the assumptions above then people will recommend using a non-parametric alternative.

**However, non-parametric tests have different null hypotheses, and are not directly interchangeable.**

For comparing a single sample to a fixed value we might use a [Wilcoxon signed-rank test](https://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test). This is *not* testing the average at all, instead it tests the more general question of whether the values tend to be close to or far away from the null value.

For comparing two samples we can use the [Mann-Whitney U test](https://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U_test). This asks the question of whether one set of numbers is systematically larger or systematically smaller than another. 

Because these questions are more general, and make less assumptions, they are also more robust. However, they can have less power to find differences.

More advanced tests also have non-parametric alternatives, e.g. the [Kruskall-Wallis](https://en.wikipedia.org/wiki/Kruskal%E2%80%93Wallis_one-way_analysis_of_variance) approach, which is quite similar to the ANOVA approach, but more robust to violations of the ANOVA assumptions.

## Correlation

Correlation measures the tendency for things to move together.

A positive correlation says that when one thing is above average the other will **tend to** also be above average (below average $\leftrightarrow$ below average). Think of temperature and watermelon consumption say.

**Correlation done not imply causation.**

We calculate the [Pearson correlation coefficients](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient) using the scaled data. Positive correlations indicate responses that tend to vary (move around) in the same direction. Negative correlations are an indication of responses going counter to each other. Correlations that are not statistically different from zero will be marked as such with a cross.

## Congeneric reliability

In this section the goal is to determine to what extent specific sets of question measure the same thing. It is based on [factor analysis](https://en.wikipedia.org/wiki/Factor_analysis) - we try to find an underlying factor explaining the responses to a group of questions as a whole. [The measures calculated here](https://en.wikipedia.org/wiki/Congeneric_reliability) are related to the famous Cronbach's $\alpha$, but considered more accurate (although alpha is also given). 

Reliability measures lower than about 0.7 suggest that multiple concepts are being measured; while measures above 0.95 suggest that you are asking exactly the same thing in different words. [If you are trying to measure a single concept reliably then a value in the target range of 0.7 to 0.95 is desired.](https://en.wikipedia.org/wiki/Cronbach%27s_alpha)

# Univariate analysis

## Frequencies of responses

We create summary tables for each nominal question and then for each scale question as if it was purely categorical.

```{r}
#| results: asis
par(mar=c(4.4,8.4,0.4,0.6))
for (i in c(nom_ques, scale_ques)) {
    currentx <- d[[i]]
    validx <- !is.na(currentx)
    x <- currentx[validx]
    y <- table(x)
    miss <- (n_obs-sum(validx))
    if (miss > 0) {    y <- c(y, `Not Answered`=miss)   }
    z <- round(y/sum(y),2)
    cat('\n\n###', info$Desc[i], '\n\n')
    w <- data.frame(Response=names(y),Frequency=c(y),Relative.Frequency=c(z), row.names = 1:length(y))
    print(kable(w))
    barplot(z*100, names.arg = wrap_strings(names(y), 28), col=viridisLite::viridis(length(y)), xlab="Percentage responded in category", main='', horiz = TRUE, las=1, cex.axis = 0.7, cex.names = 0.7, xlim = c(0,100))
    cat('\n\n---\n\n')
}
```

## Scale questions summary

We summarise the scale responses in stages.

### Quantiles

First we start by using quantiles (based on sorting in ascending order). Specifically, we note the count of responses, the minimum, the first quartile, the median, the third quartile, and the maximum.

*This is the least interesting and least relevant of all the results.* It is generally just included in an appendix for interest.

```{r}
results <- matrix(0, n_scale, 6)
for (i in 1:n_scale) {
  xx <- num_data[,i]
  xx <- xx[!is.na(xx)]
  cnt <- length(xx)
  minx <- min(xx)
  maxx <- max(xx)
  midl <- quantile(xx, c(0.25, 0.5, 0.75))
  results[i,] <- c(cnt, minx, midl, maxx)
}
colnames(results) <- c('Count','Min','Q1','Median','Q3','Max')
data.frame(Ques = info$Tag[scale_ques], results, Description = info$Desc[scale_ques]) |>
kable()
```

### Numeric summary

Then we summarise the numeric scales by their mean, standard deviation, standard error of the mean, confidence interval for the mean, expected value under the null hypothesis that the responses are purely noise, and p-value of testing that null hypothesis. 

*This is only a check that people were not neutral in their responses to the questions.* If we put the responses on a scale of 1 to 5 then 3 is a neutral response and we would like to see that people did not just answer 3 on average.

```{r}
results <- matrix(0,n_scale,7)
for (i in 1:n_scale) {
  xx <- num_data[,i]
  xx <- xx[!is.na(xx)]
  cnt <- length(xx)
  mx <- mean(xx)
  sx <- sd(xx)
  se <- sx/sqrt(cnt)
  lwr <- mx + qt(alpha/2, cnt-1)*se
  upr <- mx + qt(1 - alpha/2, cnt-1)*se
  minx <- min(xx)
  maxx <- max(xx)
  exx <- (maxx + 1)/2
  pvalx <- t.test(xx, mu=exx)$p.value
  results[i,] <- c(mx, sx, se, lwr, upr, exx, pvalx)
}
results <- cbind(results, p.adjust(results[,7]))
colnames(results) <- c('Mean', 'Std.Dev', 'SE', 'Lower', 'Upper', 'NullExp', 'pvalue', 'p_adj')
data.frame(Ques = info$Tag[scale_ques], results, Description = info$Desc[scale_ques]) |>
kable(digits = 3)
```

### Box plots

We try to plot each question group using [box plots](https://en.wikipedia.org/wiki/Box_plot).

*This is a neat visual summary of the responses, showing the variation clearly.*

```{r}
#| results: asis
par(mar=c(2.5,12,0.4,1))
for (g in seq_along(groups)) {
  current <- which(info$Group[scale_ques] == groups[g])
  cat('\n\n### Question group ', groups[g], '\n\n', sep = '')
  d[scale_ques[current]] |> boxplot(names = wrap_strings(info$Desc[scale_ques[current]], 40),
                                   col=viridisLite::viridis(length(current)), xlab="", main='',
                                   horizontal = TRUE, las=1, cex.axis = 0.7, cex.names = 0.7)
  cat('\n\n---\n\n')
}
```

## Tick Multiple Options Sections

Here we analyse those questions where people could tick multiple options. These must be analysed separately as they must be handled in groups. Note that percentages here should add up to more than 100%.

*This type of question does not lend itself to any statistical analysis and its use is generally discouraged.*

```{r}
#| results: asis
par(mar=c(2.5,15,0.4,1))
for (i in multi_ques) {
  cat('\n\n###', info$Desc[i], '\n\n')
  d[[i]] |> strsplit(";") |> unlist() |> table() -> y
  z <- round(y/sum(y),2)
  w <- data.frame(Response=names(y),Frequency=c(y),Relative.Frequency=c(z), row.names = 1:length(y))
  print(kable(w))
  pos <- barplot(z*100, names.arg = wrap_strings(names(y), 44), 
                 col=viridisLite::viridis(length(y)), xlab="", main='', horiz = TRUE, 
                 las=1, cex.axis = 0.6, cex.names = 0.6, xlim = c(0,100))
  for (j in seq_along(pos)) {
    text(z[j]*100, pos[j], paste0(round(z[j]*100,0),"%"), pos=4)
  }
}
```

## Text questions

The text questions are summarised with a crude word cloud. This is not a substitute for a formal thematic analysis, and is mostly for interest.

*Word clouds can be enlightening if differentiated based on a binary explanatory variable.*

```{r}
#| results: asis
#| warning: false
#| fig-height: 4
#| out-height: "4in"
#| fig-width: 4
#| out-width: "4in"
par(mar=rep(0.2,4))
for (i in text_ques) {
  cat('\n\n###', info$Desc[i], '\n\n')
  comments <- cleantext(d[[i]])
  comments <- comments[!is.na(comments)]
  if (length(comments) > 0) {
    merged_com <- paste(comments, collapse=" ")
    tidiedtext <- tokenizers::tokenize_words(merged_com, simplify=TRUE)
    freq <- table(tidiedtext)
    freqfilt <- freq[nchar(names(freq))>4]
    wordcloud::wordcloud(words = names(freqfilt), freq = freqfilt, min.freq = 2, max.words=20, random.order=FALSE, rot.per=0.35, colors=RColorBrewer::brewer.pal(5, "Dark2"))
  }
  cat('\n\n---\n\n')
}
```


# Cross tabulations

We test whether the responses to each question are dependent on the response to select key questions.

## Nominal versus Nominal

Here we do tests of independence on the frequency tables, similar to the classical Chi-square test of independence, but using simulated p-values to account for the imbalance in groups.

```{r}
#| results: asis
explanatory_factors <- c(1)
for (ex_fac in explanatory_factors) {
  cat('\n\n### Analysis of nominal variables by', info$Desc[ex_fac], '\n\n')
  fac <- d[[ex_fac]]
  v_fac <- !is.na(fac)
  chisq_ques <- c(2)
  chisq_ques |> sapply(\(i) {
    x <- d[[i]]
    v_x <- !is.na(x)
    tbl <- table(x[v_fac & v_x], fac[v_fac & v_x])
    pval <- chisq.test(tbl, simulate.p.value = TRUE)$p.value
    class(tbl) <- 'matrix'
    cat('\n\n####', info$Desc[i], '\n\n')
    tbl |> kable() |> print()
    pval
  }) -> pvals
  cat('\n\n#### Summary\n\n')
  data.frame(Ques = info$Tag[chisq_ques], 
             raw_p_value = pvals, 
             adj_p_value = p.adjust(pvals), 
             Description = info$Desc[chisq_ques]) |> 
    kable(digits = 3) |> print()
  cat('\n\n---\n\n')
}
```


## Scale versus Nominal

Here we would use a Kruskal-Wallis, ANOVA, or similar analysis.

```{r}
#| results: asis
par(mar=c(4.4,7.4,0.4,0.6))
explanatory_factors <- c(1,2)
for (ex_fac in explanatory_factors) {
  cat('\n\n### Analysis of scale variables by', info$Desc[ex_fac], '\n\n')
  fac <- d[[ex_fac]]
  v_fac <- !is.na(fac)
  scale_ques |> sapply(\(i) {
    x <- d[[i]]
    v_x <- !is.na(x)
    y <- x[v_fac & v_x]
    w <- fac[v_fac & v_x]
    boxplot(y ~ w, names = wrap_strings(levels(w), 28), 
            horizontal = TRUE, main = '', col = UFScolours, 
            xlab = info$Desc[i], ylab = '', las = 1, cex.axis = 0.7, 
            cex.names = 0.7)
    kruskal.test(y ~ w)$p.value
  }) -> pvals
  cat('\n\n#### Summary\n\n')
  data.frame(Ques = info$Tag[scale_ques], 
             raw_p_value = pvals, 
             adj_p_value = p.adjust(pvals), 
             Description = info$Desc[scale_ques]) |> 
    kable(digits = 3) |> print()
  cat('\n\n---\n\n')
}
```

## Correlations

```{r}
#| results: asis
#| warning: false
#| fig-height: 7
#| out-height: "7in"
#| fig-width: 7
#| out-width: "7in"
seq_len(n_grp) |> sapply(\(g) {
  cat('\n\n###', groups[g], '\n\n')
  cor_data <- num_data[,which(info$Group[scale_ques] == groups[g])]
  colnames(cor_data) <- info$Tag[which(info$Group[scale_ques] == groups[g])]
  cor_data |> corrplotfunction()
  cat('\n\n---\n\n')
}) -> nul

```

## Congeneric reliability

```{r}
#| results: asis
for (g in 1:n_grp) {
  cat('\n\n###', groups[g], '\n\n')
  cor_data <- num_data[,which(info$Group[scale_ques] == groups[g])]
  colnames(cor_data) <- info$Tag[which(info$Group[scale_ques] == groups[g])]
  cor_data |> alpha() -> alpha_results
  alpha_results$total |> list2DF() |> kable(digits = 3) |> print()
  cat('\n\n---\n\n')
}
```

Reliability calculations assume that all questions are asked in the same direction. Should this not be the case then the researcher should indicate on the metadata spreadsheet (Info) which questions need to be flipped (using a column of indicators) before doing this analysis.

# Index construction

Sometimes it is meaningful to collapse sections into a reduced number of latent dimensions. When a number of responses is collapsed to a single dimension this is often called an index. 

*In this survey the groups have low reliability. Thus, we will not attempt to collapse into reduced dimensions, nor analyse such dimensions or underlying factors at this time.* 

There are a few methods by which this might be done in general. One approach is to use a factor analysis, another is to use [principal component analysis](https://en.wikipedia.org/wiki/Principal_component_analysis). PCA seeks to find a linear combination of responses that captures as much as possible of the variation in as few dimensions as possible. 

Often it makes sense to use an average, or weighted average with weights chosen manually. This is particularly useful when there are weights suggested by theory or external data. 

*Sometimes you will see the use of simple sums to construct indices. This should be avoided as even a single missing value, or a single group with a different number of questions, will render all results not interpretable.*




---

# Referencing the software

In case it is useful, I give the referencing guides from the software itself.

```{r}
#| results: asis
citation()
```

