Dataset

To demonstrate the Likert analysis workflow, I’ll use the MACH-IV dataset, which is a freely available download. It draws on answers from Christie, R. & Geis, F. (1970) “Studies in Machiavellianism”. NY: Academic Press.), which is a tool that is used to measure Machiavellianism personality traits. The test consists of twenty Likert items. Each item is a statement that you must indicate how accurate it would be when applied to you. The data are presented in wide form with each row as a single observation.

• Answers: 1 = strongly disagree, 2 = disagree, 3 = neutral, 4 = agree, 5 = strongly agree (0 = they skipped). For the puposes of this analysis, we will consider the responses to be ordinal.
• Age: Entered into a text field.
• Gender: Chosen from a drop down list (1 = male, 2 = female, 3 = other; 0 = no answer)
• Questions:
  • Q1. Never tell anyone the real reason you did something unless it is useful to do so.
  • Q2. The best way to handle people is to tell them what they want to hear.
  • Q3. One should take action only when sure it is morally right.
  • Q4. Most people are basically good and kind.
  • Q5. It is safest to assume that all people have a vicious streak and it will come out when they are given a chance.
  • Q6. Honesty is the best policy in all cases.
  • Q7. There is no excuse for lying to someone else.
  • Q8. Generally speaking, people won’t work hard unless they’re forced to do so.
  • Q9. All in all, it is better to be humble and honest than to be important and dishonest.
  • Q10. When you ask someone to do something for you, it is best to give the real reasons for wanting it rather than giving reasons which carry more weight.
  • Q11. Most people who get ahead in the world lead clean, moral lives.
  • Q12. Anyone who completely trusts anyone else is asking for trouble.
  • Q13. The biggest difference between most criminals and other people is that the criminals are stupid enough to get caught.
  • Q14. Most people are brave.
  • Q15. It is wise to flatter important people.
  • Q16. It is possible to be good in all respects.
  • Q17. P.T. Barnum was wrong when he said that there’s a sucker born every minute.
  • Q18. It is hard to get ahead without cutting corners here and there.
  • Q19. People suffering from incurable diseases should have the choice of being put painlessly to death.
  • Q20. Most people forget more easily the death of their parents than the loss of their property.

Perhaps our research question will be to determine whether men and women demonstrate a signficant difference in Machiavellianism personality traits. We will examine each question individually, and then evaluate total score, treating that measure as interval data. Let’s start by downloading the dataset:

url <- "http://personality-testing.info/_rawdata/MACH2.zip"  # web address of dataset

if (!file.exists("./MACH2.zip")) {
    download(url, dest = "MACH2.zip", mode = "wb")  # download zip file
    unzip("MACH2.zip", exdir = "./")  # unzip datafiles into a folder /MACH2
}

ds <- "/MACH2/data.csv"  # the data set of interest
cb <- "/MACH2/codebook.txt"  # the code to the questions

dt <- data.frame(read_csv(paste(wd, ds, sep = "")))  # import data into a data.table

dtQ <- data.table(read_lines(paste(wd, cb, sep = ""))[3:22])  # reading the lines of the questions from the text file

dtQ[, `:=`(c("question", "question_text"), tstrsplit(V1, "[0-9]. ", fixed = TRUE))]  # splits the 'QX. from the Questions.' at the first '.' after a digit.

colnames(dt)[1:20] <- dtQ$question_text[1:20]  # colnames to full questions

posQ <- colnames(dt)[c(1:2, 5, 8, 12:13, 15, 17:20)]  # postive questions
negQ <- colnames(dt)[c(3:4, 6:7, 9:11, 14, 16)]  # negative questions

Explore data and remove outliers

Before we begin the analysis, let’s explore the data and remove outliers in the demographic and total score columns. The scores reflect a minimum of \(20\) and maximum of \(100\) points. When we assign gender as a factor with two levels (male and female), those respondents that did not answer will be removed when we remove NAs. Obviously there are issues with age, with a range of -9 to 999999. We will limit the respondents from 10 to 100 years old. In a similar manner, the range for elapsed time is very large. We will limit seconds_elapsed from 60 seconds (3 seconds per question) to 600 (30 seconds per question). Finally, once we limit respondents to those who answered every question, create factors and levels for each question, most of the outlier questions will be removed.

First, let’s take a look at gender first:

Number of Respondents By Gender

dt[, 1:20] <- lapply(dt[, 1:20], factor, levels = c(1, 2, 3, 4, 5), labels = c("strongly disagree", 
    "disagree", "neutral", "agree", "strongly agree"))  # change columns to factors

dt$gender <- factor(dt$gender, levels = c(1, 2), labels = c("male", "female"))  # change to factor

plot_ly(data = as.data.frame(dt), x = ~gender)

There are quite a few more male respondents than female. Now let’s take a look at age:

Distribution of Age

dt <- dt[dt$age < 110 & dt$age >= 10, ]  # remove observatons with ages outside of 18 - 100. 

