Table of contents

Direct instruction and Workshop educational approaches are equally effective in increasing dengue knowledge in 10 years old children

-

CONICET-UNAHUR, Buenos Aires, 1688, Argentina

FMS, Puerto Iguazú, 3370, Argentina

UNA, CEMIC-CONICET, Buenos Aires, 1431, Argentina

FMS, Buenos Aires, 1061, Argentina

CONICET-FMS, Buenos Aires, 1061, Argentina

FMS, Buenos Aires, 1061, Argentina

# Fresh start for project and directories

rm(list = ls())
unlink("Tables/*")
unlink("Figures/png/*")

## Packages

library(tidyverse)
library(lme4)
library(lmerTest)
library(flextable)
library(emmeans)
library(officer)
library(car)
library(asbio)
library(ggpubr)
library(openxlsx)
library(flexlsx)

## Setting table propperties

sect_properties <- prop_section(
  page_size = page_size(
    orient = "landscape",
    width = 8.3, height = 11.7
  ),
  type = "continuous",
  page_margins = page_mar()
)

0. Import database

source("data_import.R")

0.1. Selecting cases and variables

All cases with less than 5 correct answers are removed from analyses.

base_completos = base_source %>%
  filter(Score >= 5) %>%
  ungroup() %>%
  select(Case, Time, Group, Score, Proportion,
         RT, Baseline.Score, Baseline.Proportion, Baseline.RT)

1. Table 2. Descriptive statistics for each study group and assessment time

  Descriptive = base_completos %>% 
  group_by(Group, Time)%>% 
  summarise(across(c(Score), list(n = ~n(),              # Values calculation
                                      Mean = ~mean(.x, na.rm = T),
                                      SD = ~sd(.x, na.rm = T)))) %>% 
  ungroup() %>% 
                                                                              # formatting table
  pivot_longer(
    !Group:Time,
    names_to = c("Variable", "stat"),
    names_sep = "_",
    values_to = "values"
  ) %>% 
  pivot_wider(
    names_from = c("Time", "stat"),
    values_from = "values"
  ) %>%
   
  mutate(across(where(is.numeric), ~round(.x, digits = 3))) %>% 
    flextable() %>% 
  separate_header(opts = c("span-top", "bottom-vspan"),
  split = "[_\\.]") %>% 
  merge_v(j = "Group") %>% 
  fix_border_issues(part = "all") %>% 
  set_caption(caption= as_paragraph(as_chunk("Table 2",
                props = fp_text_default(bold=TRUE, color="black",
                                        font.size = 12,))), align_with_table = T)


Descriptive

Table 2

Group	Variable	Baseline			T1			T2
		n	Mean	SD	n	Mean	SD	n	Mean	SD
Direct instruction	Score	120	9.267	2.028	124	11.952	2.237	101	11.000	2.236
Abbreviated workshop	Score	112	9.196	1.888	116	12.224	2.191	98	10.806	2.345
Workshop	Score	143	9.007	1.915	146	12.253	2.283	124	11.355	2.457

Descriptive %>%                                                                 #Saving table to "Tables" folder
  save_as_docx(path = "Tables/Table_02-Descriptive.docx", pr_section = sect_properties)



Descriptive_excel <- openxlsx2::wb_workbook()$add_worksheet("Table 2") %>% 
  wb_add_flextable("Table 2", Descriptive)

Descriptive_excel$save("Tables/Table_02-Descriptive.xlsx")

2. Baseline analyses

Baseline scores are compared between groups in order to establish if all participants share the same initial level of knowledge.

2.1. ANOVA

baseline_data = base_completos %>%                                              # Preparing data
  na.omit() %>% 
  filter(Time == "Baseline")


baseline_model = lm(data=baseline_data, Score~Group)                            # Fitting model


baseline_model_table = summary(baseline_model)$coefficients %>% 
  as_tibble(rownames = NA) %>% 
  mutate(Coefficient =c("Intercept",
                        "Workshop", "Full_Workshop"),
         `CI [.25 %]` = as.vector(confint(baseline_model)[-c(1,2),1]),
         `CI [.75 %]` = as.vector(confint(baseline_model)[-c(1,2),2])) %>% 

  mutate(`Pr(>|t|)` = round(x = `Pr(>|t|)`, digits = 4)) %>% 
  mutate(`Pr(>|t|)` = ifelse(`Pr(>|t|)` == 0, "< .001", `Pr(>|t|)`)) %>% 
  mutate(across(where(is.numeric), ~round(.x, digits = 3))) %>% 
  select(Coefficient, `Estimate`,`Std. Error`, `t value`, `Pr(>|t|)`, everything()) %>% 
  flextable() %>% 
  set_caption(caption= as_paragraph(as_chunk("Baseline analyses (score)",
                props = fp_text_default(bold=TRUE, color="black",
                                        font.size = 12,))), align_with_table = T)


baseline_model_table

Baseline analyses (score)

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)	CI [.25 %]	CI [.75 %]
Intercept	9.267	0.177	52.212	< .001	-0.733	0.214
Workshop	-0.070	0.255	-0.275	0.7835	-0.733	0.214
Full_Workshop	-0.260	0.241	-1.079	0.2813	-0.733	0.214

2.2. ANOVA model assumptions

### extracting standardized residuals
baseline_model_residuals= resid(baseline_model, type = "pearson")

# variance homogeneity between groups

homogeneity_model = leveneTest(baseline_model_residuals,baseline_data$Group) %>% 
  broom::tidy() %>%
  
  mutate(across(where(is.numeric), ~round(.x, digits = 4))) %>% 
  flextable() %>% 
    set_caption(caption= as_paragraph(as_chunk("Homogenity of variance between groups - Baseline analyses (Score)",
                props = fp_text_default(bold=TRUE, color="black",
                                        font.size = 12,))), align_with_table = T)

homogeneity_model

Homogenity of variance between groups - Baseline analyses (Score)

statistic	p.value	df	df.residual
0.2035	0.816	2	372

# normality of residuals

ggplot(data.frame(baseline_model_residuals), aes(sample = baseline_model_residuals)) +
  stat_qq() +
  labs(title = "Q-Q Plot of Residuals", x = "Theoretical Quantiles", y = "Sample Quantiles") +
  theme_minimal()+
  ggtitle("Normality of residuals checks")

