Employee turnover in supermarkets: statistical consulting for an undergraduate thesis

From the Google Forms export to the thesis tables, with every decision explained.

Fictional consulting case for an undergraduate thesis: from the survey to the thesis tables, with the whole analysis in R.

Fictional data. The student, the supermarkets and the 230 responses do not exist. The data were simulated in R with programmed patterns, for demonstration purposes, and were deliberately given problems that are common in online surveys. The results show the working method; they are not a finding about real supermarkets.

Download the raw responses (CSV) · Download the question dictionary (CSV)

The questionnaire, the CSV files and the comments in the R code are in Portuguese, the case’s original language.

library(tidyverse)
library(gtsummary)
library(patchwork)
library(knitr)

# A análise roda uma vez, em scripts/analise.R, e grava os resultados que as
# versões em português, espanhol e inglês desta página leem.
caso <- "../../../projetos/caso-tcc-rotatividade"
res <- readRDS(file.path(caso, "resultados/analise.rds"))
source(file.path(caso, "scripts/rotulos.R"))
source(file.path(caso, "scripts/apresentacao.R"))
tx <- rotulos("en")
preparar(tx)

# Os trechos de código abaixo são os de scripts/analise.R e scripts/apresentacao.R.
read_chunk(file.path(caso, "scripts/analise.R"))
read_chunk(file.path(caso, "scripts/apresentacao.R"))

b <- setNames(res$resultado$beta, res$resultado$termo)
pv <- setNames(res$resultado$p, res$resultado$termo)
alfa <- setNames(res$conf$alfa, res$conf$bloco)
teste_f <- function(h) paste0("F(", h[["gl1"]], ", ", h[["gl2"]], ") = ", tx$num(h[["F"]]), ", ", tx$p_txt(h[["p"]]))

The request

A Business Administration student was finishing an undergraduate thesis titled Leadership, recognition and turnover intention among supermarket employees in central-western Minas Gerais. He ran a Google Forms survey in stores across the region: profile questions plus 20 statements rated from 1 (strongly disagree) to 5 (strongly agree), in five blocks.

Block	Items	Sample statement
Leadership	5	“My supervisor treats me with respect.”
Recognition	4	“I get praise for my work.”
Perceived pay	4	“My pay is fair for the work I do.”
Work-life balance	4	“I can plan my days off in advance.”
Turnover intention	3	“I plan to look for another job in the next 12 months.”

The project had four hypotheses:

• H1. Better-rated leadership reduces the intention to leave.
• H2. More recognition reduces the intention to leave.
• H3. Better-rated pay reduces the intention to leave, with a smaller effect than H1 and H2.
• H4. Better work-life balance reduces the intention to leave.

He had the spreadsheet exported from Forms and his advisor’s instruction to “do the statistics part.” He needed to know what to run, in what order, how to check whether he could trust the results and how to report everything in the thesis.

What I did

1. Cleaning

bruto <- read_csv("dados/respostas_brutas.csv", show_col_types = FALSE)
dic <- read_csv("dados/dicionario.csv", show_col_types = FALSE)

# O Forms exporta o texto inteiro da pergunta como nome da coluna;
# o dicionário troca cada um por um código curto (LID1, REC2...).
names(bruto) <- dic$codigo[match(names(bruto), dic$coluna)]
itens <- filter(dic, construto != "Perfil") |>
  mutate(enunciado = str_match(coluna, "\\[(.*)\\]$")[, 2])
cod <- itens$codigo
blocos <- split(cod, substr(cod, 1, 3))

brancos <- function(df) rowSums(is.na(select(df, all_of(cod))))
iguais <- function(df) apply(select(df, all_of(cod)), 1, sd, na.rm = TRUE) == 0

passo1 <- filter(bruto, aceite == "Sim")
passo2 <- filter(passo1, brancos(passo1) <= 4)      # até 20% dos 20 itens em branco
passo3 <- filter(passo2, !iguais(passo2))           # mesma resposta nos 20 itens

excluidas <- c(sem_aceite = nrow(bruto) - nrow(passo1), brancos = nrow(passo1) - nrow(passo2),
               iguais = nrow(passo2) - nrow(passo3))
celulas_brancas <- sum(is.na(select(passo3, all_of(cod))))

Of the 230 responses received, 212 went on to the analysis. Each exclusion has a reason that can be written into the thesis methods section:

Table 1: From the Forms export to the valid responses

Step	Excluded	Remaining
Responses received		230
No consent to participate	6	224
More than 20% of items blank	7	217
Same answer on every item	5	212
Valid responses		212

Someone who gives the same answer to everything, including the reversed leadership item (“My supervisor rarely gives me feedback”), is saying at the same time that the supervisor does and does not give feedback: a sign they did not read the statements.

The valid responses still had 87 scattered blank items (2.1% of cells). The rule was simple and is written into the report: each block’s score is the mean of the items answered, as long as no more than 1 item in the block is blank. Two people left 2 Work-life balance items blank and have no score for that block.

Age group was a free-text field, and everyone wrote it their own way. A text rule merged the variants (the middle column shows the original data, in Portuguese):

padronizar_faixa <- function(x) {
  case_when(str_detect(x, "18") ~ "18 a 24",
            str_detect(x, "25") ~ "25 a 34",
            str_detect(x, "35") ~ "35 a 44",
            str_detect(x, "45") ~ "45 ou mais")
}
grafias <- passo3 |>
  mutate(padronizada = padronizar_faixa(faixa_etaria)) |>
  group_by(padronizada) |>
  summarise(escritas = list(sort(unique(faixa_etaria))), pessoas = n())

# Base final: faixa padronizada, item invertido recodificado (1 vira 5, 2 vira 4...).
d <- passo3 |>
  mutate(faixa_etaria = padronizar_faixa(faixa_etaria),
         LID3 = 6L - LID3,
         turno_grupo = if_else(turno == "Noite", "Noturno", "Diurno"),
         tempo_casa = factor(tempo_casa, c("Até 1 ano", "De 1 a 3 anos", "De 3 a 5 anos", "Mais de 5 anos")))
stopifnot(!anyNA(d$faixa_etaria))

Table 2: Age group: variants found and final category

Final category	As written	People
18 to 24	“18-24”, “18 - 24 anos”, “18 a 24”, “18 a 24 anos”, “de 18 a 24”	58
25 to 34	“25-34”, “25 a 34”, “25 á 34”, “25 a 34 anos”, “De 25 a 34”	83
35 to 44	“35-44”, “35 a 44”, “35 A 44”, “35 a 44 anos”	48
45 or older	“45 anos ou mais”, “45 ou mais”, “45+”, “Mais de 45”	23

2. Sample profile

tabela_perfil <- function(res, tx) {
  res$d |>
    transmute(sexo = traduzir(sexo, tx$sexo),
              faixa = traduzir(faixa_etaria, tx$faixa),
              escolaridade = traduzir(escolaridade, tx$escolaridade, ordem = TRUE),
              funcao = traduzir(funcao, tx$funcao),
              turno = traduzir(turno, tx$turno, ordem = TRUE),
              tempo = traduzir(tempo_casa, tx$tempo, ordem = TRUE)) |>
    setNames(tx$perfil[c("sexo", "faixa", "escolaridade", "funcao", "turno", "tempo")]) |>
    tbl_summary(statistic = all_categorical() ~ "{n} ({p}%)") |>
    modify_header(label = tx$perfil[["titulo"]], stat_0 = tx$perfil[["n"]]) |>
    modify_footnote(everything() ~ NA) |>
    as_kable()
}

Table 3: Profile of the valid responses

Characteristic	n (%)
Sex
Female	112 (53%)
Male	94 (44%)
Prefer not to say	6 (2.8%)
Age group
18 to 24	58 (27%)
25 to 34	83 (39%)
35 to 44	48 (23%)
45 or older	23 (11%)
Education
Elementary school	26 (12%)
High school	130 (61%)
Some college	31 (15%)
College degree	25 (12%)
Role
Butcher and bakery	33 (16%)
Customer service	33 (16%)
Warehouse	33 (16%)
Cashier	55 (26%)
Stock clerk	58 (27%)
Shift
Morning	74 (35%)
Afternoon	76 (36%)
Night	62 (29%)
Tenure
Up to 1 year	67 (32%)
1 to 3 years	68 (32%)
3 to 5 years	39 (18%)
More than 5 years	38 (18%)

3. Scale reliability

Before adding statements up into a score, you have to show that they measure the same thing. I used two indicators: Cronbach’s alpha, which the advisor and the thesis committee expect to see, and McDonald’s omega, which is better suited when the statements do not all carry the same weight in the construct. Above .70 is the usual threshold.

