---
title: "Dyscalculia, General Intelligence, and Mathematics Grades in Hungary"
output: 
  html_document:
    toc: yes
    toc_float: yes
    df_print: paged
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
library(pacman); p_load(readxl, tidyr, dplyr, tidyverse, psych, ggplot2, cowplot, sandwich, lmtest, ggcorrplot, lavaan, semPlot, lavaanPlot)

HCReg <- function(Formula, Data, HC = "HC5"){
  require(sandwich); require(lmtest)
  BaseModel <- lm(Formula, Data)
  VC <- vcovHC(BaseModel, type = HC)
  coeftest_results <- coeftest(BaseModel, vcov = VC)
  model_summary <- summary(BaseModel)
  r_squared <- model_summary$r.squared
  adj_r_squared <- model_summary$adj.r.squared
  Output <- list(
    coeftest_results = coeftest_results,
    r_squared = r_squared,
    adj_r_squared = adj_r_squared)
  return(Output)}

data <- read_xlsx("SESIncludedData.xlsx")

data$MSZMT_total <- rowSums(data[,62:101]) #MARS Total Score

Transformed <- data %>% gather(Subtest, Score, dpv_tajek:dpv_praktik)
```

# Rationale

The Diszkalkulia Pedagogiai Vizsgalata or DPV is a dyscalculia diagnosis tool used in Hungary. It features 14 subtests which amount to parcels of different numbers of items, in addition to one sumscore for the whole battery (dpv_ossz). Dyscalculia is defined as abnormally low mathematical ability for a person's chronological age and level of intellectual functioning. Despite this, the DPV may be a measure of general intelligence, and this can be assessed by determining the latent correlation with the *g* factor from a gold-standard intelligence assessment like the WISC.

## DPV Test Details

1. Tajek, quality assessment. Binary/boolean indicating if a person has the correct bodily left-right orientation, the correct examiner left-right orientation, understands relations and left-right orientation on a plane, understands relational vocabulary in tenses in time, and relations in time, and understands relationships between numbers and understands relational grammar in numbers. 

2.1. Szaml, counting in 10,000s. Five-point scale earned based on a participant's ability to correctly name numbers, ability to count correctly in ones, tenths, and hundredths in ascending and descending order, ability to correctly place ones, tenths, and hundredths in ascending and descending order, ability to identify numbering ruel mistakes for ones, tenths, and hundredths in ascending and descending order, and whether the participant stuttered or slowed their tempo by ones, tenths, or hundredths when going between numbers, in both ascending and descending order. 

2.2. Szameml, number memory. Binary/boolean indicating if a person's number memory appears weak per the examiner's qualitative judgment. This scale is binary. 

2.3. Szamegyz, number name-digit matching in the range of 100,000. Six-point scale earned based on a participant's ability to avoid making errors when reading and writing five-digit numbers, and to avoid using zero as a placeholder when describing numbers or reading out their digits. 

2.4. Rel, quantitative relations between numerical neighbors in the range of 10,000. Four-point scale earned based on a participant's ability to correctly verbally and visually interpret manufactured quantitative relations, identify the content of a relational signal, and determine whether the neighbors of a given number are correct. 

2.5. Helyiert, designation of a local value in the 10,000 range. 

3.1. Potlas, substitution, decomposition, addition, and subtraction in numbers of 10 and 20, following arithmetic rules. One hundred-point scale with which points are earned by indicating correct understandings of operations, use of equipment, ability to mathematically abstract, to count in tens, and to switch between different operations without error. 

3.2. Muvel, addition and subtraction in 100s, 1000s, analogies, and approximate calculation. Forty-three-point scale with which points are earned by indicating correct readings of specific actions, identifying operation signals and performing with them, properly using devices, correctly interpreting counts, verbally executing operations, performing approximate calculations like rounding and estimating, and avoiding disorientation when switching between operations. 

3.3. Szorzoszt, multiplication and inclusion. Nine-point scale with which points are earned by interpreting operations correctly, and interpreting the relationship between multiplication and inclusion correctly. 

3.4. Irasbeli, written actions. Thirty-eight-point scale with which points are earned by correctly recording what written operations are, reading them out, identifying operation signals correctly, processing written transactions involving operations, using devices correctly, properly interpreting placeholders during operations, properly interpreting the concept of the remainder, retaining the remainder, obtaining a correct answer, and switching to different operations correctly. 

4. Szoveg, text tasks. 27-point scale with which points are earned by properly state mathematical questions and do so without tools or taking additional actions due to errors. 

5. Logszab, mathematical-logical rules. Nine-point scale with which points are earned by correctly interpreting the terms decreasing and increasing, and by understanding and sticking to defined rules which are used to solve a question that the participant is also scored on. 

6. Muvsorr, higher level mathematical concepts and the order of operations. Twenty-two-point scale with which points are earned by correctly interpreting fractions, reading and writing fractions, performing operations with fractions, relating decimal numbers to fractiosn, reading and writing decimal fractions, identifying operation signals with fractions, understanding when operations are wrong with fractions, processing written transactions with fractions, using devices with fractions, understanding placeholders with decimal fractions, retaining remainders with fractions, switching operations when appropriate, and interpreting negative numbers tied to specific examples like debt and thermometers correctly.

7. Praktik, practical knowledge. Binary/boolean qualitatively assessing a participant's ability to identify units of measurement, recall conversion numbers when converting between them, select appropriate operations during unit conversions, and interpret the concepts of districts and areas correctly. 

# Analyses

## Variable Distributions

```{r, fig.width = 12, fig.height = 8, warning = F}
OrdNames <- list(
  "dpv_helyiert" = "Helyiert",
  "dpv_irasbeli" = "Irasbeli", 
  "dpv_logszab"  = "Logszab", 
  "dpv_muvel"    = "Muvel", 
  "dpv_muvsorr"  = "Muvsorr", 
  "dpv_potlas"   = "Potlas", 
  "dpv_praktik"  = "Praktik", 
  "dpv_rel"      = "Rel", 
  "dpv_szamegyz" = "Szamegyz", 
  "dpv_szameml"  = "Szameml", 
  "dpv_szaml"    = "Szaml", 
  "dpv_szorzoszt"= "Szorzoszt", 
  "dpv_szoveg"   = "Szoveg", 
  "dpv_tajek"    = "Tajek")

