03-Results.Rmd

---
output: html_document
editor_options: 
  chunk_output_type: console
---

# Results

In this chapter, the relevant findings from every step performed are analyzed following the analytical strategy indicated in the previous chapter. First, latent class analysis for each country included in the sample selected is performed in order to identify the optimal number of classes for each country separately. Secondly, a global analysis is performed to identify how many latent classes are identified including all different countries, with this information it is possible to establish the number of classes that can be compared across countries. Finally, a multigroup latent class analysis is performed considering all the previous information. The multigroup analysis is constructed in multiple steps, the most restricted model until the less restricted model is evaluated. As a final step, a confirmatory latent class analysis is performed using some theoretical hypotheses that were defined based on the previous results.  

This procedure is performed for the two scales that were used to create the Students' endorsement of equal rights and opportunities indicators separately.


```{r}
ISC_lvRlca <- data_model %>% 
  dplyr::select(all_of(Id), all_of(sample), 
                all_of(Scales) 
                #all_of(Schl_cate), 
                #all_of(Stud_cate), 
                #all_of(Scores)
                ) 

cnt <- unique(ISC_lvRlca$IDCNTRY)
cbPalette <- brewer.pal(n = 8, name = "Dark2")

```

```{r}
load("data/MplusModels_LCA.RData")

load("data/MplusModels_ByCountry.RData")

load("data/MplusModels_MGCntry.RData")

```

For every analysis in this chapter, the model fit statistics table includes all the statistics that are retrieved by MPLUS software. Here is a brief description of the meaning behind every column that will be shown. The first table with the model fit statistics for all different models indicates first the number of classes used in the model, the total number of parameters estimated, the final and best Log-likelihood, the values for information criteria AIC, BIC, aBIC. The entropy indicated in each table corresponds to the relative entropy, where a perfect classification is 1, the table also indicates the log likelihood reduction (LL Reduction) from adding one class into the model. Two tests for model fit are indicated as well, the value of the statistic and the p-value associated with the Vuong-Lo-Mendell-Rubin likelihood ratio test (VLMR) and Lo-Mendell-Rubin adjusted LRT Test (LMR).  

Conditional response probabilities plots are included as a summary for every independent model, the x-axis includes all the items included in the model, the y axis corresponds to the values of the probability to agree with that item, these values range from 0 to 1, where values close to 1 indicate that is highly likely to agree with that item, in contrast, values close to 0 indicate unlikely to agree to that item. On the other hand, values around 0.5 indicate randomness in the response, for this reason, is not possible to indicate a clear tendency to agree or to disagree. The different classes identified in the plots are colored differently but the colors remain the same when the response pattern is similar across countries. When a new class was identified a new color was used. The sample size for each class appears at the end of the x-axis colored with the same color as the class.

## Students' endorsement of gender equality scale

As mentioned in the previous chapter, this scale is composed of 6 items, in the following tables and plots, these items were ordered from positive to negative items for an easier interpretation of the results. This ordering consider first all the items that were positive worded in the instrument *Men and women should have equal opportunities to take part in government* (GND1), *Men and women should have the same rights in every way* (GND2) and *Men and women should get equal pay when they are doing the same jobs* (GND5), followed by the three other items that are negatively worded *Women should stay out of politics (r)* (GND3), *Not many jobs available, men should have more right to a job than women (r)* (GND4), *Men are better qualified to be political leaders than women (r)* (GND6).  As mentioned before all these variables were recoded in two categories, Agree and Disagree. All 14 countries were analyzed independently and then pooled in the same dataset.  

### Analysis by country  

Multiple latent class models with 1 to 6-classes\footnote{Summary with all models can be found in Appendix, table @ref(tab:detailed1).} were performed in each country in order to evaluate the model fit of each one of them. The results are summarized in table \@ref(tab:summodelfitcntry1). In most European countries, the best model fit based on the different criteria indicated previously are by including 3 or 4 latent classes.  

For Belgium, Croatia, Denmark, Latvia, and The Netherlands there is no significant improvement in the log-likelihood from two to three latent classes. In this sense, BIC and aBIC simultaneously have the lowest values in the 3-class model.  

On the other hand, in Bulgaria, Estonia, Malta, Slovenia, and Sweden according to the statistical tests, BIC, and aBIC criteria, the best model is a 4-class model. In Finland, Italy and Lithuania models, the BIC, aBIC differ from the statistical test indicating a better fit for the 3-class model.  

Norway is the only country from the sample where the best model fit is the one with 5 latent classes according to the statistical tests and BIC and aBIC.    

It is a common tendency in all the evaluated countries that the AIC value is lower in the models with one more class than indicated by the statistical tests and BIC and aBIC statistics. That is consistent with the indication that this criterion tends to overfit the data.  

Values of Entropy are higher when the tests are significant but consistent with a better fit of the data, the lower entropy found in the 4-class model is in Latvia (73.7\%) and the highest value in Norway (96\%). The log-likelihood reduction is consistent in all countries, where having more than 3 latent classes reduces the log-likelihood around  0.2\% and 1\%.   

```{r modelfitcnty1, echo=FALSE, results='hide'}
#check order of country names
#resultsbyall %>% arrange(Country) %>% select(Country) %>% unique()

ModelfitByContry("ByCountry_GND", title = "Model fit statistics LCA by country Students' endorsement of gender equality", fontn = 10) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 10) %>% 
  column_spec(c(2,3,11), width = "3em") %>% 
  column_spec(c(4:10,12), width = "4em") %>% 
  column_spec(c(1), width = "6em") %>%
  pack_rows("Belgium (Flemish)",1,6) %>%
#  row_spec(2:3, background = "lightgray") %>% 
  row_spec(c(3), bold = TRUE) %>%
  pack_rows("Bulgaria",7,12) %>%
#  row_spec(9:10, background = "lightgray") %>% 
  row_spec(c(10), bold = TRUE) %>%
  pack_rows("Croatia",13,18) %>%
#  row_spec(14:15, background = "lightgray") %>% 
  row_spec(c(16), bold = TRUE) %>%
  pack_rows("Denmark",19,24) %>%
#  row_spec(21:22, background = "lightgray") %>% 
  row_spec(c(22), bold = TRUE) %>%
  pack_rows("Estonia",25,30) %>%
#  row_spec(27:28, background = "lightgray") %>% 
  row_spec(c(27), bold = TRUE) %>%
  pack_rows("Finland",31,36) %>%
#  row_spec(32:33, background = "lightgray") %>% 
  row_spec(c(33), bold = TRUE) %>%
  pack_rows("Italy",37,42) %>%
#  row_spec(39:40, background = "lightgray") %>% 
  row_spec(c(39), bold = TRUE) %>%
  pack_rows("Latvia",43,48) %>%
#  row_spec(45:46, background = "lightgray") %>% 
  row_spec(c(45), bold = TRUE) %>%
  pack_rows("Lithuania",49,54) %>%
#  row_spec(51:52, background = "lightgray") %>% 
  row_spec(c(51), bold = TRUE) %>%
  pack_rows("Malta",55,60) %>%
 # row_spec(57:58, background = "lightgray") %>% 
  row_spec(c(58), bold = TRUE) %>%
  pack_rows("The Netherlands",61,66) %>%
#  row_spec(62:63, background = "lightgray") %>% 
  row_spec(c(63), bold = TRUE) %>%
  pack_rows("Norway",67,72) %>%
#  row_spec(70:71, background = "lightgray") %>% 
  row_spec(c(71), bold = TRUE) %>%
  pack_rows("Slovenia",73,78) %>%
#  row_spec(75:76, background = "lightgray") %>% 
  row_spec(c(76), bold = TRUE) %>%
  pack_rows("Sweden",79,84) %>%
#  row_spec(81:82, background = "lightgray") %>% 
  row_spec(c(82), bold = TRUE) %>%
  collapse_rows(1, valign = "top") %>% 
  footnote(general = "The best loglikelihood value was not replicated for the following models: ",
           number = "Croatia, 6 classes model") %>% 
  print()
cat('\n')
cat('\n')
```

All models selected accomplish at least one or more of the criteria established for a good fit. The bivariate residuals were also analyzed, and all countries have residuals around the range of acceptable [-2 ; 2] as shown in the figure \@ref(fig:resid1cnt). There is just one value is outside the ranges in Malta with a 4-class model.  

\blandscape
```{r summodelfitcntry1}
ModelfitByContry("ByCountry_GND", 
                 title = "Best model, fit statistics individual country model Students' endorsement of gender equality", 
                 filterval = TRUE, fontn = 10) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 10) %>% 
  column_spec(c(2,3,9:10,12), width = "3em") %>% 
  column_spec(c(4:8,11,13), width = "4em") %>% 
  column_spec(c(1), width = "8em") %>% 
  footnote(general = "Best model based on the lowest value of BIC") %>% 
  print()
cat('\n')
cat('\n')
```
\elandscape


```{r resid1cnt, fig.width=6, fig.height=5, fig.cap="Bivariate standardized residuals individual country models for Students' endorsement of gender equality", fig.pos='H'}

residuals_ByCntryGND <-  data.frame(`BFL 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_BFL_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `BGR 4cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_BGR_C3cl4.out$tech10$bivar_model_fit_info$z,
                            `DNK 4cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_DNK_C3cl4.out$tech10$bivar_model_fit_info$z,
                            `EST 4cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_EST_C3cl4.out$tech10$bivar_model_fit_info$z,
                            `FIN 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_FIN_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `HRV 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_HRV_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `ITA 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_ITA_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `LTU 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_LTU_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `LVA 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_LVA_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `MLT 4cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_MLT_C3cl4.out$tech10$bivar_model_fit_info$z,
                            `NLD 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_NLD_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `NOR 5cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_NOR_C3cl5.out$tech10$bivar_model_fit_info$z,
                            `SVN 4cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_SVN_C3cl4.out$tech10$bivar_model_fit_info$z,
                            `SWE 4cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_SWE_C3cl4.out$tech10$bivar_model_fit_info$z)
residuals_ByCntryGND <- residuals_ByCntryGND %>% mutate(x = 1:nrow(residuals_ByCntryGND)) %>%
  reshape2::melt(id.vars = c("x"))

pByCntryGND <- residuals_ByCntryGND %>% ggplot() +
  geom_point(aes(x = x, y = value), size = 1) +
  geom_hline(yintercept=c(-1.96,1.96), linetype = "dashed", size = 0.9, color = "black") +
  scale_fill_grey() + theme_bw() +
  facet_wrap(variable~., scales = "free_y") +
  #ggtitle("Bivariate residuals for individual best country model fit") +
  labs(y="Standardized residuals")  +
  theme(legend.position = "none", legend.box="vertical",
        strip.text.y = element_text(size = 8),
        legend.spacing.y = unit(-0.2, 'cm'),
        title = element_text(size = 10),
        axis.title = element_text(size = 10),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 9),
        axis.text.y = element_text(size = 9)) 
pByCntryGND
cat("\n")
cat("\n")

```

In figure \@ref(fig:classes1cnt), the classes of each independent model can be identified by looking at the conditional probabilities. In the figure, the conditional probabilities to agree to each item are shown and plotted for each class estimated in each country. From all the models, two classes that are similar across countries are identified in the figure, Fully egalitarian and Competition-driven sexism, green and purple line respectively.  

A brief explanation of each class is described below.  

- **Fully egalitarian:** Most likely to agree to all items (green line).  This class can be observed in all countries. Conditional probabilities greater than 0.75 to agree, class sizes around 60% (Bulgaria) and 90% (Denmark).   
  
- **Competition-driven sexism:** Most likely to disagree with gender competitive items in favor of women (purple line). This class can be observed in all countries. Conditional probabilities greater than 0.75 to agree to positive views of gender equality and generally lower than 0.5 to agree to reversed negative views, class sizes around 3.6% (Denmark) and 22.5% (Bulgaria).  
  
- **Non-egalitarian:** Not likely to agree to any item (orange line). This class can be observed in four countries. Conditional probabilities lower to 0.5 to agree to any item, class sizes around 0.9% (Norway) and 2.6% (Italy).  

- **Reverse competition-driven sexism:** Most likely to agree to gender competitive items in favor of women (pink line) . This class can be observed in five countries. Conditional probabilities lower than 0.25 to agree to positive views of gender equality and generally greater than 0.75 to agree to reversed negative views, class sizes around 0.6% (Norway) and 1.6% (The Netherlands).  

- **Political egalitarian:** Likely to agree to politically related items (light-green line). This class can be observed in five countries. Conditional probabilities are greater than 0.75 in political equality items, class sizes around 3.2% (Belgium) and 1.4% (Estonia).    

- **Random response:** Not defined attitude (yellow line). This class can be observed in four countries. Conditional probabilities between 0.25 and 0.75 to agree all items, class sizes around 2.7% (Slovenia) and 16.8% (Bulgaria).  

The classes described before are not present in all countries, for this reason, a global model will be tested in the pooled sample considering not only the classes that are similar across more than one country, but additional classes will be added in order to absorb the remaining different classes that the global model will identify.  

In the following section, the global model will be tested using the pooled dataset.  

  
```{r classes1cnt, fig.width=6, fig.height=5, fig.cap="Classes for best individual country model for Students' endorsement of gender equality", fig.pos='H'}

classes_ByCntryGND <-  rbind(data.frame(Country = "Belgium (Flemish)", 
                                      left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_BFL_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_BFL_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param"))) %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 5)),
                            data.frame(Country = "Bulgaria",
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_BGR_C3cl4.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_BGR_C3cl4.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param"))) %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 6,
                                                              LatentClass == 3 ~ 4,
                                                              LatentClass == 4 ~ 1)),
                            data.frame(Country = "Denmark", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_DNK_C3cl4.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_DNK_C3cl4.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 5,
                                                              LatentClass == 3 ~ 1,
                                                              LatentClass == 4 ~ 3)),
                            data.frame(Country = "Estonia", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_EST_C3cl4.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_EST_C3cl4.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 5,
                                                              LatentClass == 2 ~ 2,
                                                              LatentClass == 3 ~ 1,
                                                              LatentClass == 4 ~ 3)),
                            data.frame(Country = "Finland", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_FIN_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_FIN_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Croatia", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_HRV_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_HRV_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 2,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Italy", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_ITA_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_ITA_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Lithuania", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_LTU_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_LTU_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 6)),
                            data.frame(Country = "Latvia", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_LVA_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_LVA_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 6)),
                            data.frame(Country = "Malta", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_MLT_C3cl4.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_MLT_C3cl4.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 4,
                                                              LatentClass == 2 ~ 6,
                                                              LatentClass == 3 ~ 3,
                                                              LatentClass == 4 ~ 1)),
                            data.frame(Country = "Netherlands", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_NLD_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_NLD_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 1,
                                                              LatentClass == 2 ~ 4,
                                                              LatentClass == 3 ~ 3)),
                            data.frame(Country = "Norway", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_NOR_C3cl5.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_NOR_C3cl5.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 4,
                                                              LatentClass == 3 ~ 1,
                                                              LatentClass == 4 ~ 2,
                                                              LatentClass == 5 ~ 5)),
                            data.frame(Country = "Slovenia", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_SVN_C3cl4.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_SVN_C3cl4.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 4,
                                                              LatentClass == 4 ~ 1)),
                            data.frame(Country = "Sweden", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_SWE_C3cl4.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_SWE_C3cl4.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 5,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 2,
                                                              LatentClass == 4 ~ 3)))

pByCntryGND <- classes_ByCntryGND %>% filter(category == 1) %>% plyr::rbind.fill(data.frame(param = "Size", Country = "Sweden")) %>% 
  mutate(LatentClass = factor(LatentClass),
         param = factor(param, levels = c(unique(classes_ByCntryGND$param), "Size")[c(1,2,5,3,4,6,7)])) %>% 
  ggplot() +
  geom_point(aes(x = param, y = est, group = LatentClass, color = LatentClass), size = 1) +
  geom_line(aes(param, est, group = LatentClass, linetype = LatentClass, color = LatentClass)) +
  geom_text(aes(x = param, y = est, label = scales::percent(proportion, accuracy = 0.1), color = LatentClass), size = 2, 
            nudge_x = 0.8, nudge_y = 0) +
  facet_wrap(Country ~ .) +
  scale_fill_grey() + theme_bw() +
  #ggtitle("Classes for individual country models for attitudes towards gender \nequality scale") +
  theme(legend.position = "none", 
        title = element_text(size = 9),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 7),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.text.x = element_text(angle = 90, size = 7, vjust = 0.5, hjust = 0)) +
  scale_y_continuous(breaks = c(0.25,0.5,0.75), limits = c(0,1)) +
  geom_hline(yintercept = c(0.25,0.5,0.75), color = "gray", size = 0.2) +
  labs(y="Response probabilities", linetype = "Latent Classes", color = "Latent Classes")  +
  scale_color_brewer(type = "qual", palette = "Dark2") 

pByCntryGND
cat("\n")
cat("\n")

```


