knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.0.5     ✓ stringr 1.4.0
## ✓ tidyr   1.1.2     ✓ forcats 0.5.0
## ✓ readr   1.4.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(ggplot2)
library(data.table)
## 
## Attaching package: 'data.table'
## The following object is masked from 'package:purrr':
## 
##     transpose
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
library(dabestr)
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
## 
##     set_names
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## Attaching package: 'dabestr'
## The following object is masked from 'package:rstatix':
## 
##     cohens_d
library(knitr)

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
corstarsl <- function(x,y){
  require(Hmisc)
  x <- as.matrix(x)
  R <- Hmisc::rcorr(x, type = y)$r
  p <- Hmisc::rcorr(x, type = y)$P
  ## define notions for significance levels; spacing is important.
  mystars <- ifelse(p < .001, "***", ifelse(p < .01, "** ", ifelse(p < .05, "* ", " ")))
  ## trunctuate the matrix that holds the correlations to two decimal
  R <- format(round(cbind(rep(-1.11, ncol(x)), R), 2))[,-1]
  ## build a new matrix that includes the correlations with their apropriate stars
  Rnew <- matrix(paste(R, mystars, sep=""), ncol=ncol(x))
  diag(Rnew) <- paste(diag(R), " ", sep="")
  rownames(Rnew) <- colnames(x)
  colnames(Rnew) <- paste(colnames(x), "", sep="")
  ## remove upper triangle
  Rnew <- as.matrix(Rnew)
  Rnew[upper.tri(Rnew, diag = TRUE)] <- ""
  Rnew <- as.data.frame(Rnew)
  ## remove last column and return the matrix (which is now a data frame)
  Rnew <- cbind(Rnew[1:length(Rnew)-1])
  return(Rnew)
}

overview of the dataset

lisa<-read_csv(file = '~/Downloads/BetterUp-TakeHomeExercise-MostafaSalariRad/ebs_takehome_data/jcc_raise.csv') 
lisa %>%
  glimpse()
## Rows: 1,795
## Columns: 9
## $ pid        <chr> "A2VFPLE7TSN1IL", "APNTEGPF7QZQI", "AGR032RZW1C76", "ADBDQ…
## $ sex        <chr> "m", "f", "m", "m", "m", "f", "m", "m", "f", "f", "m", "m"…
## $ age        <dbl> 35, 60, 49, 38, 43, 30, 35, 24, 26, 44, 32, 34, 36, 26, 34…
## $ jcc_1      <dbl> 2, 2, 3, 2, 1, 4, 1, 1, 2, 1, 4, 4, 4, 1, 3, 1, 1, 1, 1, 1…
## $ jcc_2      <dbl> 3, 1, 2, 3, 4, 3, 3, 3, 2, 2, 3, 1, 1, 3, 2, 1, 4, 4, 1, 4…
## $ jcc_3      <dbl> 2, 3, 1, 4, 3, 1, 2, 3, 2, 4, 2, 1, 1, 2, 2, 3, 2, 2, 4, 1…
## $ raise_3mo  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0…
## $ raise_6mo  <dbl> 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0…
## $ raise_12mo <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1…

The dataset contains N=1795 participants, identified with a key variable, pid. We have data on their sex and age. The three variables jss-1, jss-2, jss-3 denote how participants rated the attitude-toward-work vignettes, corresponding to view of work as a job, a career, and a calling respectively. Lastly, there is data on whether participants have received a raise in the past 3, 6, or 9 months (variable names: raise_3mo, raise_6mo, raise_9mo).

Missing values:

lisa %>%
    filter_at(.vars = names(.), any_vars(is.na(.)))
## # A tibble: 5 x 9
##   pid           sex     age jcc_1 jcc_2 jcc_3 raise_3mo raise_6mo raise_12mo
##   <chr>         <chr> <dbl> <dbl> <dbl> <dbl>     <dbl>     <dbl>      <dbl>
## 1 <NA>          f         1    19    NA     4         0         0          1
## 2 ZZZT25WTBJC11 m        30    NA     2     4         0        NA         NA
## 3 WBUT25WTBJC09 f        10     1     2    NA         1         0          1
## 4 VINN25WTBJC00 m        28    NA    NA    NA        NA        NA         NA
## 5 COUS25WTBJB0B <NA>     51     1     2     4         0         0          1

We have 5 rows with missing values. Four are missing at least one value related to the main dependent variable (the vignette ratings) and one is missing a value for sex. We remove these five participants for now but we need to check later to examine why they data is missing.

lisa %>%
    drop_na() -> lisa

Duplicates & invalid values

Identify duplicates and values that are different than the intended type, or out of range.