2.3. Kruskall wallis Test

Given that normality of standardized residuals was not met for the ANOVA model, Kruskal Wallis test for non-parametric distributions is held.

kruskal_baseline = kruskal.test(Score ~ Group, baseline_data) %>% 
  broom::tidy() %>% 
  select(statistic, p.value) %>% 
  rename(Statistic = statistic,
         `P value` = p.value) %>% 
  mutate(`P value` = round(x = `P value`, digits = 5)) %>% 
  mutate(`P value` = ifelse(`P value` == 0, "< .001", `P value`)) %>% 
  mutate(across(where(is.numeric), ~round(.x, digits = 3))) %>% 
  flextable() %>% 
  set_caption(caption= as_paragraph(as_chunk("Non parametric Baseline score analyses",
                props = fp_text_default(bold=TRUE, color="black",
                                        font.size = 12,))), align_with_table = F)

kruskal_baseline

Non parametric Baseline score analyses

Statistic	P value
1.816	0.403

3. Main Analysis

Linear mixed models analysis is held in order to study score differences between groups and time assessments. This is the main analysis.

3.1. Table 3. Generalized Mixed Model results and model fitness values

data.model <- base_completos %>%                             # Preparing data
  filter(Time != "Baseline")

model =                                                      # Model fitting
  lmer( Score ~ Time + Group + Baseline.Score + 
          Baseline.Score:Time + Time:Group + 
          (1|Case),
        data= data.model,
        control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5))) 



### validate model 
rsquared = summary(lm(predict(model)~data.model$Score))$r.squared 
AIC = AIC(model)
BIC = BIC(model)


model_table = summary(model)$coefficients %>% 
  as_tibble(rownames = NA) %>% 
  mutate(Coefficient =c("Intercept","Time",
                        "Workshop", "Full_Workshop",
                        "Baseline",
                        "Baseline : Time",
                        "Time : Workshop",
                        "Time : Full_Workshop"),
         `CI [.25 %]` = as.vector(confint(model)[-c(1,2),1]),
         `CI [.75 %]` = as.vector(confint(model)[-c(1,2),2])) %>% 
    mutate(`R squared` = c(round(rsquared, digits =3), rep("",7))) %>%
  
  mutate(AIC = c(round(AIC, digits =3), rep("",7))) %>%
  mutate(BIC = c(round(BIC, digits =3), rep("",7))) %>%
  mutate(`Pr(>|t|)` = round(x = `Pr(>|t|)`, digits = 4)) %>% 
  mutate(`Pr(>|t|)` = ifelse(`Pr(>|t|)` == 0, "< .001", `Pr(>|t|)`)) %>% 
  mutate(across(where(is.numeric), ~round(.x, digits = 3))) %>% 
  select(Coefficient, `Estimate`,`Std. Error`, df, `t value`, `Pr(>|t|)`, everything()) %>% 
  flextable() %>% 
  set_caption(caption= as_paragraph(as_chunk("Linear Mixed Models coefficients estimates (Score)",
                props = fp_text_default(bold=TRUE, color="black",
                                        font.size = 12,))), align_with_table = F)


model_table

Linear Mixed Models coefficients estimates (Score)

Coefficient	Estimate	Std. Error	df	t value	Pr(>|t|)	CI [.25 %]	CI [.75 %]	R squared	AIC	BIC
Intercept	9.128	0.480	590.700	19.025	< .001	8.191	10.065	0.768	3065.714	3111.353
Time	-1.145	0.544	362.117	-2.105	0.036	-2.207	-0.081
Workshop	0.282	0.282	582.891	0.998	0.3185	-0.269	0.833
Full_Workshop	0.340	0.267	581.769	1.273	0.2034	-0.181	0.861
Baseline	0.314	0.049	593.050	6.392	< .001	0.218	0.410
Baseline : Time	0.015	0.056	363.174	0.263	0.7924	-0.094	0.123
Time : Workshop	-0.348	0.309	343.978	-1.129	0.2596	-0.952	0.254
Time : Full_Workshop	0.204	0.292	343.191	0.701	0.4837	-0.366	0.774

model_table %>%
  save_as_docx(path = "Tables/Table03.docx", pr_section = sect_properties)

model_table_excel <- openxlsx2::wb_workbook()$add_worksheet("Table 3") %>% 
  wb_add_flextable("Table 3", model_table)

model_table_excel$save("Tables/Table03.xlsx")

3.2. Checking Model assumptions

### extracting residuals and fitted values
residuals <- residuals(model)
fitted_values <- fitted(model)

# Heterocedasticity

p1 <- ggplot(data.frame(fitted_values, residuals), aes(fitted_values, residuals)) +
  geom_point(alpha = 0.6) +
  labs(title = "Residuals vs Fitted Values", x = "Fitted Values", y = "Residuals") +
  theme_minimal()+
  ggtitle("Heterocedasticity checks")


# normality of residuals
p2 <- ggplot(data.frame(residuals), aes(sample = residuals)) +
  stat_qq() +
  labs(title = "Q-Q Plot of Residuals", x = "Theoretical Quantiles", y = "Sample Quantiles") +
  theme_minimal()+
  ggtitle("Normality of residuals checks")




# normality of random effects over Intercept, by grouping variable (Case)
exp1_model_reff = ranef(model)$Case$`(Intercept)`

p3 <- ggplot(data.frame(exp1_model_reff), aes(sample = exp1_model_reff)) +
  stat_qq() +
  labs(title = "Q-Q Plot of Residuals", x = "Theoretical Quantiles", y = "Sample Quantiles") +
  theme_minimal()+
  ggtitle("Normality of random effects checks")


ggpubr::ggarrange(p1, p2, p3)

4. Figure 1

exp1_cols= c("#aeaeae","#9D7AD2","#99ccff","#ffae3b") # group colors