```{r resid1cnt2, fig.width=6, fig.height=5, fig.cap="Bivariate model fit standardized residuals", fig.pos='H', eval=FALSE}

residuals_ByCntryGND <-  data.frame(`BFL 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_BFL_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `BGR 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_BGR_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `DNK 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_DNK_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `EST 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_EST_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `FIN 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_FIN_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `HRV 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_HRV_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `ITA 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_ITA_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `LTU 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_LTU_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `LVA 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_LVA_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `MLT 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_MLT_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `NLD 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_NLD_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `NOR 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_NOR_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `SVN 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_SVN_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `SWE 3cl` = ByCountry_GND$data.MplusModels.ByCountry.GNDlca_SWE_C3cl3.out$tech10$bivar_model_fit_info$z)
residuals_ByCntryGND <- residuals_ByCntryGND %>% mutate(x = 1:nrow(residuals_ByCntryGND)) %>%
  reshape2::melt(id.vars = c("x"))

pByCntryGND <- residuals_ByCntryGND %>% ggplot() +
  geom_point(aes(x = x, y = value), size = 1) +
  geom_hline(yintercept=c(-1.96,1.96), linetype = "dashed", size = 0.9, color = "black") +
  scale_fill_grey() + theme_bw() +
  facet_wrap(variable~., scales = "free_y") +
  ggtitle("Residuals by country models") +
  labs(y="Standardized residuals")  +
  theme(legend.position = "none", legend.box="vertical",
        strip.text.y = element_text(size = 8),
        legend.spacing.y = unit(-0.2, 'cm'),
        title = element_text(size = 9),
        axis.title = element_text(size = 8),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 8),
        axis.text.y = element_text(size = 8)) 
pByCntryGND
cat("\n")
cat("\n")

```


```{r classes1cnt2, fig.width=6, fig.height=5, fig.cap="By country 3-Classes model", fig.pos='H', eval=FALSE}

classes_ByCntryGND <-  rbind(data.frame(Country = "Belgium (Flemish)", 
                                      left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_BFL_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_BFL_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param"))) %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 2)),
                            data.frame(Country = "Bulgaria",
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_BGR_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_BGR_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param"))) %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 3)),
                            data.frame(Country = "Denmark", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_DNK_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_DNK_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 1,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 2)),
                            data.frame(Country = "Estonia", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_EST_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_EST_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Finland", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_FIN_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_FIN_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Croatia", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_HRV_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_HRV_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 2,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Italy", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_ITA_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_ITA_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Lithuania", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_LTU_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_LTU_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 2)),
                            data.frame(Country = "Latvia", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_LVA_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_LVA_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 2)),
                            data.frame(Country = "Malta", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_MLT_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_MLT_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Netherlands", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_NLD_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_NLD_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 1,
                                                              LatentClass == 2 ~ 2,
                                                              LatentClass == 3 ~ 3)),
                            data.frame(Country = "Norway", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_NOR_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_NOR_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 1,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 2)),
                            data.frame(Country = "Slovenia", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_SVN_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_SVN_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 1,
                                                              LatentClass == 2 ~ 2,
                                                              LatentClass == 3 ~ 3)),
                            data.frame(Country = "Sweden", 
                                       left_join(ByCountry_GND$data.MplusModels.ByCountry.GNDlca_SWE_C3cl3.out$parameters$probability.scale,
                                                ByCountry_GND$data.MplusModels.ByCountry.GNDlca_SWE_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "GND6") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 3)))

pByCntryGND <- classes_ByCntryGND %>% filter(category == 1) %>% 
  mutate(LatentClass = factor(LatentClass),
         param = factor(param, levels = unique(classes_ByCntryGND$param)[c(1,2,5,3,4,6)])) %>% 
  ggplot() +
  geom_point(aes(x = param, y = est, group = LatentClass, color = LatentClass), size = 1) +
  geom_line(aes(param, est, group = LatentClass, linetype = LatentClass, color = LatentClass)) +
  geom_text(aes(x = param, y = est, label = scales::percent(proportion, accuracy = 0.1), color = LatentClass), size = 2, 
            nudge_x = -0.15, nudge_y = 0.1) +
  facet_wrap(Country ~ .) +
  scale_fill_grey() + theme_bw() +
  ggtitle("Classes by country for Students' endorsement of gender equality") +
  theme(legend.position = "none", 
        title = element_text(size = 9),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 7),
        axis.text.x = element_text(angle = 90, size = 7, vjust = 0.5, hjust = 0)) +
  scale_y_continuous(breaks = c(0.25,0.5,0.75)) +
  labs(y="Response probabilities", linetype = "Latent Classes", color = "Latent Classes")  +
  scale_color_brewer(type = "qual", palette = "Dark2") 
pByCntryGND
cat("\n")
cat("\n")

```


### General model  

Table \@ref(tab:modelfitlca1) shows the results of each model using the pooled sample with all the countries. Models with 1 to 7 classes were computed and the model fit statistics were summarized in the table.  

The model that includes a single class has the largest AIC (192,838), BIC (192,891), and ABIC (192,872) values for the pooled sample, indicating that this model fits data worse than all other models. In addition, the P-values for the VLMR test, and LMR in the 2-class model are all < 0.0001; this means that both tests reject the single-class model in favor of a model with at least two latent classes. In other words, there exists heterogeneity in the target population regarding attitudes towards gender equality.  

On the other hand, the LMR LR and VLMR tests for the 6-class model are not statistically significant (P > 0.05). That is, the two tests are in favor of at most 5 classes.  

In contrast, BIC and aBIC values are all smaller in the 5-class model than those in the 6-class model; thus, consider that the models with more than 5 classes are not preferred. AIC values reach the lowest value in the 7-class model but based on the previous results this criterion will not be considered in this case due to the tendency to overfit the data.  

The relative entropy given by Mplus software decreases when including more than 4 classes and increases again with the 6-class model; this would suggest that a model with at least 6 class or 4 class is preferred.  

Together with the percentage of reduction in the log-likelihood value, that indicates that by adding two classes to the model the log-likelihood is reduced by 13.3%, this reduction is only increased by 1.2% if the model is a 3-class model and finally this value is reduced close to 0 if more than 5 classes are included.  

Now, the preferred model must be either the 5-class or the 6-class model. Considering the residuals of each model, in figure \@ref(fig:resid1) all values are around -1.96 and 1.96. But based on the parsimony principle a 4-class model can be considered as well if just one value of the residuals is outside the acceptable range.  

Theoretically, we tend to determine that the 4-class LCA model is the preferred model. We will show later that the classes identified by the 4-class model are more interpretable and representative than the rest of the models. And in particular that two classes can be compared across countries.    


\blandscape  
```{r modelfitlca1, echo=FALSE}

Modelfit("lcaGND", title = "Model fit statistics LCA Students' endorsement of gender equality", fontn = 11) %>% 
  pack_rows("All countries",1,4) %>% 
  #row_spec(4:5, background = "lightgray") %>% 
  row_spec(c(4,5), bold = TRUE) %>% 
  print()
cat('\n')
cat('\n')
```
\elandscape

```{r resid1, fig.width=5, fig.height=3, fig.cap="Bivariate model fit standardized residuals global model for Students' endorsement of gender equality", fig.pos='H'}

residuals_GND <-  data.frame(cl1 = lcaGND$GND_lca_C3cl1.out$tech10$bivar_model_fit_info$z,
                            cl2 = lcaGND$GND_lca_C3cl2.out$tech10$bivar_model_fit_info$z,
                            cl3 = lcaGND$GND_lca_C3cl3.out$tech10$bivar_model_fit_info$z,
                            cl4 = lcaGND$GND_lca_C3cl4.out$tech10$bivar_model_fit_info$z,
                            cl5 = lcaGND$GND_lca_C3cl5.out$tech10$bivar_model_fit_info$z,
                            cl6 = lcaGND$GND_lca_C3cl6.out$tech10$bivar_model_fit_info$z)
residuals_GND <- residuals_GND %>% mutate(x = 1:nrow(residuals_GND)) %>% 
  reshape2::melt(id.vars = c("x"))

pGND <- residuals_GND %>% ggplot() + 
  geom_point(aes(x = x, y = value), size = 1) +
  geom_hline(yintercept=c(-1.96,1.96), linetype = "dashed", size = 0.9, color = "black") +
  scale_fill_grey() + theme_bw() +
  facet_wrap(variable~., scales = "free_y") +
  #ggtitle("Standardized bivariate residuals for \nAttitudes towards gender equality models") +
  labs(x = "Parameters", y="Standardized residuals", color = "Number of classes")  +
  theme(legend.position = "none",
        strip.text.y = element_text(size = 8),
        legend.spacing.y = unit(-0.2, 'cm'),
        title = element_text(size = 9),
        axis.title = element_text(size = 8), 
        axis.text.x = element_blank(), 
        axis.text.y = element_text(size = 8)) +
  scale_color_brewer(type = "qual", palette = "Dark2")
pGND
cat("\n")
cat("\n")

```

```{r, echo=FALSE, results='hide'}
#----GND 3 groups----
classes3GND <- c("Fully egalitarian",
                "Competition- driven sexism",
                "Random response")
orden3GND <- c(2,3,1)
lcaGND_C3cl3 <- lcaGND$GND_lca_C3cl3.out$parameters$probability.scale %>% 
  rename_with(~ c("Class", "value")[which(c("LatentClass", "est") == .x)], .cols = c("LatentClass", "est")) %>% 
  mutate_at( c("param", "category", "Class"), ~ as.factor(.x)) %>% 
  mutate(Class = factor(Class, levels = orden3GND, labels = classes3GND))

counts3GND <- full_join(lcaGND$GND_lca_C3cl3.out$class_counts$modelEstimated,
                        lcaGND$GND_lca_C3cl3.out$class_counts$mostLikely,by = c("class"))
counts3GND  %>% 
  mutate(class = factor(class, levels = orden3GND, labels = classes3GND),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes 3-class Gender equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",3)), row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))

sizelca3_GND <- lcaGND$GND_lca_C3cl3.out$class_counts$modelEstimated %>% dplyr::select(-count) %>% 
  rename_with(~ c("Gender", "Class")[which(c("proportion", "class") == .x)], .cols = c("proportion", "class")) %>% 
  mutate(Class = factor(Class, levels = orden3GND, labels = classes3GND)) %>% 
  reshape2::melt(id.vars = c("Class"), variable.name = "Group") %>% 
  dplyr::arrange(Group) %>% 
  dplyr::group_by(Group) %>%
  dplyr::mutate(countT= sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(Class) %>%
  dplyr::mutate(per=value/countT) %>% 
  dplyr::select(Group, Class, per) 

VarClass(lcaGND_C3cl3) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.8, "Myblue", 
                                                 ifelse(value < 0.8 & value > 0.49, "Mygreen","Myred")))) %>% 
  reshape2::dcast(param ~ Class) %>% 
  kbl(caption = paste0("Probabilities to agree each item 3-class Gender equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",3)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:4, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
# graphclass(lcaGND_C3cl3, nclass = 3, orden = c(1,2,5,3,4,6), title = "LCA Gender equality with 3 classes", mg = FALSE)
# cat("\n")
# cat("\n")

#HighProb(lcaGND_C3cl3, sizelca3_GND,  orden = c(1,2,5,3,4,6), title = "Response categories probabilities and class size for\n #3-classes Gender equality model") 
cat("\n")
cat("\n")

#----GND 4 groups----
classes4GND <- c("Fully egalitarian",
                "Competition- driven sexism",
                "Non-egalitarian",
                "Political egalitarian")
orden4GND <- c(2,4,3,1)
lcaGND_C3cl4 <- lcaGND$GND_lca_C3cl4.out$parameters$probability.scale %>% 
  rename_with(~ c("Class", "value")[which(c("LatentClass", "est") == .x)], .cols = c("LatentClass", "est")) %>% 
  mutate_at( c("param", "category", "Class"), ~ as.factor(.x)) %>% 
  mutate(Class = factor(Class, levels = orden4GND, labels = classes4GND))

counts4GND <- full_join(lcaGND$GND_lca_C3cl4.out$class_counts$modelEstimated,
                        lcaGND$GND_lca_C3cl4.out$class_counts$mostLikely,by = c("class"))
counts4GND  %>% 
  mutate(class = factor(class, levels = orden4GND, labels = classes4GND),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes 4-class Gender equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))

sizelca4_GND <- lcaGND$GND_lca_C3cl4.out$class_counts$modelEstimated %>% dplyr::select(-count) %>% 
  rename_with(~ c("Gender", "Class")[which(c("proportion", "class") == .x)], .cols = c("proportion", "class")) %>% 
  mutate(Class = factor(Class, levels = orden4GND, labels = classes4GND)) %>% 
  reshape2::melt(id.vars = c("Class"), variable.name = "Group") %>% 
  dplyr::arrange(Group) %>% 
  dplyr::group_by(Group) %>%
  dplyr::mutate(countT= sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(Class) %>%
  dplyr::mutate(per=value/countT) %>% 
  dplyr::select(Group, Class, per) 

VarClass(lcaGND_C3cl4) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.8, "Myblue", 
                                                 ifelse(value < 0.8 & value > 0.49, "Mygreen","Myred")))) %>% 
  reshape2::dcast(param ~ Class) %>% 
  kbl(caption = paste0("Probabilities to agree each item 4-class Gender equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:5, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
# graphclass(lcaGND_C3cl4, nclass = 4, orden = c(1,2,5,3,4,6), title = "LCA Gender equality with 4 classes", mg = FALSE)
# cat("\n")
# cat("\n")

#HighProb(lcaGND_C3cl4, sizelca4_GND,  orden = c(1,2,5,3,4,6), title = "Response categories probabilities and class size for\n #4-classes Gender equality model") 
#cat("\n")
#cat("\n")

#----GND 5 groups----
classes5GND <- c("Fully egalitarian",
                "Competition- driven sexism",
                "Non-egalitarian",
                "Political egalitarian",
                "Reverse competition-driven sexism")
orden5GND <- c(2,3,5,4,1)
lcaGND_C3cl5 <- lcaGND$GND_lca_C3cl5.out$parameters$probability.scale %>% 
  rename_with(~ c("Class", "value")[which(c("LatentClass", "est") == .x)], .cols = c("LatentClass", "est")) %>% 
  mutate_at( c("param", "category", "Class"), ~ as.factor(.x)) %>% 
  mutate(Class = factor(Class, levels = orden5GND, labels = classes5GND))

counts5GND <- full_join(lcaGND$GND_lca_C3cl5.out$class_counts$modelEstimated,
                        lcaGND$GND_lca_C3cl5.out$class_counts$mostLikely,by = c("class"))
counts5GND  %>% 
  mutate(class = factor(class, levels = orden5GND, labels = classes5GND),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes 5-class Gender equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))

sizelca5_GND <- lcaGND$GND_lca_C3cl5.out$class_counts$modelEstimated %>% dplyr::select(-count) %>% 
  rename_with(~ c("Gender", "Class")[which(c("proportion", "class") == .x)], .cols = c("proportion", "class")) %>% 
  mutate(Class = factor(Class, levels = orden5GND, labels = classes5GND)) %>% 
  reshape2::melt(id.vars = c("Class"), variable.name = "Group") %>% 
  dplyr::arrange(Group) %>% 
  dplyr::group_by(Group) %>%
  dplyr::mutate(countT= sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(Class) %>%
  dplyr::mutate(per=value/countT) %>% 
  dplyr::select(Group, Class, per) 

VarClass(lcaGND_C3cl5) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.8, "Myblue", 
                                                 ifelse(value < 0.8 & value > 0.49, "Mygreen","Myred")))) %>% 
  reshape2::dcast(param ~ Class) %>% 
  kbl(caption = paste0("Probabilities to agree each item 5-class Gender equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",5)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:6, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
# graphclass(lcaGND_C3cl5, nclass = 5, orden = c(1,2,5,3,4,6), title = "LCA Gender equality with 5 classes", mg = FALSE)
# cat("\n")
# cat("\n")

#HighProb(lcaGND_C3cl5, sizelca5_GND,  orden = c(1,2,5,3,4,6), title = "Response categories probabilities and class size for\n #5-classes Gender equality model") 
#cat("\n")
#cat("\n")


#----GND 6 groups----
classes6GND <- c("Fully egalitarian",
                "Competition- driven sexism",
                "Non-egalitarian",
                "Political egalitarian",
                "Reverse competition sexism",
                "Pro-women pay/job")
orden6GND <- c(4,5,2,6,1,3)
lcaGND_C3cl6 <- lcaGND$GND_lca_C3cl6.out$parameters$probability.scale %>% 
  rename_with(~ c("Class", "value")[which(c("LatentClass", "est") == .x)], .cols = c("LatentClass", "est")) %>% 
  mutate_at( c("param", "category", "Class"), ~ as.factor(.x)) %>% 
  mutate(Class = factor(Class, levels = orden6GND, labels = classes6GND))

counts6GND <- full_join(lcaGND$GND_lca_C3cl6.out$class_counts$modelEstimated,
                        lcaGND$GND_lca_C3cl6.out$class_counts$mostLikely,by = c("class"))
counts6GND  %>% 
  mutate(class = factor(class, levels = orden6GND, labels = classes6GND),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes 6-class Gender equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))