confiab <- function(cols) {
  x <- select(d, all_of(cols))
  c(alfa = suppressWarnings(psych::alpha(x)$total$raw_alpha),
    omega = suppressMessages(suppressWarnings(psych::omega(x, nfactors = 1, plot = FALSE)$omega.tot)))
}
rem_sem4 <- c("REM1", "REM2", "REM3")
conf <- bind_rows(
  lapply(names(blocos), function(b) c(bloco = b, itens = length(blocos[[b]]), confiab(blocos[[b]]))),
  c(bloco = "REM_sem4", itens = 3, confiab(rem_sem4))
) |>
  mutate(across(c(itens, alfa, omega), as.numeric)) |>
  arrange(match(bloco, c("LID", "REC", "REM", "REM_sem4", "EQU", "INT")))

# Correlação de cada item de Remuneração com o total dos outros três.
r_rem <- psych::alpha(select(d, all_of(blocos$REM)))$item.stats$r.drop
names(r_rem) <- blocos$REM

Table 4: Cronbach’s alpha and omega for each block

Block	Items	Alpha	Omega
Leadership	5	.83	.83
Recognition	4	.85	.85
Perceived pay	4	.68	.73
Perceived pay, without REM4	3	.82	.83
Work-life balance	4	.81	.81
Turnover intention	3	.84	.84

Four blocks passed comfortably. Perceived pay came in at .68 because of one item: “I receive my transportation allowance on time” (REM4) barely moves with the other three (a correlation of .09 with their total, against .58 or more for the rest). That makes sense: it measures whether a benefit arrives on time, not whether the pay feels fair. Without it, alpha rises to .82. REM4 was dropped from the Perceived pay score, and the thesis records that decision and the reason for it.

4. Exploratory factor analysis

Factor analysis checks whether the 20 statements group into the five blocks the questionnaire intended. Because the answers are ordinal (1 to 5), I used polychoric correlation, which treats each answer as a band on an underlying continuous scale and does not assume that the distance between “1” and “2” equals the distance between “4” and “5”.

x_itens <- select(d, all_of(cod))
pc <- suppressWarnings(polychoric(x_itens)$rho)
kmo <- KMO(pc)$MSA
bartlett <- cortest.bartlett(pc, n = nrow(d))
fe <- suppressWarnings(fa(x_itens, nfactors = 5, cor = "poly", rotate = "oblimin", fm = "minres"))

cargas <- unclass(fe$loadings)
# Dá a cada fator o código do bloco cujos itens carregam mais nele.
dominante <- sapply(seq_len(ncol(cargas)), function(f) {
  names(which.max(tapply(abs(cargas[, f]), substr(rownames(cargas), 1, 3), mean)))
})
colnames(cargas) <- dominante
cargas <- cargas[, c("LID", "REC", "REM", "EQU", "INT")]
carga_rem4 <- max(abs(cargas["REM4", ]))

The KMO was .84 (above .80 is considered very good) and Bartlett’s test gave χ²(190) = 2,198.0, p < .001: the correlations between the statements are strong enough for factor analysis. With five factors and an oblique rotation (the constructs are allowed to correlate), each statement loaded on its own block, with one exception we already knew about: REM4 did not load on any factor (highest loading 0.14), which confirms the reliability decision.

Table 5: Factor loadings (polychoric correlation, oblimin rotation); loadings below 0.30 omitted

Item	Lead.	Recog.	Pay	Bal.	Int.
LID1	0.76
LID2	0.63
LID3	0.84
LID4	0.66
LID5	0.59
REC1		0.81
REC2		0.70
REC3		0.79
REC4		0.73
REM1			0.83
REM2			0.83
REM3			0.77
REM4
EQU1				0.73
EQU2				0.86
EQU3				0.62
EQU4				0.72
INT1					0.74
INT2					0.83
INT3					0.79

5. Scores and group comparisons

escore <- function(cols) {
  m <- as.matrix(select(d, all_of(cols)))
  if_else(rowSums(is.na(m)) <= 1, rowMeans(m, na.rm = TRUE), NA_real_)
}
d <- d |>
  mutate(LID = escore(blocos$LID), REC = escore(blocos$REC), REM = escore(rem_sem4),
         EQU = escore(blocos$EQU), INT = escore(blocos$INT))

