---
title: "RESPONDENTS’ PREFERRED SURVEY TOPICS: MEASUREMENT AND PREVALENCE"
subtitle: "Replication Files: Study 1"
author: "Tobias Gummer, Saskia Bartholomäus, and Bernd Weiss"
date: "19.01.2026"
output:
  html_document:
    theme: united
    toc: true
    toc_depth: 4
    toc_float:
      collapsed: true
      smooth_scroll: true
---

```{r loading dataset, message=FALSE, warning=FALSE, include=FALSE}
# Loading Dataset
load("data/subset.RData")
```

***
# Measurement Properties

Topic Popularity
 
```{r topic popularity, message=FALSE, warning=FALSE}
subset_tp <- select(subset, tp_nature, tp_crisis, tp_values, tp_economy, tp_well, tp_migration, tp_work, tp_pol, tp_media )
subset_tp <- subset_tp[complete.cases(subset_tp), ]

sumtable(subset_tp, vars=c("tp_well", "tp_values", "tp_nature", "tp_crisis", 
  "tp_work", "tp_economy", "tp_pol", "tp_media", "tp_migration"),
  labels = c("Well-Being", "Values", "Nature", "Crisis", "Work", "Economy", 
  "Politics", "Media", "Migration"), out = "kable", add.median = TRUE, digits = 3) %>% 
  kableExtra::kable_styling(bootstrap_options = "basic")
```
 
Test of Comparisons 
 
```{r}
topics <- c("tp_nature", "tp_crisis", "tp_values", "tp_economy", "tp_well", "tp_migration", "tp_work", "tp_pol", "tp_media")
tt_results <- list()

for (i in 1:(length(topics) - 1)) {
    for (j in (i + 1):length(topics)) {
      var1 <- topics[i]
      var2 <- topics[j]

     tab <- map_df(list(t.test(subset_tp[[var1]], subset_tp[[var2]], paired = TRUE, alternative = "two.sided")), tidy)
     tt_results[[paste(var1, "_vs_", var2, sep = "")]] <- tab 
    }
}

ttresults_df <- do.call(rbind, tt_results)
  ttresults_df$Comparison <- gsub("_vs_", " vs. ", rownames(ttresults_df))
  ttresults_df <- ttresults_df[c("Comparison", "estimate", "statistic", "p.value", "conf.low", "conf.high", "method")]

ttresults_df[c("Comparison", "estimate", "statistic", "p.value", "conf.low", "conf.high", "method")] %>%
  kbl(align = "llllll", col.names = c("Comparison", "Difference in Means", "t", "p.value", "conf.low", "conf.high", "Method"),
      digits = 3, caption = "Difference of Topic Popularity") %>%
  kable_styling() %>% 
  kableExtra::scroll_box(width = "100%", height = "300px")
```
```{r}
chi_results <- list()
sim_results <- list()

for (i in 1:(length(topics) - 1)) {
    for (j in (i + 1):length(topics)) {
      var1 <- topics[i]
      var2 <- topics[j]

     tab <- map_df(list(chisq.test(subset_tp[[var1]], subset_tp[[var2]])), tidy)
    
     tab2 <- map_df(list(chisq.test(subset_tp[[var1]], subset_tp[[var2]], simulate.p.value = TRUE)), tidy)
   
     chi_results[[paste(var1, "_vs_", var2, sep = "")]] <- tab 
     sim_results[[paste(var1, "_vs_", var2, sep = "")]] <- tab2
    }
}

chiresults_df <- do.call(rbind, chi_results)
  chiresults_df$Comparison <- gsub("_vs_", " vs. ", rownames(chiresults_df))
  chiresults_df <- chiresults_df[c("statistic", "p.value", "parameter", "method")]

chiresults_df[c("statistic", "p.value", "parameter", "method")] %>%
  kbl(align = "llllll", col.names = c("Chisq", "p.value", "df", "Method"), digits = 3, caption = "Difference of Topic Popularity") %>%
  kable_styling()  %>% 
  kableExtra::scroll_box(width = "100%", height = "300px")

simresults_df <- do.call(rbind, sim_results)
  simresults_df$Comparison <- gsub("_vs_", " vs. ", rownames(simresults_df))
  simresults_df <- simresults_df[c("statistic", "p.value", "parameter", "method")]

simresults_df[c("statistic", "p.value", "parameter", "method")] %>%
  kbl(align = "llllll", col.names = c("Chisq", "p.value", "df", "Method"), digits = 3, caption = "Difference of Topic Popularity") %>%
  kable_styling()  %>% 
  kableExtra::scroll_box(width = "100%", height = "300px")
```
 