sizelca6_GND <- lcaGND$GND_lca_C3cl6.out$class_counts$modelEstimated %>% dplyr::select(-count) %>% 
  rename_with(~ c("Gender", "Class")[which(c("proportion", "class") == .x)], .cols = c("proportion", "class")) %>% 
  mutate(Class = factor(Class, levels = orden6GND, labels = classes6GND)) %>% 
  reshape2::melt(id.vars = c("Class"), variable.name = "Group") %>% 
  dplyr::arrange(Group) %>% 
  dplyr::group_by(Group) %>%
  dplyr::mutate(countT= sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(Class) %>%
  dplyr::mutate(per=value/countT) %>% 
  dplyr::select(Group, Class, per) 

VarClass(lcaGND_C3cl6) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.8, "Myblue", 
                                                 ifelse(value < 0.8 & value > 0.49, "Mygreen","Myred")))) %>% 
  reshape2::dcast(param ~ Class) %>% 
  kbl(caption = paste0("Probabilities to agree each item 6-class Gender equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",6)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:7, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
# graphclass(lcaGND_C3cl6, nclass = 6, orden = c(1,2,5,3,4,6), title = "LCA Gender equality with 6 classes", mg = FALSE)
# cat("\n")
# cat("\n")

#HighProb(lcaGND_C3cl6, sizelca6_GND,  orden = c(1,2,5,3,4,6), title = "Response categories probabilities and class size for\n #6-classes Gender equality model") 
cat("\n")
cat("\n")
```


Three, four, five and six class models were investigated profoundly. It is difficult to choose the best model fit without a full analysis. There are some patterns that can be clearly identified in the models as can be seen in figure \@ref(fig:compare1). Class 1 with class sizes of 81.5%, 79.2%, 77% and 78.6% in each model respectively, and estimated probabilities to agree for this latent class higher than 0.92 for all six items *Men and women should have equal opportunities to take part in government*, *Men and women should have the same rights in every way*, *Men and women should get equal pay when they are doing the same job*, *Women should stay out of politics*, *Not many jobs available, men should have more right to a job than women* and *Men are better qualified to be political leaders than women* the **Fully egalitarian** group.  

The second class (Class 2 in Figure \@ref(fig:compare1)) identified in all the models, called **Competition-driven sexism**. For this class, the estimated probabilities to agree to the first 3 items are close or higher than 0.9 in all models. For the last three items, the estimated probabilities of agreement are not higher than 0.5 in all models. The class size differs in all four models, 11.3%, 11.5%, 12.8% and 8.6% in the 3, 4, 5, 6-class model respectively.  

The third class (Class 3) that can be seen with a similar pattern in all the models, it is called **Non-egalitarian**, this class appears in the 4-class model. The pattern of this class is basically showing lower estimated conditional probabilities to agree with any of these statements, no greater than 0.4, except for one item *Men and women should get equal pay when they are doing the same jobs* with an estimated probability to agree no higher than 0.55. The estimated sizes for this class are 7.7%, 7% and 5.8% in each model respectively.     

The fourth, fifth and sixth classes identified in the models differ in all the countries, nevertheless, one class appears to be consistent in the 5-class model, where this class called **Reverse competition-driven sexism** has opposite conditional probabilities compared to the second class identified previously.  

```{r compare1, fig.width=6, fig.height=6, fig.cap="Comparative conditional probabilities to agree in 3 to 6 latent class global models for Students' endorsement of gender equality", fig.pos='H'}
g0 <- ClassGraph(lcaGND_C3cl3, sizelca3_GND, orden = c(1,2,5,3,4,6), title = "3-classes", selected = c(1,3,6))
g1 <- ClassGraph(lcaGND_C3cl4, sizelca4_GND, orden = c(1,2,5,3,4,6), title = "4-classes", selected = c(1,3,2,5))
g2 <- ClassGraph(lcaGND_C3cl5, sizelca5_GND, orden = c(1,2,5,3,4,6), title = "5-classes", selected = c(1,3,2,5,4))
g3 <- ClassGraph(lcaGND_C3cl6, sizelca6_GND, orden = c(1,2,5,3,4,6), title = "6-classes", selected = c(1,3,2,7,4,5))

g <- arrangeGrob(g0,g1,g2,g3, ncol = 2)#, top = "Conditional probabilities to agree to Students' endorsement of gender equality")
grid.draw(g)
cat("\n")
cat("\n")

```

The main two classes in the solutions with five and six classes do not strongly differ from other models, and the remaining classes are not informative at all or with a sample size very small, using this as a criterion, one can prefer a four-class solution.  

In table \@ref(tab:bestfit11) the conditional probabilities to agree with the fourth latent class model are shown. These values are very close to 1 for the first class, Fully egalitarian. Similar values are obtained for the positive items in the second class Competition driven-sexism, meanwhile, for third item GND3, conditional probabilities are close to 0.5, this would mean a random response, but the last two items have lower conditional probabilities close to 0.2, which would mean not likely to agree to the statements. Table \@ref(tab:bestfit12) indicates the counts and proportions using the model estimated and the most likely probabilities.  


```{r bestfit11}
#----GND 4 groups----
classes4GND <- c("Fully egalitarian",
                "Competition- driven sexism",
                "Non-egalitarian",
                "Political egalitarian")
orden4GND <- c(2,4,3,1)
lcaGND_C3cl4 <- lcaGND$GND_lca_C3cl4.out$parameters$probability.scale %>% 
  rename_with(~ c("Class", "value")[which(c("LatentClass", "est") == .x)], 
              .cols = c("LatentClass", "est")) %>% 
  mutate_at( c("param", "category", "Class"), ~ as.factor(.x)) %>% 
  mutate(Class = factor(Class, levels = orden4GND, labels = classes4GND))

counts4GND <- full_join(lcaGND$GND_lca_C3cl4.out$class_counts$modelEstimated,
                        lcaGND$GND_lca_C3cl4.out$class_counts$mostLikely,
                        by = c("class"))

lcaGND_C3cl4$orden = rep(c(1,2,4,5,3,6), each = 2) 
VarClass(lcaGND_C3cl4) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(orden, param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.75, "Myblue", 
                  ifelse(value < 0.75 & value >= 0.25, "Mygreen","Myred")))) %>% 
  reshape2::dcast(orden + param ~ Class) %>% arrange(orden) %>% select(-orden) %>% 
  kbl(caption = paste0("Probabilities to agree each item in the four-class global model for Students' endorsement of gender equality"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), 
      row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), 
                font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:5, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
```

```{r bestfit12}
counts4GND  %>% 
  mutate(class = factor(class, levels = orden4GND, labels = classes4GND),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% 
  arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes four-class global model for Students' endorsement of gender equality"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), 
      row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), 
                font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))

sizelca4_GND <- lcaGND$GND_lca_C3cl4.out$class_counts$modelEstimated %>% 
  dplyr::select(-count) %>% 
  rename_with(~ c("Gender", "Class")[which(c("proportion", "class") == .x)], 
              .cols = c("proportion", "class")) %>% 
  mutate(Class = factor(Class, levels = orden4GND, labels = classes4GND)) %>% 
  reshape2::melt(id.vars = c("Class"), variable.name = "Group") %>% 
  dplyr::arrange(Group) %>% 
  dplyr::group_by(Group) %>%
  dplyr::mutate(countT= sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(Class) %>%
  dplyr::mutate(per=value/countT) %>% 
  dplyr::select(Group, Class, per) 
```

```{r higprob1, eval=FALSE}
HighProb(lcaGND_C3cl4, sizelca4_GND,  orden = c(1,2,5,3,4,6), 
         title = "Response categories probabilities and class size for four-\nclasses global model for Students' endorsement of gender equality") 
cat("\n")
cat("\n")

```

```{r , eval=FALSE}

prob_gender <- read.table("data/MplusModels/LCA/GND_Prob_C3cl4.dat", header = FALSE)
colnames(prob_gender) <- c("GND1", "GND2", "GND3", "GND4", "GND5", "GND6", "GND_CPROB1", "GND_CPROB2", "GND_CPROB3", "GND_CPROB4", 
                           "GND_C", "ws", "IDSTUD", "id_s", "id_j")

prob_eth <- read.table("data/MplusModels/LCA/ETH_Prob_C3cl4.dat", header = FALSE)
colnames(prob_eth) <- c("ETH1", "ETH2", "ETH3", "ETH4", "ETH5", "ETH_CPROB1", "ETH_CPROB2", "ETH_CPROB3", "ETH_CPROB4", 
                        "ETH_C", "ws", "IDSTUD", "id_s", "id_j")

ISC_lca_prob <- full_join(prob_gender, prob_eth, by = c("id_s", "id_j", "IDSTUD")) %>% 
  select(-GND1, -GND2, -GND3, -GND4, -GND5, -GND6, -ETH1, -ETH2, -ETH3, -ETH4, -ETH5, -ws.y, -ws.x)

ISC_lca_prob2 <- ISC_lvRlca %>% left_join(ISC_lca_prob, 
                                          by = c("id_s", "id_j", "IDSTUD")) %>% 
  select(-IMM1, -IMM2, -IMM3, -IMM4, -IMM5)

classes4GND <- c("Fully egalitarian",
                "Competition- driven sexism",
                "Non-egalitarian",
                "Political egalitarian")
orden4GND <- c(2,4,3,1)
ISC_lca_prob2$GND_C <- factor(ISC_lca_prob2$GND_C, levels = orden4GND, labels = classes4GND)

classes4ETH <- c("Fully egalitarian",
                "Political non-egalitarian",
                "Non-egalitarian",
                "Employment non-egalitarian")
orden4ETH <- c(2,3,1,4)
ISC_lca_prob2$ETH_C <- factor(ISC_lca_prob2$ETH_C, levels = orden4ETH, labels = classes4ETH)


#head(ISC_lvRlca)
#head(ISC_lca_prob2)

#table(ISC_lca_prob$GND_C, ISC_lca_prob$ETH_C)
#table(ISC_lca_prob2$GND_C, ISC_lca_prob2$ETH_C)

#write_sav(ISC_lca_prob2, "data/MplusModels/LCA/Profiles_probabilities.sav")
#str(ISC_lca_prob2)


#sav <- read_sav("data/MplusModels/LCA/Profiles_probabilities.sav")
```
 

### Country comparability  

```{r eval=FALSE}
load("data/MplusModels_CntCov.RData")

Modelfit("CntCovGND", title = "Model fit statistics Gender equality model with Country as covariate") %>% 
  pack_rows("Country numerical covariate",1,3) %>% 
  row_spec(2, background = "lightgray") %>% 
  row_spec(c(2), bold = TRUE) %>% 
  footnote(general = "The best loglikelihood value was not replicated in 6-class model. The solution may not be trustworthy due to local maxima.") %>% 
  print()
cat('\n')
cat('\n')
```

```{r, eval=FALSE}
#----GND 5 groups
CntCovclasses5GND <- c("Fully egalitarian",
                "Strong competition- driven sexism",
                "Equal fully egalitarian",
                "Not involved",
                "Soft competition- driven sexism")
CntCovorden5GND <- c(2,3,1,4,5)
CntCovGND_C3cl5 <- CntCovGND$GND_CntCov_C3cl5.out$parameters$probability.scale %>% 
  rename_with(~ c("Class", "value")[which(c("LatentClass", "est") == .x)], .cols = c("LatentClass", "est")) %>% 
  mutate_at( c("param", "category", "Class"), ~ as.factor(.x)) %>% 
  mutate(Class = factor(Class, levels = CntCovorden5GND, labels = CntCovclasses5GND))

CntCovcounts5GND <- full_join(CntCovGND$GND_CntCov_C3cl5.out$class_counts$modelEstimated,
                        CntCovGND$GND_CntCov_C3cl5.out$class_counts$mostLikely,by = c("class"))
CntCovcounts5GND  %>% 
  mutate(class = factor(class, levels = CntCovorden5GND, labels = CntCovclasses5GND),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes 5-class Gender equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))