TestLabeler <- function(variable, value){
  return(OrdNames[value])}

p <- ggplot(Transformed) + 
  geom_histogram(aes(Score, color = Subtest, fill = Subtest), alpha = .5, binwidth = 1) + 
  ggplot2::theme_bw() + 
  theme(legend.position = "none") + 
  facet_wrap(~ Subtest, scales = "free",
             labeller = TestLabeler) + 
  labs(
    x = NULL,
    y = NULL)

  ggdraw(p) + draw_label(x = .76, y = .15, "Subtests of the DPV")
```

## Correlations

First, some particulars of interest

```{r}
cor(data$dpv_ossz, data$WISC_TtlQ)        #DPV Total Score - FSIQ
cor(data$dpv_ossz, data$osztalyzat)       #DPV Total Score - Grades
cor(data$dpv_ossz, data$MSZMT_total)      #DPV Total - MARS Total
cor(data$dpv_ossz, as.numeric(data$SES))  #DPV Total - SES
cor(data$WISC_TtlQ, data$osztalyzat)      #FSIQ - Grades
cor(data$WISC_TtlQ, data$MSZMT_total)     #FSIQ - MARS Total
cor(data$WISC_TtlQ, as.numeric(data$SES)) #FSIQ - SES
```

and now, the DPV correlation matrix:

```{r,fig.width = 10, fig.height = 10}
VarNames <- list("Tajek", "Szaml", "Szameml", "Szamegyz", "Rel", "Helyiert", "Potlas", "Muvel", "Szorzoszt", "Irasbeli", "Szoveg", "Logszab", "Muvsorr", "Praktik")