By Political Interest 
 
```{r topic popularity by group, message=FALSE, warning=FALSE}
m1 <- mean(subset$tp_pol[subset$i_polinterest == 1], na.rm = T)
m2 <- mean(subset$tp_pol[subset$i_polinterest == 5], na.rm = T) 
tab <- rbind(m1, m2)
tab <- data.frame(`Topic Interest in Politics` = round(c(m1, m2), 2))
rownames(tab) <- c("low political interest", "high political interest")

tab %>%  kbl(align = "llllll", col.names = c("Topic Interest in Politics")) %>%
  kable_styling()
```
 
***
# Principal Component Analysis
 
Identifying Number of Components
 
```{r exploratory factor analysis topic interest, message=FALSE, warning=FALSE, fig.width=10, fig.height=4, out.width='100%', out.height='auto'}
i_topic <- select(subset, tp_nature, tp_crisis, tp_values, tp_economy, tp_well, tp_migration, tp_work, tp_pol, tp_media)
i_topic <- i_topic[complete.cases(i_topic), ]
  
scree_plot(i_topic, method="pc")
```  
 
Table Principal Component Analysis
 
```{r fit pca, message=FALSE, warning=FALSE}
fit_topic <- principal(i_topic, nfactors = 2, rotate = "varimax")
results_pca <- as.data.frame.matrix(fit_topic[["loadings"]])
results_pca <- cbind(results_pca, as.data.frame(fit_topic[["uniquenesses"]]))
results_pca <- cbind(results_pca, as.data.frame(fit_topic[["complexity"]]))
results_pca <- cbind(rownames(results_pca), data.frame(results_pca, row.names=NULL))
colnames(results_pca)<- c("Topic Interest", "Factor Loadings 1", "Factor Loadings 2", "Uniqueness", "Complexity")
results_pca[, 2:5] <- round(results_pca[, 2:5], 2)

n_row <- data.frame("Topic Interest" = "N",
  "Factor Loadings 1" = fit_topic$n.obs, "Factor Loadings 2" = NA, 
  "Uniqueness" = NA, "Complexity" = NA)
eigenvalues <- fit_topic$values[1:2]
eigenvalues_row <- data.frame("Topic Interest" = "Eigenvalues",
  "Factor Loadings 1" = round(eigenvalues[1], 2),
  "Factor Loadings 2" = round(eigenvalues[2], 2),
  "Uniqueness" = NA, "Complexity" = NA)

colnames(eigenvalues_row) <- colnames(results_pca)
colnames(n_row) <- colnames(results_pca)
table_pca <- rbind(results_pca, eigenvalues_row, n_row)
table_pca[table_pca[["Topic Interest"]] == "N", "Factor Loadings 1"] <- as.character(as.integer(table_pca[table_pca[["Topic Interest"]] == "N", "Factor Loadings 1"]))
table_pca[is.na(table_pca)] <- ""

vars <- c("tp_well", "tp_values", "tp_nature", "tp_crisis", "tp_work", "tp_economy", "tp_pol", "tp_media", "tp_migration")
labels <- c("Well-Being", "Values", "Nature", "Crisis", "Work", "Economy", "Politics", "Media", "Migration")

var_labels <- setNames(labels, vars)
table_pca <- table_pca %>%
  mutate(`Topic Interest` = ifelse(`Topic Interest` %in% names(var_labels), 
  var_labels[`Topic Interest`], `Topic Interest`))
```  
 
Table Principal Component Analysis
 