#pids
lisa %>% 
    group_by(pid) %>% 
    filter(n()>1)
## # A tibble: 2 x 9
## # Groups:   pid [1]
##   pid           sex     age jcc_1 jcc_2 jcc_3 raise_3mo raise_6mo raise_12mo
##   <chr>         <chr> <dbl> <dbl> <dbl> <dbl>     <dbl>     <dbl>      <dbl>
## 1 A5ET25WTBJC0J mm       44     1     2     4         0         0          1
## 2 A5ET25WTBJC0J m        35     1     2     4         0         0          1

There’s one duplicated pid. It’s value for all columns is the same except for sex and age. For age, it’s 44 and 35, and for sex, it’s m or mm. The latter might be a typo and can be fixed but the former is impossible to determine. I’m removing this participant’s data. Again, it’s important to go back to the survey and sample an identify the reason this pid was assigned twice.

lisa %>% 
    group_by(pid) %>% 
    filter(!n()>1) %>% 
    ungroup() -> lisa
#sex
lisa %>% 
    group_by(sex) %>% 
    count()
## # A tibble: 2 x 2
## # Groups:   sex [2]
##   sex       n
##   <chr> <int>
## 1 f       807
## 2 m       981

The values for sex are correctly coded.

#age 
range(lisa$age)
## [1] 18 98

Age range is also between 18-98 which is valid. It’d be interesting to check what people above retirement age or in their 90s are doing. Are they working? How do they construe the concept of work? Do they need to work, or is it more of a hobby? etc.

lisa  %>%
  summarise(mean = mean(age),
            sd = sd(age),
            count = n(),
            se=sd/sqrt(count),
            lower_ci = mean-1.96*se, 
            upper_ci = mean+1.96*se)
## # A tibble: 1 x 6
##    mean    sd count    se lower_ci upper_ci
##   <dbl> <dbl> <int> <dbl>    <dbl>    <dbl>
## 1  37.8  10.2  1788 0.240     37.3     38.3

To summarize, the final cleaned sample consists of N=807 females, and N=981 males. The average participant Age is about M=37.8 years old (95%CI[37.3,38.3]).

Figure below shows the distribution of age for different gender groups. In this sample, females are slightly older than males (38.6 vs. 37.1, t = 3.16, p = 0.0016).

lisa %>% 
    filter(!is.na(sex)) %>% 
    t_test(age~sex)
## # A tibble: 1 x 8
##   .y.   group1 group2    n1    n2 statistic    df       p
## * <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>   <dbl>
## 1 age   f      m        807   981      3.13 1665. 0.00177
lisa %>% 
    group_by(sex) %>% 
    summarise(mean.line=mean(age, na.rm=T)) -> mean.age.by.sex

lisa %>% 
  filter(!is.na(sex)) %>% 
  ggplot(aes(age, col=sex)) +
  geom_freqpoly(position="dodge") +
  theme_minimal() +
geom_vline(data = mean.age.by.sex, aes(xintercept=mean.line, color=sex), linetype=2)

Do most people see work as a job, career, or calling?

The table below shows the descriptive statistics of the vignette ratings. Responses were collected on a 1-4 scale where the scores were coded as follows: 1: ‘Not at all like me’, 2: ‘A little like me’, 3: ‘Somewhat like me’, and 4: Very much like me’.

# make things easier
lisa %>% 
    rename(Job=jcc_1,Career=jcc_2,Calling=jcc_3) -> lisa


lisa %>% 
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>% 
    group_by(vignette) %>% 
    get_summary_stats(value) %>% 
    knitr::kable()
vignette variable n min max median q1 q3 iqr mad mean sd se ci
Calling value 1788 1 4 2 1 3 2 1.483 2.140 1.099 0.026 0.051
Career value 1788 1 4 3 2 3 1 1.483 2.691 0.983 0.023 0.046
Job value 1788 1 4 2 1 3 2 1.483 2.298 1.137 0.027 0.053

Majority of responses fall between the second (a little like me) and third choices (somewhat like me). It’s possible that people avoided the extremes of the scales. But there is variability within this range with respect to vignette type. The highest rating is given to the vignette that views work as a career, the lowest goes to the view that one’s work is a calling, and the view of work as a job falls in between.

lisa %>% 
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>% 
    convert_as_factor(pid, vignette) %>% 
    anova_test(dv = value, wid = pid, within = vignette) #-> res.aov.1