Likert scale scores tend to be skewed and full of ties, so I compared the intention to leave across groups with non-parametric tests, always reporting the effect size, which tells you whether a difference is large or merely statistically detectable.

mw <- wilcox_test(d, INT ~ turno_grupo)
mw_r <- wilcox_effsize(d, INT ~ turno_grupo)
kw <- kruskal_test(d, INT ~ tempo_casa)
kw_e <- kruskal_effsize(d, INT ~ tempo_casa)
dunn <- dunn_test(d, INT ~ tempo_casa, p.adjust.method = "holm")
dunn_novato <- filter(dunn, group1 == "Até 1 ano", p.adj < 0.05)

resumo_grupos <- function(var) {
  d |> group_by(grupo = as.character(.data[[var]])) |>
    summarise(n = n(), mediana = median(INT), q1 = quantile(INT, .25), q3 = quantile(INT, .75))
}
grupos <- bind_rows(resumo_grupos("turno_grupo"), resumo_grupos("tempo_casa"))

• Shift. Night-shift workers have a higher intention to leave (Mann-Whitney, p < .001, r = .25, a small effect).
• Tenure. The groups differ (Kruskal-Wallis, H(3) = 10.1, p = .018, ε² = .034, a small effect). In the pairwise comparisons (Dunn, with Holm correction), the difference that holds up is between “up to 1 year” and “1 to 3 years”; the other comparisons were not significant.

Table 6: Turnover intention (score from 1 to 5) by shift and by tenure

Group	n	Median	Interquartile range
Day	150	2.67	2.33 to 3.67
Night	62	3.58	2.67 to 4.00
Up to 1 year	67	3.33	2.67 to 4.00
1 to 3 years	68	2.67	2.00 to 3.38
3 to 5 years	39	3.00	2.42 to 3.83
More than 5 years	38	3.00	2.08 to 3.67

6. Correlations

ct <- corr.test(select(d, LID, REC, REM, EQU, INT), method = "spearman", adjust = "none")
rho <- ct$r
p_int <- ct$p["INT", c("LID", "REC", "REM", "EQU")]

Spearman correlation (appropriate for ordinal scores) shows all four blocks associated with the intention to leave, all in the expected direction (p = .002 or lower). The blocks are also correlated with one another, which matters for the next step.

Table 7: Spearman correlations between scores (N ranges from 210 to 212 depending on blanks)

Block	Lead.	Recog.	Pay	Bal.	Int.
Leadership
Recognition	.48
Perceived pay	.13	.27
Work-life balance	.34	.35	.14
Turnover intention	−.49	−.51	−.31	−.22

7. Multiple regression

Correlation looks at one factor at a time. To test the hypotheses, the question is different: how much does each factor matter when the others are in the same model? I used multiple linear regression with the four scores standardized (so the effects are comparable) and three controls: night shift, up to 1 year of tenure and sex.

dm <- d |>
  transmute(INT = INT, LID = LID, REC = REC, REM = REM, EQU = EQU,
            noturno = as.numeric(turno == "Noite"),
            novato = as.numeric(tempo_casa == "Até 1 ano"),
            feminino = as.numeric(sexo == "Feminino")) |>
  drop_na() |>
  mutate(across(c(INT, LID, REC, REM, EQU), \(v) as.numeric(scale(v))))
modelo <- lm(INT ~ LID + REC + REM + EQU + noturno + novato + feminino, data = dm)

# Pressupostos
vifs <- vif(modelo)
bp <- bptest(modelo)
sw <- shapiro.test(residuals(modelo))
cook_max <- max(cooks.distance(modelo))

# Erros-padrão robustos (HC3): protegem os testes se a variância dos resíduos não for constante.
v_hc3 <- vcovHC(modelo, type = "HC3")
rob <- coeftest(modelo, vcov. = v_hc3)
ic <- coefci(modelo, vcov. = v_hc3)
resultado <- tibble(termo = rownames(rob), beta = rob[, 1], ep = rob[, 2], p = rob[, 4],
                    li = ic[, 1], ls = ic[, 2]) |>
  filter(termo != "(Intercept)")