```{r pca table, message=FALSE, warning=FALSE}
knitr::kable(table_pca)
```
 
Plot Principal Component Analysis
 
```{r pca plot, message=FALSE, warning=FALSE, fig.width=10, fig.height=4, out.width='100%', out.height='auto'}
g_pca <- ggplot(results_pca, aes(x = `Factor Loadings 2`, y = `Factor Loadings 1`)) +
  labs(x = "Factor Loadings 1", y = "Factor Loadings 2") +
  geom_point(aes(colour = `Factor Loadings 2` < 0.5 & `Factor Loadings 1` > 0.5), show.legend = FALSE) +
  geom_text(aes(label = c("nature and environment", "current crisis", "personality and values", "economy and society", "satisfaction and well-being", "flight and migration", "work and leisure", "political attiudes and behavior", "media and social networks")), 
       hjust = 0, nudge_x = 0.02, size = 3, family = "Times New Roman", colour = "black") + 
  theme_light() + custom_theme + scale_color_grey() +
  scale_y_continuous(limits = c(0, 1.)) +
  scale_x_continuous(limits = c(0, 1.)) + 
  geom_hline(yintercept = 0.5) + 
  geom_vline(xintercept = 0.5) 
 
g_pca
```
```{r save plot, message=FALSE, warning=FALSE, include=FALSE}
 ggsave(plot = g_pca, filename = "graphs/g_pca.pdf", device = cairo_pdf,
  width = 18, height = 10, units = "cm")
```

# Reasearch Question 1. Analyses

```{r loading dataset, message=FALSE, warning=FALSE, include=FALSE}
# load dataset 
load("data/subset.RData")
load("data/long.RData")
```
 
```{r crosstab, echo=TRUE, results='hide'}
crosstab <- CrossTable(subset$matrix_mean, prop.r = TRUE, prop.c = FALSE, prop.t = FALSE, prop.chisq = FALSE)
```
```{r rq1}
t <- as.data.frame(crosstab[["t"]])
t <- t %>% mutate(across(everything(), ~ as.character(.x)))
prop <- as.data.frame(crosstab[["prop.row"]])
prop <- prop %>% mutate(across(everything(), ~ .x * 100)) # <-- Multiply by 100
prop <- prop %>% mutate(across(everything(), ~ formatC(.x, format = "f", digits = 2)))
tab <- rbind(t, prop)
rownames(tab) <- c("N", "prop.row")

tab %>% kbl(align = "llllll", col.names = c("0", "0.333", "0.667", "1"), caption = "Share of Correct Assessments") %>%
  kable_styling()
```

# Reasearch Question 2. Analyses

Regression Diagnostics
 
```{r rd, message=FALSE, warning=FALSE, fig.height=10, fig.width=16}
# main model for regression diagnostics
rmatrix <- glmer(matrix ~ vb_explain + tp_interest + s_interesting + i_education + age_group + i_gender + i_east + (1 | lfdn), 
  data = long, family = "binomial", nAGQ=0)
```
```{r subset, message=FALSE, warning=FALSE}
subset_long <- select(long, matrix, vb_explain, tp_interest, s_interesting, i_education, i_gender, age_group, i_east, lfdn)
subset_long <- subset_long[complete.cases(subset_long), ]
```
 
Regression Model
 
```{r regression models, message=FALSE, warning=FALSE}
r_main <- glmer(matrix ~ vb_explain + (1 | lfdn), 
  data = subset_long, family = "binomial", nAGQ=0)

r_edu <- glmer(matrix ~ vb_explain + i_education +(1 | lfdn), 
  data = subset_long, family = "binomial", nAGQ=0)

r_socio <- glmer(matrix ~ vb_explain + i_education + age_group + i_gender + i_east +  (1 | lfdn), 
  data = subset_long, family = "binomial", nAGQ=0)
  
r_comp <- glmer(matrix ~ vb_explain + i_education + age_group + i_gender + i_east + tp_interest + s_interesting +  (1 | lfdn), 
  data = subset_long, family = "binomial", nAGQ=0)
```
 
Comparisons Test
 