Figure1 = base_completos%>%  
  group_by(Time, Group) %>% 
  summarise(Score_mean = mean(Score),
            SE = sqrt(sum((Score-mean(Score))^2/(length(Score)-1)))/sqrt(length(Score))) %>% 
  ggplot(aes(x=Time, y=Score_mean, group=Group, color=Group, fill= Group))+
  
  scale_fill_manual(values=exp1_cols)+
  scale_color_manual(values=exp1_cols)+
  scale_shape_manual(values = c(21, 25, 22,24)) +

  geom_line(lty=2,lwd=1, position=position_dodge(width = 0.2), alpha= .7)+
  geom_point(aes(shape= Group), 
             position=position_dodge(width = 0.2),size= 5)+
  geom_errorbar(aes(ymin=Score_mean-SE, ymax=Score_mean+SE),
                position=position_dodge(width = 0.2),width=0,lwd=1.75)+
  labs(y="Score",x="Time") +
  theme_bw(base_size = 15) +
  theme(legend.position=c(0.7,0.15)
        # legend.text = element_text(size = 10)
        )


Figure1

# 
# Figure1%>%
#   ggexport(filename = "Figures/png/Figure1.png", width = 5.6, height = 8.6,res = 600,)

Figure1 %>% 
  ggsave(filename = "Figures/png/Figure1.png", device = "png", width = 5.6, height = 5.6, units = "in", dpi = 600)
Figure1 %>% 
  ggsave(filename = "Figures/tif/Figure1.tif", width = 5.6, height = 5.6, units = "in", dpi = 300)

5. Table S3. Linear Mixed Models between groups comparisons

#Between groups

btw_comparisons = emmeans(model, revpairwise ~ Group|Time)$contrasts %>% 
  summary(infer = TRUE) %>% 
  as_tibble() %>% 
  rename(Contrast = contrast,
         Estimate = estimate,
         `Lower CL` = lower.CL,
         `Upper CL` = upper.CL,
         `T ratio` = t.ratio,
         `P value` = p.value) %>% 
  mutate(Contrast = rep(c("Abb. Workshop - Direct instruction", 
                      "Workshop - Direct Instruction", 
                      "Workshop - Abb. Workshop"), 2)) %>% 
  mutate(`P value` = round(x = `P value`, digits = 4)) %>% 
  mutate(`P value` = ifelse(`P value` == 0, "< .001", `P value`)) %>% 
  mutate(across(where(is.numeric), ~round(.x, digits = 2))) %>% 
  flextable() %>% 
  set_caption(caption= as_paragraph(as_chunk("Table S3. Linear Mixed Models between groups comparisons ",
                props = fp_text_default(bold=TRUE, color="black",
                                        font.size = 12,))), align_with_table = F)


btw_comparisons

Table S3. Linear Mixed Models between groups comparisons

Contrast	Time	Estimate	SE	df	Lower CL	Upper CL	T ratio	P value
Abb. Workshop - Direct instruction	T1	0.28	0.28	585.04	-0.38	0.95	1.00	0.58
Workshop - Direct Instruction	T1	0.34	0.27	583.94	-0.29	0.97	1.27	0.41
Workshop - Abb. Workshop	T1	0.06	0.27	581.37	-0.58	0.70	0.21	0.98
Abb. Workshop - Direct instruction	T2	-0.07	0.30	639.84	-0.78	0.65	-0.22	0.97
Workshop - Direct Instruction	T2	0.54	0.29	639.26	-0.13	1.22	1.89	0.14
Workshop - Abb. Workshop	T2	0.61	0.29	631.86	-0.07	1.29	2.10	0.09

btw_comparisons %>%
  save_as_docx(path = "Tables/TableS3-Pairwise_between_LMM_Score.docx", pr_section = sect_properties)

btw_comparisons_excel <- openxlsx2::wb_workbook()$add_worksheet("Table S3") %>% 
  wb_add_flextable("Table S3", btw_comparisons)

btw_comparisons_excel$save("Tables/TableS3-Pairwise_between_LMM_Score.xlsx")

6. Table S4. Linear Mixed Models within groups comparisons

### Within comparisons

wtn_comparisons = emmeans(model, revpairwise ~ Time|Group)$contrasts %>% 
  summary(infer = TRUE) %>% 
  as_tibble() %>% 
  rename(Contrast = contrast,
         Estimate = estimate,
         `Lower CL` = lower.CL,
         `Upper CL` = upper.CL,
         `T ratio` = t.ratio,
         `P value` = p.value) %>% 
  mutate(Group = c("Direct instruction", 
                      "Abb. Workshop", 
                      "Workshop")) %>% 
  mutate(`P value` = round(x = `P value`, digits = 4)) %>% 
  mutate(`P value` = ifelse(`P value` == 0, "< .001", `P value`)) %>% 
  mutate(across(where(is.numeric), ~round(.x, digits = 2))) %>% 
  flextable() %>% 
  set_caption(caption= as_paragraph(as_chunk("Linear Mixed Models within groups comparisons (Score)",
                props = fp_text_default(bold=TRUE, color="black",
                                        font.size = 12,))), align_with_table = F)


wtn_comparisons

Linear Mixed Models within groups comparisons (Score)

Contrast	Group	Estimate	SE	df	Lower CL	Upper CL	T ratio	P value
T2 - T1	Direct instruction	-1.01	0.22	353.05	-1.44	-0.59	-4.68	< .001
T2 - T1	Abb. Workshop	-1.36	0.22	343.11	-1.80	-0.93	-6.20	< .001
T2 - T1	Workshop	-0.81	0.20	339.69	-1.19	-0.43	-4.15	< .001

wtn_comparisons %>%
  save_as_docx(path = "Tables/TableS4-Pairwise_within_LMM_Score.docx", pr_section = sect_properties)


wtn_comparisons_excel <- openxlsx2::wb_workbook()$add_worksheet("Table S4") %>% 
  wb_add_flextable("Table S4", wtn_comparisons)

wtn_comparisons_excel$save("Tables/TableS4-Pairwise_within_LMM_Score.xlsx")

7. Delta analysis

# Data preparation

base_T1 = base_completos %>% 
  filter(Time != "T2") %>% 
  group_by(Case, Group) %>%
  summarise(Baseline.vs.T1 = diff(Score)) %>% 
  ungroup()


base_T2 = base_completos %>% 
  filter(Time != "T1") %>% 
  group_by(Case, Group) %>%
  summarise(Baseline.vs.T2 = diff(Score)) %>% 
  ungroup()


base_deltas = left_join(base_T1, base_T2)

# ANOVA


# Descriptive 

descriptivos_deltas = base_deltas %>% 
  group_by(Group)%>% 
  summarise(n = n(),
            across(c(Baseline.vs.T1, Baseline.vs.T2), list(Mean = ~mean(.,na.rm = T),
                                                 SD = ~sd(.,na.rm = T)))) 