ggcorrplot(round(cor(data[42:55]), 4),
           hc.order = T, 
           type = "lower",
           outline.col = "white",
           ggtheme = ggplot2::theme_bw(),
           show.legend = F,
           lab = T,
           title = "DPV Subtest Intercorrelations",
           colors = c("orangered", "gold", "steelblue"),
           show.diag = F) +
  theme(
    plot.title = element_text(hjust = .5)) +
  scale_y_discrete(labels = VarNames) +
  scale_x_discrete(labels = VarNames)
```

And now with the WAIS

```{r,fig.width = 10, fig.height = 10}
VarNames <- list("Tajek", "Szaml", "Szameml", "Szamegyz", "Rel", "Helyiert", "Potlas", "Muvel", "Szorzoszt", "Irasbeli", "Szoveg", "Logszab", "Muvsorr", "Praktik", "Szokincs", "Altmegert", "Mozaik", "Kepifog", "Matrix", "Szamterj", "Betuszam", "Kodolas", "Szimbol")

ggcorrplot(round(cor(data[c(42:55, 15, 17, 19, 23, 25, 27, 31, 33, 37, 39)]), 4),
           hc.order = F, 
           type = "lower",
           outline.col = "white",
           ggtheme = ggplot2::theme_bw(),
           show.legend = F,
           lab = T,
           title = "DPV and WISC Subtest Intercorrelations",
           colors = c("orangered", "gold", "steelblue"),
           show.diag = F) +
  theme(
    plot.title = element_text(hjust = .5)) +
  scale_y_discrete(labels = VarNames) +
  scale_x_discrete(labels = VarNames)
```

### EFA

There is no available theoretical factor structure for the DPV, so exploratory results follow. Since the DPV is supposed to measure one thing, a single-factor structure will be defaulted to unless another can be found. Results from both will be presented regardless.

```{r, fig.width = 8, fig.height = 8, warning = F}
fa.parallel(data[42:55])
fa.parallel(sqrt(data[42:55]))
```

```{r}
print(fa(data[42:55], nfactors = 2), cut = .3, sort = T)
print(fa(data[42:55], nfactors = 3), cut = .3, sort = T)
print(fa(data[42:55], nfactors = 4), cut = .3, sort = T)
```

The 4-factor models fit best, but none of these look good at all.

```{r}
ModBi4 <- '
g =~ dpv_tajek + dpv_szaml + dpv_szameml + dpv_szamegyz + dpv_rel + dpv_helyiert + dpv_potlas + dpv_muvel + dpv_szorzoszt + dpv_irasbeli + dpv_szoveg + dpv_logszab + dpv_muvsorr + dpv_praktik

F1 =~ dpv_muvsorr + dpv_praktik
F2 =~ dpv_rel + dpv_irasbeli
F3 =~ dpv_szaml + dpv_tajek'

ModGF4 <- '
F1 =~ dpv_praktik + dpv_muvsorr + dpv_logszab + dpv_szoveg + dpv_szameml + dpv_muvel + dpv_szorzoszt
F2 =~ dpv_szoveg + dpv_rel + dpv_szamegyz + dpv_muvel + dpv_szorzoszt
F3 =~ dpv_logszab + dpv_tajek + dpv_szaml + dpv_helyiert
F4 =~ dpv_irasbeli + dpv_szorzoszt

dpv_szamegyz ~~ 0*dpv_szamegyz
dpv_szorzoszt ~~ 0*dpv_szorzoszt'