```{r regression comp, message=FALSE, warning=FALSE}
c <- avg_comparisons(r_edu, variables = list(vb_explain = c(0, 1)), comparison = "difference")

c1 <- avg_comparisons(r_comp, variables = list(s_interesting = c(1, 5)), comparison = "difference")
c2 <- avg_comparisons(r_comp, variables = list(tp_interest = c(1, 7)), comparison = "difference")
c3 <- avg_comparisons(r_comp, variables = list(i_education = c(1, 3)), comparison = "difference")

c_matrix <- rbind(c, c1, c2, c3)
c_matrix[c("term", "contrast", "estimate", "std.error", "statistic", "p.value")] %>% 
    kbl(align = "llll", digits = 3) %>%
    kable_styling()
```
 
## Regression Tables
 
```{r table, message=FALSE, warning=FALSE}
tab_model(r_main, r_edu, r_socio, r_comp, 
  show.est = TRUE, show.se = TRUE, show.ci = FALSE, show.aic = TRUE, collapse.se = TRUE, linebreak = TRUE,  p.style = "numeric", show.reflvl = TRUE,
  pred.labels = c("Intercept", "Topic Explanation", "Education: Intermediate", "Education: High", "31-48", "49-63", "63+", "Gender: Female", "Germany: East", "Topic Interest", "Survey Interest" ),
  dv.labels = c("Main Effect", "Controlled Effect I", "Controlled Effect II", "Competence Effect"), 
  title = "Mixed Effects Logistic Regression Model: Correct Assessment of a Survey Questions' Topic",
  CSS = list(css.thead = 'border-top: 1px solid;', css.summary= 'border-bottom: 1px solid;', css.table = 'width: 100%;'))
```

## Regression Graphs
 
```{r education graph, message=FALSE, warning=FALSE, include=FALSE}
p_edu <- ggpredict(r_comp, terms = c("i_education"), ci.lvl = 0.95)

g_edu <- ggplot(p_edu, aes(x = x, y = predicted)) + 
  labs(x = labels$i_education, y = NULL, title = NULL) +
  geom_point() +  geom_pointrange(aes(ymin=conf.low, ymax=conf.high)) +
  theme_light() + custom_theme + 
  scale_y_continuous(limits = c(0.0, 1.0)) +
  scale_x_discrete(breaks = c(1, 2, 3), labels = c("low","intermediate","high")) + 
  theme(strip.text = element_blank(), plot.margin = margin(t = 5, r = 5, b = 5, l = 30))
```

```{r education graph, message=FALSE, warning=FALSE, include=FALSE}
p_tp <- ggpredict(r_comp, terms = c("tp_interest"), ci.lvl = 0.95)

g_tp <- ggplot(p_tp, aes(x = x, y = predicted)) + 
  labs(x = "Topic Interest", y = NULL, title = NULL) +
  geom_point() + geom_pointrange(aes(ymin=conf.low, ymax=conf.high)) +
  theme_light() + custom_theme + 
  scale_y_continuous(limits = c(0.0, 1.0)) +
  scale_x_continuous(breaks = c(1, 2, 3, 4, 5, 6, 7), labels = c("low","","","","","","high")) + 
  theme(strip.text = element_blank(), plot.margin = margin(t = 5, r = 5, b = 5, l = 30))
```

```{r education graph, message=FALSE, warning=FALSE, include=FALSE}
p_s <- ggpredict(r_comp, terms = c("s_interesting"), ci.lvl = 0.95)

g_s <- ggplot(p_s, aes(x = x, y = predicted)) + 
  labs(x = "Survey Interest", y = NULL, title = NULL) +
  geom_point() + geom_pointrange(aes(ymin=conf.low, ymax=conf.high)) +
  theme_light() + custom_theme + 
  scale_y_continuous(limits = c(0.0, 1.0)) +
  scale_x_continuous(breaks = c(1, 2, 3, 4, 5), labels = c("low","","","","high")) + 
  theme(strip.text = element_blank(), plot.margin = margin(t = 5, r = 5, b = 5, l = 30))
```