sizeCntCov5_GND <- CntCovGND$GND_CntCov_C3cl5.out$class_counts$modelEstimated %>% dplyr::select(-count) %>% 
  rename_with(~ c("Gender", "Class")[which(c("proportion", "class") == .x)], .cols = c("proportion", "class")) %>% 
  mutate(Class = factor(Class, levels = CntCovorden5GND, labels = CntCovclasses5GND)) %>% 
  reshape2::melt(id.vars = c("Class"), variable.name = "Group") %>% 
  dplyr::arrange(Group) %>% 
  dplyr::group_by(Group) %>%
  dplyr::mutate(countT= sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(Class) %>%
  dplyr::mutate(per=value/countT) %>% 
  dplyr::select(Group, Class, per) 

VarClass(CntCovGND_C3cl5) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.8, "Myblue", 
                                                 ifelse(value < 0.8 & value > 0.49, "Mygreen","Myred")))) %>% 
  reshape2::dcast(param ~ Class) %>% 
  kbl(caption = paste0("Probabilities to agree each item 5-class Gender equality model with country as covariate"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",5)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:6, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
# graphclass(CntCovGND_C3cl5, nclass = 5, orden = c(1,2,5,3,4,6), title = "LCA Gender equality with 5 classes with country as covariate", mg = FALSE)
# cat("\n")
# cat("\n")

HighProb(CntCovGND_C3cl5, sizeCntCov5_GND,  orden = c(1,2,5,3,4,6), title = "Response categories probabilities and class size for\n 5-classes Gender equality model with country as covariate") 
cat("\n")
cat("\n")

```

To evaluate the country comparability, the classes that were found in the independent model were identified to later check how many of them could be tested for comparability using a multigroup latent class model.  

The different classes identified according to each global model is summarized next, indicating in which countries the same class is present.  

- Global three-class model:  
  1. Fully egalitarian - ALL COUNTRIES
  2. Competition-driven sexism - ALL COUNTRIES
  3. Random response - BGR, LVA, LTU, MLT

- Global four-class model:  
  1. Fully egalitarian - ALL COUNTRIES
  2. Competition-driven sexism - ALL COUNTRIES
  3. Non-egalitarian - HRV, DNK, EST, FIN, ITA, NOR, SLV, SWE
  4. Political egalitarian - BFL, DNK, EST, NOR, SWE

- Global five-class model:  
  1. Fully egalitarian - ALL COUNTRIES
  2. Competition-driven sexism - ALL COUNTRIES
  3. Non-egalitarian - HRV, DNK, EST, FIN, ITA, NOR, SLV, SWE
  4. Political egalitarian - BFL, DNK, EST, NOR, SWE
  5. Reverse competition-driven sexism - BGR, MLT, NLD, NOR, SLV

- Global six-class model:  
  1. Fully egalitarian - ALL COUNTRIES
  2. Competition-driven sexism - ALL COUNTRIES
  3. Non-egalitarian - HRV, DNK, EST, FIN, ITA, NOR, SLV, SWE
  4. Political egalitarian - BFL, DNK, EST, NOR, SWE
  5. Reverse competition-driven sexism - BGR, MLT, NLD, NOR, SLV
  6. Pro-women pay/job - Not defined in individual country models  

With three classes, the random response class is not very interpretable. With six classes, a new no-identified class appears, which is not interpretable. With five classes, the reverse competition-driven sexism class is present in five countries but with class sizes lower than 1%, which would be not representative. With four classes, the four classes are identified across countries and two of them are present in all countries, which is the best model for comparability. 


#### Country multigroup analysis 

In table \@ref(tab:mgmodelfit1) different models with multigroup analysis are tested, first the more restricted model is evaluated, complete homogeneity. In this model, all conditional and unconditional probabilities are fixed to be equal across the groups.  

Then, the partial homogeneity is tested where only the conditional probabilities are constrained to be equal across the groups, and the class sizes are estimated freely. The second approach of partial homogeneity was also tested, but not included in the table, where only the conditional probabilities for the two common classes identified are constrained across groups, and the remaining are freely estimated along with the unconditional probabilities. Finally, the complete heterogeneous model is tested, where not only the unconditional probabilities are estimated freely but all the conditional probabilities as well.
In the last two models the best log-likelihood is not replicated, this means that the solution may not be trustworthy due to local maxima. These results cannot be considered valid.  

Just by looking at the valid results, the partial homogeneity where all conditional probabilities are constrained to be equal across groups shows a better fit compared to the more restricted model, the complete homogeneity. With it is valid to indicate that the 4 classes identified do not share the same unconditional probabilities (class sizes) across the groups, but the conditional probabilities can be considered as equal in all groups.  


<!-- \blandscape   -->
```{r mgmodelfit1}
ModelfitMGCntry("MGCntry_GND", title = "Multigroup model fit statistics global model with four-classes with Students' endorsement of gender equality", 
                filterval=c(4), fontn = 10) %>% 
  pack_rows("Four-class model", 1, 3) %>% 
  pack_rows("Complete homogeneity", 1, 1) %>% 
  pack_rows("Partial homogeneity", 2, 2) %>% 
  pack_rows("Complete heterogeneity", 3, 3) %>%    
  #pack_rows("Five-class model", 4, 6) %>% 
  #pack_rows("Complete homogeneity", 4, 4) %>% 
  #pack_rows("Partial homogeneity", 5, 5) %>% 
  #pack_rows("Complete heterogeneity",6, 6) %>%    
  row_spec(c(2), bold = TRUE) %>% 
  footnote(general = "The best log-likelihood value was not replicated for the following models:",
           number = c("4-class Complete heterogeneity model."
                      )) %>% 
  print() 
cat('\n') 
cat('\n') 
```

<!-- \elandscape   -->


Figure \@ref(fig:MGchom1) indicates the values for the patterns with the conditional probabilities fixed in all countries, but also the unconditional probabilities are constrained to be equal in all groups. Here can be observed that the patterns are like the ones identified in the independent models and the global model as well. But this would force the classes sizes to be 79.2% for the Fully egalitarian, 11.5% for the Competition-driven sexism, 7.7% for the Political egalitarian and 1.7% for the Non-egalitarian for all the countries which is not optimal.  

```{r MGchom1, fig.width=3, fig.height=2, fig.cap="Conditional probabilities to agree in a 4-class complete homogeneous multigroup model for Students' endorsement of gender equality scale", fig.pos='H'}
MGCntyCHo_GND_C3cl4 <- data.frame(full_join(MGCntry_GND$GND_MGCnty_C3cl4_3CHom.out$parameters$probability.scale %>% 
                                              mutate(class = as.numeric(substr(LatentClass,3,3))), 
                                            MGCntry_GND$GND_MGCnty_C3cl4_3CHom.out$class_counts$modelEstimated %>% 
                                              filter(variable == "C") %>% 
                                              select(class, proportion), by = c("class"))) %>% 
  mutate(ClassN = case_when(class == 1 ~ 4, class == 2 ~ 3, class == 3 ~ 2, class == 4 ~ 1)#,
         #ClassN = factor(ClassN, levels = 1:4, labels = classes4HetMG)
         )
cHomMGCntryGND <- MGCntyCHo_GND_C3cl4 %>% filter(category == 1) %>% 
  mutate(ClassN = factor(ClassN),
         param = factor(param, levels = unique(MGCntyCHo_GND_C3cl4$param)[c(1,2,5,3,4,6)]),
         proportion = ifelse(param == "GND6", proportion, NA)) %>% 
  ggplot() +
  geom_point(aes(x = param, y = est, group = ClassN, color = ClassN), size = 1) +
  geom_line(aes(param, est, group = ClassN, linetype = ClassN, color = ClassN)) +
  geom_text(aes(x = param, y = est, label = scales::percent(proportion, accuracy = 0.1), color = ClassN), size = 2, 
            nudge_x = -0.15, nudge_y = 0.1) +
  scale_fill_grey() + theme_bw() +
  #ggtitle("Classes in a country multigroup model with \ncomplete homogeneous classes") +
  theme(legend.position = "none", 
        title = element_text(size = 9),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 7),
        axis.text.x = element_text(angle = 90, size = 7, vjust = 0.5, hjust = 0)) +
  scale_y_continuous(breaks = c(0.25,0.5,0.75)) +
  labs(y="Response probabilities", linetype = "Latent Classes", color = "Latent Classes")  +
  scale_colour_manual(values=cbPalette[c(3,5,2,1)]) +
  scale_linetype_manual(values=c(3,5,2,1))
#scale_color_brewer(type = "qual", palette = "Dark2") 
cHomMGCntryGND
cat("\n")
cat("\n")
```


```{r}
MGCntyPHo_GND_C3cl4 <- data.frame(full_join(MGCntry_GND$GND_MGCnty_C3cl4_2PHom.out$parameters$probability.scale %>% mutate(C = as.numeric(substr(LatentClass,3,3))), 
                                            MGCntry_GND$GND_MGCnty_C3cl4_2PHom.out$class_counts$modelEstimated.patterns %>% 
                                              mutate(LatentClass = paste0(G,".",C)) %>% 
                                              group_by(G) %>% mutate(countT= sum(count, na.rm = TRUE)) %>% ungroup() %>% 
                                              mutate(proportion=round(count/countT*100, 1)) %>% 
                                              select(G, LatentClass, C, proportion), by = c("C"))) %>% 
  mutate(Group = G,
         ClassN = case_when(C == 1 ~ 3, C == 2 ~ 1, C == 3 ~ 2, C == 4 ~ 4),
         Group = as.character(factor(Group, levels = 1:length(cnt), labels = cnt))#,
         #ClassN = factor(ClassN, levels = 1:4, labels = classes4HetMG)
         )
```

The proportions based on a global model cannot be extrapolated to a country level, because only partial invariance was obtained those values are not applicable to each country. When considering the partial invariance results (indicated in Table \@ref(tab:classGND)) it is possible to identify that even though Fully egalitarian is the predominant class regarding size, in countries such as Bulgaria, Latvia, and Lithuania most predominant class is the Political egalitarian with around 45% of the population. Denmark, Finland, Norway and Sweden obtained the larger number of members in the Fully egalitarian class based on the model estimated counts. In contrast, Bulgaria and Lithuania has the larger number of members in the Non-egalitarian class compared to the rest of the European countries.   

```{r classGND}
MGCntyPHo_GND_C3cl4 %>% 
  select(Group, ClassN, proportion) %>% 
  unique() %>% 
  mutate(ClassN = factor(ClassN, levels = 1:4, labels = classes4GND)) %>% 
  reshape2::dcast(Group ~ ClassN) %>% 
  kbl(caption = paste0("Class sizes for partial homogeneity model for Students' endorsement of gender equality scale"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:5, width = "5em") %>%  
  print()
```


```{r MGphom1, fig.width=6, fig.height=6, fig.cap="Conditional probabilities to agree in a 4-class partial homogeneous multigroup model for Students' endorsement of gender equality scale", fig.pos='H', eval=FALSE}

pHomMGCntryGND <- MGCntyPHo_GND_C3cl4 %>% filter(category == 1) %>% plyr::rbind.fill(data.frame(param = "Size", Group = "Sweden")) %>% 
  mutate(ClassN = factor(ClassN),
         param = factor(param, levels = c(unique(MGCntyPHo_GND_C3cl4$param), "Size")[c(1,2,5,3,4,6,7)]), 
         proportion = ifelse(param == "GND6", proportion, NA)) %>% 
  ggplot() +
  geom_point(aes(x = param, y = est, group = ClassN, color = ClassN), size = 1) +
  geom_line(aes(param, est, group = ClassN, linetype = ClassN, color = ClassN)) +
  geom_text(aes(x = param, y = est, label = scales::percent(proportion/100, accuracy = 0.1), color = ClassN), size = 2, 
            nudge_x = 0.8, nudge_y = 0) +
  facet_wrap(Group ~ .) +
  scale_fill_grey() + theme_bw() +
  #ggtitle("Classes in a 4-class country multigroup partially homogeneous model, \nStudents' endorsement of gender equality") +
  theme(legend.position = "none", 
        title = element_text(size = 9),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 7),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.text.x = element_text(angle = 90, size = 7, vjust = 0.5, hjust = 0)) +
  scale_y_continuous(breaks = c(0.25,0.5,0.75), limits = c(0,1)) +
  geom_hline(yintercept = c(0.25,0.5,0.75), color = "gray", size = 0.2) +
  labs(y="Response probabilities", linetype = "Latent Classes", color = "Latent Classes")  +
  scale_colour_manual(values=cbPalette[c(1,2,3,5)]) +
  scale_linetype_manual(values=c(1,2,3,5))
#scale_color_brewer(type = "qual", palette = "Dark2") 
pHomMGCntryGND
cat("\n")
cat("\n")
```

```{r, fig.width=6, fig.height=5, fig.cap="Conditional probabilities to agree in a partial homogeneous multigroup model with constrains, \nGender equalyty scale", fig.pos='H', eval=FALSE}
MGCnty4PH_GND_C3cl4 <- data.frame(left_join(MGCntry_GND$GND_MGCnty_C3cl4_2.2PHom.out$parameters$probability.scale, 
                                            MGCntry_GND$GND_MGCnty_C3cl4_2.2PHom.out$class_counts$modelEstimated.patterns %>% 
                                              mutate(LatentClass = paste0(G,".",C)) %>% 
                                              group_by(G) %>% mutate(countT= sum(count, na.rm = TRUE)) %>% ungroup() %>% 
                                              mutate(proportion=round(count/countT*100, 1)) %>% 
                                              select(LatentClass, proportion), by = c("LatentClass"))) %>% 
  mutate(Group = sub('\\..*', '', LatentClass),
         Class = sub('*.{1,2}\\.', '', LatentClass),
         ClassN = case_when(LatentClass == 1.1 ~ 1, LatentClass == 1.2 ~ 2, LatentClass == 1.3 ~ 3, LatentClass == 1.4 ~ 5,
                            LatentClass == 2.1 ~ 1, LatentClass == 2.2 ~ 2, LatentClass == 2.3 ~ 3, LatentClass == 2.4 ~ 5,
                            LatentClass == 3.1 ~ 1, LatentClass == 3.2 ~ 2, LatentClass == 3.3 ~ 3, LatentClass == 3.4 ~ 5,
                            LatentClass == 4.1 ~ 1, LatentClass == 4.2 ~ 2, LatentClass == 4.3 ~ 3, LatentClass == 4.4 ~ 5,
                            LatentClass == 5.1 ~ 1, LatentClass == 5.2 ~ 2, LatentClass == 5.3 ~ 3, LatentClass == 5.4 ~ 5,
                            LatentClass == 6.1 ~ 1, LatentClass == 6.2 ~ 2, LatentClass == 6.3 ~ 3, LatentClass == 6.4 ~ 5,
                            LatentClass == 7.1 ~ 1, LatentClass == 7.2 ~ 2, LatentClass == 7.3 ~ 3, LatentClass == 7.4 ~ 5,
                            LatentClass == 8.1 ~ 1, LatentClass == 8.2 ~ 2, LatentClass == 8.3 ~ 3, LatentClass == 8.4 ~ 5,
                            LatentClass == 9.1 ~ 1, LatentClass == 9.2 ~ 2, LatentClass == 9.3 ~ 3, LatentClass == 9.4 ~ 5,
                            LatentClass == 10.1 ~ 1, LatentClass == 10.2 ~ 2, LatentClass == 10.3 ~ 4, LatentClass == 10.4 ~ 5,
                            LatentClass == 11.1 ~ 1, LatentClass == 11.2 ~ 2, LatentClass == 11.3 ~ 4, LatentClass == 11.4 ~ 5,
                            LatentClass == 12.1 ~ 1, LatentClass == 12.2 ~ 2, LatentClass == 12.3 ~ 4, LatentClass == 12.4 ~ 5,
                            LatentClass == 13.1 ~ 1, LatentClass == 13.2 ~ 2, LatentClass == 13.3 ~ 4, LatentClass == 13.4 ~ 5,
                            LatentClass == 14.1 ~ 1, LatentClass == 14.2 ~ 2, LatentClass == 14.3 ~ 4, LatentClass == 14.4 ~ 5),
         Group = as.character(factor(Group, levels = 1:length(cnt), labels = cnt))#,
         #ClassN = factor(ClassN, levels = 1:4, labels = classes4HetMG)
         )

MGCnty4PH_GND_C3cl4 %>% arrange(Group) %>% select(Group, param, est, ClassN, category) %>% VarClass(orden = c(1,2,5,3,4,6)) %>% 
  bind_rows(MGCnty4PH_GND_C3cl4 %>% select(Group, proportion, ClassN) %>% unique() %>% mutate(param = "Class Size (\\%)", category = "1") %>% 
              rename(est = proportion)) %>% 
  group_by(ClassN, Group, param) %>% filter(category == 1) %>%
  select(param, ClassN, Group, est) %>%
  mutate(value = cell_spec(est, color = ifelse(param == "Class Size (\\%)", "black", ifelse(est > 0.8, "Myblue", "Myred")))) %>%
  reshape2::dcast(Group + param ~  ClassN) %>% 
  mutate(param = factor(param, levels = unique(param)[c(1,2,3,6,4,5,7)])) %>% arrange(Group,param) %>% 
  kbl(caption = "Probabilities to Agree in a 5-class country multigroup model with 2 constrained classes,\nStudents' endorsement of gender equality", 
      booktabs = TRUE, longtable = TRUE,
      align = "llrrrrr", row.names = FALSE, digits = 3, escape = FALSE#, col.names = c("Item", classes2EU, classes2LA)
      ) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 10) %>%
  column_spec(1, width = "5em") %>%
  column_spec(2, width = "16em") %>%
  column_spec(3:7, width = "3em") %>%
  collapse_rows(1, valign = "top") %>%
  print()
cat('\n')
cat('\n')

MGCnty4PH_GND_C3cl4 %>% filter(category == 1) %>% plyr::rbind.fill(data.frame(param = "Size", Country = "Sweden")) %>% 
  mutate(ClassN = factor(ClassN),
         param = factor(param, levels = c(unique(MGCnty4PH_GND_C3cl4$param), "Size")[c(1,2,5,3,4,6,7)]),
         proportion = ifelse(param == "GND6", proportion, NA)) %>% 
  ggplot() +
  geom_point(aes(x = param, y = est, group = ClassN, color = ClassN), size = 1) +
  geom_line(aes(param, est, group = ClassN, linetype = ClassN, color = ClassN)) +
  geom_text(aes(x = param, y = est, label = scales::percent(proportion/100, accuracy = 0.1), color = ClassN), size = 2, 
            nudge_x = 0.8, nudge_y = 0.0) +
  facet_wrap(Group ~ .) +
  scale_fill_grey() + theme_bw() +
  ggtitle("Classes in a 5-class multigroup model with 2 constrained classes,\nStudents' endorsement of gender equality") +
  theme(legend.position = "none", 
        title = element_text(size = 9),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 7),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.text.x = element_text(angle = 90, size = 7, vjust = 0.5, hjust = 0)) +
  scale_y_continuous(breaks = c(0.25,0.5,0.75), limits = c(0,1)) +
  geom_hline(yintercept = c(0.25,0.5,0.75), color = "gray", size = 0.2) +
  labs(y="Response probabilities", linetype = "Latent Classes", color = "Latent Classes")  +
  scale_colour_manual(values=cbPalette[c(5,2,1,3,4)]) +
  scale_linetype_manual(values=c(5,2,1,3,4))
cat("\n")
cat("\n")
```


### Confirmatory Latent Class Analysis  

The confirmatory model was performed by establishing some constraints based on the results obtained in the previous analysis. For the Students' endorsement of gender equality, the hypothesis tested is that students that would be classified into the Fully egalitarian class would highly agree (totally) with that *Men and women should have equal opportunities to take part in government* which is the first item in the scale. Also, they would have the same probability to agree (in the same level) to both items *Men and Women should have the same rights in every way* and *Women should stay out of politics (r)* that are the next items in the questionnaire.    

The confirmatory approach also test that the students that would be classified in the second group Competition-driven sexism class would be equally likely to agree to the first item in the questionnaire *Men and women should have equal opportunities to take part in government* as the members of the first class agree with the second item and third item. And that these students would not have a clear attitude (either agree or disagree) with the item *Women should stay out of politics (r)* which means that they will tend to give a random response to this item.  

This can be tested by setting the conditional probabilities for the first item in the Fully egalitarian class to 1. The probabilities to agree to the second and third item in the first class and the first item in the Competition-driven sexism class to be equal. The third hypothesis can be tested by setting the conditional probability of the third item in the Competition-driven sexism class to 0.5. The rest of the conditional probabilities are estimated freely as can be seen in table \@ref(tab:probconf1). Values for the thresholds, class sizes and the probabilities to agree to each item can be found in the appendix table \@ref(tab:detailed2).  

In table \@ref(tab:confm1) the model fit statistics of this model indicate that they do not differ considerably from the exploratory approach where the log-likelihood, BIC, and aBIC have higher values, only the AIC value is better in the exploratory model.  


```{r cfa1, results='hide'}
orden4GND <- c(1,2,3,4)
ConfUnstlcaGND <- ConflcaGND$GND_Conflca_C3cl4.out$parameters$unstandardized %>% 
  rename_with(~ c("Parameter", "Classo", "value")[which(c("param", "LatentClass", "est") == .x)], .cols = c("param", "LatentClass", "est")) %>% 
  mutate_at( c("Parameter", "paramHeader", "Classo"), ~ as.factor(.x)) %>%  
  mutate(Class = case_when(Parameter == "C#1" ~ 1,
                        Parameter == "C#2" ~ 2,
                        Parameter == "C#3" ~ 3,
                        TRUE ~ as.numeric(Classo)),
         Parameter = case_when(Parameter %in% c("C#1", "C#2", "C#3") ~ "Means",
                        TRUE ~ as.character(Parameter))) %>% 
  mutate(Class = factor(Class, levels = orden4GND, labels = classes4GND)) %>% 
  group_by(Class, Parameter) %>% 
  select(Parameter, Class, value) 