# Main model

base_deltas %>%
  select(Baseline.vs.T1, Baseline.vs.T2) %>%
  map(~ aov(.x ~ Group, data = base_deltas)) %>%
  map_dfr(~ broom::tidy(.), .id = 'source') %>%
  filter(term != "Residuals") %>% 
  rename(F.statistic = statistic) %>% 
  select(source,F.statistic, p.value) %>% 
  pivot_wider(names_from = source, values_from = c("F.statistic", "p.value"),
              names_glue = "{source}_{.value}") %>% 
  mutate(across(where(is.numeric),
                ~ifelse(.x < 0.001 & .x != 0, "< .001", .x)))-> anova_deltas

## Final table

tabla_deltas = cbind(descriptivos_deltas,
                               anova_deltas) %>% 
  
  
  select(Group:Baseline.vs.T1_SD, Baseline.vs.T1_F.statistic, Baseline.vs.T1_p.value,
         Baseline.vs.T2_Mean: Baseline.vs.T2_SD, Baseline.vs.T2_F.statistic, Baseline.vs.T2_p.value) %>%
  mutate(across(where(is.numeric), ~round(.x,3))) %>%

  flextable() %>%
  separate_header(split = "[_\\_]") %>%
    merge_v(j = c("Baseline.vs.T1_p.value", "Baseline.vs.T1_F.statistic",
                  "Baseline.vs.T2_p.value", "Baseline.vs.T2_F.statistic")) %>%
  fix_border_issues(part = "all")

tabla_deltas

Group	n	Baseline.vs.T1				Baseline.vs.T2
		Mean	SD	F.statistic	p.value	Mean	SD	F.statistic	p.value
Direct instruction	117	2.752	2.435	1.349	0.261	1.806	2.276	2.893	0.057
Abbreviated workshop	112	3.116	2.485			1.747	2.701
Workshop	143	3.231	2.300			2.484	2.585

8. Question-by-question analysis

8.1. Descriptive

Frecuencias_conceptos = base_long_concept %>% 
    group_by(Group, Time, Pregunta) %>%
    summarise(Frecuencia = sum(Respuesta, na.rm = T),
              Porcentaje = (Frecuencia/n())*100,
              n = n()) %>% 
    mutate(Time = factor(Time, levels = c("T2", "T1", "Baseline"))) %>% 
  ungroup()


  Descriptive = Frecuencias_conceptos %>% 
  select(Time, Group, Pregunta, Porcentaje) %>% 
    pivot_wider(id_cols = c("Pregunta", "Group"),
    names_from = c( "Time"),
    values_from = "Porcentaje"
  ) %>%
    arrange(desc(Pregunta)) %>% 
   
  mutate(across(where(is.numeric), ~round(.x, digits = 3))) %>% 
    flextable() %>% 
  merge_v(j = "Pregunta") %>% 
     width(j = 1, width = 10) %>% 
  valign(j = 1, valign = "top") %>% 
  fix_border_issues(part = "all") %>% 
  set_caption(caption= as_paragraph(as_chunk("Descriptive statistics",
                props = fp_text_default(bold=TRUE, color="black",
                                        font.size = 12,))), align_with_table = T)



Descriptive

Descriptive statistics

Pregunta	Group	Baseline	T1	T2
To avoid mosquitoes, do you have to throw away containers that collect stagnant water?	Direct instruction	62.992	76.800	80.583
	Abb. Workshop	64.957	81.197	74.000
	Workshop	65.753	80.137	80.645
To avoid mosquito bites, do you have to wear short-sleeved clothing?	Direct instruction	77.165	79.200	84.466
	Abb. Workshop	80.342	78.632	76.000
	Workshop	76.712	82.877	88.710
Is the aedes aegypti black with white spots and stripes on its body and legs?	Direct instruction	63.780	89.600	84.466
	Abb. Workshop	61.538	91.453	86.000
	Workshop	58.219	91.096	92.742
If you think you may have dengue, should you take an aspirin?	Direct instruction	45.669	87.200	77.670
	Abb. Workshop	50.427	81.197	69.000
	Workshop	44.521	82.192	70.968
If a container collects water and it cannot be emptied, should we cover it so that mosquitoes do not breed there?	Direct instruction	85.827	87.200	84.466
	Abb. Workshop	79.487	88.889	82.000
	Workshop	82.192	89.726	85.484
Does the use of repellents help avoid mosquito bites?	Direct instruction	88.976	90.400	94.175
	Abb. Workshop	85.470	91.453	87.000
	Workshop	92.466	93.836	90.323
Does the mosquito that transmits dengue bite at noon?	Direct instruction	40.945	55.200	39.806
	Abb. Workshop	35.897	60.684	42.000
	Workshop	35.616	55.479	50.000
Does the mosquito that spreads dengue fever breed more in the jungle than in the city?	Direct instruction	37.795	65.600	54.369
	Abb. Workshop	29.060	67.521	59.000
	Workshop	28.767	76.712	65.323
Does the aedes aegypti mosquito lay its eggs in streams?	Direct instruction	26.772	78.400	49.515
	Abb. Workshop	28.205	76.923	51.000
	Workshop	24.658	81.507	54.839
Does dengue cause fever and body aches?	Direct instruction	82.677	92.000	88.350
	Abb. Workshop	87.179	94.017	87.000
	Workshop	91.096	90.411	88.710
Do mosquitoes bite more on the face and back than on the legs and arms?	Direct instruction	37.795	53.600	38.835
	Abb. Workshop	38.462	58.974	52.000
	Workshop	38.356	67.808	54.032
Can dengue kill you?	Direct instruction	54.331	84.000	57.282
	Abb. Workshop	62.393	88.034	76.000
	Workshop	71.918	92.466	80.645
Besides mosquitoes, can people also spread dengue fever?	Direct instruction	43.307	84.000	71.845
	Abb. Workshop	41.880	75.214	62.000
	Workshop	32.192	59.589	53.226
Before flying, is the mosquito a larva that lives in water?	Direct instruction	63.780	85.600	89.320
	Abb. Workshop	66.667	91.453	76.000
	Workshop	60.274	90.411	93.548