dt <- na.omit(dt)  # remove respondents with NAs

plot_ly(data = as.data.frame(dt), x = ~age)  # histogram of age

Most of the respondents are fairly young. Next, lets take a look at seconds_elapsed:

Distrubtion of Seconds Elapsed

dt <- dt[dt$seconds_elapsed < 1200 & dt$seconds_elapsed >= 60, ]  # remove observatons outside of 60 - 1200 seconds. 

plot_ly(data = as.data.frame(dt), x = ~seconds_elapsed)  # histogram of seconds_elapsed

The vast majority of repondents spend very less than 4 minutes on then questionnaire. Finally, let’s take a look at the distribution of scores:

Distrubtion of Scores (Machiavellianism Index)

plot_ly(data = as.data.frame(dt), x = ~score)  # histogram of age

The total scores are normally distributed. The scores are a composite of all the answers. The range is 20-100. The test

pt <- tabular(gender ~ (age + seconds_elapsed + score) * (mean + sd), data = dt)

pander(pt)

Descriptive Statistics

According to Stevens, descriptive statistics for ordinal data should be median. But first, one simple method to visualize the percentages by gender is to use the likert function in the likert package. These figures allow you to group by a factor, and provides color coordinated horizontal bar charts to compare the percent of answers.

dt_pos <- dt[posQ]  # select only the 11 positive questions according to the codetext

dt_pos$gender <- dt$gender  # create a standalone data.table for this analysis.

mach_l1 <- likert::likert(dt_pos[, 1:4, drop = FALSE], grouping = dt_pos$gender)  # view the first 4 (1:4 questions) and group them by gender
likert.bar.plot(mach_l1)

Similarly, you can look at the fiirst 4 negative questions. These visual graphics take up quite a bit of room.

dt_neg <- dt[negQ]
dt_neg$gender <- dt$gender
mach_l2 <- likert::likert(dt_neg[, 1:4, drop = FALSE], grouping = dt_neg$gender)
likert.bar.plot(mach_l2)

In the following examples, I’ve selected the first 2 questions from males to show the summary tables that sjtPlot provides with the sjt.frq function.

library(arsenal)
# Controls used for the arsenal package / demographics tables
mycontrols <- tableby.control(test = TRUE, total = TRUE, digits = 1, numeric.test = "anova", 
    cat.test = "chisq", ordered.test = "trend", numeric.stats = c("N", "meansd", 
        "medianrange", "q1q3", "iqr"), cat.stats = c("countpct"), stats.labels = list(N = "Number", 
        meansd = "Mean (SD)", medianrange = "Median (Range)", q1q3 = "Q1, Q3", iqr = "IQR"))

dtt <- dt
colnames(dtt)[c(1:20, 22)] <- c("Q1", "Q2", "Q3", "Q4", "Q5", "Q6", "Q7", "Q8", "Q9", 
    "Q10", "Q11", "Q12", "Q13", "Q14", "Q15", "Q16", "Q17", "Q18", "Q19", "Q20", 
    "Sex")
tab1 <- tableby(Sex ~ Q1 + Q2, data = dtt, control = mycontrols)
t1 <- "Q1 = Never tell anyone the real reason you did something unless it is useful to do so. Q2 = The best way to handle people is to tell them what they want to hear. N = Number, Q1 = First Quartile, Q3 = Third Quartile, SD = Standard Deviation."
summary(tab1, text = F, pfootnote = T, title = t1)

1. Pearson’s Chi-squared test

rm(tab1, td1, t1)

# sjt.frq(dt[dt$gender == 'male',c(1:2)], emph.md = TRUE, show.summary = FALSE,
# emph.quart = TRUE, no.output = FALSE) # the sjtPlot provides quite a few
# options! emph.md emphasizes the median, emph.quart draws a line are the
# quartiles.  If you show.summary = TRUE, it provides mean, sigma and other
# values in the bottom of the table.

sjt.grpmean(var.cnt = score, var.grp = gender, digits = 1, no.output = FALSE)

## Error in sjt.grpmean(var.cnt = score, var.grp = gender, digits = 1, no.output = FALSE): could not find function "sjt.grpmean"

wt <- wilcox.test(as.numeric(dt_wc[, 1]) ~ as.factor(gender), data = dt_wc, alternative = "two.sided", 
    exact = FALSE)

print(wt)

Wilcoxon rank sum test with continuity correction

data: as.numeric(dt_wc[, 1]) by as.factor(gender) W = 18789658, p-value < 2.2e-16 alternative hypothesis: true location shift is not equal to 0