ajuste <- summary(modelo)
# H3 diz que o efeito da remuneração é menor: testa se os coeficientes são iguais.
igual_lid <- linearHypothesis(modelo, "REM = LID", vcov. = v_hc3)
igual_rec <- linearHypothesis(modelo, "REM = REC", vcov. = v_hc3)
teste <- function(h) c(F = h$F[2], gl1 = h$Df[2], gl2 = h$Res.Df[2], p = h$`Pr(>F)`[2])
menor_confirmado <- igual_lid$`Pr(>F)`[2] < 0.05 && igual_rec$`Pr(>F)`[2] < 0.05
# O texto das páginas foi escrito para o caso em que a diferença não se confirma.
stopifnot(!menor_confirmado)

# H3 tem duas partes: reduz (testada pelo β) e com efeito menor (testada acima).
hipoteses <- tibble(termo = c("LID", "REC", "REM", "EQU")) |>
  left_join(resultado, by = "termo") |>
  mutate(decisao = case_when(!(beta < 0 & p < 0.05) ~ "nao",
                             termo == "REM" & !menor_confirmado ~ "parcial",
                             TRUE ~ "sim"))

Before reading the coefficients, I checked the assumptions:

• Multicollinearity: highest VIF of 1.47, well below the usual limit of 5. The blocks are correlated, but not to the point of being confounded.
• Residual variance: Breusch-Pagan test, p = .713. No strong sign of non-constant variance; even so, the tests use robust standard errors (HC3), which cost little and protect the result.
• Normality of residuals: Shapiro-Wilk, p = .049. With more than 200 cases, a slight departure does not compromise the coefficient tests.
• Influential cases: highest Cook’s distance of 0.04, far below 1. No single response changes the result.

The model explains 43% of the variation in the intention to leave (adjusted R² = .41, N = 210).

H3 says the effect of pay is smaller than those of leadership and recognition. Looking at the βs is not enough: you have to test whether the difference between them is real. The test of equality of coefficients (with the same HC3 correction) did not confirm the difference: F(1, 202) = 3.02, p = .084 against leadership and F(1, 202) = 1.92, p = .167 against recognition.

What the data showed

grafico_likert <- function(res, tx) {
  enunciados <- if (is.null(tx$enunciados)) setNames(res$itens$enunciado, res$itens$codigo) else tx$enunciados
  cod <- res$itens$codigo
  inv <- if_else(res$itens$invertido == "sim", paste0(" (", tx$invertido, ")"), "")
  rotulos_itens <- setNames(paste0(cod, inv, ": ", enunciados[cod]), cod)
  cores_likert <- setNames(c("#6a4fa3", "#b9aee0", "#d9dbe1", "#9dbde8", cor[["azul"]]), tx$likert)

  # Barras divergentes: a resposta neutra fica metade de cada lado do zero.
  lik <- res$respostas |>
    pivot_longer(everything(), names_to = "codigo", values_to = "resposta") |>
    drop_na() |>
    count(codigo, resposta) |>
    group_by(codigo) |>
    complete(resposta = 1:5, fill = list(n = 0)) |>
    arrange(resposta, .by_group = TRUE) |>
    mutate(p = n / sum(n),
           inicio = -(p[1] + p[2] + p[3] / 2) + lag(cumsum(p), default = 0),
           fim = inicio + p, meio = (inicio + fim) / 2) |>
    ungroup() |>
    mutate(nivel = factor(tx$likert[resposta], tx$likert),
           bloco = factor(tx$blocos[substr(codigo, 1, 3)], tx$blocos),
           item = factor(rotulos_itens[codigo], rev(rotulos_itens)))
  totais <- lik |>
    group_by(item, bloco) |>
    summarise(discorda = sum(p[resposta <= 2]), concorda = sum(p[resposta >= 4]), .groups = "drop")
  rotulados <- filter(lik, p >= 0.12)

  ggplot(lik, aes(y = item)) +
    geom_segment(aes(x = inicio, xend = fim, yend = item, colour = nivel), linewidth = 7) +
    geom_vline(xintercept = 0, colour = "white", linewidth = 0.6) +
    geom_text(data = rotulados, aes(x = meio, label = tx$pct(p)), size = 3, family = "Montserrat",
              colour = if_else(rotulados$resposta %in% c(1, 5), "white", cor[["profundo"]])) +
    geom_text(data = totais, aes(x = -0.82, label = tx$pct(discorda)), hjust = 1, size = 2.8,
              family = "Montserrat", colour = cor[["cinza"]]) +
    geom_text(data = totais, aes(x = 0.82, label = tx$pct(concorda)), hjust = 0, size = 2.8,
              family = "Montserrat", colour = cor[["cinza"]]) +
    facet_grid(rows = vars(bloco), scales = "free_y", space = "free_y", labeller = label_wrap_gen(14)) +
    scale_colour_manual(values = cores_likert, drop = FALSE) +
    scale_x_continuous(limits = c(-0.92, 0.92), breaks = c(-0.5, 0, 0.5), labels = \(x) tx$pct(abs(x))) +
    scale_y_discrete(labels = label_wrap_gen(40)) +
    guides(colour = guide_legend(nrow = 2, byrow = TRUE, override.aes = list(linewidth = 4))) +
    labs(x = NULL, y = NULL, colour = NULL) +
    theme(strip.text.y = element_text(angle = 0, hjust = 0, face = "bold"),
          axis.text.y = element_text(size = 10, colour = cor[["profundo"]], lineheight = 0.95),
          legend.text = element_text(size = 10), legend.justification = "left", legend.location = "plot",
          panel.grid.major.y = element_blank(), panel.spacing.y = unit(10, "pt"))
}