After it rains, do we have to eliminate the containers with water in them to avoid mosquito breeding?	Direct instruction	80.315	91.200	89.320
	Abb. Workshop	78.632	89.744	85.000
	Workshop	83.562	91.096	86.290

##8.1. GLMM coefficients for each question

table_list = list()

Questions_vector = unique(base_long_concept$Pregunta)

for (i in 1:length(Questions_vector)) {
  
  baseline_response_df = base_long_concept %>% 
    filter(Time == "Baseline") %>% 
    select(Caso, Pregunta, Group, Respuesta) %>% 
    rename(Baseline_response = "Respuesta")
  
  data.model <- base_long_concept %>%                             # Preparing data
    filter(Time != "Baseline") %>%
    filter(Pregunta == Questions_vector[i]) %>% 
    left_join(baseline_response_df) %>% 
    filter(!is.na(Respuesta))%>% 
    filter(!is.na(Baseline_response))
  
  
  model =                                                      # Model fitting
    glmer( Respuesta ~ Time + Group + Baseline_response + 
            Baseline_response:Time + Time:Group + 
            (1|Caso),
          data= data.model,
          family=binomial,
          control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5))) 
  
  
  # summary(model)
  
  ### validate model 
  rsquared = summary(lm(predict(model)~data.model$Respuesta))$r.squared 
  AIC = AIC(model)
  BIC = BIC(model)
  
  
table_list[[i]] = summary(model)$coefficients %>% 
    as_tibble(rownames = NA) %>% 
    mutate(Coefficient =c("Intercept","Time",
                          "Workshop", "Full_Workshop",
                          "Baseline",
                          "Baseline : Time",
                          "Time : Workshop",
                          "Time : Full_Workshop")
           ) %>%
    mutate(`R squared` = c(round(rsquared, digits =3), rep("",7))) %>%
    mutate(AIC = c(round(AIC, digits =3), rep("",7))) %>%
    mutate(BIC = c(round(BIC, digits =3), rep("",7))) %>%
    add_row(Coefficient = Questions_vector[i], .before = 1)                     ## Adding question text to table
  

}

table_questions_models = reduce(table_list, full_join) %>% 

  mutate(`Pr(>|z|)` = round(x = `Pr(>|z|)`, digits = 4)) %>%
  mutate(`Pr(>|z|)` = ifelse(`Pr(>|z|)` < 0.001, "< .001", `Pr(>|z|)`)) %>% 
  mutate(across(where(is.numeric), ~round(.x, digits = 3))) %>% 
  select(Coefficient, `Estimate`,`Std. Error`, `z value`, `Pr(>|z|)`, everything()) 


table_questions_models_flex = table_questions_models%>% 
    flextable() %>%
    merge_h_range(i = seq(1, nrow(table_questions_models), by = 9),
                  j1 = 1, j2 = ncol(table_questions_models)) %>%
    set_caption(caption= as_paragraph(as_chunk(paste("Generalized Linear Mixed Models coefficients estimates for each question", Questions_vector[i]),
                  props = fp_text_default(bold=TRUE, color="black",
                                          font.size = 12,))), align_with_table = F) 
#   
table_questions_models_flex

Generalized Linear Mixed Models coefficients estimates for each question Does the mosquito that transmits dengue bite at noon?