```{r graph matrix, message=FALSE, warning=FALSE, fig.width=10, out.width='100%', out.height='auto'}
graph_matrix = ggarrange( g_edu, g_tp, g_s, ncol = 1, nrow = 3) 

graph_matrix = annotate_figure(graph_matrix,
  top = text_grob("Predicted Values\nCorrect Assessment of a Survey Questions' Topic\nwith 95% Confidence Intervals", size = 10, family = "Times New Roman", face = "plain"))

graph_matrix
```
```{r save graph mot, message=FALSE, warning=FALSE, include=FALSE}
 ggsave(plot = graph_matrix, filename = "graphs/graph_matrix.pdf", device = cairo_pdf,
  width = 16, height = 16, units = "cm")
```

# Reasearch Question 3. Analyses
 
## Correct Assessment 

t-test for Difference in Means

```{r matrix, message=FALSE, warning=FALSE}
bartlett.test(matrix_mean ~ vb_explain, data = subset)

tab <- map_df(list(t.test(subset$matrix_mean ~ subset$vb_explain, var.equal = TRUE, alternative = "two.sided")), tidy)

tab[c("estimate", "statistic", "p.value", "conf.low", "conf.high", "method")] %>% 
  kbl(align = "llllll", col.names = c("Difference in Means", "t", "p.value", "conf.low", "conf.high", "Method"), digits = 3) %>%
  kable_styling()

t.test.cluster(long$matrix, cluster = long$lfdn, group = long$vb_explain)
```
 
## Topic Preferences
 
t-test for Difference in Means

```{r topic popularity by group, message=FALSE, warning=FALSE}
vb_results <- list()

for (i in topics) {
  
    bartlett_result <- bartlett.test(subset[[i]] ~ subset$vb_explain)
    p_value_bartlett <- bartlett_result$p.value
    var_equal <- ifelse(p_value_bartlett > 0.05, TRUE, FALSE)
    
    t_test_result <- t.test(subset[[i]] ~ subset$vb_explain, alternative = "two.sided", var.equal = var_equal)
    
    t_test_tab <- broom::tidy(t_test_result)
    t_test_tab$Comparison <- paste(i, "_by_", "vb_explain", sep = "")
    vb_results[[paste(i, "_by_", "vb_explain", sep = "")]] <- rbind(t_test_tab)
}

vbresults_df <- do.call(rbind, vb_results)
vbresults_df$Comparison <- gsub("_by_", " by ", rownames(vbresults_df))
vbresults_df <- vbresults_df[c("Comparison", "estimate", "statistic", "p.value", "conf.low", "conf.high", "method")]

vbresults_df[c("Comparison", "estimate", "statistic", "p.value", "conf.low", "conf.high", "method")] %>% 
  kbl(align = "llllll", col.names = c("Comparison", "Difference in Means", "t", "p.value", "conf.low", "conf.high", "Method"), digits = 3, 
  caption = "Topic Popularity by Experimental Group") %>%
  kable_styling() %>% 
  kableExtra::scroll_box(width = "100%", height = "300px")
```

## Acquiescence
 
t-test for Difference in Means

```{r acquiescence, message=FALSE, warning=FALSE}
bartlett.test(acquiescence ~ vb_explain, data = subset)
  
tab <- map_df(list(t.test(subset$acquiescence ~ subset$vb_explain, var.equal = TRUE, alternative = "two.sided")), tidy)

tab[c("estimate", "statistic", "p.value", "conf.low", "conf.high", "method")] %>% 
  kbl(align = "llllll", col.names = c("Difference in Means", "t", "p.value", "conf.low", "conf.high", "Method"), digits = 3) %>%
    kable_styling()

rm(tab)
```

Proportion Test

```{r acquiescence, message=FALSE, warning=FALSE}
prop_test_data <- subset %>%
  group_by(vb_explain) %>%
  dplyr::summarise( mean_acquiescence = mean(acquiescence, na.rm = TRUE), n = n()) %>%
  ungroup()

prop_test_data <- prop_test_data %>%  mutate(count = mean_acquiescence * n)
prop_test <- prop.test(prop_test_data$count, prop_test_data$n)
tab <- tidy(prop_test)

tab %>%  select(estimate1, estimate2, statistic, p.value, conf.low, conf.high, method) %>%
  kbl(align = "llllll", digits = 3) %>%
  kable_styling()

rm(tab)
```