## ANOVA Table (type III tests)
## 
## $ANOVA
##     Effect DFn  DFd      F        p p<.05   ges
## 1 vignette   2 3574 94.355 1.17e-40     * 0.044
## 
## $`Mauchly's Test for Sphericity`
##     Effect   W        p p<.05
## 1 vignette 0.9 1.74e-41     *
## 
## $`Sphericity Corrections`
##     Effect   GGe        DF[GG]    p[GG] p[GG]<.05  HFe        DF[HF]    p[HF]
## 1 vignette 0.909 1.82, 3249.79 3.09e-37         * 0.91 1.82, 3252.92 2.86e-37
##   p[HF]<.05
## 1         *

An analysis of variance to comparing the three ratings showed a significant effect of vignette type (F(1.82,3249.79) = 94.355, p < 0.0001, \(\eta^{2}{g}\) = 0.045). The plot below shows the distribution of response to each vignette.

lisa %>% 
    pivot_longer(cols = c(Job,Career,Calling), names_to = "vignette") %>% 
    group_by(vignette) %>% 
    count(value) %>% 
    mutate(Freq = n/sum(n)) -> descriptive.props
    
lisa %>% 
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>% 
    mutate(value.factor = factor(5-value, labels = c("Very much like me", "Somewhat like me", "A little like me", "Not at all like me"))) %>% 
    convert_as_factor(pid, vignette) %>% 
    ggplot(aes(x = factor(vignette, level=c("Job", "Career", "Calling")), fill=value.factor)) + 
    geom_bar() +
    ylab("count") +
    theme_classic() +
    scale_fill_discrete(name = "response") + 
    theme(axis.text.x.bottom = element_text(size=12), axis.text.y.left = element_text(size=12), legend.text = element_text(size=12))+
    xlab("vignette")

Table below shows the pairwise comparisons between different vignettes. The biggest difference is between Career and Calling followed by the difference between Career and Job. The Job vignette is also slightly rated higher than the Calling vignette.

lisa %>% 
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>% 
    convert_as_factor(pid, vignette) %>% 
    pairwise_t_test(
        value ~ vignette, paired = TRUE,
        p.adjust.method = "bonferroni") %>% 
    knitr::kable()
.y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
value Calling Career 1788 1788 -14.535120 1787 0.000000 0.000 ****
value Calling Job 1788 1788 -3.332316 1787 0.000879 0.003 **
value Career Job 1788 1788 10.369666 1787 0.000000 0.000 ****

Do gender differences exist in the characterization of work as a job, career, or calling?

lisa %>% 
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>% 
    convert_as_factor(pid, vignette) %>% 
    anova_test(dv = value, wid = pid, between = sex, within = vignette) #-> res.aov.2
## ANOVA Table (type III tests)
## 
## $ANOVA
##         Effect DFn  DFd      F        p p<.05      ges
## 1          sex   1 1786  0.040 8.42e-01       2.73e-06
## 2     vignette   2 3572 91.479 1.80e-39     * 4.30e-02
## 3 sex:vignette   2 3572  2.699 6.70e-02       1.00e-03
## 
## $`Mauchly's Test for Sphericity`
##         Effect     W       p p<.05
## 1     vignette 0.901 2.4e-41     *
## 2 sex:vignette 0.901 2.4e-41     *
## 
## $`Sphericity Corrections`
##         Effect  GGe        DF[GG]    p[GG] p[GG]<.05  HFe       DF[HF]   p[HF]
## 1     vignette 0.91 1.82, 3248.77 3.66e-36         * 0.91 1.82, 3251.9 3.4e-36
## 2 sex:vignette 0.91 1.82, 3248.77 7.30e-02           0.91 1.82, 3251.9 7.3e-02
##   p[HF]<.05
## 1         *
## 2

A mixed ANOVA testing the Vignette X Sex interaction showed a main effect of vignette (F(1.82,3248.77) = 91.479, p < 0.0001, \(\eta^{2}{g}\) = 0.043) and no effect of sex (F < 1).

The interaction was also not significant (p < 0.6). But given the small p-value, I looked at the pairwise comparisons (People could object to this, but we’ll address it further down). Results showed that female participants viewed their works as a Calling than male participants. Table below shows these contrasts.

lisa %>% 
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>% 
    convert_as_factor(pid, vignette) %>% 
    group_by(vignette) %>% 
    pairwise_t_test(
        value ~ sex,
        p.adjust.method = "bonferroni") %>% 
    knitr::kable()
vignette .y. group1 group2 n1 n2 p p.signif p.adj p.adj.signif
Calling value f m 807 981 0.0283 * 0.0283 *
Career value f m 807 981 0.3240 ns 0.3240 ns
Job value f m 807 981 0.2850 ns 0.2850 ns

Here we can see the distribution of response split by participant gender:

lisa %>% 
  pivot_longer(cols = c(Job,Career,Calling), names_to = "vignette") %>% 
  group_by(vignette, sex) %>% 
  count(value) %>% 
  mutate(Freq = n/sum(n)) -> descriptive.props


descriptive.props %>% 
    mutate(value.factor = factor(5-value, labels = c("Very much like me", "Somewhat like me", "A little like me", "Not at all like me"))) %>% 
    ggplot(aes(x = sex, y = Freq, fill=value.factor)) + 
    geom_bar(stat="identity", position="dodge") +
    ylab("percentage") +
    theme_classic() +
    scale_fill_discrete(name = "response") + 
    theme(axis.text.x.bottom = element_text(size=12), axis.text.y.left = element_text(size=12), legend.text = element_text(size=12)) + 
    facet_grid(.~vignette)+
    coord_flip()

The plot below shows the bias-corrected and accelerated confidence intervals from 5000 bootstrap resamples for the differences between female and male participants in their ratings of the three vignettes. The gender difference in ratings of the Job and Career vignettes overlaps with 0 (Job: Mdiff=-0.0577 [95CI -0.166; 0.0478], Career: Mdiff=-0.0461 [95CI -0.139; 0.0494]) but not for the Calling vignette (Mdiff=0.115 [95CI 0.0119; 0.218]).

lisa %>% 
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>% 
    mutate(group= paste0(sex,".", vignette)) %>% 
    dabest(x = group, y = value, paired = F, idx = list(c("m.Job","f.Job"), c("m.Career","f.Career"), c("m.Calling","f.Calling"))) %>%  
    mean_diff() -> gender.vignette.comp
plot(gender.vignette.comp)

Is there a relationship between dominant work attitudes and the likelihood of having received a recent raise?

lisa %>% 
  count(raise_3mo, raise_6mo, raise_12mo) %>% 
  kable()
raise_3mo raise_6mo raise_12mo n
0 0 0 695
0 0 1 556
0 1 0 312
1 0 0 225

Checking the distribution of raises shows that nobody has received more than one raise.

df.reverese<-lisa
lisa %>% 
    select(starts_with("raise")) %>% 
    transmute_all(function(x) ifelse(x == 1, deparse(substitute(x)), NA)) %>% 
    unite(col = "raises", na.rm = T) %>% 
    mutate(raises = fct_recode(raises)) %>% 
    bind_cols(lisa,.) -> lisa

# recode and reorder levels
levels(lisa$raises)<-levels(lisa$raises)[c(1,2,4,3)]


lisa %>% 
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>% 
    convert_as_factor(pid, vignette) %>%
    group_by(vignette, raises) %>% 
    ggplot(aes(y = value, x = as.numeric(raises), col=vignette)) +
    theme_classic() + geom_smooth(method="lm") +
    coord_cartesian(xlim = c(1,4.5), ylim = c(1,4))+
    xlab(label = "raises") +
    scale_x_continuous(breaks=c(1:4), labels=c("None", "1 year", "6 months", "3 months")) +
    ylab("attitudes toward work") +
    theme(legend.title = element_blank(), legend.position = c(0.2,0.8), axis.text.x.bottom=element_text(size=12)) -> p.left


lisa %>%
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>%
    convert_as_factor(pid, vignette) %>%
    group_by(vignette, raises) %>%
    get_summary_stats(value) %>%
    ggplot(aes(raises, mean, fill=vignette)) +
    geom_bar(stat = "identity", position = "dodge") +
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se), position=position_dodge(0.9), width=0.2) +
    theme_classic() +
    ylab("attitudes toward work (SE)") +
    scale_x_discrete(labels=c("None", "1 year","6 months","3 months")) +
    theme(legend.position = "none", axis.text.x.bottom=element_text(size=12)) +
    coord_cartesian(ylim=c(1,4)) +
    annotate("text", x=1, y=3.2, label= "*", size=14) -> p.right

multiplot(p.left, p.right, cols = 2)

The left panel shows the linear relationship between the raises and attitudes toward work. The right panel shows that a job view of work is significantly rated higher than a calling view of work in people who did not get a raise, but not others.

lisa %>%     
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>% 
    plyr::dlply(.variables = c("vignette"), function(df)return(cor.test(df$value, as.numeric(df$raises))))