Figure 1: Distribution of answers to the 20 items, as given (LID3 not yet reverse-coded)

• Pay is the lowest-rated area. REM4 stands apart from the other three items in its block, the same pattern that the reliability analysis and the factor analysis pointed to.
• The reversed item works: LID3 shows the opposite pattern to the other leadership items, as it should.

grafico_grupos <- function(res, tx) {
  d <- res$d |>
    mutate(turno_grupo = traduzir(turno_grupo, tx$turno_grupo, ordem = TRUE),
           tempo_casa = traduzir(tempo_casa, tx$tempo, ordem = TRUE))
  caixa <- function(var, destaque, titulo, subtitulo, eixo_y = tx$g_eixo_int) {
    ggplot(d, aes(x = .data[[var]], y = INT)) +
      geom_boxplot(aes(fill = .data[[var]] == destaque), width = 0.55, outlier.shape = NA,
                   colour = cor[["profundo"]], linewidth = 0.4, alpha = 0.85) +
      geom_jitter(width = 0.14, height = 0.06, size = 1.3, alpha = 0.35, colour = cor[["profundo"]]) +
      scale_fill_manual(values = c(`TRUE` = cor[["azul"]], `FALSE` = cor[["lavanda"]]), guide = "none") +
      scale_y_continuous(limits = c(0.85, 5.15), breaks = 1:5) +
      scale_x_discrete(labels = label_wrap_gen(10)) +
      labs(title = titulo, subtitle = subtitulo, x = NULL, y = eixo_y) +
      theme(plot.title = element_text(face = "bold", size = 12),
            plot.subtitle = element_text(size = 10, colour = cor[["cinza"]]),
            panel.grid.major.x = element_blank())
  }
  set.seed(1)  # mesmo espalhamento dos pontos nas três línguas
  caixa("turno_grupo", tx$turno_grupo[["Noturno"]], tx$g_turno,
        paste0("Mann-Whitney, ", tx$p_txt(res$mw[["p"]]), "; r = ", tx$num1(res$mw[["r"]]))) +
    caixa("tempo_casa", tx$tempo[["Até 1 ano"]], tx$g_tempo,
          paste0("Kruskal-Wallis, ", tx$p_txt(res$kw[["p"]]), "; ε² = ", tx$num1(res$kw[["eps2"]], 3)),
          eixo_y = NULL) +
    plot_layout(widths = c(1, 2))
}

Figure 2: Turnover intention by shift and by tenure

• Nights and the first year matter. People who work nights and people with up to 1 year of tenure think more about leaving, and both effects remain in the regression model, even with the other factors controlled for.

