This webpage contains outputs and R codes used for the two-level multilevel modeling in Simmons, Shin, and Hart (2024). Data analysis scripts have been posted through an online data repository, accessible at https://osf.io/gqpfk.

Data Analysis

Preprocess

Load libraries.

suppressPackageStartupMessages({
    library(tidyverse)
    library(dplyr)
    library(fastDummies)
    library(nlme)
    library(lmeInfo)
    library(clubSandwich)
    library(plotly)
    library(ggplot2)
})

load("output.RData")

Preprocess data.

data = read.csv("data.csv") 

data$NCES_district_locale <- gsub("-.*", "", data$NCES_district_locale)

data <- dummy_cols(data, select_columns = "NCES_district_locale", remove_selected_columns = FALSE)
colnames(data) <-  gsub("^NCES_district_locale_", "", colnames(data))

data <- data %>% 
    mutate(prevalence_rate = round((dyslexia_n / total_n) * 100, 2)) %>%
    mutate(year_2017 = year - 2017) %>%
    mutate(slope_2017 = year_2017) %>%
    mutate(year_2018 = year - 2018) %>%
    mutate(knot_2018 = as.integer(I(year >= 2018 ))) %>%
    mutate(slope_change_2018 = year_2018 * knot_2018) %>%
    mutate(year_2019 = year - 2019) %>%
    mutate(knot_2019 = as.integer(I(year >= 2019 ))) %>%
    mutate(slope_change_2019 = year_2019 * knot_2019) %>%
    mutate(year_2020 = year - 2020) %>%
    mutate(knot_2020 = as.integer(I(year >= 2020 ))) %>%
    mutate(slope_change_2020 = year_2020 * knot_2020) %>%
    mutate(year_2021 = year - 2021) %>%
    mutate(knot_2021 = as.integer(I(year >= 2021 ))) %>%
    mutate(slope_change_2021 = year_2021 * knot_2021) %>%
    mutate(year_2022 = year - 2022) %>%
    mutate(knot_2022 = as.integer(I(year >= 2022 ))) %>%
    mutate(slope_change_2022 = year_2022 * knot_2022)

Visualize

Null Model

Check residuals for normality

Model.0 <- lme(
    prevalence_rate ~ 1, 
    data = data,
    method = "REML",
    random = ~ 1 | district,
    control = list(
        maxIter = 100, 
        msMaxIter = 100, 
        tolerance = 1e-3, 
        opt = "optim",
        optimMethod = "BFGS"
    )
)

Cluster-robust variance estimation

options(width = 100)
Model.0.vcov <- vcovCR(Model.0, type = "CR2") # Calculate variance-covariance matrix
Model.0.crve <-
  coef_test(Model.0, vcov = Model.0.vcov , test = "Satterthwaite") # Test regression coefficients 
Model.0.crve

##        Coef. Estimate   SE t-stat d.f. (Satt) p-val (Satt) Sig.
##  (Intercept)      258 19.5   13.3         810       <0.001  ***

Variance components

VarCorr(Model.0, which = "var-cov")

## district = pdLogChol(1) 
##             Variance  StdDev  
## (Intercept) 302620.93 550.1099
## Residual     30090.81 173.4670

Standard errors of variances

sqrt(diag(varcomp_vcov(Model.0)))

## Tau.district.var((Intercept))                      sigma_sq 
##                    15251.2076                      610.0463

Intraclass correlation coefficient (ICC)

total_var <- as.numeric(VarCorr(Model.0)[[2]]) +
  as.numeric(VarCorr(Model.0)[[1]]) 

ICC.L2 = as.numeric(VarCorr(Model.0)[[1]]) / total_var
ICC.L2

## [1] 0.9095589

Model 1

Legislative trends in the prevalence rate of students identified as having characteristics of dyslexia in Texas school districts.

Model.1 <- lme(
  prevalence_rate ~ slope_2017 + slope_change_2018 + slope_change_2019 + 
    slope_change_2020 + slope_change_2021 + slope_change_2022, 
  data = data,
  method = "REML",
  random = ~ 1 | district,
  correlation = corAR1(value = 0.2, form = ~ year | district),
  control = list(
    maxIter = 100, 
    msMaxIter = 100, 
    tolerance = 1e-3, 
    opt = "optim",
    optimMethod = "BFGS"
  )
)

Cluster-robust variance estimation

options(width = 100)
Model.1.vcov <- vcovCR(Model.1, type = "CR2") # Calculate variance-covariance matrix
Model.1.crve <-
  coef_test(Model.1, vcov = Model.1.vcov , test = "Satterthwaite") # Test regression coefficients 
Model.1.crve

##              Coef. Estimate     SE t-stat d.f. (Satt) p-val (Satt) Sig.
##        (Intercept)   4.0977 0.0826 49.580         810       <0.001  ***
##         slope_2017   0.2922 0.0291 10.028         810       <0.001  ***
##  slope_change_2018   0.2211 0.0481  4.601         810       <0.001  ***
##  slope_change_2019   0.1055 0.0477  2.213         810       0.0272    *
##  slope_change_2020  -0.1063 0.0497 -2.139         810       0.0327    *
##  slope_change_2021  -0.0105 0.0503 -0.209         810       0.8346     
##  slope_change_2022  -0.0130 0.0514 -0.254         810       0.7996

Variance components

VarCorr(Model.1, which = "var-cov")

## district = pdLogChol(1) 
##             Variance  StdDev   
## (Intercept) 0.1091445 0.3303703
## Residual    6.8870020 2.6243098