## $Calling
## 
##  Pearson's product-moment correlation
## 
## data:  df$value and as.numeric(df$raises)
## t = -0.11007, df = 1786, p-value = 0.9124
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.04895584  0.04375802
## sample estimates:
##          cor 
## -0.002604506 
## 
## 
## $Career
## 
##  Pearson's product-moment correlation
## 
## data:  df$value and as.numeric(df$raises)
## t = 2.4945, df = 1786, p-value = 0.0127
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.01259937 0.10499265
## sample estimates:
##       cor 
## 0.0589222 
## 
## 
## $Job
## 
##  Pearson's product-moment correlation
## 
## data:  df$value and as.numeric(df$raises)
## t = -2.0274, df = 1786, p-value = 0.04277
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.094065675 -0.001563614
## sample estimates:
##         cor 
## -0.04791738 
## 
## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   vignette
## 1  Calling
## 2   Career
## 3      Job

We can see the correlations between the three vignettes rating and raises below.

lisa %>% 
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>% 
    convert_as_factor(pid, raises) %>% 
    group_by(raises) %>% 
    pairwise_t_test(
        value ~ vignette,
        p.adjust.method = "bonferroni") %>% 
    knitr::kable()
raises .y. group1 group2 n1 n2 p p.signif p.adj p.adj.signif
value Calling Career 695 695 0.00e+00 **** 0.00e+00 ****
value Calling Job 695 695 5.98e-05 **** 1.79e-04 ***
value Career Job 695 695 3.00e-06 **** 9.20e-06 ****
raise_12mo value Calling Career 556 556 0.00e+00 **** 0.00e+00 ****
raise_12mo value Calling Job 556 556 1.61e-01 ns 4.83e-01 ns
raise_12mo value Career Job 556 556 0.00e+00 **** 0.00e+00 ****
raise_6mo value Calling Career 225 225 0.00e+00 **** 0.00e+00 ****
raise_6mo value Calling Job 225 225 4.13e-02 * 1.24e-01 ns
raise_6mo value Career Job 225 225 5.60e-06 **** 1.68e-05 ****
raise_3mo value Calling Career 312 312 0.00e+00 **** 0.00e+00 ****
raise_3mo value Calling Job 312 312 4.17e-01 ns 1.00e+00 ns
raise_3mo value Career Job 312 312 0.00e+00 **** 0.00e+00 ****

The above pairwise comparisons show a pattern: Only in thso who did not receive a raises is work considered significant more as a job than a calling. In those who did receive raises, this difference is not significant.

Below, I created a dichotomous raise variable (None vs. Some) to look at how the two groups differ in attitudes toward work. The difference between None vs. Some are not significant but still in the expected direction.

# create a dual factor
lisa$raises.dual<-lisa$raises
levels(lisa$raises.dual)<-c("None", "Some","Some","Some")

lisa %>% 
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>% 
    anova_test(dv = value, wid = pid, between = raises.dual, within = vignette) #-> res.aov.3
## ANOVA Table (type III tests)
## 
## $ANOVA
##                 Effect DFn  DFd      F        p p<.05      ges
## 1          raises.dual   1 1786  0.006 9.39e-01       4.06e-07
## 2             vignette   2 3572 86.279 2.55e-37     * 4.10e-02
## 3 raises.dual:vignette   2 3572  2.646 7.10e-02       1.00e-03
## 
## $`Mauchly's Test for Sphericity`
##                 Effect   W        p p<.05
## 1             vignette 0.9 1.24e-41     *
## 2 raises.dual:vignette 0.9 1.24e-41     *
## 
## $`Sphericity Corrections`
##                 Effect   GGe        DF[GG]    p[GG] p[GG]<.05  HFe
## 1             vignette 0.909 1.82, 3246.81 3.47e-34         * 0.91
## 2 raises.dual:vignette 0.909 1.82, 3246.81 7.60e-02           0.91
##          DF[HF]    p[HF] p[HF]<.05
## 1 1.82, 3249.94 3.24e-34         *
## 2 1.82, 3249.94 7.60e-02
lisa %>% 
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>% 
    mutate(group= paste0(raises.dual,".", vignette)) %>% 
    dabest(x = group, y = value, paired = F, idx = list(c("None.Job","Some.Job"), c("None.Career","Some.Career"), c("None.Calling","Some.Calling"))) %>% 
    mean_diff() -> raises.vignette.comp


lisa %>% 
    pivot_longer(cols = c(Job, Career, Calling), names_to = "vignette") %>% 
    convert_as_factor(pid, vignette) %>% 
    group_by(vignette) %>% 
    pairwise_t_test(
        value ~ raises.dual,
        p.adjust.method = "bonferroni") %>% 
    knitr::kable()
vignette .y. group1 group2 n1 n2 p p.signif p.adj p.adj.signif
Calling value None Some 695 1093 0.6410 ns 0.6410 ns
Career value None Some 695 1093 0.0751 ns 0.0751 ns
Job value None Some 695 1093 0.0558 ns 0.0558 ns 
plot(raises.vignette.comp)