Coefficient	Estimate	Std. Error	z value	Pr(>|z|)	R squared	AIC	BIC
Does the use of repellents help avoid mosquito bites?
Intercept	4.587	1.661	2.761	0.0058	0.603	325.107	365.976
Time	2.432	1.559	1.560	0.1187
Workshop	0.881	1.392	0.633	0.5267
Full_Workshop	1.242	1.407	0.883	0.3774
Baseline	4.033	1.366	2.951	0.0032
Baseline : Time	-2.615	1.365	-1.916	0.0554
Time : Workshop	-1.694	1.416	-1.196	0.2317
Time : Full_Workshop	-1.883	1.436	-1.311	0.1899
Does the mosquito that spreads dengue fever breed more in the jungle than in the city?
Intercept	0.785	0.269	2.916	0.0035	0.522	860.557	901.269
Time	-0.712	0.360	-1.979	0.0478
Workshop	0.259	0.342	0.758	0.4482
Full_Workshop	0.660	0.335	1.971	0.0487
Baseline	0.010	0.293	0.036	0.9717
Baseline : Time	0.422	0.393	1.073	0.2833
Time : Workshop	0.121	0.460	0.264	0.792
Time : Full_Workshop	-0.095	0.443	-0.214	0.8302
Can dengue kill you?
Intercept	1.997	0.470	4.246	< .001	0.493	605.595	646.333
Time	-1.998	0.539	-3.709	< .001
Workshop	0.337	0.480	0.703	0.4822
Full_Workshop	0.990	0.513	1.931	0.0535
Baseline	0.600	0.417	1.437	0.1507
Baseline : Time	0.206	0.502	0.411	0.681
Time : Workshop	1.019	0.598	1.702	0.0887
Time : Full_Workshop	0.564	0.606	0.930	0.3526
To avoid mosquitoes, do you have to throw away containers that collect stagnant water?
Intercept	0.726	0.318	2.282	0.0225	0.35	662.544	703.15
Time	0.791	0.470	1.684	0.0922
Workshop	0.446	0.379	1.176	0.2395
Full_Workshop	0.294	0.354	0.830	0.4066
Baseline	1.070	0.312	3.428	< .001
Baseline : Time	-0.515	0.431	-1.195	0.232
Time : Workshop	-0.942	0.537	-1.755	0.0793
Time : Full_Workshop	-0.481	0.512	-0.940	0.347
To avoid mosquito bites, do you have to wear short-sleeved clothing?
Intercept	0.715	0.475	1.506	0.132	0.534	598.998	639.88
Time	0.662	0.593	1.115	0.2649
Workshop	-0.080	0.465	-0.171	0.8643
Full_Workshop	0.466	0.457	1.020	0.3079
Baseline	1.865	0.499	3.733	< .001
Baseline : Time	-0.096	0.560	-0.171	0.8639
Time : Workshop	-0.713	0.609	-1.170	0.242
Time : Full_Workshop	0.246	0.639	0.385	0.7006
Do mosquitoes bite more on the face and back than on the legs and arms?
Intercept	-0.265	0.253	-1.048	0.2945	0.435	866.469	907.048
Time	-0.843	0.367	-2.298	0.0216
Workshop	0.398	0.328	1.211	0.2259
Full_Workshop	0.886	0.324	2.737	0.0062
Baseline	0.991	0.283	3.506	< .001
Baseline : Time	0.601	0.383	1.568	0.1168
Time : Workshop	0.414	0.461	0.898	0.3693
Time : Full_Workshop	-0.089	0.442	-0.201	0.841
Does the aedes aegypti mosquito lay its eggs in streams?
Intercept	1.708	0.319	5.349	< .001	0.487	821.045	861.901
Time	-1.628	0.384	-4.238	< .001
Workshop	-0.018	0.387	-0.048	0.962
Full_Workshop	0.176	0.371	0.475	0.635
Baseline	-0.104	0.344	-0.302	0.7626
Baseline : Time	-0.031	0.434	-0.072	0.9426
Time : Workshop	0.079	0.491	0.161	0.8718
Time : Full_Workshop	0.141	0.468	0.300	0.764
If a container collects water and it cannot be emptied, should we cover it so that mosquitoes do not breed there?
Intercept	4.938	1.212	4.074	< .001	0.585	464.52	505.153
Time	1.226	1.151	1.065	0.2867
Workshop	1.057	0.983	1.075	0.2823
Full_Workshop	0.892	0.898	0.993	0.3207
Baseline	1.670	0.938	1.780	0.075
Baseline : Time	-1.110	1.019	-1.089	0.276
Time : Workshop	-1.310	1.054	-1.243	0.2138
Time : Full_Workshop	-0.893	0.965	-0.925	0.3548
After it rains, do we have to eliminate the containers with water in them to avoid mosquito breeding?
Intercept	8.594	7.149	1.202	0.2293	0.549	252.707	292.958
Time	10.679	9.100	1.174	0.2406
Workshop	3.159	6.886	0.459	0.6465
Full_Workshop	2.758	6.604	0.418	0.6763
Baseline	1.687	6.698	0.252	0.8011
Baseline : Time	-1.760	5.753	-0.306	0.7596
Time : Workshop	-13.205	8.152	-1.620	0.1053
Time : Full_Workshop	-12.358	8.061	-1.533	0.1253
Besides mosquitoes, can people also spread dengue fever?
Intercept	1.843	0.374	4.928	< .001	0.457	751.541	791.972
Time	-0.976	0.441	-2.212	0.027
Workshop	-0.831	0.440	-1.888	0.0591
Full_Workshop	-1.697	0.422	-4.023	< .001
Baseline	1.416	0.360	3.927	< .001
Baseline : Time	-0.237	0.438	-0.541	0.5885
Time : Workshop	0.268	0.550	0.487	0.6263
Time : Full_Workshop	0.581	0.513	1.131	0.258
Before flying, is the mosquito a larva that lives in water?
Intercept	6.822	0.002	4,376.678	< .001	0.529	399.404	439.929
Time	2.203	1.087	2.026	0.0427
Workshop	2.481	0.002	1,592.007	< .001
Full_Workshop	0.446	0.669	0.668	0.5044
Baseline	0.284	0.002	182.051	< .001
Baseline : Time	1.395	0.918	1.519	0.1286
Time : Workshop	-6.307	1.132	-5.574	< .001
Time : Full_Workshop	-0.496	1.393	-0.356	0.7218
If you think you may have dengue, should you take an aspirin?
Intercept	2.739	0.509	5.379	< .001	0.583	646.412	686.91
Time	-0.900	0.504	-1.787	0.0739
Workshop	-0.508	0.498	-1.019	0.3081
Full_Workshop	-0.642	0.467	-1.375	0.1691
Baseline	0.252	0.380	0.663	0.5075
Baseline : Time	0.534	0.474	1.126	0.26
Time : Workshop	-0.355	0.618	-0.574	0.5659
Time : Full_Workshop	-0.059	0.578	-0.102	0.9187
Is the aedes aegypti black with white spots and stripes on its body and legs?
Intercept	6.573	0.003	2,099.183	< .001	0.544	358.356	398.718
Time	-1.073	0.003	-341.884	< .001
Workshop	0.057	0.849	0.067	0.9468
Full_Workshop	0.033	0.003	10.643	< .001
Baseline	1.961	0.003	626.222	< .001
Baseline : Time	-0.814	0.003	-259.469	< .001
Time : Workshop	0.685	0.832	0.823	0.4105
Time : Full_Workshop	1.310	0.003	417.392	< .001
Does dengue cause fever and body aches?
Intercept	7.838	1.810	4.330	< .001	0.624	318.988	359.486
Time	-1.194	1.542	-0.774	0.4388
Workshop	0.484	1.607	0.301	0.7632
Full_Workshop	-0.992	1.310	-0.757	0.4489
Baseline	1.245	1.583	0.786	0.4317
Baseline : Time	0.158	1.486	0.106	0.9153
Time : Workshop	-0.543	1.529	-0.355	0.7227
Time : Full_Workshop	1.033	1.327	0.778	0.4365
Does the mosquito that transmits dengue bite at noon?
Intercept	0.131	0.246	0.532	0.5949	0.386	888.916	929.167
Time	-0.637	0.359	-1.775	0.0759
Workshop	0.223	0.314	0.710	0.4777
Full_Workshop	-0.065	0.293	-0.221	0.8253
Baseline	0.543	0.253	2.150	0.0315
Baseline : Time	-0.115	0.353	-0.327	0.7438
Time : Workshop	-0.083	0.446	-0.187	0.852
Time : Full_Workshop	0.419	0.419	1.001	0.3166

# 
# 
# table_questions_models_flex %>%
#   save_as_docx(path = "Tables/TableS4-GLMM_questions.docx", pr_section = sect_properties)
# 

8.2. Table S5. Linear Mixed Models between groups comparisons for significant questions

significant_questions <- c("Do mosquitoes bite more on the face and back than on the legs and arms?",
                           "Besides mosquitoes, can people also spread dengue fever?",
                           "Before flying, is the mosquito a larva that lives in water?",
                           "Is the aedes aegypti black with white spots and stripes on its body and legs?")

table_list = list()
  