ConfUnstlcaGND$orden <- c(rep(c(1,2,4,5,3,6),4), rep(7,3))
ConfUnstlcaGND %>% 
  reshape2::dcast(orden + Parameter ~ Class) %>% arrange(orden) %>% select(-orden) %>% 
  kbl(caption = paste0("Thresholds 4-class Confirmatory LCA Students' endorsement of gender equality"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 3, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "4em") %>%  
  column_spec(2:5, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()


SummConflcaGND <- ConflcaGND$GND_Conflca_C3cl4.out$parameters$probability.scale %>% 
  rename_with(~ c("Class", "value")[which(c("LatentClass", "est") == .x)], .cols = c("LatentClass", "est")) %>% 
  mutate_at( c("param", "category", "Class"), ~ as.factor(.x)) %>% 
  mutate(Class = factor(Class, levels = orden4GND, labels = classes4GND))

ConfcountsGND <- full_join(ConflcaGND$GND_Conflca_C3cl4.out$class_counts$modelEstimated,
                        ConflcaGND$GND_Conflca_C3cl4.out$class_counts$mostLikely,by = c("class"))
ConfcountsGND  %>% 
  mutate(class = factor(class, levels = orden4GND, labels = classes4GND),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes 4-class Confirmatory LCA Students' endorsement of gender equality"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))

sizeConflca_GND <- ConflcaGND$GND_Conflca_C3cl4.out$class_counts$modelEstimated %>% dplyr::select(-count) %>% 
  rename_with(~ c("Gender", "Class")[which(c("proportion", "class") == .x)], .cols = c("proportion", "class")) %>% 
  mutate(Class = factor(Class, levels = orden4GND, labels = classes4GND)) %>% 
  reshape2::melt(id.vars = c("Class"), variable.name = "Group") %>% 
  dplyr::arrange(Group) %>% 
  dplyr::group_by(Group) %>%
  dplyr::mutate(countT= sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(Class) %>%
  dplyr::mutate(per=value/countT) %>% 
  dplyr::select(Group, Class, per) 

```

```{r probconf1}

SummConflcaGND$orden = rep(c(1,2,4,5,3,6), each = 2)
VarClass(SummConflcaGND) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(orden, param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.8, "Myblue", 
                                                 ifelse(value < 0.8 & value > 0.49, "Mygreen","Myred")))) %>% 
  reshape2::dcast(orden + param ~ Class) %>% arrange(orden) %>% select(-orden) %>% 
  kbl(caption = paste0("Probabilities to agree each item in a 4-class Confirmatory LCA for Students' endorsement of gender equality scale"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:5, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
# graphclass(ConflcaGND, nclass = 4, orden = c(1,2,5,3,4,6), title = "Confirmatory LCA Gender equality with 4 classes", mg = FALSE)
# cat("\n")
# cat("\n")
```

This result suggests that it can be plausible to test some other similar hypotheses in order to identify straightforward patterns based on the theoretical background of the items analysed that could help to interpret them even more.  


```{r confm1}

FinalCompGender <- list("GND_lca_C3cl4.out" = lcaGND$GND_lca_C3cl4.out, "GND_Conflca_C3cl4.out"=ConflcaGND$GND_Conflca_C3cl4.out) 
ModelfitConf("FinalCompGender", title = "Model fit statistics 4-class Confirmatory LCA for Students' endorsement of gender equality scale") %>% 
  pack_rows("All countries",1,2) %>% 
  print()
cat('\n')
cat('\n')
```


```{r CFAm1, eval=FALSE}
HighProb(SummConflcaGND, sizeConflca_GND,  orden = c(1,2,5,3,4,6)#, title = "Conditional probabilities and class size for\n 4-classes Confirmatory LCA model Gender equality scale"
         ) 
cat("\n")
cat("\n")
```

\newpage  

## Students' endorsement of equal rights for all ethnic/racial groups scale  

This scale is composed of 5 items, in the following results, these items were ordered in the output for an easier interpretation of the results. This ordering considers first *All ethnic and racial groups should have equal chance to get a good education* (ETH1), *All ethnic and racial groups should have an equal chance to get good jobs* (ETH2), *Members of all ethnic and racial groups should have same rights and responsibilities* (ETH5), and *Schools should teach students to respect members of all ethnic and racial groups* (ETH3), followed by *Members of all ethnic and racial groups should be encouraged to run in elections* (ETH4). As mentioned before all these variables were recorded in two categories, as Agree and Disagree. All 14 countries were analyzed independently and then pooled in the same dataset.    

### Analysis by country  

A latent class analysis with 1 to 6-class models was performed in each country to evaluate the model fit of each one of them\footnote{Model fit statistics for each model can be found in the Appendix @ref(tab:detailed3)}. The results are summarized in table \@ref(tab:summodelfitcntry2). In most European countries, the best model fit based on the different criteria indicated previously are by including 3 or 4 latent classes. 

For Belgium, Bulgaria, Estonia, Italy, Lithuania, Latvia, Slovenia, The Netherlands, Norway, Slovenia and Sweden, according to the statistical tests, BIC, and aBIC criteria, the best model is a 3-class model. 

On the other hand, Denmark and Malta have a better model fit in a 4-class model, consistently between statistical tests and BIC criteria.  

In Croatia models, tests indicate that a 2-class model is better for their data, even though BIC indicates a 3-class model to have the lowest value.

Norway is the only country from the sample that the best model fit is the one with 5 latent classes according to the statistical tests and BIC and aBIC.

It is a common tendency in all the evaluated countries the AIC value is lower in the models with one more class than indicated by the statistical tests and BIC and aBIC. This is consistent with the indication that this criterion tends to overfit the data. Values of Entropy are higher when the tests are significant, but consistent with a better fit of the data the lowest entropy found in the 3-class model is in Belgium (60.5%) and the highest value in Sweden (90.2%). The log-likelihood reduction is consistent in all countries, where having more than 3 latent classes reduces the log-likelihood around 0.1% and 0.6%.

The bivariate residuals were also analyzed, and all countries have residuals around the range of acceptable [-2 ; 2] as shown in the figure \@ref(fig:resid2cnt). There are just two values outside the ranges in Italy with a 3-class model. 


```{r modelfitcnty2, echo=FALSE, results='hide'}

ModelfitByContry("ByCountry_ETH", title = "Model fit statistics LCA by country Students' endorsement of equal rights for all ethnic/racial groups scale") %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 10) %>% 
  column_spec(c(2,3,11), width = "3em") %>% 
  column_spec(c(4:10,12), width = "4em") %>% 
  column_spec(c(1), width = "6em") %>%
  pack_rows("Belgium (Flemish)",1,6) %>%
#  row_spec(2:3, background = "lightgray") %>% 
  row_spec(c(3), bold = TRUE) %>%
  pack_rows("Bulgaria",7,12) %>%
#  row_spec(9:10, background = "lightgray") %>% 
  row_spec(c(10), bold = TRUE) %>%
  pack_rows("Croatia",13,18) %>%
#  row_spec(14:15, background = "lightgray") %>% 
  row_spec(c(16), bold = TRUE) %>%
  pack_rows("Denmark",19,24) %>%
#  row_spec(21:22, background = "lightgray") %>% 
  row_spec(c(22), bold = TRUE) %>%
  pack_rows("Estonia",25,30) %>%
#  row_spec(27:28, background = "lightgray") %>% 
  row_spec(c(27), bold = TRUE) %>%
  pack_rows("Finland",31,36) %>%
#  row_spec(32:33, background = "lightgray") %>% 
  row_spec(c(33), bold = TRUE) %>%
  pack_rows("Italy",37,42) %>%
#  row_spec(39:40, background = "lightgray") %>% 
  row_spec(c(39), bold = TRUE) %>%
  pack_rows("Latvia",43,48) %>%
#  row_spec(45:46, background = "lightgray") %>% 
  row_spec(c(45), bold = TRUE) %>%
  pack_rows("Lithuania",49,54) %>%
#  row_spec(51:52, background = "lightgray") %>% 
  row_spec(c(51), bold = TRUE) %>%
  pack_rows("Malta",55,60) %>%
 # row_spec(57:58, background = "lightgray") %>% 
  row_spec(c(58), bold = TRUE) %>%
  pack_rows("The Netherlands",61,66) %>%
#  row_spec(62:63, background = "lightgray") %>% 
  row_spec(c(63), bold = TRUE) %>%
  pack_rows("Norway",67,72) %>%
#  row_spec(70:71, background = "lightgray") %>% 
  row_spec(c(71), bold = TRUE) %>%
  pack_rows("Slovenia",73,78) %>%
#  row_spec(75:76, background = "lightgray") %>% 
  row_spec(c(76), bold = TRUE) %>%
  pack_rows("Sweden",79,84) %>%
#  row_spec(81:82, background = "lightgray") %>% 
  row_spec(c(82), bold = TRUE) %>%
  collapse_rows(1, valign = "top") %>% 
  footnote(general = "The best loglikelihood value was not replicated for the following models",
           number = c("Denmark - 6-classes complete heterogeneity; ")) %>% 
  print()
cat('\n')
cat('\n')
```

\blandscape
```{r summodelfitcntry2}
ModelfitByContry("ByCountry_ETH", 
                 title = "Best model, fit statistics individual country model Students' endorsement of \nequal rights for all ethnic/racial groups", 
                 filterval = TRUE, fontn = 10) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 10) %>% 
  column_spec(c(2,3,9:10,12), width = "3em") %>% 
  column_spec(c(4:8,11,13), width = "4em") %>% 
  column_spec(c(1), width = "8em") %>% 
  footnote(general = "Best model based on the lowest value of BIC") %>% 
  print()
cat('\n')
cat('\n')
```
\elandscape


```{r resid2cnt, fig.width=6, fig.height=5, fig.cap="Bivariate standardized residuals for individual country models for Students' endorsement of equal rights for all ethnic/racial groups scale", fig.pos='H'}

residuals_ByCntryETH <-  data.frame(`BFL 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_BFL_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `BGR 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_BGR_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `DNK 4cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_DNK_C3cl4.out$tech10$bivar_model_fit_info$z,
                            `EST 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_EST_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `FIN 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_FIN_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `HRV 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_HRV_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `ITA 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_ITA_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `LTU 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_LTU_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `LVA 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_LVA_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `MLT 4cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_MLT_C3cl4.out$tech10$bivar_model_fit_info$z,
                            `NLD 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_NLD_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `NOR 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_NOR_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `SVN 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_SVN_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `SWE 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_SWE_C3cl3.out$tech10$bivar_model_fit_info$z)
residuals_ByCntryETH <- residuals_ByCntryETH %>% mutate(x = 1:nrow(residuals_ByCntryETH)) %>%
  reshape2::melt(id.vars = c("x"))

pByCntryETH <- residuals_ByCntryETH %>% ggplot() +
  geom_point(aes(x = x, y = value), size = 1) +
  geom_hline(yintercept=c(-1.96,1.96), linetype = "dashed", size = 0.9, color = "black") +
  scale_fill_grey() + theme_bw() +
  facet_wrap(variable~., scales = "free_y") +
  #ggtitle("Bivariate residuals by independent best country model fit") +
  labs(y="Standardized residuals")  +
  theme(legend.position = "none", legend.box="vertical",
        strip.text.y = element_text(size = 8),
        legend.spacing.y = unit(-0.2, 'cm'),
        title = element_text(size = 9),
        axis.title = element_text(size = 10),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 9),
        axis.text.y = element_text(size = 9)) 
pByCntryETH
cat("\n")
cat("\n")

```

In figure \@ref(fig:classes2cnt) the classes of each independent model can be identified. In the figure, the conditional probabilities for agreement to each item are shown and plotted for each of the classes modeled in each country. Here can be identified two classes that are similar in all the models, the green and purple lines.

- **Fully egalitarian:** Most likely to agree to all items in the scale  (green line). This class can be observed in all countries. Conditional probabilities greater than 0.7 to agree, class sizes around 61.8% (Latvia) and 90% (Sweden)

- **Political non-egalitarian:** Most likely to agree to all items but a random answer to *All ethnic and racial groups should be encouraged to run in elections (ETH4)* item (orange line). This class can be observed in all countries. Conditional probabilities to agree higher than 0.5 in all items but to the political item (< 0.5), class sizes around 7.6% (Denmark) and 36% (Latvia).

- **Non-egalitarian:** Not likely to agree to any item in the scale (purple line). This class can be observed in all countries. Conditional probabilities lower than 0.5 to agree all items, class sizes around 1.4% (Lithuania) and 5.4% (Bulgaria).

- **Country specific class:** (pink line). This class can be observed in two countries.  
  - Employment non-egalitarian: Not likely to agree to *All ethnic and racial groups should have an equal chance to get good jobs (ETH2)* item. Class size 8.3% (Malta)
  - Strong political non-egalitarian: Most likely to agree to most items in the scale but not likely to agree to *All ethnic and racial groups should be encouraged to run in elections (ETH4)* item. Class size 21.3% (Denmark).

```{r classes2cnt, fig.width=6, fig.height=5, fig.cap="Classes for best individual country model for Students' endorsement of equal rights for all ethnic/racial groups", fig.pos='H'}

classes_ByCntryETH <-  rbind(data.frame(Country = "Belgium (Flemish)", 
                                      left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_BFL_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_BFL_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param"))) %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Bulgaria",
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_BGR_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_BGR_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param"))) %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 2,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Denmark", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_DNK_C3cl4.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_DNK_C3cl4.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 4,
                                                              LatentClass == 3 ~ 1,
                                                              LatentClass == 4 ~ 2)),
                            data.frame(Country = "Estonia", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_EST_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_EST_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 1,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 2)),
                            data.frame(Country = "Finland", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_FIN_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_FIN_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 2,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Croatia", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_HRV_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_HRV_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 2)),
                            data.frame(Country = "Italy", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_ITA_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_ITA_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Lithuania", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_LTU_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_LTU_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Latvia", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_LVA_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_LVA_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Malta", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_MLT_C3cl4.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_MLT_C3cl4.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 4,
                                                              LatentClass == 3 ~ 2,
                                                              LatentClass == 4 ~ 1)),
                            data.frame(Country = "Netherlands", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_NLD_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_NLD_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 2)),
                            data.frame(Country = "Norway", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_NOR_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_NOR_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 1,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 2)),
                            data.frame(Country = "Slovenia", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_SVN_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_SVN_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Sweden", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_SWE_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_SWE_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 2)))

pByCntryETH <- classes_ByCntryETH %>% filter(category == 1) %>% plyr::rbind.fill(data.frame(param = "Size", Country = "Sweden")) %>% 
  mutate(LatentClass = factor(LatentClass),
         param = factor(param, levels = c(unique(classes_ByCntryETH$param), "Size")[c(1,2,5,3,4,6)])) %>% 
  ggplot() +
  geom_point(aes(x = param, y = est, group = LatentClass, color = LatentClass), size = 1) +
  geom_line(aes(param, est, group = LatentClass, linetype = LatentClass, color = LatentClass)) +
  geom_text(aes(x = param, y = est, label = scales::percent(proportion, accuracy = 0.1), color = LatentClass), size = 2, 
            nudge_x = 0.8, nudge_y = 0) +
  facet_wrap(Country ~ .) +
  scale_fill_grey() + theme_bw() +
  #ggtitle("Classes by country for Students' endorsement of equal rights for all ethnic/racial groups scale") +
  theme(legend.position = "none", 
        title = element_text(size = 9),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 7),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.text.x = element_text(angle = 90, size = 7, vjust = 0.5, hjust = 0)) +
  scale_y_continuous(breaks = c(0.25,0.5,0.75), limits = c(0,1)) +
  labs(y="Response probabilities", linetype = "Latent Classes", color = "Latent Classes")  +
  scale_color_brewer(type = "qual", palette = "Dark2") 
pByCntryETH
cat("\n")
cat("\n")

```

```{r resid2cnt2, fig.width=6, fig.height=5, fig.cap="Bivariate model fit standardized residuals", fig.pos='H', eval=FALSE}

