---
title: "Questionnaires and Surveys: Analyses with R"
author: "Martin Schweinberger"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output:
  bookdown::html_document2:
    includes:
      in_header: GoogleAnalytics.html
bibliography: bibliography.bib
link-citations: yes
---


```{r uq1, echo=F, fig.cap="", message=FALSE, warning=FALSE, out.width='100%'}
knitr::include_graphics("https://slcladal.github.io/images/uq1.jpg")
```

# Introduction

This tutorial offers some advice on what to consider when creating surveys and questionnaires, provides tips on visualizing survey data, and exemplifies how survey and questionnaire data can be analyzed. As this tutorial is introductory, issues relating to what software to use when creating a survey (e.g. SurveyMonkey, Qualtrics, GoogleForms, etc.) or how to program questionnaires or online experiments in Java or R are not discussed. 

```{r diff, echo=FALSE, out.width= "15%", out.extra='style="float:right; padding:10px"'}
knitr::include_graphics("https://slcladal.github.io/images/gy_chili.jpg")
```

This tutorial is aimed at beginners and intermediate users of R with the aim of showcasing how to visualize and analyze survey and questionnaire data using R. The aim is not to provide a fully-fledged analysis but rather to show and exemplify selected useful methods associated with surveys and questionnaires. 


<div class="warning" style='padding:0.1em; background-color:#f2f2f2; color:#51247a'>
<span>
<p style='margin-top:1em; text-align:center'>
The entire R Notebook for the tutorial can be downloaded [**here**](https://slcladal.github.io/content/surveys.Rmd).  If you want to render the R Notebook on your machine, i.e. knitting the document to html or a pdf, you need to make sure that you have R and RStudio installed and you also need to download the [**bibliography file**](https://slcladal.github.io/content/bibliography.bib) and store it in the same folder where you store the Rmd file. <br></p>
<p style='margin-left:1em;'>
</p></span>
</div>

<br>


**How to use the R Notebook for this tutorial**

As all calculations and visualizations in this tutorial rely on R, it is necessary to install R and RStudio. If these programs (or, in the case of R, environments) are not already installed on your machine, please search for them in your favorite search engine and add the term "download". Open any of the first few links and follow the installation instructions (they are easy to follow, do not require any specifications, and are pretty much self-explanatory).

To follow this tutorial interactively (by using the R Notebook), follow the instructions listed below.

1. Create a folder somewhere on your computer
2. Download the [R Notebook](https://slcladal.github.io/content/surveys.Rmd) and save it in the folder you have just created
3. Open R Studio
4. Click on "File" in the upper left corner of the R Studio interface
5. Click on "New Project..."
6. Select „Existing Directory“
7. Browse to the folder you have just created and click on “Create New Project”
8. Now click on "Files" above the lower right panel
9. Click on the file "surveys.Rmd"
  + The Markdown file of this tutorial should now be open in the upper left panel of R Studio. To execute the code blocks used in this session, simply click on the green arrows in the top right corner of the code boxes.
  + To render a PDF of this tutorial, simply click on "Knit" above the upper left panel in R Studio.

# Design Basics

A survey is a research method for gathering information based on a sample of people. Questionnaires are a research instrument and typically represent a part of a survey, i.e. that part where participants are asked to answer a set of questions. [@brown2001surveys, 6] defines questionnaires as "any written instruments that present respondents with a series of  questions or statements to which they are to react either by writing out their answers or selecting among existing answers."  

Questionnaires elicit three types of data:

* Factual 
* Behavioral
* Attitudinal

While factual and behavioral questions are about what the respondent is and does, attitudinal questions tap into what the respondent thinks or feels.

The advantages of surveys are that they 
* offer a relative cheap, quick, and effective way to collect (targeted) data from a comparatively large set of people; and 
* that they can be distributed or carried out in various formats (face-to-face, by telephone, by computer or via social media, or by postal service).

Disadvantages of questionnaires are that they are prone to providing unreliable or unnatural data. Data gathered via surveys can be unreliable due to the social desirability bias which is the tendency of respondents to answer questions in a manner that will be viewed favorably by others. Thus, the data that surveys provide may not necessarily be representative of actual natural behavior. 

Questionnaires and surveys are widely used in language research and thus one of the most common research designs. In this section, we will discuss what needs to kept in mind when designing questionnaires and surveys, what pieces of software or platforms one can use, options for visualizing questionnaire and survey data, statistical methods that are used to evaluate questionnaire and survey data (reliability), and which statistical methods are used in analyzing the data.

# Things to consider

Here are some rules to consider during the creation of questionnaires and before  any survey is distributed.

* Surveys should not be longer than they have to be while they have to be long enough to collect all the data that are needed. It is crucial that a questionnaire collects all necessary data (including socio-demographic details). 

* The language should be simple and easy to understand - this means that jargon should be avoided. Also, leading questions and value judgement on the side of the creators of questionnaires should be avoided to prevent social desirability bias.

* Before distributing a questionnaire, it should be piloted. Piloting is essential to check if respondents understand the questions as intended, and to check how long it takes to answer the questions. Also, the people who are involved in the piloting should be allowed to provide feedback to avoid errors. 

* When questions go beyond simply collecting socio-demographic details and if the data contains test and filler items, the order of questions (within blocks) should be quasi-randomized. Quasi-randomization means that test items are not asked in direct succession and that they do not appear as first or last items. Quasi-randomization helps to avoid fatigue effects or results that are caused by the ordering of questions. However, the questions should still follow an internally consistent logic so that related questions appear in the same block. Also, more specific questions should be asked after more general questions.

* Questions must be unambiguous and they cannot ask multiple aspects at once. Take, for instance, the following question "Do you consider UQ to be a good university with respect to teaching and research?" If the respondent answers positively, then no issues arise but if the answer is negative, i.e. "No", then we do not know if the respondent thinks that UQ is not a good university with respect to teaching OR with respect to research OR both! In such cases, questions should be split:

  + Do you consider UQ to be a good university with respect to teaching?

  + Do you consider UQ to be a good university with respect to research?"

* To check if respondents are concentrated, read the questions carefully, and answering truthfully, it is useful to include reverse questions, i.e. questions that have the opposite polarity. Reverse questions allow to check if respondents only answers "very satisfied" or "completely agree" without respect to the content of the question. Giving the same answer to questions which have opposite propositions would indicate that respondents do not read questions carefully or do not answer truthfully.

* If questions are not open or unstructured, i.e. if different options to answer to a question are provided, it is crucial that the options are fine-grained enough so that the data that is collected allows us to answer the research question that we want to investigate. In this context, the scaling of answer options is important. Scales reflect different types of answering options and they come in three basic forms: nominal and categorical, ordinal, or numeric. 

  + *Nominal and categorical scales*: 
Nominal and categorical scales only list the membership of a particular class. Nominal scales offer exactly two options (yes/no or on/off), while categorical scales offer several options (e.g. the state in which someone was born).

  + *Ordinal scales*: 
With *ordinal scales* it is possible to rank the values, but the distances between the ranks can not be exactly quantified. An example of an ordinal scales is the ranking in a 100-meter run. The 2nd in a 100-meter run did not go twice as fast as the 4th. It is often the case that ordinal variables consist of integer, positive numbers (1, 2, 3, 4, etc.). In the context of surveys, ordinal scales are the most important as all Likert scales (after the psychologist Rensis Likert) are ordinal scales. The levels of the typical five-level Likert item could be: *Strongly disagree (1)*, *Disagree (2) *, *Neither agree nor disagree (3)*, *Agree (4)*, and *Strongly agree (5)*. As such, the Likert scale is a bipolar scale that can be balanced, if there is an uneven number of options with the center option being neutral, or unbalanced, if there are an even number of options which forces respondents to express a preferences for wither of the two poles (this is called a  "forced choice" method.

  + *(True) Numeric scales*: 
There are two basic types of numeric scales: *interval-scales* and *ratio-scales*. For *interval scales*, the differences between levels are significant, but not the relationship between levels. For instance, 20 degree Celsius is not twice as hot as 10 degree Celsius. For *ratio-scales* both the differences and the relationships between the levels are significant (e.g. the times in a 100-meter dash: 10 is exactly twice as high as 5 and half as much as 20).

Of these scales, numeric is the most informative and questionnaires should always aim to extract the most detailed information without becoming to long.

# Visualizing survey data

Just as the data that is provided by surveys and questionnaires can take various forms, there are numerous ways to display survey data. In the following, we will have a look at some of the most common or useful ways in which survey and questionnaire data can be visualized. However, before we can begin, we need to set up our R session as shown below.

**Preparation and session set up**

To run the scripts shown below without errors, certain *packages* need to be installed from an R *library*. Before turning to the code below, please install the packages by running the code below this paragraph. If you have already installed the packages mentioned below, then you can skip ahead ignore this section. To install the necessary packages, simply run the following code - it may take some time (between 1 and 5 minutes to install all of the libraries so you do not need to worry if it takes some time).

```{r prep_01, echo=T, eval = F, message=FALSE, warning=FALSE}
# install packages
install.packages("knitr")
install.packages("lattice")
install.packages("tidyverse")
install.packages("likert")
install.packages("MASS")
install.packages("psych")
install.packages("viridis")
install.packages("ggplot2")
install.packages("here")
install.packages("flextable")
install.packages("devtools")
install.packages("GPArotation")
# devtools::install_github("matherion/userfriendlyscience", dependencies=T)
install.packages("ufs")
# install klippy for copy-to-clipboard button in code chunks
install.packages("remotes")
remotes::install_github("rlesur/klippy")
```

You can now activate the packages by running the code chunk below. 

```{r prep_02, message=FALSE, warning=FALSE}
# set options
options(stringsAsFactors = F)         # no automatic data transformation
options("scipen" = 100, "digits" = 4) # suppress math annotation
# install packages
library(knitr)
library(lattice)
library(tidyverse)
library(likert)
library(MASS)
library(psych)
library(viridis)
library(ggplot2)
library(here)
library(flextable)
library(devtools)
# library(userfriendlyscience)
# activate klippy for copy-to-clipboard button
klippy::klippy()
```

Once you have installed R, RStudio, and have also initiated the session by executing the code shown above, you are good to go.

## Line graphs for Likert-scaled data{-}

A special case of line graphs is used when dealing with Likert-scaled variables (we will talk about Likert scales in more detail below). In such cases, the line graph displays the density of cumulative frequencies of responses. The difference between the cumulative frequencies of responses displays differences in preferences. We will only focus on how to create such graphs using the `ggplot` environment here as it has an in-build function (`ecdf`) which is designed to handle such data.

In a first step, we load a data set (ldat) which contains Likert-scaled variables. This data set represents fictitious rating of students from courses about how satisfied they were with their learning experience. The response to the Likert item is numeric so that *strongly disagree/very dissatisfied* would get the lowest (1) and *strongly agree/very satisfied* the highest numeric value (5).

```{r line_0a1, message=FALSE, warning=FALSE}
# define color vectors
clrs3 <- c("firebrick4",  "gray70", "darkblue")
clrs5 <- c("firebrick4", "firebrick1", "gray70", "blue", "darkblue")
# load data
ldat <- base::readRDS(url("https://slcladal.github.io/data/lid.rda", "rb"))
```

Let's briefly inspect the `ldat` data set.

```{r line_01b, echo = F}
ldat %>%
  as.data.frame() %>%
  head(10) %>%
  flextable() %>%
  flextable::set_table_properties(width = .5, layout = "autofit") %>%
  flextable::theme_zebra() %>%
  flextable::fontsize(size = 12) %>%
  flextable::fontsize(size = 12, part = "header") %>%
  flextable::align_text_col(align = "center") %>%
  flextable::set_caption(caption = "First 10 rows of the ldat data set.")  %>%
  flextable::border_outer()
```

The ldat data set has only two columns: a column labeled *Course* which has three levels (*German*, *Japanese*, and *Chinese*) and a column labeled *Satisfaction* which contains values from 1 to 5 which represent values ranging from *very dissatisfied* to *very satisfied*. Now that we have data resembling a Likert-scaled item from a questionnaire, we will display the data in a cumulative line graph.

```{r line_03, echo=T, message=FALSE, warning=FALSE}
# create cumulative density plot
ldat %>%
  ggplot(aes(x = Satisfaction, color = Course)) + 
  geom_step(aes(y = ..y..), stat = "ecdf") +
  labs(y = "Cumulative Density") + 
  scale_x_discrete(limits = c("1","2","3","4","5"), 
                   breaks = c(1,2,3,4,5),
                   labels=c("very dissatisfied", "dissatisfied", 
                            "neutral", "satisfied", "very satisfied")) + 
  scale_colour_manual(values = clrs3) + 
  theme_bw() 
```


The satisfaction of the German course was the lowest as the red line shows the highest density (frequency of responses) of *very dissatisfied* and *dissatisfied* ratings. The students in our fictitious data set were most satisfied with the Chinese course as the blue line is the lowest for *very dissatisfied* and "dissatisfied" ratings while the difference between the courses shrinks for "satisfied" and *very satisfied*. The Japanese language course is in-between the German and the Chinese course.  

## Pie charts{-}

Most commonly, the data for visualization comes from tables of absolute frequencies associated with a categorical or nominal variable. The default way to visualize such frequency tables are pie charts and bar plots. In a first step, we modify the data to get counts and percentages. 

```{r pie_01, echo = T, eval=T, message=FALSE, warning=FALSE}
# create bar plot data
bdat <- ldat %>%
  dplyr::group_by(Satisfaction) %>%
  dplyr::summarise(Frequency = n()) %>%
  dplyr::mutate(Percent = round(Frequency/sum(Frequency)*100, 1)) %>%
  # order the levels of Satisfaction manually so that the order is not alphabetical
  dplyr::mutate(Satisfaction = factor(Satisfaction, 
                                      levels = 1:5,
                                      labels = c("very dissatisfied",
                                                 "dissatisfied", 
                                                 "neutral", 
                                                 "satisfied", 
                                                 "very satisfied")))
```


Let's briefly inspect the new data set.

```{r pie_02, echo = F}
bdat %>%
  as.data.frame() %>%
  flextable() %>%
  flextable::set_table_properties(width = .5, layout = "autofit") %>%
  flextable::theme_zebra() %>%
  flextable::fontsize(size = 12) %>%
  flextable::fontsize(size = 12, part = "header") %>%
  flextable::align_text_col(align = "center") %>%
  flextable::set_caption(caption = "Overview of the bardata data set.")  %>%
  flextable::border_outer()
```

Before creating bar plots, we will briefly turn to pie charts because pie charts are very common despite suffering from certain shortcomings. Consider the following example which highlights some of the issues that arise when using pie charts.

```{r pie_03, echo = T, results = 'asis', message=FALSE, warning=FALSE}
# create pie chart
bdat %>%
  ggplot(aes("", Percent, fill = Satisfaction)) + 
  geom_bar(stat="identity", width=1, color = "white") +
  coord_polar("y", start=0) +
  scale_fill_manual(values = clrs5) +
  theme_void()
```


If the slices of the pie chart are not labelled, it is difficult to see which slices are smaller or bigger compared to other slices. This problem can easily be avoided when using a bar plot instead. This issue can be avoided by adding labels to pie charts. The labeling of pie charts is, however, somewhat tedious as the positioning is tricky. Below is an example for adding labels without specification.

```{r pie_05, echo = T, results = 'asis', message=FALSE, warning=FALSE}
# create pie chart
bdat %>%
  ggplot(aes("", Percent, fill = Satisfaction)) +
  geom_bar(stat="identity", width=1, color = "white") +
  coord_polar("y", start=0) +
  scale_fill_manual(values = clrs5) +
  theme_void() +
  geom_text(aes(y = Percent, label = Percent), color = "white", size=6)
```

To place the labels where they make sense, we will add another variable to the data called "Position".

```{r pie_07, echo = T, results = 'asis', message=FALSE, warning=FALSE}
pdat <- bdat %>%
  dplyr::arrange(desc(Satisfaction)) %>%
  dplyr::mutate(Position = cumsum(Percent)- 0.5*Percent)
```

Let's briefly inspect the new data set.

```{r pie_08, echo = F}
pdat %>%
  as.data.frame() %>%
  flextable() %>%
  flextable::set_table_properties(width = .5, layout = "autofit") %>%
  flextable::theme_zebra() %>%
  flextable::fontsize(size = 12) %>%
  flextable::fontsize(size = 12, part = "header") %>%
  flextable::align_text_col(align = "center") %>%
  flextable::set_caption(caption = "Overview of the new data.")  %>%
  flextable::border_outer()
```

Now that we have specified the position, we can include it into the pie chart.

```{r pie_09, echo = T, results = 'asis', message=FALSE, warning=FALSE}
# create pie chart
pdat %>%
  ggplot(aes("", Percent, fill = Satisfaction)) + 
  geom_bar(stat="identity", width=1, color = "white") +
  coord_polar("y", start=0) +
  scale_fill_manual(values = clrs5) +
  theme_void() +
  geom_text(aes(y = Position, label = Percent), color = "white", size=6)
```


We will now create separate pie charts for each course. In a first step, we create a data set that does not only contain the Satisfaction levels and their frequency but also the course.

```{r pie_11, message=FALSE, warning=FALSE}
# create grouped pie data
gldat <- ldat %>%
  dplyr::group_by(Course, Satisfaction) %>%
  dplyr::summarise(Frequency = n()) %>%
  dplyr::mutate(Percent = round(Frequency/sum(Frequency)*100, 1),
                Satisfaction = factor(Satisfaction, 
                                      levels = 1:5,
                                      labels = c("very dissatisfied",
                                                 "dissatisfied", 
                                                 "neutral", 
                                                 "satisfied", 
                                                 "very satisfied"))) %>%
  dplyr::arrange(desc(Satisfaction)) %>%
  dplyr::mutate(Position = cumsum(Percent)- 0.5*Percent)
```


Let's briefly inspect the new data set.

```{r pie_12, echo = F}
gldat %>%
  as.data.frame() %>%
  flextable() %>%
  flextable::set_table_properties(width = .5, layout = "autofit") %>%
  flextable::theme_zebra() %>%
  flextable::fontsize(size = 12) %>%
  flextable::fontsize(size = 12, part = "header") %>%
  flextable::align_text_col(align = "center") %>%
  flextable::set_caption(caption = "Overview of the grouped piedata.")  %>%
  flextable::border_outer()
```

Now that we have created the data, we can plot separate pie charts for each course.

```{r pie_15, echo = T, results = 'asis', message=FALSE, warning=FALSE}
# create pie chart
gldat %>%
  ggplot(aes("", Percent, fill = Satisfaction)) + 
  facet_wrap(~Course) +
  geom_bar(stat="identity", width=1, color = "white") +
  coord_polar("y", start=0) +
  scale_fill_manual(values = clrs5) +
  theme_void() +
  geom_text(aes(y = Position, label = Percent), color = "white", size=4)
```

## Bar plots{-}

Like pie charts, bar plot display frequency information across categorical variable levels.

```{r bar_01, echo = T, results='hide', message=FALSE, warning=FALSE}
# bar plot
bdat %>%
  ggplot(aes(Satisfaction, Percent, fill = Satisfaction)) +
  # determine type of plot
  geom_bar(stat="identity") +          
  # use black & white theme
  theme_bw() +                         
  # add and define text
  geom_text(aes(y = Percent-5, label = Percent), color = "white", size=3) + 
  # add colors
  scale_fill_manual(values = clrs5) +
  # suppress legend
  theme(legend.position="none")
```


Compared with the pie chart, it is much easier to grasp the relative size and order of the percentage values which shows that pie charts are unfit to show relationships between elements in a graph and, as a general rule of thumb, should be avoided.

Bar plots can be grouped which adds another layer of information that is particularly useful when dealing with frequency counts across multiple categorical variables. But before we can create grouped bar plots, we need to create an appropriate data set. 

```{r bar_03, echo = T, results='hide', message=FALSE, warning=FALSE}
# create bar plot data
gldat <- ldat %>%
  dplyr::group_by(Course, Satisfaction) %>%
  dplyr::summarise(Frequency = n()) %>%
  dplyr::mutate(Percent = round(Frequency/sum(Frequency)*100, 1)) %>%
  dplyr::mutate(Satisfaction = factor(Satisfaction, 
                                      levels = 1:5,
                                      labels = c("very dissatisfied",
                                                 "dissatisfied", 
                                                 "neutral", 
                                                 "satisfied", 
                                                 "very satisfied")))
```


Let's briefly inspect the data set.

```{r bar_04, echo = F}
gldat %>%
  as.data.frame() %>%
  flextable() %>%
  flextable::set_table_properties(width = .5, layout = "autofit") %>%
  flextable::theme_zebra() %>%
  flextable::fontsize(size = 12) %>%
  flextable::fontsize(size = 12, part = "header") %>%
  flextable::align_text_col(align = "center") %>%
  flextable::set_caption(caption = "Overview of the grouped bardata.")  %>%
  flextable::border_outer()
```

We have now added *Course* as an additional categorical variable and will include Course as the "fill" argument in our bar plot. To group the bars, we use the command "position=position_dodge()".

```{r bar_05, echo = T, results='hide', message=FALSE, warning=FALSE}
# bar plot
gldat %>%
  ggplot(aes(Satisfaction, Frequency, fill = Course)) + 
  geom_bar(stat="identity", position = position_dodge()) +
  # define colors
  scale_fill_manual(values = clrs3) + 
  # add text
  geom_text(aes(label=Frequency), vjust=1.6, color="white", 
            # define text position and size
            position = position_dodge(0.9),  size=3.5) + 
  theme_bw()                         
```

Bar plots are particularly useful when visualizing data obtained through Likert items. As this is a very common issue that empirical researchers face. There are two basic ways to display Likert items using bar plots: grouped bar plots and more elaborate scaled bar plots.

Although we have seen above how to create grouped bar plots, we will repeat it here with the language course example used above when we used cumulative density line graphs to visualize how to display Likert data.  

In a first step, we recreate the data set which we have used above. The data set consists of a Likert-scaled variable (Satisfaction) which represents rating of students from three courses about how satisfied they were with their language-learning course. The response to the Likert item is numeric so that "strongly disagree/very dissatisfied" would get the lowest and "strongly agree/very satisfied" the highest numeric value. 

Again, we can also plot separate bar graphs for each class by specifying "facets".

```{r bar_13, echo=T, message=FALSE, warning=FALSE}
# create grouped bar plot
gldat %>%
  ggplot(aes(Satisfaction, Frequency, 
             fill = Satisfaction, 
             color = Satisfaction)) +
  facet_grid(~Course) +
  geom_bar(stat="identity", position=position_dodge()) +
  geom_line() +
  # define colors
  scale_fill_manual(values=clrs5) +
  scale_color_manual(values=clrs5) +
  # add text and define color
  geom_text(aes(label=Frequency), vjust=1.6, color="white", 
            # define text position and size
            position = position_dodge(0.9),  size=3.5) +     
  theme_bw() +
  theme(axis.text.x=element_blank())
```

Another and very interesting way to display such data is by using the Likert package. In a first step, we need to activate the package, clean the data, and extract a subset for the data visualization example.

One aspect that is different to previous visualizations is that, when using the Likert package, we need to transform the data into a "likert" object (which is, however, very easy and is done by using the "likert()" function as shown below).

```{r likert1}
sdat  <- base::readRDS(url("https://slcladal.github.io/data/sdd.rda", "rb"))
```

```{r likert2, echo = F}
sdat %>%
  as.data.frame() %>%
  head(10) %>%
  flextable() %>%
  flextable::set_table_properties(width = .5, layout = "autofit") %>%
  flextable::theme_zebra() %>%
  flextable::fontsize(size = 12) %>%
  flextable::fontsize(size = 12, part = "header") %>%
  flextable::align_text_col(align = "center") %>%
  flextable::set_caption(caption = "First 10 rows of the survey data.")  %>%
  flextable::border_outer()
```

As you can see, we need to clean and adapt the column names. To do this, we will 

* add an identifier which shows which question we are dealing with (e.g. Q 1: question text)
* remove the dots between words with spaces
* add a question mark at the end of questions
* remove superfluous white spaces

```{r likert7}
# clean column names
colnames(sdat)[3:ncol(sdat)] <- paste0("Q ", str_pad(1:10, 2, "left", "0"), ": ", colnames(sdat)[3:ncol(sdat)]) %>%
  stringr::str_replace_all("\\.", " ") %>%
  stringr::str_squish() %>%
  stringr::str_replace_all("$", "?")
# inspect column names
colnames(sdat)
```  

Now, that we have nice column names, we will replace the numeric values (1 to 5) with labels ranging from *disagree* to *agree* and convert our data into a data frame.

```{r likert9}
lbs <- c("disagree", "somewhat disagree", "neither agree nor disagree",  "somewhat agree", "agree")
survey <- sdat %>%
  dplyr::mutate_if(is.character, factor) %>%
  dplyr::mutate_if(is.numeric, factor, levels = 1:5, labels = lbs) %>%
  drop_na() %>%
  as.data.frame()
```

```{r likert11, echo = F}
survey %>%
  as.data.frame() %>%
  head(10) %>%
  flextable() %>%
  flextable::set_table_properties(width = .5, layout = "autofit") %>%
  flextable::theme_zebra() %>%
  flextable::fontsize(size = 12) %>%
  flextable::fontsize(size = 12, part = "header") %>%
  flextable::align_text_col(align = "center") %>%
  flextable::set_caption(caption = "First 10 rows of the edited survey data.")  %>%
  flextable::border_outer()
```

Now, we can use the `plot` and the `likert` function to visualize the survey data.

```{r likert13}
plot(likert(survey[,3:12]), ordered = F, wrap= 60)
```

To save this plot, you can use the `save_plot` function from the `cowplot` package as shown below.

```{r likert15, eval = F}
survey_p1 <- plot(likert(survey[,3:12]), ordered = F, wrap= 60)
# save plot
cowplot::save_plot(here("images", "stu_p1.png"), # where to save the plot
                   survey_p1,        # object to plot
                   base_asp = 1.5,  # ratio of space fro questions vs space for plot
                   base_height = 8) # size! higher for smaller font size

```


An additional and very helpful feature is that the `likert` package enables grouping the data as shown below. The display columns 3 to 8 and use column 1 for grouping.

```{r likert17, fig.height = 10, fig.width = 8, fig.align = "center"}
# create plot
plot(likert(survey[,3:8], grouping = survey[,1]))
```


# Useful statistics

This section introduces some statistical measures or tests that useful when dealing with survey data. We will begin with measures of reliability (Cronbach's $\alpha$), then move on to methods for merging variables (Factor analysis and principle component analysis), and finally to ordinal regression which tests which variables correlate with a certain outcome.

## Evaluating the reliability of questions{-}

### Cronbach's Alpha{-}

Oftentimes several questions in one questionnaire aim to tap into the same cognitive concept or attitude or whatever we are interested in. The answers to these related questions should be internally consistent, i.e. the responses should correlate strongly and positively. 

Cronbach's $\alpha$ [@cronbach1951alpha] is measure of internal consistency or reliability that provides information on how strongly the responses to a set of questions correlate. The formula for Cronbach's $\alpha$ is shown below (N: number of items, $\bar c$: average inter-item co-variance among items, $\bar v$: average variance).

$\alpha = \frac{N*\bar c}{\bar v + (N-1)\bar c}$

If the values for Cronbach's $\alpha$ are low (below .7), then this indicates that the questions are not internally consistent (and do not tap into the same concept) or that the questions are not uni-dimensional (as they should be). 

While Cronbach's $\alpha$ is the most frequently used measures of reliability (probably because it is conceptually simple and can be computed very easily), it underestimates the reliability of a test and overestimates the first factor saturation. This can be a problem is the data is *lumpy*. Thus, various other measures of reliability have been proposed. Also,Cronbach's $\alpha$ assumes that scale items are repeated measurements, an assumption that is often violated.

An alternative reliability measure that takes the amount of variance per item into account and thus performs better when dealing with *lumpy* data (although it is still affected by *lumpiness*) is Guttman's Lambda 6 (G6) [@guttman1945lambda]. In contrast to Cronbach's $\alpha$, G6 is mostly used to evaluate the reliability of individual test items though. This means that it provides information about how well individual questions reflect the concept that they aim to tap into. 

Probably the best measures of reliability are $\omega$ (*omega*) measures. Hierarchical $\omega$ provides more appropriate estimates of the general factor saturation while total $\omega$ is a better estimate of the reliability of the total test compared to both Cronbach's $\alpha$ and G6 [@revelle2009coefficients].

**Calculating Cronbach's alpha in R**

We will now calculate Cronbach's $\alpha$ in R. In a first step, we activate the "psych" package and load as well as inspect the data. 

```{r cron_01, message=FALSE, warning=FALSE}
# load data
surveydata <- base::readRDS(url("https://slcladal.github.io/data/sud.rda", "rb"))
```


```{r cron_02a, echo = F}
surveydata %>%
  as.data.frame() %>%
  flextable() %>%
  flextable::set_table_properties(width = .5, layout = "autofit") %>%
  flextable::theme_zebra() %>%
  flextable::fontsize(size = 12) %>%
  flextable::fontsize(size = 12, part = "header") %>%
  flextable::align_text_col(align = "center") %>%
  flextable::set_caption(caption = "Overview of the survey data.")  %>%
  flextable::border_outer()
```


The inspection of the data shows that the responses of participants represent the rows and that the questions represent columns. The column names show that we have 15 questions and that the first five questions aim to test how outgoing respondents are. To check if the first five questions reliably test "outgoingness" (or "extraversion"), we calculate Cronbach's alpha for these five questions.

Thus, we use the "alpha()" function and provide the questions that tap into the concept we want to assess. In addition to Cronbach's $\alpha$, the "alpha()" function also reports Guttman's lambda_6 which is an alternative measure for reliability. This is an advantage because Cronbach's $\alpha$ underestimates the reliability of a test and overestimates the first factor saturation.

```{r cron_02b, message=FALSE, warning=FALSE}
# calculate cronbach's alpha
Cronbach <- psych::alpha(surveydata[c("Q01_Outgoing",	
                   "Q02_Outgoing",	
                   "Q03_Outgoing",	
                   "Q04_Outgoing",	
                   "Q05_Outgoing")], check.keys=F)
# inspect results
Cronbach
```

The output of the "alpha()" function  is rather extensive and we will only interpret selected output here. 

The value under alpha is Cronbach's $\alpha$ and it should be above 0.7. The values to its left and right are the lower and upper bound of its confidence interval. The values in the column with the header "G6" show how well each question represents the concept it aims to reflect. Low values indicate that the question does not reflect the underlying concept while high values (.7 and higher) indicate that the question captures that concept well (or to an acceptable degree).

### Omega{-}

The omega ($\omega$) coefficient is also a reliability measure of internal consistency. $\omega$ represents an estimate of the general factor saturation of a test that was proposed by McDonald. [@zinbarg2005cronbach] compare McDonald's Omega to Cronbach's $\alpha$ and Revelle's $\beta$. They conclude that omega is the best estimate [@zinbarg2006estimating].

A very handy way to calculate McDonald's $\omega$ is to use the `scaleReliability()` function from the `userfriendlyscience` package (which also provides Cronbach's $\alpha$ and the Greatest Lower Bound (GLB) estimate which is also a very good and innovative measure of reliability) [see also @peters2014alpha]. 

```{r omeg_02, message=FALSE, warning=FALSE}
# activate package
library(ufs)
# extract reliability measures
reliability <- ufs::scaleStructure(surveydata[c("Q01_Outgoing",	
                                                "Q02_Outgoing",	
                                                "Q03_Outgoing",	
                                                "Q04_Outgoing",	
                                                "Q05_Outgoing")])
# inspect results
print(reliability)
```

## Factor analysis{-}

When dealing with many variables it is often the case that several variables are related and represent a common, underlying factor. To find such underlying factors, we can use a factor analysis.

Factor analysis is a method that allows to find commonalities or structure in data. This is particularly useful when dealing with many variables. Factors can be considered hidden latent variables or driving forces that affect or underlie several variables at once.

This becomes particularly apparent when considering socio-demographic variables as behaviors are not only dependent on single variables, e.g., economic status, but on the interaction of several additional variables such as education level, marital status, number of children, etc. All of these variables can be combined into a single factor (or hidden latent variable).

```{r fa_01, message=FALSE, warning=FALSE}
# remove respondent
surveydata <- surveydata %>% 
  dplyr::select(-Respondent)
factoranalysis <- factanal(surveydata, 3, rotation="varimax")
print(factoranalysis, digits=2, cutoff=.2, sort=TRUE)
```

The results of a factor analysis can be visualized so that questions which reflect the same underlying factor are grouped together.

```{r fa_04, message=FALSE, warning=FALSE}
# plot factor 1 by factor 2
load <- factoranalysis$loadings[,1:2]
# set up plot
plot(load, type="n", xlim = c(-1.5, 1.5)) 
# add variable names
text(load,
     # define labels
     labels=names(surveydata),
     # define font size 
     # (smaller than default = values smaller than 1)
     cex=.7)  
```

The plot shows that the questions form groups which indicates that the questions do a rather good job at reflecting the concepts that they aim to tap into. The only problematic question is question 10 (Q10) which aimed to tap into the intelligence of respondents but appears not to correlate strongly with the other questions that aim to extract information about the respondents intelligence. In such cases, it makes sense, to remove a question (in this case Q10) from the survey as it does not appear to reflect what we wanted it to.

## Principle component analysis{-}

Principal component analysis is used when several questions or variables reflect a common factor and they should be combined into a single variable, e.g. during  the statistical analysis of the data. Thus, principal component analysis can be used to collapse different variables (or questions) into one. 

Imagine you have measured lengths of sentences in different ways (in words, syllables, characters, time it takes to pronounce, etc.). You could combine all these different measures of length by applying a PCA to those measures and using the first principal component as a single proxy for all these different measures. 

```{r pca_03, message=FALSE, warning=FALSE}
# entering raw data and extracting PCs  from the correlation matrix
PrincipalComponents <- princomp(surveydata[c("Q01_Outgoing",	
                   "Q02_Outgoing",	
                   "Q03_Outgoing",	
                   "Q04_Outgoing",	
                   "Q05_Outgoing")], cor=TRUE)
summary(PrincipalComponents) # print variance accounted for
```

The output shows that the first component (Comp.1) explains 91.58 percent of the variance. This shows that we only lose 8.42 percent of the variance if we use this component as a proxy for "outgoingness" if we use the collapsed component rather than the five individual items.

```{r pca_05, message=FALSE, warning=FALSE}
loadings(PrincipalComponents) # pc loadings
```

We now check if the five questions that are intended to tap into "outgoingness" represent one (and not more) underlying factors. Do check this, we create a scree plot. 

```{r pca_07, message=FALSE, warning=FALSE}
plot(PrincipalComponents,type="lines") # scree plot
```

The scree plot shown above indicates that we only need a single component to explain the variance as there is a steep decline from the first to the second component. This confirms that the questions that tap into "outgoingness" represent one (and not more) underlying factors.

```{r pca_09, message=FALSE, warning=FALSE}
PrincipalComponents$scores # the principal components
```

You could now replace the five items which tap into "outgoingness" with the single first component shown in the table above.

## Ordinal Regression{-}

Ordinal regression is very similar to multiple linear regression but takes an ordinal dependent variable [@agresti2010analysis]. For this reason, ordinal regression is one of the key methods in analyzing Likert data. 

To see how an ordinal regression is implemented in R, we load and inspect the "ordinaldata" data set. The data set consists of 400 observations of students that were either educated at this school (Internal = 1) or not (Internal = 0). Some of the students have been abroad (Exchange = 1) while other have not (Exchange = 0). In addition, the data set contains the students' final score of a language test (FinalScore) and the dependent variable which the recommendation of a committee for an additional, very prestigious program. The recommendation has three levels ("very likely", "somewhat likely", and "unlikely") and reflects the committees’ assessment of whether the student is likely to succeed in the program.

```{r orr1, message=FALSE, warning=FALSE}
# load data
ordata <- base::readRDS(url("https://slcladal.github.io/data/oda.rda", "rb"))
# inspect data
str(ordata)
```

In a first step, we need to re-level the ordinal variable to represent an ordinal factor (or a progression from "unlikely" over "somewhat likely" to "very likely". And we will also factorize Internal and Exchange to make it easier to interpret the output later on.


```{r orr2, message=FALSE, warning=FALSE}
# relevel data
ordata <- ordata %>%
dplyr::mutate(Recommend = factor(Recommend, 
                           levels=c("unlikely", "somewhat likely", "very likely"),
                           labels=c("unlikely",  "somewhat likely",  "very likely"))) %>%
  dplyr::mutate(Exchange = ifelse(Exchange == 1, "Exchange", "NoExchange")) %>%
  dplyr::mutate(Internal = ifelse(Internal == 1, "Internal", "External"))
```

Now that the dependent variable is re-leveled, we check the distribution of the variable levels by tabulating the data. To get a better understanding of the data we create frequency tables across variables rather than viewing the variables in isolation.

```{r orr3, message=FALSE, warning=FALSE}
## three way cross tabs (xtabs) and flatten the table
ftable(xtabs(~ Exchange + Recommend + Internal, data = ordata))
```

We also check the mean and standard deviation of the final score as final score is a numeric variable and cannot be tabulated (unless we convert it to a factor).

```{r orr4, message=FALSE, warning=FALSE}

summary(ordata$FinalScore); sd(ordata$FinalScore)
```

The lowest score is 1.9 and the highest score is a 4.0 with a mean of approximately 3. Finally, we inspect the distributions graphically.

```{r orr5, message=FALSE, warning=FALSE}
# visualize data
ordata %>%
  ggplot(aes(x = Recommend, y = FinalScore)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) +
  facet_grid(Exchange ~ Internal, margins = TRUE) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
```

We see that we have only few students that have taken part in an exchange program and there are also only few internal students overall. With respect to recommendations, only few students are considered to very likely succeed in the program. We can now start with the modelling by using the "polr" function. To make things easier for us, we will only consider the main effects here as this tutorial only aims to how to implement an ordinal regression but not how it should be done in a proper study - then, the model fitting and diagnostic procedures would have to be performed accurately, of course. 

```{r orr6, message=FALSE, warning=FALSE}
## fit ordered logit model and store results 'm'
m <- polr(Recommend ~ Internal + Exchange + FinalScore, data = ordata, Hess=TRUE)
## view a summary of the model
summary(m)
```

The results show that having studied here at this school increases the chances of receiving a positive recommendation but that having been on an exchange has a negative but insignificant effect on the recommendation. The final score also correlates positively with a positive recommendation but not as much as having studied here.

```{r orr7, message=FALSE, warning=FALSE}
## store table
(ctable <- coef(summary(m)))
```

As the regression report does not provide p-values, we have to calculate them separately (after having calculated them, we add them to the coefficient table).

```{r orr8, message=FALSE, warning=FALSE}
## calculate and store p values
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
## combined table
(ctable <- cbind(ctable, "p value" = p))
```
As predicted, Exchange does not have a significant effect but FinalScore and Internal both correlate significantly with the likelihood of receiving a positive recommendation.

```{r orr9, message=FALSE, warning=FALSE}
# extract profiled confidence intervals
ci <- confint(m)
# calculate odds ratios and combine them with profiled CIs
exp(cbind(OR = coef(m), ci))
```

The odds ratios show that internal students are 2.85 or 285 percent more likely compared to non-internal students to receive positive evaluations and that a 1-point increase in the test score lead to a 1.85 or 185 percent increase in the chances of receiving a positive recommendation. The effect of an exchange is slightly negative but, as we have seen above, not significant.

In a final step, we will visualize the results of the ordinal regression model. To do that, we need to reformat the data and add the predictions.

```{r orr11, message=FALSE, warning=FALSE}
# extract predictions
predictions <- predict(m, data = ordata, type = "prob")
# add predictions to the data
newordata <- cbind(ordata, predictions)
# rename columns
colnames(newordata)[6:7] <- c("somewhat_likely", "very_likely")
# reformat data
newordata <- newordata %>%
  dplyr::select(-Recommend) %>%
  tidyr::gather(Recommendation, Probability, unlikely:very_likely)  %>%
  dplyr::mutate(Recommendation = factor(Recommendation, 
                                        levels = c("unlikely",
                                                   "somewhat_likely",
                                                   "very_likely")))
```

```{r orr11b, message=FALSE, warning=FALSE}
newordata %>%
  as.data.frame() %>%
  head(10) %>%
  flextable() %>%
  flextable::set_table_properties(width = .5, layout = "autofit") %>%
  flextable::theme_zebra() %>%
  flextable::fontsize(size = 12) %>%
  flextable::fontsize(size = 12, part = "header") %>%
  flextable::align_text_col(align = "center") %>%
  flextable::set_caption(caption = "First 10 rows of the newordata.")  %>%
  flextable::border_outer()
```



We can now visualize the predictions of the model.

```{r orr13, message=FALSE, warning=FALSE}
# bar plot
newordata %>%
  ggplot(aes(x = FinalScore, Probability,
             color = Recommendation, 
             group = Recommendation)) + 
  facet_grid(Exchange~Internal) +
  geom_smooth() +  
  # define colors
  scale_fill_manual(values = clrs3) +
  scale_color_manual(values = clrs3) +
  theme_bw()  
```

For more information about regression modeling, model fitting, and model diagnostics, please have a look at the tutorial on [fixed-effects regressions](/regression.html#1_Fixed_Effects_Regression).


# Citation & Session Info {-}

Schweinberger, Martin. 2020. *Questionnaires and Surveys: Analyses with R*. Brisbane: The University of Queensland. url: https://slcladal.github.io/surveys.html (Version 2020.12.11).

```
@manual{schweinberger2020survey,
  author = {Schweinberger, Martin},
  title = {Questionnaires and Surveys: Analyses with R},
  note = {https://slcladal.github.io/survey.html},
  year = {2020},
  organization = "The University of Queensland, Australia. School of Languages and Cultures},
  address = {Brisbane},
  edition = {2020/12/11}
}
```


```{r fin}
sessionInfo()
```


# References{-}

***

[Back to top](#introduction)

[Back to HOME](https://slcladal.github.io/index.html)

***