Regression Diagnostics
 
```{r main, message=FALSE, warning=FALSE}
# main model for regression diagnostics
racqu <- lm(acquiescence ~ vb_explain + i_education, data= subset)
```
```{r rd test, message=FALSE, warning=FALSE}
# Breusch-Pagan test
lmtest::bptest(racqu)  
```
```{r rd pol vote cases, message=FALSE, warning=FALSE}
# identifying influential cases
cooksd <- cooks.distance(racqu)
influential_acqu <- as.numeric(names(cooksd)[(cooksd > 4*mean(cooksd, na.rm=T))])
subset_acq <- subset[-c(influential_acqu), ]
subset_acq <- select(subset_acq, acquiescence, vb_explain, i_education)
subset_acq <- subset_acq[complete.cases(subset_acq), ]
```
 
Regression Model
 
```{r main, message=FALSE, warning=FALSE}
racqu <- lm_robust(acquiescence ~ vb_explain + i_education, data= subset)
```
 
Comparisons Test
 
```{r regression comp, message=FALSE, warning=FALSE}
c <- avg_comparisons(racqu, variables = list(vb_explain = c(0, 1)), comparison = "difference")

c[c("term", "contrast", "estimate", "std.error", "statistic", "p.value")] %>% 
    kbl(align = "llll", digits = 3) %>%
    kable_styling()
```

## Extreme Responding
 
t-test for Difference in Means

```{r extreme, message=FALSE, warning=FALSE}
bartlett.test(extreme ~ vb_explain, data = subset)
  
tab <- map_df(list(t.test(subset$extreme ~ subset$vb_explain, var.equal = TRUE, alternative = "two.sided")), tidy)

tab[c("estimate", "statistic", "p.value", "conf.low", "conf.high", "method")] %>% 
  kbl(align = "llllll", col.names = c("Difference in Means", "t", "p.value", "conf.low", "conf.high", "Method"), digits = 3) %>%
    kable_styling()

rm(tab)
```

Proportion Test

```{r extreme, message=FALSE, warning=FALSE}
prop_test_data <- subset %>%
  group_by(vb_explain) %>%
  dplyr::summarise(mean_extreme = mean(extreme, na.rm = TRUE), n = n()) %>%
  ungroup()

prop_test_data <- prop_test_data %>% mutate(count = mean_extreme * n)
prop_test <- prop.test(prop_test_data$count, prop_test_data$n)
tab <- tidy(prop_test)

tab %>% select(estimate1, estimate2, statistic, p.value, conf.low, conf.high, method) %>%
  kbl(align = "llllll", digits = 3) %>%
  kable_styling()

rm(tab)
```

Regression Diagnostics
 
```{r main, message=FALSE, warning=FALSE}
# main model for regression diagnostics
rext <- lm(extreme ~ vb_explain + i_education, data= subset)
```
```{r rd test, message=FALSE, warning=FALSE}
# Breusch-Pagan test
lmtest::bptest(rext)  
```
```{r rd pol vote cases, message=FALSE, warning=FALSE}
# identifying influential cases
cooksd <- cooks.distance(rext)
influential_ext <- as.numeric(names(cooksd)[(cooksd > 4*mean(cooksd, na.rm=T))])
subset_ext <- subset[-c(influential_ext), ]
subset_ext <- select(subset_ext, extreme, vb_explain, i_education)
subset_ext <- subset_ext[complete.cases(subset_ext), ]
```
 
Regression Models
 
```{r main, message=FALSE, warning=FALSE}
rext <- lm(extreme ~ vb_explain + i_education, data= subset)
```
 
Comparisons Test
 