residuals_ByCntryETH <-  data.frame(`BFL 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_BFL_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `BGR 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_BGR_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `DNK 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_DNK_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `EST 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_EST_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `FIN 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_FIN_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `HRV 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_HRV_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `ITA 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_ITA_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `LTU 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_LTU_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `LVA 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_LVA_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `MLT 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_MLT_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `NLD 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_NLD_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `NOR 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_NOR_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `SVN 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_SVN_C3cl3.out$tech10$bivar_model_fit_info$z,
                            `SWE 3cl` = ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_SWE_C3cl3.out$tech10$bivar_model_fit_info$z)
residuals_ByCntryETH <- residuals_ByCntryETH %>% mutate(x = 1:nrow(residuals_ByCntryETH)) %>%
  reshape2::melt(id.vars = c("x"))

pByCntryETH <- residuals_ByCntryETH %>% ggplot() +
  geom_point(aes(x = x, y = value), size = 1) +
  geom_hline(yintercept=c(-1.96,1.96), linetype = "dashed", size = 0.9, color = "black") +
  scale_fill_grey() + theme_bw() +
  facet_wrap(variable~., scales = "free_y") +
  ggtitle("Residuals by country models") +
  labs(y="Standardized residuals")  +
  theme(legend.position = "none", legend.box="vertical",
        strip.text.y = element_text(size = 8),
        legend.spacing.y = unit(-0.2, 'cm'),
        title = element_text(size = 9),
        axis.title = element_text(size = 8),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 8),
        axis.text.y = element_text(size = 8)) 
pByCntryETH
cat("\n")
cat("\n")

```

```{r classes2cnt2, fig.width=6, fig.height=5, fig.cap="By country 3-Classes model", fig.pos='H', eval=FALSE}

classes_ByCntryETH <-  rbind(data.frame(Country = "Belgium (Flemish)", 
                                      left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_BFL_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_BFL_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param"))) %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Bulgaria",
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_BGR_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_BGR_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param"))) %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 2,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Denmark", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_DNK_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_DNK_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 3)),
                            data.frame(Country = "Estonia", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_EST_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_EST_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 1,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 2)),
                            data.frame(Country = "Finland", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_FIN_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_FIN_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 2,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Croatia", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_HRV_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_HRV_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 2)),
                            data.frame(Country = "Italy", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_ITA_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_ITA_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Lithuania", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_LTU_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_LTU_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Latvia", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_LVA_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_LVA_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Malta", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_MLT_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_MLT_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 2,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Netherlands", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_NLD_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_NLD_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 2)),
                            data.frame(Country = "Norway", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_NOR_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_NOR_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 1,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 2)),
                            data.frame(Country = "Slovenia", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_SVN_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_SVN_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 2,
                                                              LatentClass == 2 ~ 3,
                                                              LatentClass == 3 ~ 1)),
                            data.frame(Country = "Sweden", 
                                       left_join(ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_SWE_C3cl3.out$parameters$probability.scale,
                                                ByCountry_ETH$data.MplusModels.ByCountry.ETHlca_SWE_C3cl3.out$class_counts$modelEstimated %>%
                                                  mutate(LatentClass = as.character(class), param = "ETH4") %>% 
                                                  select(LatentClass, param, proportion), by = c("LatentClass", "param")))  %>% 
                               mutate(LatentClass = case_when(LatentClass == 1 ~ 3,
                                                              LatentClass == 2 ~ 1,
                                                              LatentClass == 3 ~ 2)))

pByCntryETH <- classes_ByCntryETH %>% filter(category == 1) %>% 
  mutate(LatentClass = factor(LatentClass),
         param = factor(param, levels = unique(classes_ByCntryETH$param)[c(1,2,3,5,4)])) %>% 
  ggplot() +
  geom_point(aes(x = param, y = est, group = LatentClass, color = LatentClass), size = 1) +
  geom_line(aes(param, est, group = LatentClass, linetype = LatentClass, color = LatentClass)) +
  geom_text(aes(x = param, y = est, label = scales::percent(proportion, accuracy = 0.1), color = LatentClass), size = 2, 
            nudge_x = -0.15, nudge_y = 0.1) +
  facet_wrap(Country ~ .) +
  scale_fill_grey() + theme_bw() +
  ggtitle("Classes by country for Students' endorsement of equal rights for all ethnic/racial groups scale") +
  theme(legend.position = "none", 
        title = element_text(size = 9),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 7),
        axis.text.x = element_text(angle = 90, size = 7, vjust = 0.5, hjust = 0)) +
  scale_y_continuous(breaks = c(0.25,0.5,0.75)) +
  labs(y="Response probabilities", linetype = "Latent Classes", color = "Latent Classes")  +
  scale_color_brewer(type = "qual", palette = "Dark2") 
pByCntryETH
cat("\n")
cat("\n")

```


### General model

The model with a single class has the largest AIC (159379), BIC (159423), and aBIC (159407) values for the European countries model in Table \@ref(tab:modelfitlca2), indicating that this model fits data worse than all other models. In addition, the P-values of the VLMR test, and LMR in the 2-class model are all < 0.0001; this means that both tests reject the single-class model in favor of a model with at least two latent classes. In other words, there exists heterogeneity in the target population regarding attitudes towards gender equality.  

In the 6-class model, the LMR LR and VLMR are not statistically significant (P > 0.05). That is, the two tests are in favor of at most 5 classes. In contrast, AIC, BIC and aBIC values are all smaller in the 5-class model than those in the 6-class model; thus, consider that the models with more than 5 classes are not preferred. The relative entropy given by Mplus software decreases when including more than 4 classes and increases again with the 6-class model; this would suggest that a model with at least 6 class or 4 class is preferred. Together with the percentage of reduction in the log-likelihood value, that indicates that by adding two classes to the model the log-likelihood is reduced by 12.1%, this reduction is only increased by 1.5% if the model is a 3-class model and finally this value is reduced close to 0 if more than 5 classes are included.  

Now, the preferred model must be either the 4-class or higher model considering the residuals of each model in figure \@ref(fig:resid2), where all values are around -1.96 and 1.96. Theoretically, we tend to determine that the 4-class LCA model is the preferred model. We will show later that the classes identified by the 4-class model are more interpretable and representative than the rest of the models. And in particular that 3-classes can be compared across countries.  


\blandscape  
```{r modelfitlca2, echo=FALSE}
Modelfit("lcaETH", title = "Model fit statistics LCA Students' endorsement of equal rights for all ethnic/racial groups scale", fontn = 11) %>%
  pack_rows("All countries",1,7) %>%
  #row_spec(3:4, background = "lightgray") %>% 
  row_spec(c(3,4), bold = TRUE) %>%
  print()
cat('\n')
cat('\n')
```
\elandscape

```{r resid2, fig.width=5, fig.height=3, fig.cap="Bivariate model fit standardized residuals global model for Students' endorsement of equal rights for all ethnic/racial groups scale", fig.pos='H'}

residuals_ETH <-  data.frame(cl1 = lcaETH$ETH_lca_C3cl1.out$tech10$bivar_model_fit_info$z,
                            cl2 = lcaETH$ETH_lca_C3cl2.out$tech10$bivar_model_fit_info$z,
                            cl3 = lcaETH$ETH_lca_C3cl3.out$tech10$bivar_model_fit_info$z,
                            cl4 = lcaETH$ETH_lca_C3cl4.out$tech10$bivar_model_fit_info$z,
                            cl5 = lcaETH$ETH_lca_C3cl5.out$tech10$bivar_model_fit_info$z,
                            cl6 = lcaETH$ETH_lca_C3cl6.out$tech10$bivar_model_fit_info$z)
residuals_ETH <- residuals_ETH %>% mutate(x = 1:nrow(residuals_ETH)) %>%
  reshape2::melt(id.vars = c("x"))

pETH <- residuals_ETH %>% ggplot() +
  geom_point(aes(x = x, y = value)) +
  geom_hline(yintercept=c(-1.96,1.96), linetype = "dashed", size = 0.9, color = "black") +
  scale_fill_grey() + theme_bw() +
  facet_wrap(variable~., scales = "free_y") +
  #ggtitle("Standardized bivariate residuals for Attitudes towards ethnic and race equal rights models") +
  labs(x = "Parameters", y="Standardized residuals", color = "Number of classes")  +
  theme(legend.position = "none", legend.box="vertical",
        strip.text.y = element_text(size = 8),
        legend.spacing.y = unit(-0.2, 'cm'),
        title = element_text(size = 9),
        axis.title = element_text(size = 8),
        axis.text.x = element_blank(),
        axis.text.y = element_text(size = 8)) +
  scale_color_brewer(type = "qual", palette = "Dark2")
pETH
cat("\n")
cat("\n")

```

```{r, echo=FALSE, fig.width=6, fig.height=5, results='hide'}
#----ETH 3 groups----
classes3ETH <- c("Fully egalitarian",
                "Political non-egalitarian",
                "Not egalitarian")
orden3ETH <- c(2,3,1)
lcaETH_C3cl3 <- lcaETH$ETH_lca_C3cl3.out$parameters$probability.scale %>% 
  rename_with(~ c("Class", "value")[which(c("LatentClass", "est") == .x)], .cols = c("LatentClass", "est")) %>% 
  mutate_at( c("param", "category", "Class"), ~ as.factor(.x)) %>% 
  mutate(Class = factor(Class, levels = orden3ETH, labels = classes3ETH))

counts3ETH <- full_join(lcaETH$ETH_lca_C3cl3.out$class_counts$modelEstimated,
                        lcaETH$ETH_lca_C3cl3.out$class_counts$mostLikely,by = c("class"))
counts3ETH  %>% 
  mutate(class = factor(class, levels = orden3ETH, labels = classes3ETH),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes 3-class Ethnic and race equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",3)), row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))


sizelca3_ETH <- lcaETH$ETH_lca_C3cl3.out$class_counts$modelEstimated %>% dplyr::select(-count) %>% 
  rename_with(~ c("Ethnic", "Class")[which(c("proportion", "class") == .x)], .cols = c("proportion", "class")) %>% 
  mutate(Class = factor(Class, levels = orden3ETH, labels = classes3ETH)) %>% 
  reshape2::melt(id.vars = c("Class"), variable.name = "Group") %>% 
  dplyr::arrange(Group) %>% 
  dplyr::group_by(Group) %>%
  dplyr::mutate(countT= sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(Class) %>%
  dplyr::mutate(per=value/countT) %>% 
  dplyr::select(Group, Class, per) 

VarClass(lcaETH_C3cl3) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.8, "Myblue", 
                                                 ifelse(value < 0.8 & value > 0.49, "Mygreen","Myred")))) %>% 
  reshape2::dcast(param ~ Class) %>% 
  kbl(caption = paste0("Probabilities to agree each item 3-class Ethnic and race equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",3)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:4, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
# graphclass(lcaETH_C3cl3, nclass = 3, orden = c(1,2,3,5,4), 
#            title = "LCA Ethnic and race equality with 3 classes", mg = FALSE)
# cat("\n")
# cat("\n")

#HighProb(lcaETH_C3cl3, sizelca3_ETH, orden = c(1,2,3,5,4),  
#         title = "Response categories probabilities and class size for\n 3-classes Ethnic and race equality #model")
cat("\n")
cat("\n")

#----ETH 4 groups----
classes4ETH <- c("Fully egalitarian",
                "Political non-egalitarian",
                "Non-egalitarian",
                "Employment non-egalitarian")
orden4ETH <- c(2,3,1,4)
lcaETH_C3cl4 <- lcaETH$ETH_lca_C3cl4.out$parameters$probability.scale %>% 
  rename_with(~ c("Class", "value")[which(c("LatentClass", "est") == .x)], .cols = c("LatentClass", "est")) %>% 
  mutate_at( c("param", "category", "Class"), ~ as.factor(.x)) %>% 
  mutate(Class = factor(Class, levels = orden4ETH, labels = classes4ETH))

counts4ETH <- full_join(lcaETH$ETH_lca_C3cl4.out$class_counts$modelEstimated,
                        lcaETH$ETH_lca_C3cl4.out$class_counts$mostLikely,by = c("class"))
counts4ETH  %>% 
  mutate(class = factor(class, levels = orden4ETH, labels = classes4ETH),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes 4-class Ethnic and race equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))


sizelca4_ETH <- lcaETH$ETH_lca_C3cl4.out$class_counts$modelEstimated %>% dplyr::select(-count) %>% 
  rename_with(~ c("Ethnic", "Class")[which(c("proportion", "class") == .x)], .cols = c("proportion", "class")) %>% 
  mutate(Class = factor(Class, levels = orden4ETH, labels = classes4ETH)) %>% 
  reshape2::melt(id.vars = c("Class"), variable.name = "Group") %>% 
  dplyr::arrange(Group) %>% 
  dplyr::group_by(Group) %>%
  dplyr::mutate(countT= sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(Class) %>%
  dplyr::mutate(per=value/countT) %>% 
  dplyr::select(Group, Class, per) 

VarClass(lcaETH_C3cl4) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.8, "Myblue", 
                                                 ifelse(value < 0.8 & value > 0.49, "Mygreen","Myred")))) %>% 
  reshape2::dcast(param ~ Class) %>% 
  kbl(caption = paste0("Probabilities to agree each item 4-class Ethnic and race equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:5, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
# graphclass(lcaETH_C3cl4, nclass = 4, orden = c(1,2,3,5,4), 
#            title = "LCA Ethnic and race equality with 4 classes", mg = FALSE)
# cat("\n")
# cat("\n")

#HighProb(lcaETH_C3cl4, sizelca4_ETH, orden = c(1,2,5,3,4),  
#         title = "Response categories probabilities and class size for\n 4-classes Ethnic and race equality model")
cat("\n")
cat("\n")

#----ETH 5 groups----
classes5ETH <- c("Fully egalitarian",
                "Political non-egalitarian",
                "Non-egalitarian",
                "Employment non-egalitarian",
                "Random response")
orden5ETH <- c(5,4,2,1,3)
lcaETH_C3cl5 <- lcaETH$ETH_lca_C3cl5.out$parameters$probability.scale %>% 
  rename_with(~ c("Class", "value")[which(c("LatentClass", "est") == .x)], .cols = c("LatentClass", "est")) %>% 
  mutate_at( c("param", "category", "Class"), ~ as.factor(.x)) %>% 
  mutate(Class = factor(Class, levels = orden5ETH, labels = classes5ETH))

counts5ETH <- full_join(lcaETH$ETH_lca_C3cl5.out$class_counts$modelEstimated,
                    lcaETH$ETH_lca_C3cl5.out$class_counts$mostLikely,by = c("class"))
counts5ETH  %>% 
  mutate(class = factor(class, levels = orden5ETH, labels = classes5ETH),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes 5-class Ethnic and race equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))


sizelca5_ETH <- lcaETH$ETH_lca_C3cl5.out$class_counts$modelEstimated %>% dplyr::select(-count) %>% 
  rename_with(~ c("Ethnic", "Class")[which(c("proportion", "class") == .x)], .cols = c("proportion", "class")) %>% 
  mutate(Class = factor(Class, levels = orden5ETH, labels = classes5ETH)) %>% 
  reshape2::melt(id.vars = c("Class"), variable.name = "Group") %>% 
  dplyr::arrange(Group) %>% 
  dplyr::group_by(Group) %>%
  dplyr::mutate(countT= sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(Class) %>%
  dplyr::mutate(per=value/countT) %>% 
  dplyr::select(Group, Class, per) 

VarClass(lcaETH_C3cl5) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.8, "Myblue", 
                                                 ifelse(value < 0.8 & value > 0.49, "Mygreen","Myred")))) %>% 
  reshape2::dcast(param ~ Class) %>% 
  kbl(caption = paste0("Probabilities to agree each item 5-class Ethnic and race equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",5)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:6, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
# graphclass(lcaETH_C3cl5, nclass = 5, orden = c(1,2,3,5,4), 
#            title = "LCA Ethnic and race equality with 5 classes", mg = FALSE)
# cat("\n")
# cat("\n")

#HighProb(lcaETH_C3cl5, sizelca5_ETH, orden = c(1,2,3,5,4),  
#         title = "Response categories probabilities and class size for\n 5-classes Ethnic and race equality model")
cat("\n")
cat("\n")

#----ETH 6 groups----
classes6ETH <- c("Fully egalitarian",
                "Political non-egalitarian",
                "Non-egalitarian",
                "Employment non-egalitarian",
                "Respect non-egalitarian",
                "Random response")
orden6ETH <- c(1,3,2,4,6,5)
lcaETH_C3cl6 <- lcaETH$ETH_lca_C3cl6.out$parameters$probability.scale %>% 
  rename_with(~ c("Class", "value")[which(c("LatentClass", "est") == .x)], .cols = c("LatentClass", "est")) %>% 
  mutate_at( c("param", "category", "Class"), ~ as.factor(.x)) %>% 
  mutate(Class = factor(Class, levels = orden6ETH, labels = classes6ETH))

counts6ETH <- full_join(lcaETH$ETH_lca_C3cl6.out$class_counts$modelEstimated,
                    lcaETH$ETH_lca_C3cl6.out$class_counts$mostLikely,by = c("class"))
counts6ETH  %>% 
  mutate(class = factor(class, levels = orden6ETH, labels = classes6ETH),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes 6-class Ethnic and race equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))


sizelca6_ETH <- lcaETH$ETH_lca_C3cl6.out$class_counts$modelEstimated %>% dplyr::select(-count) %>% 
  rename_with(~ c("Ethnic", "Class")[which(c("proportion", "class") == .x)], .cols = c("proportion", "class")) %>% 
  mutate(Class = factor(Class, levels = orden6ETH, labels = classes6ETH)) %>% 
  reshape2::melt(id.vars = c("Class"), variable.name = "Group") %>% 
  dplyr::arrange(Group) %>% 
  dplyr::group_by(Group) %>%
  dplyr::mutate(countT= sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(Class) %>%
  dplyr::mutate(per=value/countT) %>% 
  dplyr::select(Group, Class, per) 

VarClass(lcaETH_C3cl6) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.8, "Myblue", 
                                                 ifelse(value < 0.8 & value > 0.49, "Mygreen","Myred")))) %>% 
  reshape2::dcast(param ~ Class) %>% 
  kbl(caption = paste0("Probabilities to agree each item 6-class Ethnic and race equality model"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",6)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:7, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
# graphclass(lcaETH_C3cl6, nclass = 6, orden = c(1,2,3,5,4), 
#            title = "LCA Ethnic and race equality with 6 classes", mg = FALSE)
# cat("\n")
# cat("\n")

#HighProb(lcaETH_C3cl6, sizelca6_ETH, orden = c(1,2,3,5,4),  
#         title = "Response categories probabilities and class size for\n 6-classes Ethnic and race equality model")
cat("\n")
cat("\n")
```

Similarly to the previous scale, the three, four, five and six class models were investigated profoundly. Is not easy to choose the best model fit without doing a full analysis. There are some patterns that can be clearly identified in all the models, Class 1 with sizes of 75.9%, 75.4%, 74% and 78.8% in each model respectively, the estimated probabilities to agree for this latent class, the **Fully egalitarian** group, for all four first items are higher than 0.99 and 0.83 for the item All ethnic and racial groups should be encouraged to run in elections.  

The second class called **Political non-egalitarian** (Class 2 in figure \@ref(fig:compare2)) is identified in all the models and have estimated probabilities to agree to the first 2 items higher than 0.93. For the next two items, the estimated probabilities to agree are around 0.66 and 0.75 in all models and for the last item the probability decreases to 0.5. The class size differs in all four models, 21.7%, 16%, 15% and 8.4% in the 3, 4, 5, 6-class model respectively.  

The third class that can be seen with a similar pattern in all the models is called **Non-egalitarian**, this class appears from the 3-class model. The pattern of this class is basically showing lower estimated conditional probabilities to agree to any of these statements, no greater than 0.13. The estimated sizes for this class are 2.4%, 6.5%, 7.1% and 7% in each model respectively.  

The fifth and sixth classes identified in the models differ in all the countries, nevertheless, one class appears to be consistent in the 5-class model, where this class called **Employment non-egalitarian** has low conditional probabilities to agree (0.2) to the item *All ethnic and racial groups should have an equal chance to get good jobs*.  


```{r compare2, fig.width=6, fig.height=6, fig.cap="Comparative conditional probabilities to agree in 3 to 6 latent class global models for Students' endorsement of equal rights for all ethnic/racial groups", fig.pos='H'}

g0 <- ClassGraph(lcaETH_C3cl3, sizelca3_ETH, orden = c(1,2,5,3,4), title = "3-classes", selected = c(1,2,3))
g1 <- ClassGraph(lcaETH_C3cl4, sizelca4_ETH, orden = c(1,2,5,3,4), title = "4-classes", selected = c(1,2,3,4))
g2 <- ClassGraph(lcaETH_C3cl5, sizelca5_ETH, orden = c(1,2,5,3,4), title = "5-classes", selected = c(1,2,3,4,5))
g3 <- ClassGraph(lcaETH_C3cl6, sizelca6_ETH, orden = c(1,2,5,3,4), title = "6-classes", selected = c(1,2,3,4,5,6))

g<-arrangeGrob(g0,g1,g2,g3, ncol = 2)#, top = "Conditional probabilities to agree to attitudes towards \nethnic and race equality scale")
grid.draw(g)
cat("\n")
cat("\n")

```

The main three classes in the solutions with three, four and five classes do not strongly differ from other models, and the remaining classes are not informative at all or with very small size; using this as a criterion, one can prefer a four-class solution. In table \@ref(tab:bestfit21) the conditional probabilities to agree are shown. These values are very close to 1 in the first class, Fully egalitarian. Similar values are obtained for the first two items in the second class Political non-egalitarian; the next three items start decreasing the conditional probability to agree from 0.76 to 0.46. Class sizes shown in \@ref(tab:bestfit22) indicate that proportions of unconditional probabilities even though are not exactly the same, have similar values among the Model estimated and Most likely classification.  


```{r bestfit21}
lcaETH_C3cl4$orden = rep(c(1,2,4,5,3), each = 2)
VarClass(lcaETH_C3cl4) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(orden, param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.75, "Myblue", 
                                                 ifelse(value < 0.75 & value >= 0.25, "Mygreen","Myred")))) %>% 
  reshape2::dcast(orden + param ~ Class) %>% arrange(orden) %>% select(-orden) %>% 
  kbl(caption = paste0("Probabilities to agree each item four-class global model for Students' endorsement of equal rights for all ethnic/racial groups scale"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:5, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
# graphclass(lcaETH_C3cl4, nclass = 4, orden = c(1,2,3,5,4), 
#            title = "LCA Ethnic and race equality with 4 classes", mg = FALSE)
# cat("\n")
# cat("\n")
```


```{r bestfit22}
counts4ETH  %>% 
  mutate(class = factor(class, levels = orden4ETH, labels = classes4ETH),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes four-class global model for Students' endorsement of equal rights for all ethnic/racial groups scale"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))
```

```{r higprob2, eval=FALSE}

HighProb(lcaETH_C3cl4, sizelca4_ETH, orden = c(1,2,5,3,4))#, title = "Response categories probabilities and class size for\n 4-classes Ethnic and race equality model")
cat("\n")
cat("\n")
```


### Country comparability  

To evaluate the country comparability, the classes that were found in the independent models were identified to later check how many of them could be tested for comparability using a multigroup latent class model.

- 3-class model:
  1. Fully egalitarian: ALL COUNTRIES
  2. Political non-egalitarian: ALL COUNTRIES
  3. Non-egalitarian: ALL COUNTRIES

- 4-class model:
  1. Fully egalitarian: ALL COUNTRIES
  2. Political non-egalitarian: ALL COUNTRIES
  3. Non-egalitarian: ALL COUNTRIES
  4. Employment non-egalitarian: MLT

\newpage  

- 5-class model:
  1. Fully egalitarian: ALL COUNTRIES
  2. Political non-egalitarian: ALL COUNTRIES
  3. Non-egalitarian: ALL COUNTRIES
  4. Employment non-egalitarian: MLT
  5. Random response: Not identified in individual country models

With 3 classes, all classes are very interpretable. With 5 classes, a random response class is identified, which is not interpretable. With 4 classes, the Employment non-egalitarian class is present in just one country, which is not representative. 

With a 4-classes model, three main classes are identified across countries. All the classes are present in all countries which means that is the best model for comparability. One remaining class can be freely estimated that varies in each country and/or with a class size of 0.


```{r, eval = FALSE}
load("data/MplusModels_CntCov.RData")

Modelfit("CntCovETH", title = "Model fit statistics Ethnic and race equal rights model with Country as covariate") %>% 
  pack_rows("Country numerical covariate",1,3) %>% 
  row_spec(2, background = "lightgray") %>% 
  row_spec(c(2), bold = TRUE) %>% 
  footnote(general = "The best loglikelihood value was not replicated in 6-class model. The solution may not be trustworthy due to local maxima.") %>% 
  print()
cat('\n')
cat('\n')
```

```{r, eval = FALSE}
#----ETH 5 groups
CntCovclasses5ETH <- c("Fully egalitarian",
                "No equal political rights",
                "Not egalitarian",
                "No good jobs",
                "Not involved")
CntCovorden5ETH <- c(2,3,1,4,5)
CntCovETH_C3cl5 <- CntCovETH$ETH_CntCov_C3cl5.out$parameters$probability.scale %>% 
  rename_with(~ c("Class", "value")[which(c("LatentClass", "est") == .x)], .cols = c("LatentClass", "est")) %>% 
  mutate_at( c("param", "category", "Class"), ~ as.factor(.x)) %>% 
  mutate(Class = factor(Class, levels = CntCovorden5ETH, labels = CntCovclasses5ETH))

CntCovcounts5ETH <- full_join(CntCovETH$ETH_CntCov_C3cl5.out$class_counts$modelEstimated,
                        CntCovETH$ETH_CntCov_C3cl5.out$class_counts$mostLikely,by = c("class"))
CntCovcounts5ETH  %>% 
  mutate(class = factor(class, levels = CntCovorden5ETH, labels = CntCovclasses5ETH),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes 5-class Ethnic and race equal rights model covariate"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))

sizeCntCov5_ETH <- CntCovETH$ETH_CntCov_C3cl5.out$class_counts$modelEstimated %>% dplyr::select(-count) %>% 
  rename_with(~ c("Gender", "Class")[which(c("proportion", "class") == .x)], .cols = c("proportion", "class")) %>% 
  mutate(Class = factor(Class, levels = CntCovorden5ETH, labels = CntCovclasses5ETH)) %>% 
  reshape2::melt(id.vars = c("Class"), variable.name = "Group") %>% 
  dplyr::arrange(Group) %>% 
  dplyr::group_by(Group) %>%
  dplyr::mutate(countT= sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(Class) %>%
  dplyr::mutate(per=value/countT) %>% 
  dplyr::select(Group, Class, per) 

VarClass(CntCovETH_C3cl5) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.8, "Myblue", 
                                                 ifelse(value < 0.8 & value > 0.49, "Mygreen","Myred")))) %>% 
  reshape2::dcast(param ~ Class) %>% 
  kbl(caption = paste0("Probabilities to agree each item 5-class Ethnic and race equal rights model with country as covariate"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",5)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:6, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
# graphclass(CntCovETH_C3cl5, nclass = 5, orden = c(1,2,5,3,4,6), title = "LCA Ethnic and race equal rights with 5 classes with country as covariate", mg = FALSE)
# cat("\n")
# cat("\n")

HighProb(CntCovETH_C3cl5, sizeCntCov5_ETH,  orden = c(1,2,5,3,4,6), title = "Response categories probabilities and class size for\n 5-classes Ethnic and race equal rights model with country as covariate") 
cat("\n")
cat("\n")

```


#### Country multigroup analysis 

In table \@ref(tab:mgmodelfit2) different models with multigroup analysis are tested, first the more restricted model is evaluated, the complete homogeneity. In this model, all conditional and unconditional probabilities are fixed to be equal across the groups. Then, the partial homogeneity is tested where only the conditional probabilities are constrained to be equal across the groups, and the class sizes are estimated freely. 

The second approach of partial homogeneity is tested too but not included in the table, where only the conditional probabilities for the three common classes identified are constrained across groups, and the remaining are freely estimated along with the unconditional probabilities. 

Finally, the complete heterogeneous model is tested, where not only the unconditional probabilities are estimated freely but all the conditional probabilities as well. In the last two models the best log-likelihood is not replicated, this means that the solution may not be trustworthy due to local maxima. These results cannot be considered valid.

<!-- \blandscape   -->

```{r mgmodelfit2}

ModelfitMGCntry("MGCntry_ETH", title = "Multigroup model fit statistics, global model with four-classes for Students' endorsement of equal rights for all ethnic/racial groups scale",
                filterval=c(4), fontn = 10) %>%
  pack_rows("Four-class model", 1, 3) %>%
  pack_rows("Complete homogeneity", 1, 1) %>%
  pack_rows("Partial homogeneity", 2, 2) %>%
  pack_rows("Complete heterogeneity", 3, 3) %>%
  #pack_rows("Five-class model", 4, 5) %>%
  #pack_rows("Complete homogeneity", 4, 4) %>%
  #pack_rows("Partial homogeneity", 5, 5) %>%
  #pack_rows("Complete heterogeneity",6, 6) %>%
  row_spec(c(2), bold = TRUE) %>%
  footnote(general = "The best log-likelihood value was not replicated for the following models:",
           number = c(" 4-class Complete heterogeneity model"#, "The model estimation did not terminate normally due to an insufficient number of e steps: 5-class Complete heterogeneity; "
                      )) %>%
  print()
cat('\n')
cat('\n')


```

<!-- \elandscape   -->

Just by looking at the valid results, the partial homogeneity where all conditional probabilities are constrained to be equal across groups shows a better fit compared to the more restricted model, with complete homogeneity. It is valid to indicate that the 4 classes identified do not share the same unconditional probabilities (class sizes) across the groups, but the conditional probabilities can be considered as equal in all groups.


```{r MGchom2, fig.width=3, fig.height=2, fig.cap="Conditional probabilities to agree in a 4-class complete homogeneous multigroup model for Students' endorsement of equal rights for all ethnic/racial groups scale", fig.pos='H'}
MGCntyCHo_ETH_C3cl4 <- data.frame(full_join(MGCntry_ETH$ETH_MGCnty_C3cl4_3CHom.out$parameters$probability.scale %>% 
                                              mutate(class = as.numeric(substr(LatentClass,3,3))), 
                                            MGCntry_ETH$ETH_MGCnty_C3cl4_3CHom.out$class_counts$modelEstimated %>% filter(variable == "C") %>% 
                                              select(class, proportion), by = c("class"))) %>% 
  mutate(ClassN = case_when(class == 1 ~ 2, class == 2 ~ 3, class == 3 ~ 4, class == 4 ~ 1)#,
         #ClassN = factor(ClassN, levels = 1:4, labels = classes4HetMG)
         )
cHomMGCntryETH <- MGCntyCHo_ETH_C3cl4 %>% filter(category == 1) %>% 
  mutate(ClassN = factor(ClassN),
         param = factor(param, levels = unique(MGCntyCHo_ETH_C3cl4$param)[c(1,2,5,3,4,6)]),
         proportion = ifelse(param == "ETH4", proportion, NA)) %>% 
  ggplot() +
  geom_point(aes(x = param, y = est, group = ClassN, color = ClassN), size = 1) +
  geom_line(aes(param, est, group = ClassN, linetype = ClassN, color = ClassN)) +
  geom_text(aes(x = param, y = est, label = scales::percent(proportion, accuracy = 0.1), color = ClassN), size = 2, 
            nudge_x = -0.15, nudge_y = 0.1) +
  scale_fill_grey() + theme_bw() +
  #ggtitle("Classes in a complete homogeneous multigroup model, \nEthnic and race equal rights scale") +
  theme(legend.position = "none", 
        title = element_text(size = 9),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 7),
        axis.text.x = element_text(angle = 90, size = 7, vjust = 0.5, hjust = 0)) +
  scale_y_continuous(breaks = c(0.25,0.5,0.75)) +
  labs(y="Response probabilities", linetype = "Latent Classes", color = "Latent Classes")  +
  scale_colour_manual(values=cbPalette[c(1,2,3,4)]) +
  scale_linetype_manual(values=c(1,2,3,4))
#scale_color_brewer(type = "qual", palette = "Dark2") 
cHomMGCntryETH
cat("\n")
cat("\n")
```


Similarly to the gender equality scale, it was not possible to establish complete invariance between the countries for the ethnic and race equality scale, the Fully egalitarian class concentrate the 75.4% of the population when looking at the global model (figure \@ref(fig:MGchom2)), Political non-egalitarian class includes 16% of the population,  6.5% of the population belongs to Employment non-egalitarian class and 2.2% to the Non-egalitarian class. 


```{r}
MGCntyPHo_ETH_C3cl4 <- data.frame(full_join(MGCntry_ETH$ETH_MGCnty_C3cl4_2PHom.out$parameters$probability.scale %>% mutate(C = as.numeric(substr(LatentClass,3,3))), 
                                            MGCntry_ETH$ETH_MGCnty_C3cl4_2PHom.out$class_counts$modelEstimated.patterns %>% 
                                              mutate(LatentClass = paste0(G,".",C)) %>% 
                                              group_by(G) %>% mutate(countT= sum(count, na.rm = TRUE)) %>% ungroup() %>% 
                                              mutate(proportion=round(count/countT*100, 1)) %>% select(G, LatentClass, C, proportion), by = c("C"))) %>% 
  mutate(Group = G,
         ClassN = case_when(C == 1 ~ 4, C == 2 ~ 2, C == 3 ~ 1, C == 4 ~ 3),
         Group = as.character(factor(Group, levels = 1:length(cnt), labels = cnt))#,
         #ClassN = factor(ClassN, levels = 1:4, labels = classes4HetMG)
         )
```

```{r classETH}
MGCntyPHo_ETH_C3cl4 %>% 
  select(Group, ClassN, proportion) %>% 
  unique() %>% 
  mutate(ClassN = factor(ClassN, levels = 1:4, labels = classes4ETH)) %>% 
  reshape2::dcast(Group ~ ClassN) %>% 
  kbl(caption = paste0("Class sizes for partial homogeneity model for Students' endorsement of equal rights for all ethnic/racial groups scale"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:5, width = "5em") %>%  
  print()
```

These proportions vary when looking at the partial invariance (indicated in Table \@ref(tab:classETH)) where only Finland, Norway and Sweden class sizes include more than the 70% of the population for the Fully egalitarian class. The size of Political non-egalitarian class includes in average 26% of the population, more specifically in Belgium (Flemish), Bulgaria and Latvia sizes are higher than 37.9% of the population. Similar to the previous class, the number of members in the Employment non-egalitarian class differs from the global model to the partial invariance model, these numbers go from 5.2% in Finland to 22.5% in Malta. Finally, the Non-egalitarian class could be the only one that maintains the same proportion between the global and the country-specific size around 2.2%, with the exception of Bulgaria that includes 5.4% of the population. 

```{r MGphom2, fig.width=6, fig.height=6, fig.cap="Conditional probabilities to agree in a 4-class partial homogeneous multigroup model for Students' endorsement of equal rights for all ethnic/racial groups scale", fig.pos='H', eval=FALSE}

pHomMGCntryETH <- MGCntyPHo_ETH_C3cl4 %>% filter(category == 1) %>% plyr::rbind.fill(data.frame(param = "Size", Group = "Sweden")) %>% 
  mutate(ClassN = factor(ClassN),
         param = factor(param, levels = c(unique(MGCntyPHo_ETH_C3cl4$param), "Size")[c(1,2,5,3,4,6,7)]), 
         proportion = ifelse(param == "ETH4", proportion, NA)) %>% 
  ggplot() +
  geom_point(aes(x = param, y = est, group = ClassN, color = ClassN), size = 1) +
  geom_line(aes(param, est, group = ClassN, linetype = ClassN, color = ClassN)) +
  geom_text(aes(x = param, y = est, label = scales::percent(proportion/100, accuracy = 0.1), color = ClassN), size = 2, 
            nudge_x = 0.8, nudge_y = 0) +
  facet_wrap(Group ~ .) +
  scale_fill_grey() + theme_bw() +
  #ggtitle("Classes in a 4-class country multigroup partially homogeneous model, \nStudents' endorsement of equal rights for all ethnic/racial groups scale") +
  theme(legend.position = "none", 
        title = element_text(size = 9),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 7),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.text.x = element_text(angle = 90, size = 7, vjust = 0.5, hjust = 0)) +
  scale_y_continuous(breaks = c(0.25,0.5,0.75), limits = c(0,1)) +
  geom_hline(yintercept = c(0.25,0.5,0.75), color = "gray", size = 0.2) +
  labs(y="Response probabilities", linetype = "Latent Classes", color = "Latent Classes")  +
  scale_colour_manual(values=cbPalette[c(1,2,3,6)]) +
  scale_linetype_manual(values=c(1,2,3,6))
#scale_color_brewer(type = "qual", palette = "Dark2") 
pHomMGCntryETH
cat("\n")
cat("\n")
```


```{r, fig.width=6, fig.height=5, fig.cap="Conditional probabilities to agree in a partial homogeneous multigroup model with constrains, \nEthnic and race equal rights scale", fig.pos='H', eval=FALSE}
MGCnty4PH_ETH_C3cl4 <- data.frame(left_join(MGCntry_ETH$ETH_MGCnty_C3cl4_2.2PHom.out$parameters$probability.scale, 
                                            MGCntry_ETH$ETH_MGCnty_C3cl4_2.2PHom.out$class_counts$modelEstimated.patterns %>% 
                                              mutate(LatentClass = paste0(G,".",C)) %>% 
                                              group_by(G) %>% mutate(countT= sum(count, na.rm = TRUE)) %>% ungroup() %>% 
                                              mutate(proportion=round(count/countT*100, 1)) 
                                            %>% select(LatentClass, proportion), by = c("LatentClass"))) %>% 
  mutate(Group = sub('\\..*', '', LatentClass),
         Class = sub('*.{1,2}\\.', '', LatentClass),
         ClassN = case_when(LatentClass == 1.1 ~ 1, LatentClass == 1.2 ~ 2, LatentClass == 1.3 ~ 3, LatentClass == 1.4 ~ 5,
                            LatentClass == 2.1 ~ 1, LatentClass == 2.2 ~ 2, LatentClass == 2.3 ~ 3, LatentClass == 2.4 ~ 5,
                            LatentClass == 3.1 ~ 1, LatentClass == 3.2 ~ 2, LatentClass == 3.3 ~ 3, LatentClass == 3.4 ~ 5,
                            LatentClass == 4.1 ~ 1, LatentClass == 4.2 ~ 2, LatentClass == 4.3 ~ 3, LatentClass == 4.4 ~ 5,
                            LatentClass == 5.1 ~ 1, LatentClass == 5.2 ~ 2, LatentClass == 5.3 ~ 3, LatentClass == 5.4 ~ 5,
                            LatentClass == 6.1 ~ 1, LatentClass == 6.2 ~ 2, LatentClass == 6.3 ~ 3, LatentClass == 6.4 ~ 5,
                            LatentClass == 7.1 ~ 1, LatentClass == 7.2 ~ 2, LatentClass == 7.3 ~ 3, LatentClass == 7.4 ~ 5,
                            LatentClass == 8.1 ~ 1, LatentClass == 8.2 ~ 2, LatentClass == 8.3 ~ 3, LatentClass == 8.4 ~ 5,
                            LatentClass == 9.1 ~ 1, LatentClass == 9.2 ~ 2, LatentClass == 9.3 ~ 3, LatentClass == 9.4 ~ 5,
                            LatentClass == 10.1 ~ 1, LatentClass == 10.2 ~ 2, LatentClass == 10.3 ~ 4, LatentClass == 10.4 ~ 5,
                            LatentClass == 11.1 ~ 1, LatentClass == 11.2 ~ 2, LatentClass == 11.3 ~ 4, LatentClass == 11.4 ~ 5,
                            LatentClass == 12.1 ~ 1, LatentClass == 12.2 ~ 2, LatentClass == 12.3 ~ 4, LatentClass == 12.4 ~ 5,
                            LatentClass == 13.1 ~ 1, LatentClass == 13.2 ~ 2, LatentClass == 13.3 ~ 4, LatentClass == 13.4 ~ 5,
                            LatentClass == 14.1 ~ 1, LatentClass == 14.2 ~ 2, LatentClass == 14.3 ~ 4, LatentClass == 14.4 ~ 5),
         Group = as.character(factor(Group, levels = 1:length(cnt), labels = cnt))#,
         #ClassN = factor(ClassN, levels = 1:4, labels = classes4HetMG)
         )

MGCnty4PH_ETH_C3cl4 %>% arrange(Group) %>% select(Group, param, est, ClassN, category) %>% VarClass(orden = c(1,2,5,3,4,6)) %>% 
  bind_rows(MGCnty4PH_ETH_C3cl4 %>% select(Group, proportion, ClassN) %>% unique() %>% mutate(param = "Class Size (\\%)", category = "1") %>% 
              rename(est = proportion)) %>% 
  group_by(ClassN, Group, param) %>% filter(category == 1) %>%
  select(param, ClassN, Group, est) %>%
  mutate(value = cell_spec(est, color = ifelse(param == "Class Size (\\%)", "black", ifelse(est > 0.8, "Myblue", "Myred")))) %>%
  reshape2::dcast(Group + param ~  ClassN) %>% 
  mutate(param = factor(param, levels = unique(param)[c(1,2,3,6,4,5,7)])) %>% arrange(Group,param) %>% 
  kbl(caption = "Probabilities to Agree in a 5-class country multigroup model with 2 constrained classes,\nStudents' endorsement of gender equality", 
      booktabs = TRUE, longtable = TRUE,
      align = "llrrrrr", row.names = FALSE, digits = 3, escape = FALSE#, col.names = c("Item", classes2EU, classes2LA)
      ) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 10) %>%
  column_spec(1, width = "5em") %>%
  column_spec(2, width = "16em") %>%
  column_spec(3:7, width = "3em") %>%
  collapse_rows(1, valign = "top") %>%
  print()
cat('\n')
cat('\n')

MGCnty4PH_ETH_C3cl4 %>% filter(category == 1) %>% plyr::rbind.fill(data.frame(param = "Size", Country = "Sweden")) %>% 
  mutate(ClassN = factor(ClassN),
         param = factor(param, levels = c(unique(MGCnty4PH_ETH_C3cl4$param), "Size")[c(1,2,5,3,4,6,7)]),
         proportion = ifelse(param == "ETH6", proportion, NA)) %>% 
  ggplot() +
  geom_point(aes(x = param, y = est, group = ClassN, color = ClassN), size = 1) +
  geom_line(aes(param, est, group = ClassN, linetype = ClassN, color = ClassN)) +
  geom_text(aes(x = param, y = est, label = scales::percent(proportion/100, accuracy = 0.1), color = ClassN), size = 2, 
            nudge_x = 0.8, nudge_y = 0.0) +
  facet_wrap(Group ~ .) +
  scale_fill_grey() + theme_bw() +
  ggtitle("Classes in a 5-class multigroup model with 2 constrained classes,\nStudents' endorsement of gender equality") +
  theme(legend.position = "none", 
        title = element_text(size = 9),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 7),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.text.x = element_text(angle = 90, size = 7, vjust = 0.5, hjust = 0)) +
  scale_y_continuous(breaks = c(0.25,0.5,0.75), limits = c(0,1)) +
  geom_hline(yintercept = c(0.25,0.5,0.75), color = "gray", size = 0.2) +
  labs(y="Response probabilities", linetype = "Latent Classes", color = "Latent Classes")  +
  scale_colour_manual(values=cbPalette[c(5,2,1,3)]) +
  scale_linetype_manual(values=c(5,2,1,3))
cat("\n")
cat("\n")
```


### Confirmatory Latent Class Analysis

The confirmatory model was performed by establishing some constraints based on the previous research. For the Students' endorsement of equal rights for all ethnic/racial groups, two hypotheses were tested. First, that students classified in the Fully egalitarian and Political non-egalitarian classes would agree equally to first two items in the scale *All ethnic and racial groups should have equal chance to get good education* and *All ethnic and racial groups should have an equal chance to get good jobs*. 
The other hypothesis test that beside the previous equality, students that would be classified in the Non-egalitarian class would have the remaining probability to agree to all the items, which would mean to be unlikely to agree.  


```{r cfa2, results='hide'}
#----ETH 3 groups----
orden4ETH <- c(1,2,3,4)
ConfUnstlcaETH <- ConflcaETH$ETH_Conflca_C3cl4.out$parameters$unstandardized %>% 
  rename_with(~ c("Parameter", "Classo", "value")[which(c("param", "LatentClass", "est") == .x)], .cols = c("param", "LatentClass", "est")) %>% 
  mutate_at( c("Parameter", "paramHeader", "Classo"), ~ as.factor(.x)) %>%  
  mutate(Class = case_when(Parameter == "C#1" ~ 1,
                        Parameter == "C#2" ~ 2,
                        Parameter == "C#3" ~ 3,
                        TRUE ~ as.numeric(Classo)),
         Parameter = case_when(Parameter %in% c("C#1", "C#2", "C#3") ~ "Means",
                        TRUE ~ as.character(Parameter))) %>% 
  mutate(Class = factor(Class, levels = orden4ETH, labels = classes4ETH)) %>% 
  group_by(Class, Parameter) %>% 
  select(Parameter, Class, value) 

ConfUnstlcaETH$orden <- c(rep(c(1,2,4,5,3),4), rep(6,3))
ConfUnstlcaETH %>% 
  reshape2::dcast(orden + Parameter ~ Class) %>% arrange(orden) %>% select(-orden) %>% 
  kbl(caption = paste0("Thresholds 4-class Confirmatory LCA Attitude towards ethnic and race equal rights scale"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 3, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "4em") %>%  
  column_spec(2:5, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()


SummConflcaETH <- ConflcaETH$ETH_Conflca_C3cl4.out$parameters$probability.scale %>% 
  rename_with(~ c("Class", "value")[which(c("LatentClass", "est") == .x)], .cols = c("LatentClass", "est")) %>% 
  mutate_at( c("param", "category", "Class"), ~ as.factor(.x)) %>% 
  mutate(Class = factor(Class, levels = orden4ETH, labels = classes4ETH))

ConfcountsETH <- full_join(ConflcaETH$ETH_Conflca_C3cl4.out$class_counts$modelEstimated,
                        ConflcaETH$ETH_Conflca_C3cl4.out$class_counts$mostLikely,by = c("class"))
ConfcountsETH  %>% 
  mutate(class = factor(class, levels = orden4ETH, labels = classes4ETH),
         proportion.x = scales::percent(proportion.x,accuracy = 0.1),
         proportion.y = scales::percent(proportion.y,accuracy = 0.1)) %>% arrange(desc(count.y)) %>% 
  kbl(col.names = c("Class", "Counts", "Proportion", "Counts", "Proportion"),
      caption = paste0("Class sizes 4-class Students' endorsement of equal rights for all ethnic/racial groups scale"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 1, escape = TRUE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  add_header_above(c(" " = 1 , "Model estimated" = 2, "Most likely" = 2))
```


```{r probconf2}
sizeConflca_ETH <- ConflcaETH$ETH_Conflca_C3cl4.out$class_counts$modelEstimated %>% dplyr::select(-count) %>% 
  rename_with(~ c("Gender", "Class")[which(c("proportion", "class") == .x)], .cols = c("proportion", "class")) %>% 
  mutate(Class = factor(Class, levels = orden4ETH, labels = classes4ETH)) %>% 
  reshape2::melt(id.vars = c("Class"), variable.name = "Group") %>% 
  dplyr::arrange(Group) %>% 
  dplyr::group_by(Group) %>%
  dplyr::mutate(countT= sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(Class) %>%
  dplyr::mutate(per=value/countT) %>% 
  dplyr::select(Group, Class, per) 

SummConflcaETH$orden = rep(c(1,2,4,5,3), each = 2)
VarClass(SummConflcaETH) %>% group_by(Class, param) %>% 
  filter(category == 1) %>% 
  select(orden, param, Class, value) %>% 
  mutate(value = cell_spec(value, color = ifelse(value >= 0.8, "Myblue", 
                                                 ifelse(value < 0.8 & value > 0.49, "Mygreen","Myred")))) %>% 
  reshape2::dcast(orden + param ~ Class) %>% arrange(orden) %>% select(-orden) %>% 
  kbl(caption = paste0("Probabilities to agree each item 4-class Confirmatory LCA Students' endorsement of equal rights for all ethnic/racial groups scale"),
      booktabs = TRUE, longtable = TRUE, align = c("l", rep("r",4)), row.names = FALSE, digits = 3, escape = FALSE) %>%
  kable_classic_2(full_width = F) %>% 
  kable_styling(latex_options = c("repeat_header", "HOLD_position"), font_size = 9) %>% 
  column_spec(1, width = "15em") %>%  
  column_spec(2:5, width = "5em") %>%  
  collapse_rows(1, valign = "top") %>% 
  print()
# graphclass(ConflcaETH, nclass = 4, orden = c(1,2,3,5,4), title = "Confirmatory LCA Ethnic and race equal rights with 4 classes", mg = FALSE)
# cat("\n")
# cat("\n")
```

For this test, as stated in table \@ref(tab:probconf2), the conditional probabilities for the first and second latent class are equal for the first two items respectively. On the other hand, probabilities in the Non-egalitarian class are the difference between the total probability to agree (1) and the probability obtained in the Fully egalitarian class. The rest of the conditional probabilities were estimated freely. 

In table \@ref(tab:confm2) the model fit statistics of this model do not differ considerably from the exploratory approach analyzed previously but statistics such as AIC, BIC, and aBIC are still lower in the exploratory approach. This way it cannot be suggested that these restrictions apply to the population. Other hypotheses can be tested until obtaining a better model fit. 

```{r confm2}
FinalCompEthnic <- list("ETH_lca_C3cl4.out" = lcaETH$ETH_lca_C3cl4.out, 
                        "ETH_Conflca_C3cl4.out" = ConflcaETH$ETH_Conflca_C3cl4.out) 
ModelfitConf("FinalCompEthnic", title = "Model fit statistics Confirmatory LCA Students' endorsement of equal rights for all ethnic/racial groups scale") %>% 
  pack_rows("All countries",1,2) %>% 
  print()
cat('\n')
cat('\n')
```



```{r CFAm2, eval=FALSE}
HighProb(SummConflcaETH, sizeConflca_ETH,  orden = c(1,2,5,3,4))#, title = "Conditional probabilities and class size for\n 4-classes Confirmatory LCA model Ethnic and race equal rights scale") 
cat("\n")
cat("\n")
```

\clearpage