The Kruskal Wallis test is performed in R with the kruskal.test function, which uses multivariate factors, such as Education & Certification Levels (Associates, Bachelors, Masters, or Doctorates. The current study does not contain multivariate data. I randomly added unknown gender to the categories, such that it is now multivariate.

dt_kw <- dt

dt_kw$gender <- factor(dt_kw$gender, levels = c("male", "female", "unknown"))  # create a third category
dt_kw$gender[sample(nrow(dt_kw), 2000)] <- "unknown"  # label

dt_kw$`Q1. Never tell anyone the real reason you did something unless it is useful to do so.` <- as.numeric(dt_kw$`Q1. Never tell anyone the real reason you did something unless it is useful to do so.`)

kable(Summarize(`Q1. Never tell anyone the real reason you did something unless it is useful to do so.` ~ 
    gender, data = dt_kw, digit = 3)) %>% kable_styling(bootstrap_options = "basic", 
    full_width = F)

kt <- kruskal.test(as.numeric(dt_kw[, 1]) ~ as.factor(gender), data = dt_kw)

print(kt)

Kruskal-Wallis rank sum test

data: as.numeric(dt_kw[, 1]) by as.factor(gender) Kruskal-Wallis chi-squared = 272.66, df = 2, p-value < 2.2e-16

The One-Way ANOVA is a multivariate analysis only used when ordinal data (Likert data) are transformed to interval data. The data, however, still have to meet the assumptions of the test (normality, etc.). For example, the ordinal data can be ordered, and then a table is produced and the data converted to probability or frequency. For this analysis we used the score values and

For the interval data, such as scores, we can use the HH packages to conduct an ANOVA. The hov and hovPlot functions indicate:

Oneway analysis of variance makes the assumption that the variances of the groups are equal. Brown and Forsyth, 1974 present the recommended test of this assumption. The Brown and Forsyth test statistic is the F statistic resulting from an ordinary one-way analysis of variance on the absolute deviations from the median.

The “trellis” object with three panels containing boxplots for each group: The observed data “y”, the data with the median subtracted “y-med(y)”, and the absolute deviations from the median “abs(y-med(y))” The Brown and Forsyth test statistic is the F statistic resulting from an ordinary one-way analysis of variance on the data points in the third panel.

library(HH)

dt_gen <- dt

dt_gen$gender <- factor(dt_gen$gender, levels = c("male", "female", "unknown"))  # create a third category
dt_gen$gender[sample(nrow(dt_gen), 2000)] <- "unknown"  # label

hovPlot(as.numeric(score) ~ as.factor(gender), data = dt_gen)  # example test of the assumtions

library(texreg)
hv <- hov(as.numeric(score) ~ as.factor(gender), data = dt_gen)  # example parametric test
hv <- hov(as.numeric(score) ~ as.factor(gender), data = dt)
screenreg(hv)

## Error in (function (classes, fdef, mtable) : unable to find an inherited method for function 'extract' for signature '"htest"'

set.seed(101)
fac_num <- function(x) {
    if (x == "strongly disagree") {
        x = 1
    }
    if (x == "disagree") {
        x = 2
    }
    if (x == "neutral") {
        x = 3
    }
    if (x == "agree") {
        x = 4
    }
    if (x == "strongly agree") {
        x = 5
    }
}
dt_gen_num <- sapply(dt_gen[, c(1:20)], function(x) as.integer(factor(as.character(x), 
    levels = c("strongly disagree", "disagree", "neutral", "agree", "strongly agree"))))


dt_gen_num <- cbind(dt_gen_num, dt_gen[, c(21:24)])
colnames(dt_gen_num) <- colnames(dt_gen)
mat <- dt_gen_num[, c(1:20)]

myPCA <- prcomp(mat, scale. = F, center = F)
plot(myPCA)

# myPCA$rotation # loadings myPCA$x # scores
biplot(myPCA, scale = 0)

# Perform SVD
mySVD <- svd(mat)
# mySVD # the diagonal of Sigma mySVD$d is given as a vector
sigma <- matrix(0, 20, 20)  # we have 4 PCs, no need for a 5th column
diag(sigma) <- mySVD$d  # sigma is now our true sigma matrix


# Compare PCA scores with the SVD's U*Sigma
theoreticalScores <- mySVD$u %*% sigma
all(round(myPCA$x, 5) == round(theoreticalScores, 5))  # TRUE

[1] TRUE

# Compare PCA loadings with the SVD's V
all(round(myPCA$rotation, 5) == round(mySVD$v, 5))  # TRUE

[1] TRUE

# Show that mat == U*Sigma*t(V)
recoverMatSVD <- theoreticalScores %*% t(mySVD$v)
all(round(mat, 5) == round(recoverMatSVD, 5))  # TRUE

[1] TRUE

# Show that mat == scores*t(loadings)
recoverMatPCA <- myPCA$x %*% t(myPCA$rotation)
all(round(mat, 5) == round(recoverMatPCA, 5))  # TRUE

[1] TRUE

plot(myPCA$x)

library(survey)

Finally, regression techniques include Ordered logistic regression, Multinomial logistic regression, or we could collapse the levels of the Dependent variable into two levels and run binary logistic regression. Regression models are ways to develop predictive models based on the questions.

Until next time…