for (i in 1:length(significant_questions)) {

  baseline_response_df = base_long_concept %>% 
    filter(Time == "Baseline") %>% 
    select(Caso, Pregunta, Group, Respuesta) %>% 
    rename(Baseline_response = "Respuesta")
  
  data.model <- base_long_concept %>%                             # Preparing data
    filter(Time != "Baseline") %>%
    filter(Pregunta == significant_questions[i]) %>% 
    left_join(baseline_response_df) %>% 
    filter(!is.na(Respuesta))%>% 
    filter(!is.na(Baseline_response))
  
  
  model =                                                      # Model fitting
    glmer( Respuesta ~ Time + Group + Baseline_response + 
            Baseline_response:Time + Time:Group + 
            (1|Caso),
          data= data.model,
          family=binomial,
          control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5))) 

#Between groups

table_list[[i]] = emmeans(model, revpairwise ~ Group|Time, p.adjust = "holm")$contrasts %>% 
    summary(infer = TRUE) %>% 
    as_tibble() %>% 
    rename(Contrast = contrast,
           Estimate = estimate,
           `Lower CL` = asymp.LCL,
           `Upper CL` = asymp.UCL,
           `Z ratio` = z.ratio,
           `P value` = p.value) %>% 
  mutate(Contrast = rep(c("Abb. Workshop - Direct instruction", 
                      "Workshop - Direct instruction", 
                      "Workshop - Abb. Workshop"), 2)) %>% 

    add_row(Contrast = significant_questions[i], .before = 1)                     ## Adding question text to table


}

btw_commparisons = reduce(table_list, full_join) %>% 
    mutate(`P value` = round(x = `P value`, digits = 4)) %>%
    mutate(`P value` = ifelse(`P value` < 0.001, "< .001", `P value`)) %>%
    mutate(across(where(is.numeric), ~round(.x, digits = 2))) %>% 
  select(Contrast, Time, Estimate, SE, `Lower CL`, `Upper CL`, `Z ratio`, `P value`)
  
  
btw_comparsons_flex = btw_commparisons %>% 
    flextable() %>%
  
    merge_h_range(i = seq(1, nrow(btw_commparisons), by = 7), j1 = 1, j2 = ncol(btw_commparisons)) %>%
    set_caption(caption= as_paragraph(as_chunk("Table S5. Linear Mixed Models between groups comparisons for significant questions",
                  props = fp_text_default(bold=TRUE, color="black",
                                          font.size = 12,))), align_with_table = F)
btw_comparsons_flex

Table S5. Linear Mixed Models between groups comparisons for significant questions

Contrast	Time	Estimate	SE	Lower CL	Upper CL	Z ratio	P value
Do mosquitoes bite more on the face and back than on the legs and arms?
Abb. Workshop - Direct instruction	T1	0.40	0.33	-0.37	1.17	1.21	0.4466
Workshop - Direct instruction	T1	0.89	0.32	0.13	1.64	2.74	0.017
Workshop - Abb. Workshop	T1	0.49	0.32	-0.27	1.25	1.51	0.2863
Abb. Workshop - Direct instruction	T2	0.81	0.38	-0.07	1.69	2.16	0.0784
Workshop - Direct instruction	T2	0.80	0.36	-0.04	1.64	2.23	0.067
Workshop - Abb. Workshop	T2	-0.01	0.35	-0.83	0.80	-0.04	0.999
Besides mosquitoes, can people also spread dengue fever?
Abb. Workshop - Direct instruction	T1	-0.83	0.44	-1.86	0.20	-1.89	0.1422
Workshop - Direct instruction	T1	-1.70	0.42	-2.69	-0.71	-4.02	< .001
Workshop - Abb. Workshop	T1	-0.87	0.37	-1.74	0.01	-2.31	0.0543
Abb. Workshop - Direct instruction	T2	-0.56	0.41	-1.53	0.40	-1.36	0.3601
Workshop - Direct instruction	T2	-1.12	0.39	-2.04	-0.19	-2.83	0.013
Workshop - Abb. Workshop	T2	-0.55	0.38	-1.44	0.33	-1.47	0.3078
Before flying, is the mosquito a larva that lives in water?
Abb. Workshop - Direct instruction	T1	2.48	0.00	2.48	2.49	1,592.01	< .001
Workshop - Direct instruction	T1	0.45	0.67	-1.12	2.01	0.67	0.7823
Workshop - Abb. Workshop	T1	-2.04	0.67	-3.60	-0.47	-3.04	0.0066
Abb. Workshop - Direct instruction	T2	-3.83	1.13	-6.48	-1.17	-3.38	0.0021
Workshop - Direct instruction	T2	-0.05	1.46	-3.48	3.38	-0.03	0.9994
Workshop - Abb. Workshop	T2	3.78	1.26	0.83	6.73	3.00	0.0076
Is the aedes aegypti black with white spots and stripes on its body and legs?
Abb. Workshop - Direct instruction	T1	0.06	0.85	-1.93	2.05	0.07	0.9975
Workshop - Direct instruction	T1	0.03	0.00	0.03	0.04	10.64	< .001
Workshop - Abb. Workshop	T1	-0.02	0.85	-2.01	1.97	-0.03	0.9996
Abb. Workshop - Direct instruction	T2	0.74	0.77	-1.06	2.54	0.97	0.5985
Workshop - Direct instruction	T2	1.34	0.00	1.33	1.35	303.29	< .001
Workshop - Abb. Workshop	T2	0.60	0.77	-1.20	2.40	0.78	0.7135

btw_comparsons_flex %>%
  save_as_docx(path = "Tables/TableS5-Pairwise_between_questions.docx", pr_section = sect_properties)

btw_comparsons_flex_excel <- openxlsx2::wb_workbook()$add_worksheet("Table S5") %>% 
  wb_add_flextable("Table S5", btw_comparsons_flex)

btw_comparsons_flex_excel$save("Tables/TableS5-Pairwise_within_LMM_Score.xlsx")

8.3. Table S6. Linear Mixed Models within groups comparisons for significant questions

table_list = list()