Standard errors of variances

sqrt(diag(varcomp_vcov(Model.1)))

## Tau.district.var((Intercept))                    cor_params                      sigma_sq 
##                    1.77682797                    0.01780919                    1.79907633

Model 2

Prevalence rate change of students identified as having characteristics of dyslexia over seven school years moderated by school district locale (i.e., city, suburban, town, or rural)

Model.2 <- lme(
    prevalence_rate ~ slope_2017 + slope_change_2018 + slope_change_2019 + 
        slope_change_2020 + slope_change_2021 + slope_change_2022 +
        (slope_2017 + slope_change_2018 + slope_change_2019 + 
             slope_change_2020 + slope_change_2021 + slope_change_2022)*Suburb + 
        (slope_2017 + slope_change_2018 + slope_change_2019 + 
             slope_change_2020 + slope_change_2021 + slope_change_2022)*Town +
        (slope_2017 + slope_change_2018 + slope_change_2019 + 
             slope_change_2020 + slope_change_2021 + slope_change_2022)*Rural, 
    data = data,
    method = "REML",
    random = ~ 1 | district,
    correlation = corAR1(value = 0.2, form = ~ year | district),
    control = list(
        maxIter = 100, 
        msMaxIter = 100, 
        tolerance = 1e-3, 
        opt = "optim",
        optimMethod = "BFGS"
    )
)

Cluster-robust variance estimation

options(width = 100)
Model.2.vcov <- vcovCR(Model.2, type = "CR2") # Calculate variance-covariance matrix
Model.2.crve <-
    coef_test(Model.2, vcov = Model.2.vcov , test = "Satterthwaite") # Test regression coefficients 
Model.2.crve

##                     Coef. Estimate     SE t-stat d.f. (Satt) p-val (Satt) Sig.
##               (Intercept)   2.8017 0.1741 16.096          94      < 0.001  ***
##                slope_2017   0.1672 0.0589  2.839          94      0.00555   **
##         slope_change_2018   0.2982 0.1250  2.386          94      0.01905    *
##         slope_change_2019   0.1205 0.1213  0.994          94      0.32289     
##         slope_change_2020  -0.2268 0.1214 -1.868          94      0.06482    .
##         slope_change_2021   0.0829 0.0961  0.863          94      0.39009     
##         slope_change_2022   0.0255 0.1173  0.217          94      0.82849     
##                    Suburb   0.3260 0.2463  1.323         201      0.18725     
##                      Town   1.0761 0.2360  4.560         188      < 0.001  ***
##                     Rural   1.9708 0.2117  9.311         141      < 0.001  ***
##         slope_2017:Suburb   0.1196 0.0803  1.489         201      0.13807     
##  slope_change_2018:Suburb  -0.1753 0.1568 -1.118         201      0.26506     
##  slope_change_2019:Suburb   0.0647 0.1386  0.467         201      0.64115     
##  slope_change_2020:Suburb   0.0611 0.1525  0.401         201      0.68912     
##  slope_change_2021:Suburb  -0.0609 0.1652 -0.369         201      0.71264     
##  slope_change_2022:Suburb   0.0662 0.1622  0.408         201      0.68378     
##           slope_2017:Town   0.0869 0.0772  1.126         188      0.26159     
##    slope_change_2018:Town  -0.1027 0.1466 -0.701         188      0.48439     
##    slope_change_2019:Town   0.0190 0.1449  0.131         188      0.89611     
##    slope_change_2020:Town   0.1310 0.1421  0.922         188      0.35774     
##    slope_change_2021:Town  -0.1255 0.1179 -1.065         188      0.28845     
##    slope_change_2022:Town   0.0308 0.1397  0.220         188      0.82593     
##          slope_2017:Rural   0.1733 0.0763  2.272         141      0.02459    *
##   slope_change_2018:Rural  -0.0555 0.1478 -0.375         141      0.70791     
##   slope_change_2019:Rural  -0.0567 0.1453 -0.390         141      0.69690     
##   slope_change_2020:Rural   0.1603 0.1475  1.086         141      0.27916     
##   slope_change_2021:Rural  -0.1093 0.1273 -0.858         141      0.39213     
##   slope_change_2022:Rural  -0.1091 0.1453 -0.751         141      0.45397

Variance components

VarCorr(Model.2, which = "var-cov")

## district = pdLogChol(1) 
##             Variance  StdDev   
## (Intercept) 0.1076066 0.3280345
## Residual    6.1009470 2.4700095

Standard errors of variances

sqrt(diag(varcomp_vcov(Model.2)))

## Tau.district.var((Intercept))                    cor_params                      sigma_sq 
##                    1.39041715                    0.01782691                    1.41182923

Dyslexia Identification: Texas Legislative Trends in Prevalence Rate of Students by School District Locale

2024-05-24

Data Analysis

Preprocess

Load libraries.

Preprocess data.

Visualize

Null Model

Check residuals for normality

Cluster-robust variance estimation

Variance components

Standard errors of variances

Intraclass correlation coefficient (ICC)

Model 1

Legislative trends in the prevalence rate of students identified as having characteristics of dyslexia in Texas school districts.

Cluster-robust variance estimation

Variance components

Standard errors of variances

Model 2

Prevalence rate change of students identified as having characteristics of dyslexia over seven school years moderated by school district locale (i.e., city, suburban, town, or rural)

Cluster-robust variance estimation

Variance components

Standard errors of variances