```{r regression comp, message=FALSE, warning=FALSE}
c <- avg_comparisons(rext, variables = list(vb_explain = c(0, 1)), comparison = "difference")

c[c("term", "contrast", "estimate", "std.error", "statistic", "p.value")] %>% 
    kbl(align = "llll", digits = 3) %>%
    kable_styling()
```

## Midpoint Responding
 
t-test for Difference in Means

```{r midpoint, message=FALSE, warning=FALSE}
bartlett.test(midpoint ~ vb_explain, data = subset)

tab <- map_df(list(t.test(subset$midpoint ~ subset$vb_explain, var.equal = TRUE, alternative = "two.sided")), tidy)

tab[c("estimate", "statistic", "p.value", "conf.low", "conf.high", "method")] %>% 
  kbl(align = "llllll", col.names = c("Difference in Means", "t", "p.value", "conf.low", "conf.high", "Method"), digits = 3) %>%
  kable_styling()

rm(tab)
```

Proportion Test 

```{r midpoint, message=FALSE, warning=FALSE}
prop_test_data <- subset %>%
  group_by(vb_explain) %>%
  dplyr::summarise(mean_midpoint = mean(midpoint, na.rm = TRUE), n = n()) %>%
  ungroup()

prop_test_data <- prop_test_data %>% mutate(count = mean_midpoint * n)
prop_test <- prop.test(prop_test_data$count, prop_test_data$n)
tab <- tidy(prop_test)

tab %>% select(estimate1, estimate2, statistic, p.value, conf.low, conf.high, method) %>%
  kbl(align = "llllll", digits = 3) %>%
  kable_styling()

rm(tab)
```

Regression Diagnostics
 
```{r main, message=FALSE, warning=FALSE}
# main model for regression diagnostics
rmid <- lm(midpoint  ~ vb_explain + i_education, data= subset)
```
```{r rd test, message=FALSE, warning=FALSE}
# Breusch-Pagan test
lmtest::bptest(rmid)  
```
```{r rd pol vote cases, message=FALSE, warning=FALSE}
# identifying influential cases
cooksd <- cooks.distance(rmid)
influential_mid <- as.numeric(names(cooksd)[(cooksd > 4*mean(cooksd, na.rm=T))])
subset_mid <- subset[-c(influential_mid), ]
subset_mid <- select(subset_mid, midpoint, vb_explain, i_education)
subset_mid <- subset_mid[complete.cases(subset_mid), ]
```
 
Regression Models
 
```{r main, message=FALSE, warning=FALSE}
rmid <- lm(midpoint ~ vb_explain + i_education, data= subset)
```
 
Comparisons Test
 
```{r regression comp, message=FALSE, warning=FALSE}
c <- avg_comparisons(rmid, variables = list(vb_explain = c(0, 1)), comparison = "difference")

c[c("term", "contrast", "estimate", "std.error", "statistic", "p.value")] %>% 
    kbl(align = "llll", digits = 3) %>%
    kable_styling()
```

## Item Nonresponse
 
t-test for Difference in Means

```{r na, message=FALSE, warning=FALSE}
bartlett.test(na ~ vb_explain, data = subset)

tab <- map_df(list(t.test(subset$na ~ subset$vb_explain, var.equal = FALSE, alternative = "two.sided")), tidy)

tab[c("estimate", "statistic", "p.value", "conf.low", "conf.high", "method")] %>% 
  kbl(align = "llllll", col.names = c("Difference in Means", "t", "p.value", "conf.low", "conf.high", "Method"), digits = 3) %>%
  kable_styling()

rm(tab)
```

Proportion Test

```{r na, message=FALSE, warning=FALSE}
prop_test_data <- subset %>%
  group_by(vb_explain) %>%
  dplyr::summarise(mean_na = mean(na, na.rm = TRUE), n = n()) %>%
  ungroup()

prop_test_data <- prop_test_data %>% mutate(count = mean_na * n)
prop_test <- prop.test(prop_test_data$count, prop_test_data$n)
tab <- tidy(prop_test)

tab %>% select(estimate1, estimate2, statistic, p.value, conf.low, conf.high, method) %>%
  kbl(align = "llllll", digits = 3) %>%
  kable_styling()

rm(tab)
```