grafico_coef <- function(res, tx, base = 12, rotulo = 3.6) {
  res$resultado |>
    mutate(nome = factor(tx$termos[termo], rev(tx$termos[termo])),
           grupo = if_else(p < 0.05, tx$g_detectado, tx$g_sem_efeito)) |>
    ggplot(aes(x = beta, y = nome, colour = grupo)) +
    geom_vline(xintercept = 0, colour = cor[["cinza"]], linewidth = 0.5, linetype = "dashed") +
    geom_errorbar(aes(xmin = li, xmax = ls), width = 0, linewidth = 1.6) +
    geom_point(size = 3.6, shape = 21, fill = "white", stroke = 2) +
    geom_text(aes(label = tx$num(beta)), vjust = -1.1, size = rotulo, family = "Montserrat",
              colour = cor[["profundo"]], show.legend = FALSE) +
    scale_colour_manual(values = setNames(c(cor[["azul"]], cor[["medio"]]), c(tx$g_detectado, tx$g_sem_efeito))) +
    scale_x_continuous(labels = \(x) tx$num(x)) +
    labs(x = tx$g_eixo_coef, y = NULL, colour = NULL) +
    theme_minimal(base_family = "Montserrat", base_size = base) +
    theme(text = element_text(colour = cor[["profundo"]]),
          axis.text.y = element_text(colour = cor[["profundo"]]),
          panel.grid.minor = element_blank(), panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_line(colour = cor[["grade"]], linewidth = 0.4),
          plot.background = element_blank(), panel.background = element_blank(),
          legend.position = "top", legend.justification = "left")
}

# Capa da página (listagem e compartilhamento), no estilo da capa do Radar:
# fundo azul profundo, logo, título, subtítulo e uma versão pequena do gráfico
# de coeficientes. A listagem mostra só a faixa do meio da imagem (cerca de
# 2/3 da altura), então título e gráfico ficam nessa faixa.
salvar_capa <- function(res, tx, arquivo = "img/capa.png", logo) {
  claro <- c(texto = "#ffffff", sub = "#cacfee", apoio = "#9aa3c7", linha = "#34344f",
             detectado = "#9dbde8", sem_efeito = "#6b7299")
  grafico <- res$resultado |>
    mutate(nome = factor(tx$termos[termo], rev(tx$termos[termo])),
           detectado = p < 0.05) |>
    ggplot(aes(x = beta, y = nome, colour = detectado)) +
    geom_vline(xintercept = 0, colour = claro[["apoio"]], linewidth = 0.5, linetype = "dashed") +
    geom_errorbar(aes(xmin = li, xmax = ls), width = 0, linewidth = 1.8) +
    geom_point(size = 3.8, shape = 21, fill = cor[["profundo"]], stroke = 2) +
    scale_colour_manual(values = c(`TRUE` = claro[["detectado"]], `FALSE` = claro[["sem_efeito"]]), guide = "none") +
    scale_x_continuous(labels = \(x) tx$num(x)) +
    labs(title = tx$capa_grafico, x = NULL, y = NULL) +
    theme_minimal(base_family = "Montserrat", base_size = 15) +
    theme(text = element_text(colour = claro[["sub"]]),
          axis.text.y = element_text(colour = claro[["sub"]], size = 14),
          axis.text.x = element_text(colour = claro[["apoio"]], size = 12),
          plot.title = element_text(colour = claro[["apoio"]], size = 14, margin = margin(b = 10)),
          plot.title.position = "plot",
          panel.grid.minor = element_blank(), panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_line(colour = claro[["linha"]], linewidth = 0.5),
          plot.background = element_blank(), panel.background = element_blank())

  texto <- function(rotulo, y, tamanho, cor_texto, peso = "plain", altura = 1.15) {
    grid::grid.text(rotulo, x = grid::unit(72, "native"), y = grid::unit(675 - y, "native"),
                    just = c("left", "top"),
                    gp = grid::gpar(fontfamily = "Montserrat", fontface = peso, fontsize = tamanho,
                                    col = cor_texto, lineheight = altura),
                    vp = grid::viewport(xscale = c(0, 1200), yscale = c(0, 675)))
  }

  dir.create(dirname(arquivo), showWarnings = FALSE)
  ragg::agg_png(arquivo, width = 1200, height = 675, res = 72, background = cor[["profundo"]])
  on.exit(grDevices::dev.off())
  grid::grid.raster(png::readPNG(logo), x = grid::unit(72, "points"), y = grid::unit(675 - 48, "points"),
                    width = grid::unit(64, "points"), just = c("left", "top"))
  linhas <- function(x) lengths(strsplit(x, "\n"))
  y_titulo <- 175
  texto(tx$capa_titulo, y_titulo, 46, claro[["texto"]], "bold", 1.1)
  y_sub <- y_titulo + 51 * linhas(tx$capa_titulo) + 22
  texto(tx$capa_subtitulo, y_sub, 22, claro[["sub"]], altura = 1.25)
  y_frase <- y_sub + 28 * linhas(tx$capa_subtitulo) + 26
  texto(tx$capa_frase, y_frase, 17, claro[["apoio"]], altura = 1.4)
  texto(tx$capa_rodape, 630, 12, claro[["apoio"]])
  print(grafico, vp = grid::viewport(x = grid::unit(1150, "points"), y = grid::unit(675 / 2, "points"),
                                     width = grid::unit(500, "points"), height = grid::unit(410, "points"),
                                     just = c("right", "center")))
}