FitBi4 <- cfa(ModBi4, std.lv = T, meanstructure = T, orthogonal = T, data, estimator = "DWLS", ordered = c("dpv_szaml", "dpv_szamegyz", "dpv_rel", "dpv_potlas", "dpv_muvel", "dpv_szorzoszt", "dpv_irasbeli", "dpv_szoveg", "dpv_logszab", "dpv_muvsorr"), check.gradient = F, control=list(rel.tol=1e-4))
FitGF4 <- cfa(ModGF4, std.lv = T, meanstructure = T, orthogonal = F, data, estimator = "DWLS", ordered = c("dpv_szaml", "dpv_szamegyz", "dpv_rel", "dpv_potlas", "dpv_muvel", "dpv_szorzoszt", "dpv_irasbeli", "dpv_szoveg", "dpv_logszab", "dpv_muvsorr"), check.gradient = F, control=list(rel.tol=1e-4))

summary(FitBi4, stand = T, fit = T)
summary(FitGF4, stand = T, fit = T)
```

Both fail to converge and simpler models aren't much better. A single-factor model will be used.

### CFAs

WISC first, using the theoretical model that was desired in its construction.

```{r}
ModWISC <- '
wm =~ betuszam_ertek + szamterj_ertek
vc =~ szokincs_ertek + kozjel_ertek + altmegert_ertek
ps =~ kodolas_ertek + szimbol_ertek
pr =~ kepifog_ertek + mozaik_ertek + matrix_ertek
g =~ wm + vc + ps + pr'

FitWISC <- cfa(ModWISC, std.lv = T, meanstructure = T, data)

summary(FitWISC, stand = T, fit = T)
```

DPV next, single-factor, as suggested might work by the parallel analysis, however poor that appeared.

```{r, warning = F}
ModDPV <- '
dpv =~ dpv_tajek + dpv_szaml + dpv_szameml + dpv_szamegyz + dpv_rel + dpv_helyiert + dpv_potlas + dpv_muvel + dpv_szorzoszt + dpv_irasbeli + dpv_szoveg + dpv_logszab + dpv_muvsorr + dpv_praktik'

FitDPV <- sem(ModDPV, std.lv = T, meanstructure = T, data, estimator = "DWLS", ordered = c("dpv_szaml", "dpv_szamegyz", "dpv_rel", "dpv_potlas", "dpv_muvel", "dpv_szorzoszt", "dpv_irasbeli", "dpv_szoveg", "dpv_logszab", "dpv_muvsorr"))

summary(FitDPV, stand = T, fit = T)
```

And then both at once.

```{r, warning = F}
ModBoth <- '
wm =~ betuszam_ertek + szamterj_ertek
vc =~ szokincs_ertek + kozjel_ertek + altmegert_ertek
ps =~ kodolas_ertek + szimbol_ertek
pr =~ kepifog_ertek + mozaik_ertek + matrix_ertek
g =~ wm + vc + ps + pr

dpv =~ dpv_tajek + dpv_szaml + dpv_szameml + dpv_szamegyz + dpv_rel + dpv_helyiert + dpv_potlas + dpv_muvel + dpv_szorzoszt + dpv_irasbeli + dpv_szoveg + dpv_logszab + dpv_muvsorr + dpv_praktik' 

FitBoth <- cfa(ModBoth, std.lv = T, meanstructure = T, data, estimator = "DWLS")