Regression Diagnostics
 
```{r main, message=FALSE, warning=FALSE}
# main model for regression diagnostics
rna <- lm(na  ~ vb_explain + i_education, data= subset)
```
```{r rd test, message=FALSE, warning=FALSE}
# Breusch-Pagan test
lmtest::bptest(rna)  
```
```{r rd pol vote cases, message=FALSE, warning=FALSE}
# identifying influential cases
cooksd <- cooks.distance(rna)
influential_na <- as.numeric(names(cooksd)[(cooksd > 4*mean(cooksd, na.rm=T))])
subset_na <- subset[-c(influential_na), ]
subset_na <- select(subset_na, na, vb_explain, i_education)
subset_na <- subset_na[complete.cases(subset_na), ]
```
 
Regression Models
 
```{r main, message=FALSE, warning=FALSE}
rna <- lm(na ~ vb_explain + i_education, data= subset)
```
 
Comparisons Test
 
```{r regression comp, message=FALSE, warning=FALSE}
c <- avg_comparisons(rna, variables = list(vb_explain = c(0, 1)), comparison = "difference")

c[c("term", "contrast", "estimate", "std.error", "statistic", "p.value")] %>% 
    kbl(align = "llll", digits = 3) %>%
    kable_styling()
```

## Comprehensibility 
 
t-test for Difference in Means

```{r matrix, message=FALSE, warning=FALSE}
bartlett.test(i_comp ~ vb_explain, data = subset)

tab <- map_df(list(t.test(subset$i_comp ~ subset$vb_explain, var.equal = TRUE, alternative = "two.sided")), tidy)

tab[c("estimate", "statistic", "p.value", "conf.low", "conf.high", "method")] %>% 
  kbl(align = "llllll", col.names = c("Difference in Means", "t", "p.value", "conf.low", "conf.high", "Method"), digits = 3) %>%
  kable_styling()

rm(tab)
```

Regression Diagnostics
 
```{r main, message=FALSE, warning=FALSE}
# main model for regression diagnostics
rcomp <- lm(i_comp  ~ vb_explain + i_education, data= subset)
```
```{r rd test, message=FALSE, warning=FALSE}
# Breusch-Pagan test
lmtest::bptest(rcomp)  
```
```{r rd pol vote cases, message=FALSE, warning=FALSE}
# identifying influential cases
cooksd <- cooks.distance(rcomp)
influential_comp <- as.numeric(names(cooksd)[(cooksd > 4*mean(cooksd, na.rm=T))])
subset_comp <- subset[-c(influential_comp), ]
subset_comp <- select(subset_comp, i_comp, vb_explain, i_education)
subset_comp <- subset_comp[complete.cases(subset_comp), ]
```
 
Regression Models
 
```{r main, message=FALSE, warning=FALSE}
rcomp <- lm(i_comp ~ vb_explain + i_education, data= subset)
```
 
Comparisons Test
 
```{r regression comp, message=FALSE, warning=FALSE}
c <- avg_comparisons(rcomp, variables = list(vb_explain = c(0, 1)), comparison = "difference")

c[c("term", "contrast", "estimate", "std.error", "statistic", "p.value")] %>% 
    kbl(align = "llll", digits = 3) %>%
    kable_styling()
```
 
## Regression Tables

```{r table, message=FALSE, warning=FALSE}
tab_model(racqu, rext, rmid, rna, rcomp,
  show.est = TRUE, show.se = TRUE, show.ci = FALSE, show.aic = TRUE, collapse.se = TRUE, linebreak = TRUE,  p.style = "numeric", show.reflvl = TRUE,
  pred.labels = c("Intercept", "Topic Explanation", "Education: Intermediate", "Education: High"),
  dv.labels = c("Acquiescence", "Extreme Responding", "Midpoint Responding", "Item Nonresponse", "Comprehension"), 
  title = "OLS Regression Model: Data Quality by Experimental Group",
  CSS = list(css.thead = 'border-top: 1px solid;', css.summary= 'border-bottom: 1px solid;', css.table = 'width: 100%;'))
```
 