for (i in 1:length(significant_questions)) {

  baseline_response_df = base_long_concept %>% 
    filter(Time == "Baseline") %>% 
    select(Caso, Pregunta, Group, Respuesta) %>% 
    rename(Baseline_response = "Respuesta")
  
  data.model <- base_long_concept %>%                             # Preparing data
    filter(Time != "Baseline") %>%
    filter(Pregunta == significant_questions[i]) %>% 
    left_join(baseline_response_df) %>% 
    filter(!is.na(Respuesta))%>% 
    filter(!is.na(Baseline_response))
  
  
  model =                                                      # Model fitting
    glmer( Respuesta ~ Time + Group + Baseline_response + 
            Baseline_response:Time + Time:Group + 
            (1|Caso),
          data= data.model,
          family=binomial,
          control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5))) 

#Between groups

table_list[[i]] = emmeans(model, revpairwise ~ Time|Group, p.adjust = "holm")$contrasts %>% 
    summary(infer = TRUE) %>% 
    as_tibble() %>% 
    rename(Contrast = contrast,
           Estimate = estimate,
           `Lower CL` = asymp.LCL,
           `Upper CL` = asymp.UCL,
           `Z ratio` = z.ratio,
           `P value` = p.value) %>% 
  mutate(Group = c("Direct instruction", 
                      "Abb. Workshop", 
                      "Workshop")) %>% 

    add_row(Contrast = significant_questions[i], .before = 1)                     ## Adding question text to table


}

wtn_commparisons = reduce(table_list, full_join) %>% 
    mutate(`P value` = round(x = `P value`, digits = 4)) %>%
    mutate(`P value` = ifelse(`P value` == 0, "< .001", `P value`)) %>%
    mutate(across(where(is.numeric), ~round(.x, digits = 2))) %>% 
  select(Contrast, Group, Estimate, SE, `Lower CL`, `Upper CL`, `Z ratio`, `P value`)
  
  
wtn_comparsons_flex = wtn_commparisons %>% 
    flextable() %>%
  
    merge_h_range(i = seq(1, nrow(wtn_commparisons), by = 4), j1 = 1, j2 = ncol(wtn_commparisons)) %>%
    set_caption(caption= as_paragraph(as_chunk("Table S5. Linear Mixed Models within groups comparisons for significant questions",
                  props = fp_text_default(bold=TRUE, color="black",
                                          font.size = 12,))), align_with_table = F)
wtn_comparsons_flex

Table S5. Linear Mixed Models within groups comparisons for significant questions

Contrast	Group	Estimate	SE	Lower CL	Upper CL	Z ratio	P value
Do mosquitoes bite more on the face and back than on the legs and arms?
T2 - T1	Direct instruction	-0.54	0.33	-1.19	0.11	-1.64	0.1014
T2 - T1	Abb. Workshop	-0.13	0.33	-0.77	0.52	-0.39	0.6977
T2 - T1	Workshop	-0.63	0.31	-1.23	-0.03	-2.06	0.0394
Besides mosquitoes, can people also spread dengue fever?
T2 - T1	Direct instruction	-1.09	0.43	-1.93	-0.26	-2.56	0.0105
T2 - T1	Abb. Workshop	-0.83	0.38	-1.56	-0.09	-2.19	0.0282
T2 - T1	Workshop	-0.51	0.32	-1.14	0.12	-1.60	0.1094
Before flying, is the mosquito a larva that lives in water?
T2 - T1	Direct instruction	2.90	0.97	1.00	4.80	3.00	0.0027
T2 - T1	Abb. Workshop	-3.41	0.59	-4.57	-2.24	-5.73	< .001
T2 - T1	Workshop	2.40	1.00	0.45	4.36	2.41	0.0157
Is the aedes aegypti black with white spots and stripes on its body and legs?
T2 - T1	Direct instruction	-1.48	0.00	-1.49	-1.47	-459.69	< .001
T2 - T1	Abb. Workshop	-0.79	0.83	-2.43	0.84	-0.95	0.3397
T2 - T1	Workshop	-0.17	0.00	-0.18	-0.16	-44.87	< .001

wtn_comparsons_flex %>%
  save_as_docx(path = "Tables/TableS6-Pairwise_within_questions.docx", pr_section = sect_properties)


wtn_comparsons_flex_excel <- openxlsx2::wb_workbook()$add_worksheet("Table S6") %>% 
  wb_add_flextable("Table S6", wtn_comparsons_flex)

wtn_comparsons_flex_excel$save("Tables/TableS6-Pairwise_within_questions.xlsx")

8.4. Figure 2. Percentage of correct responses for questions with statistically significant differences between or within groups.

plot_list = list()

for (i in 1:length(significant_questions)) {

      plot_list[[i]] <- Frecuencias_conceptos%>%  
        filter(Pregunta == significant_questions[i]) %>% 
        mutate(Time = factor(Time, levels = c("Baseline", "T1", "T2"))) %>% 
        mutate(p = Porcentaje/100,
               # n =n(),
               SE = sqrt((p * (1 - p)) / n)*100) %>% 

        ggplot(aes(x=Time, y=Porcentaje, group=Group, color=Group, fill= Group))+

        scale_fill_manual(values=exp1_cols)+
        scale_color_manual(values=exp1_cols)+
        scale_shape_manual(values = c(21, 25, 22,24)) +

        geom_line(lty=2,lwd=1, position=position_dodge(width = 0.2), alpha= .7)+
        geom_point(aes(shape= Group),
                   position=position_dodge(width = 0.2),size= 3)+
        geom_errorbar(aes(ymin=Porcentaje-SE, ymax=Porcentaje+SE),
                      position=position_dodge(width = 0.2),width=0,lwd=1.25)+
        annotate("text", label = paste(strwrap(significant_questions[i],
                                               width = 40),
                                       collapse = "\n"),
                 x = 0.6, y = 15, size = 4, colour = "black", hjust = 0)+
        labs(y="Percentage",x="Time") +
        theme_bw(base_size = 15) +
        coord_cartesian(ylim = c(10,100))

      
      }


Figure2 = ggarrange(plotlist = plot_list, common.legend = T,
                    legend = "bottom", ncol = 2, nrow = 2,
                    labels = c("A", "B", "C", "D"))

Figure2

# Figure2%>%
#   ggexport(filename = "Figures/png/Figure2.png", width = 720, height = 720)

Figure2 %>% 
  ggsave(filename = "Figures/png/Figure2.png", width = 5.6, height = 5.6, units = "in", dpi = 600)

Figure2 %>% 
  ggsave(filename = "Figures/tif/Figure2.tif", width = 8.6, height = 8.5, units = "in", dpi = 300, bg = "white")