summary(FitBoth, stand = T, fit = T)
```

## Plots

```{r, fig.width = 10, fig.height = 10}
lavaanPlot(model = FitBoth, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"), coefs = TRUE, covs = TRUE, stars = TRUE, stand = T, digits = 1)
```

```{r, fig.width = 14, fig.height = 10, warning = F}
DPVLats <- list(
  DPV = c("dpv_tajek", "dpv_szaml", "dpv_szameml", "dpv_szamegyz", "dpv_rel", "dpv_helyiert", "dpv_potlas", "dpv_muvel", "dpv_szorzoszt", "dpv_irasbeli", "dpv_szoveg", "dpv_logszab", "dpv_muvsorr", "dpv_praktik"),
  wm = c("betuszam_ertek", "szamterj_ertek"),
  vc = c("szokincs_ertek", "kozjel_ertek", "altmegert_ertek"),
  ps = c("kodolas_ertek", "szimbol_ertek"),
  pr = c("kepifog_ertek", "mozaik_ertek", "matrix_ertek"))

semPaths(FitBoth, "model", "std", title = F, residuals = F,
         groups = "DPVLats", pastel = T, mar = c(2, 1, 3, 1),
         bifactor = c("dpv", "g"),
         layout = "tree",
         exoCov = T, intercepts = F)
```

## Robust Regressions

```{r}
HCReg(osztalyzat ~ WISC_TtlQ, data)
HCReg(osztalyzat ~ dpv_ossz, data)
HCReg(osztalyzat ~ as.numeric(SES), data) 
HCReg(osztalyzat ~ MSZMT_total, data)
HCReg(osztalyzat ~ WISC_TtlQ + dpv_ossz, data)
HCReg(osztalyzat ~ WISC_TtlQ + as.numeric(SES), data) 
HCReg(osztalyzat ~ WISC_TtlQ + MSZMT_total, data)
HCReg(osztalyzat ~ WISC_TtlQ + dpv_ossz + as.numeric(SES) + MSZMT_total, data)
```

## SAM Models

```{r}
WISCSAM <- '
g =~ betuszam_ertek + szamterj_ertek + szokincs_ertek + kozjel_ertek + altmegert_ertek + kodolas_ertek + szimbol_ertek + kepifog_ertek + mozaik_ertek + matrix_ertek

osztalyzat ~ g'

DPVSAM <- '
dpv =~ dpv_tajek + dpv_szaml + dpv_szameml + dpv_szamegyz + dpv_rel + dpv_helyiert + dpv_potlas + dpv_muvel + dpv_szorzoszt + dpv_irasbeli + dpv_szoveg + dpv_logszab + dpv_muvsorr + dpv_praktik

osztalyzat ~ dpv'

SESSAM <- '
osztalyzat ~ SES'

MARSSAM <- '
MARS =~ mszmt1 + mszmt2 + mszmt3 + mszmt4 + mszmt5 + mszmt6 + mszmt7 + mszmt8 + mszmt9 + mszmt10 + mszmt11 + mszmt12 + mszmt13 + mszmt14 + mszmt15 + mszmt16 + mszmt17 + mszmt18 + mszmt19 + mszmt20 + mszmt21 + mszmt22 + mszmt23 + mszmt24 + mszmt25 + mszmt26 + mszmt27 + mszmt28 + mszmt29 + mszmt30 + mszmt31 + mszmt32 + mszmt33 + mszmt34 + mszmt35 + mszmt36 + mszmt37 + mszmt38 + mszmt39 + mszmt40

osztalyzat ~ MARS'

WISCDPVSAM <- '
g =~ betuszam_ertek + szamterj_ertek + szokincs_ertek + kozjel_ertek + altmegert_ertek + kodolas_ertek + szimbol_ertek + kepifog_ertek + mozaik_ertek + matrix_ertek

dpv =~ dpv_tajek + dpv_szaml + dpv_szameml + dpv_szamegyz + dpv_rel + dpv_helyiert + dpv_potlas + dpv_muvel + dpv_szorzoszt + dpv_irasbeli + dpv_szoveg + dpv_logszab + dpv_muvsorr + dpv_praktik

g ~~ dpv

osztalyzat ~ g + dpv'

WISCSESSAM <- '
g =~ betuszam_ertek + szamterj_ertek + szokincs_ertek + kozjel_ertek + altmegert_ertek + kodolas_ertek + szimbol_ertek + kepifog_ertek + mozaik_ertek + matrix_ertek

g ~~ SES

osztalyzat ~ g + SES'

WISCMARSSAM <- '
g =~ betuszam_ertek + szamterj_ertek + szokincs_ertek + kozjel_ertek + altmegert_ertek + kodolas_ertek + szimbol_ertek + kepifog_ertek + mozaik_ertek + matrix_ertek

MARS =~ mszmt1 + mszmt2 + mszmt3 + mszmt4 + mszmt5 + mszmt6 + mszmt7 + mszmt8 + mszmt9 + mszmt10 + mszmt11 + mszmt12 + mszmt13 + mszmt14 + mszmt15 + mszmt16 + mszmt17 + mszmt18 + mszmt19 + mszmt20 + mszmt21 + mszmt22 + mszmt23 + mszmt24 + mszmt25 + mszmt26 + mszmt27 + mszmt28 + mszmt29 + mszmt30 + mszmt31 + mszmt32 + mszmt33 + mszmt34 + mszmt35 + mszmt36 + mszmt37 + mszmt38 + mszmt39 + mszmt40

g ~~ MARS

osztalyzat ~ g + MARS'

WISCDPVSESMARSSAM <- '
g =~ betuszam_ertek + szamterj_ertek + szokincs_ertek + kozjel_ertek + altmegert_ertek + kodolas_ertek + szimbol_ertek + kepifog_ertek + mozaik_ertek + matrix_ertek

dpv =~ dpv_tajek + dpv_szaml + dpv_szameml + dpv_szamegyz + dpv_rel + dpv_helyiert + dpv_potlas + dpv_muvel + dpv_szorzoszt + dpv_irasbeli + dpv_szoveg + dpv_logszab + dpv_muvsorr + dpv_praktik

MARS =~ mszmt1 + mszmt2 + mszmt3 + mszmt4 + mszmt5 + mszmt6 + mszmt7 + mszmt8 + mszmt9 + mszmt10 + mszmt11 + mszmt12 + mszmt13 + mszmt14 + mszmt15 + mszmt16 + mszmt17 + mszmt18 + mszmt19 + mszmt20 + mszmt21 + mszmt22 + mszmt23 + mszmt24 + mszmt25 + mszmt26 + mszmt27 + mszmt28 + mszmt29 + mszmt30 + mszmt31 + mszmt32 + mszmt33 + mszmt34 + mszmt35 + mszmt36 + mszmt37 + mszmt38 + mszmt39 + mszmt40

g ~~ dpv + SES + MARS
dpv ~~ SES + MARS
SES ~~ MARS

osztalyzat ~ g + dpv + SES + MARS'

WISCSAMFit <- sam(WISCSAM, data)
DPVSAMFit <- sam(DPVSAM, data)
SESSAMFit <- sem(SESSAM, data)
MARSSAMFit <- sam(MARSSAM, data)
WISCDPVSAMFit <- sam(WISCDPVSAM, data)
WISCSESSAMFit <- sam(WISCSESSAM, data)
WISCMARSSAMFit <- sam(WISCMARSSAM, data)
WISCDPVSESMARSSAMFit <- sam(WISCDPVSESMARSSAM, data)

summary(WISCSAMFit, stand = T, fit = T)
summary(DPVSAMFit, stand = T, fit = T)
summary(SESSAMFit, stand = T, fit = T)
summary(MARSSAMFit, stand = T, fit = T)
summary(WISCDPVSAMFit, stand = T, fit = T)
summary(WISCSESSAMFit, stand = T, fit = T)
summary(WISCMARSSAMFit, stand = T, fit = T)
summary(WISCDPVSESMARSSAMFit, stand = T, fit = T)
```

# Discussion

The DPV measures a *g* that is highly similar to the *g* measured by the WISC (*r* = .78), which is very similar to the *r* = .72 result achieved by Zaboski, Kranzler & Gage (2018). This is likely attenuated because of psychometric sampling error, since the DPV is purely concerned with mathematics and the *g* factor from the WISC is a higher-order *g* factor that is much less colored by a particular form of content.

# References

Zaboski, B. A., Kranzler, J. H., & Gage, N. A. (2018). Meta-analysis of the relationship between academic achievement and broad abilities of the Cattell-horn-Carroll theory. *Journal of School Psychology*, 71, 42--56. <https://doi.org/10.1016/j.jsp.2018.10.001>

```{r}
sessionInfo()
```