Figure 3: Effect of each factor on the intention to leave, holding the other factors constant (β and 95% confidence interval with robust standard errors)

• Leadership and recognition are the two strongest factors (β = −0.31 and −0.30). Pay also reduces the intention (β = −0.16). The estimate is smaller, but the test did not confirm the difference: with this sample, you cannot claim that pay matters less.
• Work-life balance is the most instructive result. In the simple correlation it appears associated with the intention to leave (ρ = −.22), but in the model the effect drops to zero (β = 0.04, p = .468). Work-life balance goes hand in hand with leadership and recognition: where the supervisor is good, days off tend to be planned. On its own, it “borrows” the effect of the others.

Table 8: Hypothesis tests (multiple regression, HC3 robust standard errors)

Hypothesis	β	95% CI	p	Result
H1. Leadership reduces intention	−0.31	[−0.45; −0.18]	< .001	Supported
H2. Recognition reduces intention	−0.30	[−0.44; −0.16]	< .001	Supported
H3. Pay reduces it, with a smaller effect	−0.16	[−0.27; −0.05]	.004	Partially supported
H4. Work-life balance reduces intention	0.04	[−0.07; 0.15]	.468	Not supported

How this goes into the thesis

The student receives the tables in APA format, which most programs accept, and an example of how to report the numbers. The table below is the one that goes into the results chapter:

Table 9: Multiple linear regression of turnover intention

Predictor	β	SE	95% CI	p
Leadership	−0.31	0.07	[−0.45; −0.18]	< .001
Recognition	−0.30	0.07	[−0.44; −0.16]	< .001
Perceived pay	−0.16	0.06	[−0.27; −0.05]	.004
Work-life balance	0.04	0.06	[−0.07; 0.15]	.468
Night shift	0.45	0.12	[0.22; 0.68]	< .001
Up to 1 year of tenure	0.40	0.12	[0.16; 0.64]	.001
Female	0.12	0.11	[−0.10; 0.33]	.280

Note. N = 210. Scores and dependent variable standardized; night shift, up to 1 year of tenure and female coded as 0 or 1. SE: robust standard error (HC3). R² = .43; adjusted R² = .41; F(7, 202) = 21.9, p < .001.

Sample write-up. This paragraph shows which numbers go in and in what order. The writing of the thesis, the discussion of the results and the conclusions are the student’s.

To test the hypotheses, a multiple linear regression with robust standard errors (HC3) was estimated, with turnover intention as the dependent variable. The model explained 43% of the variance (adjusted R² = .41). Perceived leadership (β = −0.31, p < .001) and recognition (β = −0.30, p < .001) were negatively associated with turnover intention, supporting H1 and H2. Perceived pay was also negatively associated with it (β = −0.16, p = .004). Although its estimate was smaller than those for leadership and recognition, tests of equality of coefficients did not indicate a significant difference (F(1, 202) = 3.02, p = .084, and F(1, 202) = 1.92, p = .167, respectively); H3 was therefore partially supported. Work-life balance was not associated with turnover intention after controlling for the other variables (β = 0.04, p = .468); H4 was not supported.

What the student receives

• The complete report in Quarto, with the code and output for every step. If the data change (for example, responses that come in later), you just run it again.
• The commented R script and the clean dataset, so the student and the advisor can check them.
• The tables in APA format, ready to paste into Word.
• A plain-language explanation of each test: why it was chosen, what the result means and what its limits are. It helps with writing the methods section and with the defense.
• A meeting to go over the results and the questions the thesis committee usually asks about the statistics.

The work is the student’s. I do the analysis and explain it; the writing, the interpretation and the conclusions are his.

Have a survey to analyze?

It could be an undergraduate thesis, a master’s thesis or an internal company survey. Send me the project and the spreadsheet of responses, and I will tell you what can be done and how.

Get in touch