Simulation of Major Analyses in ‘Minorities Are Disproportionately Underrepresented in Special Education: Longitudinal Evidence Across Five Disability Conditions’ by Morgan, Farkas, Hillemeier, Maczuga, & Cook (2015) in Educational Researcher
Abstract
This report simulates the major analyses conducted in the study “Minorities Are Disproportionately Underrepresented in Special Education: Longitudinal Evidence Across Five Disability Conditions”. The study concluded that “Minority children were consistently less likely than otherwise similar White, English-speaking children to be identified as disabled and so to receive special education services. From kindergarten entry to the end of middle school, racial- and ethnic-minority children were less likely to be identified as having (a) learning disabilities, (b) speech or language impairments, (c) intellectual disabilities, (d) health impairments, or (e) emotional disturbances.” (Morgan et al., 2015). This conclusion was made after controlling for many demographic and socioeconomic status variables as compared to uncontrolled prevalence rates. Hazard/survival modeling of K-8 longitudinal data through discrete-time logistic regression is conducted on simulated data of the relevant variables. The simulated data structure, R code, and summary outputs are provided for replication of the analyses and graphic and are presented in R Markdown/HTML formating.
Data Simulation
The following variables were simulated for the data structure that will be used to conduct the discrete-time logistic regression. Values do not correspond to those from the original study as neither a correlation matrix nor multivariate distributions were provided to simulate from and should be viewed as nonrepresentative. Many variables were sampled from the univariate distributions that were provided in the study. Standardized scores were randomly sampled from a normal distribution (0,1). Diagnosis status was uniformly sampled using the 5th root of the prevalence rates (by 8th grade) provided, to produce a similar overall rate in the simulated data. The diagnosis data was censored by student so that once a student was diagnosed, he/she would be removed from analyses for later grades.
Time invarying variables
- Ethnicity | Factors: White, Black, Hispanic, Other
- Gender | Factors: Female, Male
- Child’s age in month’s (fall Kindergarten): numeric
- Mother’s marital status (fall K) | Factor: Married, Unmarried
- Whether child was born at low birth weight | Factors: Y/N
- Mother’s age at child’s birth | Factors: Younger than 18, between 18 and 38, older than 38
- SES Quintile (fall K)| Factors: 1st, 2nd, 3rd, 4th, 5th
- Whether child has health insurance (fall K) | Factors: Y/N
- Parent interviewed in non-English | Factors: Y/N
Time varying variables measured at grades K, 1, 3, 5, 8
- Z score of averaged reading and math achievement test scores | Numeric
- Z score of Externalizing problem behaviors score | Numeric
- Z score of Self-regulatory behaviors score | Numeric
- Learning Disability Diagnosis Status | Binary numeric 0/1
- Speech or Language Impairment Diagnosis Status | Binary numeric 0/1
- Intellectual Disability Diagnosis Status | Binary numeric 0/1
- Health Impairment Disability Diagnosis Status | Binary numeric 0/1
- Emotional Disturbance Diagnosis Status | Binary numeric 0/1
set.seed(1)
data <- data.frame(matrix(ncol=0,nrow=20100))
data$id <- c(1:20100)
data$ethnicity <- as.factor(sample(c("White","Black", "Hispanic", "Other"),20100, replace=TRUE, prob = c(.596,.147,.187,.07)))
data$gender <- as.factor(sample(c("Male","Female"),20100, replace=TRUE, prob = c(.512,.488)))
data$agemonthsK <- sample(59:75,20100, replace = TRUE)
data$mothermarried <- as.factor(sample(c("Not Married","Married"),20100, replace=TRUE, prob = c(.307,.693)))
data$lowbirthweight <- as.factor(sample(c("Low Birth Weight","Not Low Birth Weight"),20100, replace=TRUE, prob = c(.04,.96)))
data$motheragebirth <- as.factor(sample(c("< 18","18-38", "> 38"),20100, replace=TRUE, prob = c(.056,.902,.042)))
data$SESquintile <- as.factor(sample(c("1st","2nd", "3rd", "4th", "5th"),20100, replace=TRUE, prob = c(.2,.2,.2,.2,.2)))
data$healthinsurance <- as.factor(sample(c("Has health insurance","Does not have health insurance"),20100, replace=TRUE, prob = c(.908,.092)))
data$interviewlanguage <- as.factor(sample(c("English","Not English"),20100, replace=TRUE, prob = c(.926,.073)))
data$testK <- rnorm(20100, 0,1); data$test1 <- rnorm(20100,0,1); data$test3 <- rnorm(20100,0,1);data$test5 <- rnorm(20100,0,1);data$test8 <- rnorm(20100,0,1)
data$extK <- rnorm(20100,0,1);data$ext1 <- rnorm(20100,0,1);data$ext3 <- rnorm(20100,0,1);data$ext5 <- rnorm(20100,0,1);data$ext8 <- rnorm(20100,0,1)
data$srK<-rnorm(20100,0,1);data$sr1<-rnorm(20100,0,1);data$sr3<-rnorm(20100,0,1);data$sr5<-rnorm(20100,0,1);data$sr8<-rnorm(20100,0,1);
data$LDdiagK <- sample(c(0,1),20100, replace=TRUE, prob = c(.987,.013))
data$LDdiag1 <- ifelse(data$LDdiagK==0,rbinom(20100,1,.013),NA)
data$LDdiag3 <- ifelse(data$LDdiag1==0,rbinom(20100,1,.013),NA)
data$LDdiag5 <- ifelse(data$LDdiag3==0,rbinom(20100,1,.013),NA)
data$LDdiag8 <- ifelse(data$LDdiag5==0,rbinom(20100,1,.013),NA)
data$SLdiagK <- sample(c(0,1),20100, replace=TRUE, prob = c(.987,.013))
data$SLdiag1 <- ifelse(data$SLdiagK==0,rbinom(20100,1,.013),NA)
data$SLdiag3 <- ifelse(data$SLdiag1==0,rbinom(20100,1,.013),NA)
data$SLdiag5 <- ifelse(data$SLdiag3==0,rbinom(20100,1,.013),NA)
data$SLdiag8 <-ifelse(data$SLdiag5==0,rbinom(20100,1,.013),NA)
data$IDdiagK <- sample(c(0,1),20100, replace=TRUE, prob = c(.9986,.0014))
data$IDdiag1 <- ifelse(data$IDdiagK==0,rbinom(20100,1,.0014),NA)
data$IDdiag3 <- ifelse(data$IDdiag1==0,rbinom(20100,1,.0014),NA)
data$IDdiag5 <- ifelse(data$IDdiag3==0,rbinom(20100,1,.0014),NA)
data$IDdiag8 <- ifelse(data$IDdiag5==0,rbinom(20100,1,.0014),NA)
data$HIdiagK <- sample(c(0,1),20100, replace=TRUE, prob = c(.998,.002))
data$HIdiag1 <- ifelse(data$HIdiagK==0,rbinom(20100,1,.002),NA)
data$HIdiag3 <- ifelse(data$HIdiag1==0,rbinom(20100,1,.002),NA)
data$HIdiag5 <- ifelse(data$HIdiag3==0,rbinom(20100,1,.002),NA)
data$HIdiag8 <- ifelse(data$HIdiag5==0,rbinom(20100,1,.002),NA)
data$EDdiagK <- sample(c(0,1),20100, replace=TRUE, prob = c(.9986,.0014))
data$EDdiag1 <- ifelse(data$EDdiagK==0,rbinom(20100,1,.0014),NA)
data$EDdiag3 <- ifelse(data$EDdiag1==0,rbinom(20100,1,.0014),NA)
data$EDdiag5 <- ifelse(data$EDdiag3==0,rbinom(20100,1,.0014),NA)
data$EDdiag8 <- ifelse(data$EDdiag5==0,rbinom(20100,1,.0014),NA)
head(data,5)## id ethnicity gender agemonthsK mothermarried lowbirthweight
## 1 1 White Female 62 Married Not Low Birth Weight
## 2 2 White Female 59 Married Not Low Birth Weight
## 3 3 White Female 72 Married Not Low Birth Weight
## 4 4 Black Female 65 Married Not Low Birth Weight
## 5 5 White Male 62 Married Not Low Birth Weight
## motheragebirth SESquintile healthinsurance interviewlanguage
## 1 18-38 5th Has health insurance English
## 2 18-38 1st Has health insurance English
## 3 18-38 1st Has health insurance English
## 4 18-38 4th Has health insurance English
## 5 18-38 1st Has health insurance English
## testK test1 test3 test5 test8 extK
## 1 0.99156556 1.0795464 -1.3397263 1.02541344 -1.3214049 0.1338242
## 2 -0.20545273 -0.7436489 -1.0153784 -0.03700107 1.0386168 -0.3815717
## 3 -0.02247959 1.0709135 1.3260162 -0.40732599 0.9601089 0.6238057
## 4 0.04536568 0.3966972 2.5649551 -1.06535994 -1.5744353 1.2061409
## 5 2.32504783 0.4447884 0.8016381 -0.96594534 -1.0870676 1.2976219
## ext1 ext3 ext5 ext8 srK sr1
## 1 -0.2996678 -1.9680970 -0.5198568 -0.86789015 1.1629558 -0.4079413
## 2 -1.0618775 1.2805163 -0.5851074 -1.30413781 -1.8611894 0.4529624
## 3 -2.0783756 -1.1739912 -0.7845906 0.01649076 -0.4802180 0.2405847
## 4 1.3450941 0.3341574 1.6805346 0.65372460 0.6639543 2.4464585
## 5 -1.3837864 -0.7691503 0.9835491 0.09686152 -0.1276034 0.1839210
## sr3 sr5 sr8 LDdiagK LDdiag1 LDdiag3 LDdiag5
## 1 -0.5230175 -0.66740193 -1.3021425 0 0 0 0
## 2 0.9814047 -0.52114598 -1.6268343 0 1 NA NA
## 3 -1.6679944 -0.57271408 -0.2105019 0 0 0 0
## 4 -0.8112780 0.22076212 -1.2651730 0 0 0 0
## 5 0.6000380 -0.05425816 0.2735865 0 0 0 0
## LDdiag8 SLdiagK SLdiag1 SLdiag3 SLdiag5 SLdiag8 IDdiagK IDdiag1 IDdiag3
## 1 0 0 0 0 0 0 0 0 0
## 2 NA 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0
## IDdiag5 IDdiag8 HIdiagK HIdiag1 HIdiag3 HIdiag5 HIdiag8 EDdiagK EDdiag1
## 1 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0
## EDdiag3 EDdiag5 EDdiag8
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0Data Restructuring
The data structure was reorganized from a wide format, where each student would have one row with all variable values, to a long format, where each student would have 5 rows corresponding to grades K, 1, 3, 5, 8 with the time varying variable values shown accordingly. A new variable called Grade was made to indicate what grade each row corresponded to. The ordering of the ethnicity factors was also rearranged so that later analyses would use White as the reference group to compare other ethnicities’ diagnosis rates to.
library(tidyr)
var <- c('test','ext','sr','LDdiag', 'SLdiag', 'IDdiag', 'HIdiag', 'EDdiag')
data2 <- reshape(data,dir='long', varying=lapply(var,grep,names(data)), v.names=var)
data2 <- data2[order(data2$id, data2$time),]
print(levels(data2$ethnicity))## [1] "Black" "Hispanic" "Other" "White"data2$ethnicity <- factor(data$ethnicity, levels(data$ethnicity)[c(4,1,2,3)])
print(levels(data2$ethnicity))## [1] "White" "Black" "Hispanic" "Other"data2$time <- as.factor(data2$time)
levels(data2$time) <- c('K','1st','3rd','5th','8th')
names(data2)[names(data2)=='time'] <- 'Grade'
head(data2,10)## id ethnicity gender agemonthsK mothermarried lowbirthweight
## 1.1 1 White Female 62 Married Not Low Birth Weight
## 1.2 1 White Female 62 Married Not Low Birth Weight
## 1.3 1 White Female 62 Married Not Low Birth Weight
## 1.4 1 Black Female 62 Married Not Low Birth Weight
## 1.5 1 White Female 62 Married Not Low Birth Weight
## 2.1 2 Black Female 59 Married Not Low Birth Weight
## 2.2 2 Other Female 59 Married Not Low Birth Weight
## 2.3 2 Hispanic Female 59 Married Not Low Birth Weight
## 2.4 2 Hispanic Female 59 Married Not Low Birth Weight
## 2.5 2 White Female 59 Married Not Low Birth Weight
## motheragebirth SESquintile healthinsurance interviewlanguage
## 1.1 18-38 5th Has health insurance English
## 1.2 18-38 5th Has health insurance English
## 1.3 18-38 5th Has health insurance English
## 1.4 18-38 5th Has health insurance English
## 1.5 18-38 5th Has health insurance English
## 2.1 18-38 1st Has health insurance English
## 2.2 18-38 1st Has health insurance English
## 2.3 18-38 1st Has health insurance English
## 2.4 18-38 1st Has health insurance English
## 2.5 18-38 1st Has health insurance English
## Grade test ext sr LDdiag SLdiag IDdiag HIdiag
## 1.1 K 0.99156556 0.1338242 1.1629558 0 0 0 0
## 1.2 1st 1.07954639 -0.2996678 -0.4079413 0 0 0 0
## 1.3 3rd -1.33972633 -1.9680970 -0.5230175 0 0 0 0
## 1.4 5th 1.02541344 -0.5198568 -0.6674019 0 0 0 0
## 1.5 8th -1.32140490 -0.8678902 -1.3021425 0 0 0 0
## 2.1 K -0.20545273 -0.3815717 -1.8611894 0 0 0 0
## 2.2 1st -0.74364893 -1.0618775 0.4529624 1 0 0 0
## 2.3 3rd -1.01537838 1.2805163 0.9814047 NA 0 0 0
## 2.4 5th -0.03700107 -0.5851074 -0.5211460 NA 0 0 0
## 2.5 8th 1.03861682 -1.3041378 -1.6268343 NA 0 0 0
## EDdiag
## 1.1 0
## 1.2 0
## 1.3 0
## 1.4 0
## 1.5 0
## 2.1 0
## 2.2 0
## 2.3 0
## 2.4 0
## 2.5 0Time Discrete Longitudinal Logistic Regression
Models
The following code creates generalized linear models using a logistic fit to predict diagnosis longitudinally across grades K, 1, 3, 5, and 8 using all other predictors. Three models are created for each disability category; one using just ethnicity and grade levels as predictors, the second using all predictors, and the third including an ethnicity x grade interaction. The first two types of models allow one to comapre how diagnosis rates may compare between minority groups and whites before and after all the other variables are controlled for. The third type of models is for later hazard and survival functions that will be plotted. The original study did not include an ethnicity x grade interaction in their table of model coefficients. The R output summary for the modeling of learning disability rates with all predictors is provided as well as a table comparing the first two types of models made for a total of 10 models. R summary coefficients are in logit units and table coefficients are provided as odds ratios.
LDjusteth <- glm(LDdiag~Grade + ethnicity, family="binomial", data=data2)
glmLD <- glm(LDdiag~Grade + ethnicity+gender+agemonthsK+mothermarried+
lowbirthweight+motheragebirth+SESquintile+healthinsurance+interviewlanguage+test+ext+sr, family="binomial", data=data2)
SLjusteth <- glm(LDdiag~Grade + ethnicity, family="binomial", data=data2)
glmSL <- glm(SLdiag~Grade + ethnicity+gender+agemonthsK+mothermarried+
lowbirthweight+motheragebirth+SESquintile+healthinsurance+interviewlanguage+test+ext+sr, family="binomial", data=data2)
IDjusteth <- glm(IDdiag~Grade + ethnicity, family="binomial", data=data2)
glmID <- glm(IDdiag~Grade + ethnicity+gender+agemonthsK+mothermarried+
lowbirthweight+motheragebirth+SESquintile+healthinsurance+interviewlanguage+test+ext+sr, family="binomial", data=data2)
HIjusteth <- glm(HIdiag~Grade + ethnicity, family="binomial", data=data2)
glmHI <- glm(HIdiag~Grade + ethnicity+gender+agemonthsK+mothermarried+
lowbirthweight+motheragebirth+SESquintile+healthinsurance+interviewlanguage+test+ext+sr, family="binomial", data=data2)
EDjusteth <- glm(EDdiag~Grade + ethnicity, family="binomial", data=data2)
glmED <- glm(EDdiag~Grade + ethnicity+gender+agemonthsK+mothermarried+
lowbirthweight+motheragebirth+SESquintile+healthinsurance+interviewlanguage+test+ext+sr, family="binomial", data=data2)
glmLD1 <- glm(LDdiag~Grade -1 + ethnicity + Grade:ethnicity +gender+agemonthsK+mothermarried+
lowbirthweight+motheragebirth+SESquintile+healthinsurance+interviewlanguage+test+ext+sr, family="binomial", data=data2)
glmSL1 <- glm(SLdiag~Grade -1 + ethnicity + Grade:ethnicity +gender+agemonthsK+mothermarried+
lowbirthweight+motheragebirth+SESquintile+healthinsurance+interviewlanguage+test+ext+sr, family="binomial", data=data2)
glmID1 <- glm(IDdiag~Grade -1 + ethnicity + Grade:ethnicity +gender+agemonthsK+mothermarried+
lowbirthweight+motheragebirth+SESquintile+healthinsurance+interviewlanguage+test+ext+sr, family="binomial", data=data2)
glmHI1 <- glm(HIdiag~Grade -1 + ethnicity + Grade:ethnicity +gender+agemonthsK+mothermarried+
lowbirthweight+motheragebirth+SESquintile+healthinsurance+interviewlanguage+test+ext+sr, family="binomial", data=data2)
glmED1 <- glm(EDdiag~Grade -1 + ethnicity + Grade:ethnicity +gender+agemonthsK+mothermarried+
lowbirthweight+motheragebirth+SESquintile+healthinsurance+interviewlanguage+test+ext+sr, family="binomial", data=data2)
summary(glmLD1)##
## Call:
## glm(formula = LDdiag ~ Grade - 1 + ethnicity + Grade:ethnicity +
## gender + agemonthsK + mothermarried + lowbirthweight + motheragebirth +
## SESquintile + healthinsurance + interviewlanguage + test +
## ext + sr, family = "binomial", data = data2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.2277 -0.1667 -0.1603 -0.1539 3.0651
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## GradeK -3.747899 0.441114 -8.496 <2e-16
## Grade1st -3.717740 0.441229 -8.426 <2e-16
## Grade3rd -3.643333 0.440740 -8.266 <2e-16
## Grade5th -3.868527 0.442311 -8.746 <2e-16
## Grade8th -3.729017 0.441508 -8.446 <2e-16
## ethnicityBlack -0.040506 0.186090 -0.218 0.8277
## ethnicityHispanic -0.097257 0.173300 -0.561 0.5747
## ethnicityOther 0.074447 0.244705 0.304 0.7610
## genderMale 0.013754 0.056577 0.243 0.8079
## agemonthsK -0.003137 0.005740 -0.547 0.5847
## mothermarriedNot Married 0.011366 0.061149 0.186 0.8525
## lowbirthweightNot Low Birth Weight -0.106839 0.137752 -0.776 0.4380
## motheragebirth> 38 -0.265137 0.186441 -1.422 0.1550
## motheragebirth18-38 -0.184103 0.118702 -1.551 0.1209
## SESquintile2nd 0.053461 0.090928 0.588 0.5566
## SESquintile3rd 0.046586 0.090942 0.512 0.6085
## SESquintile4th 0.046956 0.090827 0.517 0.6052
## SESquintile5th 0.085615 0.089692 0.955 0.3398
## healthinsuranceHas health insurance -0.181800 0.090037 -2.019 0.0435
## interviewlanguageNot English 0.015624 0.109015 0.143 0.8860
## test -0.004265 0.028204 -0.151 0.8798
## ext -0.001832 0.028351 -0.065 0.9485
## sr 0.031964 0.028295 1.130 0.2586
## Grade1st:ethnicityBlack -0.078014 0.266252 -0.293 0.7695
## Grade3rd:ethnicityBlack -0.137921 0.261034 -0.528 0.5972
## Grade5th:ethnicityBlack -0.076913 0.277634 -0.277 0.7818
## Grade8th:ethnicityBlack 0.100988 0.261006 0.387 0.6988
## Grade1st:ethnicityHispanic -0.005687 0.244016 -0.023 0.9814
## Grade3rd:ethnicityHispanic 0.054643 0.239387 0.228 0.8194
## Grade5th:ethnicityHispanic 0.303724 0.241545 1.257 0.2086
## Grade8th:ethnicityHispanic 0.051219 0.243605 0.210 0.8335
## Grade1st:ethnicityOther 0.026846 0.335438 0.080 0.9362
## Grade3rd:ethnicityOther -0.021691 0.338221 -0.064 0.9489
## Grade5th:ethnicityOther 0.366469 0.334297 1.096 0.2730
## Grade8th:ethnicityOther 0.128059 0.335984 0.381 0.7031
##
## GradeK ***
## Grade1st ***
## Grade3rd ***
## Grade5th ***
## Grade8th ***
## ethnicityBlack
## ethnicityHispanic
## ethnicityOther
## genderMale
## agemonthsK
## mothermarriedNot Married
## lowbirthweightNot Low Birth Weight
## motheragebirth> 38
## motheragebirth18-38
## SESquintile2nd
## SESquintile3rd
## SESquintile4th
## SESquintile5th
## healthinsuranceHas health insurance *
## interviewlanguageNot English
## test
## ext
## sr
## Grade1st:ethnicityBlack
## Grade3rd:ethnicityBlack
## Grade5th:ethnicityBlack
## Grade8th:ethnicityBlack
## Grade1st:ethnicityHispanic
## Grade3rd:ethnicityHispanic
## Grade5th:ethnicityHispanic
## Grade8th:ethnicityHispanic
## Grade1st:ethnicityOther
## Grade3rd:ethnicityOther
## Grade5th:ethnicityOther
## Grade8th:ethnicityOther
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 135774 on 97940 degrees of freedom
## Residual deviance: 13531 on 97905 degrees of freedom
## (2560 observations deleted due to missingness)
## AIC: 13601
##
## Number of Fisher Scoring iterations: 7library(stargazer,quietly=TRUE)
stargazer(LDjusteth,glmLD,SLjusteth,glmSL,IDjusteth,glmID,HIjusteth,glmHI,EDjusteth,glmED, type="html", column.labels=c("Learning Disability", "Speech/Language Impairment", "Intellectual Disability","Health Impairment", "Emotional Disturbances"),column.separate=c(2,2,2,2,2),dep.var.caption = "",apply.coef=exp)| LDdiag | SLdiag | IDdiag | HIdiag | EDdiag | ||||||
| Learning Disability | Speech/Language Impairment | Intellectual Disability | Health Impairment | Emotional Disturbances | ||||||
| (1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | |
| Grade1st | 1.020*** | 1.020*** | 1.020*** | 1.044*** | 0.657** | 0.657** | 1.030*** | 1.034*** | 1.073*** | 1.071*** |
| (0.089) | (0.089) | (0.089) | (0.090) | (0.295) | (0.295) | (0.228) | (0.228) | (0.263) | (0.263) | |
| Grade3rd | 1.096*** | 1.097*** | 1.096*** | 1.041*** | 1.279*** | 1.282*** | 1.135*** | 1.140*** | 1.007*** | 1.004*** |
| (0.088) | (0.088) | (0.088) | (0.090) | (0.248) | (0.248) | (0.223) | (0.223) | (0.267) | (0.268) | |
| Grade5th | 0.958*** | 0.959*** | 0.958*** | 0.835*** | 1.073*** | 1.076*** | 1.244*** | 1.246*** | 1.041*** | 1.040*** |
| (0.091) | (0.091) | (0.091) | (0.096) | (0.259) | (0.259) | (0.218) | (0.218) | (0.265) | (0.265) | |
| Grade8th | 1.053*** | 1.054*** | 1.053*** | 1.056*** | 0.729** | 0.731** | 1.143*** | 1.149*** | 0.933*** | 0.930*** |
| (0.089) | (0.089) | (0.089) | (0.091) | (0.287) | (0.287) | (0.223) | (0.223) | (0.273) | (0.273) | |
| ethnicityBlack | 0.925*** | 0.924*** | 0.925*** | 0.932*** | 0.989*** | 0.992*** | 1.054*** | 1.054*** | 0.742*** | 0.740*** |
| (0.085) | (0.085) | (0.085) | (0.086) | (0.244) | (0.244) | (0.193) | (0.193) | (0.272) | (0.272) | |
| ethnicityHispanic | 0.984*** | 0.982*** | 0.984*** | 0.932*** | 0.846*** | 0.848*** | 0.825*** | 0.823*** | 1.048*** | 1.053*** |
| (0.076) | (0.076) | (0.076) | (0.079) | (0.239) | (0.239) | (0.195) | (0.195) | (0.217) | (0.217) | |
| ethnicityOther | 1.186*** | 1.187*** | 1.186*** | 0.952*** | 0.892** | 0.892** | 0.838*** | 0.834*** | 0.962*** | 0.964*** |
| (0.104) | (0.104) | (0.104) | (0.117) | (0.351) | (0.351) | (0.291) | (0.291) | (0.334) | (0.334) | |
| genderMale | 1.014*** | 1.020*** | 1.114*** | 0.899*** | 0.791*** | |||||
| (0.057) | (0.058) | (0.172) | (0.138) | (0.169) | ||||||
| agemonthsK | 0.997*** | 0.995*** | 1.004*** | 0.993*** | 1.020*** | |||||
| (0.006) | (0.006) | (0.017) | (0.014) | (0.017) | ||||||
| mothermarriedNot Married | 1.012*** | 1.031*** | 0.772*** | 1.035*** | 0.869*** | |||||
| (0.061) | (0.063) | (0.196) | (0.149) | (0.188) | ||||||
| lowbirthweightNot Low Birth Weight | 0.898*** | 0.893*** | 1.402*** | 2.161*** | 1.932*** | |||||
| (0.138) | (0.142) | (0.508) | (0.505) | (0.584) | ||||||
| motheragebirth> 38 | 0.767*** | 0.794*** | 0.427 | 1.581*** | 3.377*** | |||||
| (0.186) | (0.201) | (0.817) | (0.450) | (0.678) | ||||||
| motheragebirth18-38 | 0.832*** | 0.941*** | 1.229*** | 1.209*** | 2.456*** | |||||
| (0.119) | (0.129) | (0.418) | (0.342) | (0.584) | ||||||
| SESquintile2nd | 1.054*** | 1.023*** | 0.948*** | 0.745*** | 0.837*** | |||||
| (0.091) | (0.095) | (0.261) | (0.218) | (0.257) | ||||||
| SESquintile3rd | 1.047*** | 1.129*** | 0.938*** | 0.761*** | 0.834*** | |||||
| (0.091) | (0.093) | (0.261) | (0.216) | (0.257) | ||||||
| SESquintile4th | 1.047*** | 0.997*** | 0.873*** | 0.895*** | 0.829*** | |||||
| (0.091) | (0.095) | (0.266) | (0.207) | (0.257) | ||||||
| SESquintile5th | 1.089*** | 1.138*** | 0.698** | 0.798*** | 0.693*** | |||||
| (0.090) | (0.092) | (0.281) | (0.212) | (0.268) | ||||||
| healthinsuranceHas health insurance | 0.834*** | 0.932*** | 1.479*** | 1.111*** | 1.135*** | |||||
| (0.090) | (0.097) | (0.345) | (0.247) | (0.302) | ||||||
| interviewlanguageNot English | 1.017*** | 0.773*** | 0.707* | 0.651** | 0.895*** | |||||
| (0.109) | (0.126) | (0.388) | (0.324) | (0.345) | ||||||
| test | 0.996*** | 0.956*** | 1.160*** | 1.093*** | 0.854*** | |||||
| (0.028) | (0.029) | (0.085) | (0.069) | (0.084) | ||||||
| ext | 0.998*** | 1.013*** | 1.007*** | 1.157*** | 0.882*** | |||||
| (0.028) | (0.029) | (0.086) | (0.069) | (0.085) | ||||||
| sr | 1.032*** | 1.068*** | 1.112*** | 1.050*** | 1.038*** | |||||
| (0.028) | (0.029) | (0.086) | (0.069) | (0.084) | ||||||
| Constant | 0.013 | 0.023 | 0.013 | 0.020 | 0.002 | 0.001 | 0.002 | 0.001 | 0.001 | 0.0001 |
| (0.067) | (0.439) | (0.067) | (0.455) | (0.197) | (1.393) | (0.171) | (1.154) | (0.200) | (1.464) | |
| Observations | 97,940 | 97,940 | 97,940 | 98,061 | 100,222 | 100,222 | 100,098 | 100,098 | 100,213 | 100,213 |
| Log Likelihood | -6,772.650 | -6,767.802 | -6,772.650 | -6,431.161 | -1,036.117 | -1,028.591 | -1,503.379 | -1,495.202 | -1,065.848 | -1,056.976 |
| Akaike Inf. Crit. | 13,561.300 | 13,581.600 | 13,561.300 | 12,908.320 | 2,088.233 | 2,103.182 | 3,022.758 | 3,036.403 | 2,147.695 | 2,159.952 |
| Note: | p<0.1; p<0.05; p<0.01 | |||||||||
Hazard and Survival Models Setup
The following code takes the models created earlier predicting diagnoses for each disability category from all other variables as well as an ethnicity x grade interaction, extracts the logit coefficient terms of the effect of grade on diagnosis by ethnicity after controlling for all other variables, and converts them into their corresponding hazard probabilities, the probability of a student within an ethnicity being diagnosed in a grade given that one was not previously diagnosed after controlling for all other variables. These probabilities are plotted across the grade levels by ethnicity to produce hazard function plots for each disability category. Survival plots are also created by plotting by grade, ethnicity, and disability category the proportion of students who are predicted to survive to the next grade level without being diagnosed based on their corresponding hazard functions controlling for all other variables.
library(effects)## Loading required package: carData## lattice theme set by effectsTheme()
## See ?effectsTheme for details.effLD <- Effect(focal.predictors=c("Grade", "ethnicity"), mod=glmLD1, multiline=TRUE, style="stacked")
effdfLD <- data.frame(effLD)
effSL <- Effect(focal.predictors=c("Grade", "ethnicity"), mod=glmSL1, multiline=TRUE, style="stacked")
effdfSL <- data.frame(effSL)
effID <- Effect(focal.predictors=c("Grade", "ethnicity"), mod=glmID1, multiline=TRUE, style="stacked")
effdfID <- data.frame(effID)
effHI <- Effect(focal.predictors=c("Grade", "ethnicity"), mod=glmHI1, multiline=TRUE, style="stacked")
effdfHI <- data.frame(effHI)
effED <- Effect(focal.predictors=c("Grade", "ethnicity"), mod=glmED1, multiline=TRUE, style="stacked")
effdfED <- data.frame(effED)
effdfLD$survival <- c(1:20)
for(i in 1:20){
if(i%%5==1){effdfLD$survival[i]=1-effdfLD$fit[i]}
else{effdfLD$survival[i]=effdfLD$survival[i-1]-(effdfLD$survival[i-1]*effdfLD$fit[i])}
}
effdfSL$survival <- c(1:20)
for(i in 1:20){
if(i%%5==1){effdfSL$survival[i]=1-effdfSL$fit[i]}
else{effdfSL$survival[i]=effdfSL$survival[i-1]-(effdfSL$survival[i-1]*effdfSL$fit[i])}
}
effdfID$survival <- c(1:20)
for(i in 1:20){
if(i%%5==1){effdfID$survival[i]=1-effdfID$fit[i]}
else{effdfID$survival[i]=effdfID$survival[i-1]-(effdfID$survival[i-1]*effdfID$fit[i])}
}
effdfHI$survival <- c(1:20)
for(i in 1:20){
if(i%%5==1){effdfHI$survival[i]=1-effdfHI$fit[i]}
else{effdfHI$survival[i]=effdfHI$survival[i-1]-(effdfHI$survival[i-1]*effdfHI$fit[i])}
}
effdfED$survival <- c(1:20)
for(i in 1:20){
if(i%%5==1){effdfED$survival[i]=1-effdfED$fit[i]}
else{effdfED$survival[i]=effdfED$survival[i-1]-(effdfED$survival[i-1]*effdfED$fit[i])}
}Hazard Function Plots
library(ggplot2)
LDplotH <- ggplot(effdfLD, aes(x=Grade,y=fit)) + geom_point(aes(col=ethnicity)) + scale_x_discrete(limits=c("K","1st","3rd","5th","8th")) +
geom_line(aes(group=ethnicity, color=ethnicity)) + labs(x="Grade", y="Hazard Function", title="Learning Disability Diagnosis")
SLplotH <- ggplot(effdfSL, aes(x=Grade,y=fit)) + geom_point(aes(col=ethnicity)) + scale_x_discrete(limits=c("K","1st","3rd","5th","8th")) +
geom_line(aes(group=ethnicity, color=ethnicity)) + labs(x="Grade", y="Hazard Function", title="Speech or Language Impairment Diagnosis")
IDplotH <- ggplot(effdfID, aes(x=Grade,y=fit)) + geom_point(aes(col=ethnicity)) + scale_x_discrete(limits=c("K","1st","3rd","5th","8th")) +
geom_line(aes(group=ethnicity, color=ethnicity)) + labs(x="Grade", y="Hazard Function", title="Intellectual Disability Diagnosis")
HIplotH <- ggplot(effdfHI, aes(x=Grade,y=fit)) + geom_point(aes(col=ethnicity)) + scale_x_discrete(limits=c("K","1st","3rd","5th","8th")) +
geom_line(aes(group=ethnicity, color=ethnicity)) + labs(x="Grade", y="Hazard Function", title="Health Impairment Diagnosis")
EDplotH <- ggplot(effdfED, aes(x=Grade,y=fit)) + geom_point(aes(col=ethnicity)) + scale_x_discrete(limits=c("K","1st","3rd","5th","8th")) +
geom_line(aes(group=ethnicity, color=ethnicity)) + labs(x="Grade", y="Hazard Function", title="Emotional Disturbance Diagnosis")
plot(LDplotH);plot(SLplotH);plot(IDplotH); plot(HIplotH);plot(EDplotH)
Survival Function Plots
LDplotS <- ggplot(effdfLD, aes(x=Grade,y=survival)) + geom_point(aes(col=ethnicity)) + scale_x_discrete(limits=c("K","1st","3rd","5th","8th")) +
geom_line(aes(group=ethnicity, color=ethnicity)) + labs(x="Grade", y="Survival Function", title="Learning Disability Diagnosis")
SLplotS <- ggplot(effdfSL, aes(x=Grade,y=survival)) + geom_point(aes(col=ethnicity)) + scale_x_discrete(limits=c("K","1st","3rd","5th","8th")) +
geom_line(aes(group=ethnicity, color=ethnicity)) + labs(x="Grade", y="Survival Function", title="Speech or Language Impairment Diagnosis")
IDplotS <- ggplot(effdfID, aes(x=Grade,y=survival)) + geom_point(aes(col=ethnicity)) + scale_x_discrete(limits=c("K","1st","3rd","5th","8th")) +
geom_line(aes(group=ethnicity, color=ethnicity)) + labs(x="Grade", y="Survival Function", title="Intellectual Disability Diagnosis")
HIplotS <- ggplot(effdfHI, aes(x=Grade,y=survival)) + geom_point(aes(col=ethnicity)) + scale_x_discrete(limits=c("K","1st","3rd","5th","8th")) +
geom_line(aes(group=ethnicity, color=ethnicity)) + labs(x="Grade", y="Survival Function", title="Health Impairment Diagnosis")
EDplotS <- ggplot(effdfED, aes(x=Grade,y=survival)) + geom_point(aes(col=ethnicity)) + scale_x_discrete(limits=c("K","1st","3rd","5th","8th")) +
geom_line(aes(group=ethnicity, color=ethnicity)) + labs(x="Grade", y="Survival Function", title="Emotional Disturbance Diagnosis")
plot(LDplotS); plot(SLplotS);plot(IDplotS);plot(HIplotS);plot(EDplotS)
Reference
Morgan, P.L., Farkas, G., Hillemeier, M.M., Mattison, R., Maczuga, S., Li, H., & Cook, M. (2015). Minorities Are Disproportionately Underrepresented in Special Education: Longitudinal Evidence Across Five Disability Conditions. Educational Researcher