Features of One Gender Schools

Demographics

Starting with one of the categorical demographic variables, there are 8 columns of interest:

• HBCU (Flag for Historically Black College or University)
• PBI (Flag for predominantly black institution)
• ANNHI (Flag for Alaska Native Hawaiian serving institution)
• TRIBAL (Flag for tribal college and university)
• AANAPII (Flag for Asian American, Native American, and Pacific Islander-serving institution)
• HSI (Flag for Hispanic-serving institution)
• NANTI (Flag for Native American non-tribal institution)
• RELAFFIL (Religious affiliation of the institution)

Let’s take a look it look at what values one of these columns has.

table(college_scorecard_completely_prepped$PBI)

## 
##    0    1 NULL 
## 5949   94    3

Like the sex-specific variables above, we have null values which don’t give us useful information. Moreover, because there are so many categories, we don’t want to create a chart for each one. Instead, we will be fusing the first 7 variables that are flagged for particular ethnicities, just like we did for the sex-specific categories.

Additionally, because the religious column has dozens of categories, creating a chart for all of them does not make sense. Instead, we are going to collapse all of the categories into religious and non-religious.

# Data Manipulations 
library(forcats) 


# Change values so that they are strings 
college_scorecard_cat_demographics_rw <- college_scorecard_completely_prepped %>%  
  mutate(HBCU = dplyr::recode(HBCU,
                             '0' = 'Not HBCU', 
                             '1' = 'HBCU', 
                             'NULL'= 'Not HBCU' )) %>%
  mutate(PBI = dplyr::recode(PBI,
                             '0' = 'Not PBI', 
                             '1' = 'PBI', 
                             'NULL'= 'Not PBI' )) %>%
  mutate(ANNHI = dplyr::recode(ANNHI,
                             '0' = 'Not ANNHI', 
                             '1' = 'ANNHI', 
                             'NULL'= 'Not ANNHI' )) %>%
  mutate(TRIBAL = dplyr::recode(TRIBAL,
                             '0' = 'Not TRIBAL', 
                             '1' = 'TRIBAL', 
                             'NULL'= 'Not TRIBAL' )) %>%
  mutate(AANAPII = dplyr::recode(AANAPII,
                             '0' = 'Not AANAPII', 
                             '1' = 'AANAPII', 
                             'NULL'= 'Not AANAPII' )) %>%
  mutate(HSI = dplyr::recode(HSI,
                             '0' = 'Not HSI', 
                             '1' = 'HSI', 
                             'NULL'= 'Not HSI' )) %>%
  mutate(NANTI = dplyr::recode(NANTI,
                             '0' = 'Not NANTI', 
                             '1' = 'NANTI', 
                             'NULL'= 'Not NANTI' )) %>%  
# Collapse categories so that we have No Affiliation and  Religious Affiliation 
  mutate(RELAFFIL = fct_other(RELAFFIL,
                              keep = "NULL"))%>%
  mutate(RELAFFIL = dplyr::recode(RELAFFIL, 
                           "NULL" = "No Affiliation", 
                           "Other" = "Religious Affiliation"))

# Fuse categories into one variable 
college_scorecard_cat_demographics_fused <-  unite(college_scorecard_cat_demographics_rw, HBCU, PBI, ANNHI, TRIBAL, AANAPII,HSI, NANTI, col = "DEMOGRAPHICS", sep = "/")

table(college_scorecard_cat_demographics_fused$DEMOGRAPHICS)

## 
##     HBCU/Not PBI/Not ANNHI/Not TRIBAL/Not AANAPII/Not HSI/Not NANTI 
##                                                                  99 
##         Not HBCU/Not PBI/ANNHI/Not TRIBAL/AANAPII/Not HSI/Not NANTI 
##                                                                  17 
##             Not HBCU/Not PBI/ANNHI/Not TRIBAL/Not AANAPII/HSI/NANTI 
##                                                                   1 
##         Not HBCU/Not PBI/ANNHI/Not TRIBAL/Not AANAPII/Not HSI/NANTI 
##                                                                   9 
##     Not HBCU/Not PBI/ANNHI/Not TRIBAL/Not AANAPII/Not HSI/Not NANTI 
##                                                                   2 
##             Not HBCU/Not PBI/ANNHI/TRIBAL/AANAPII/Not HSI/Not NANTI 
##                                                                   1 
##         Not HBCU/Not PBI/ANNHI/TRIBAL/Not AANAPII/Not HSI/Not NANTI 
##                                                                   5 
##         Not HBCU/Not PBI/Not ANNHI/Not TRIBAL/AANAPII/HSI/Not NANTI 
##                                                                  80 
##     Not HBCU/Not PBI/Not ANNHI/Not TRIBAL/AANAPII/Not HSI/Not NANTI 
##                                                                  60 
##         Not HBCU/Not PBI/Not ANNHI/Not TRIBAL/Not AANAPII/HSI/NANTI 
##                                                                   1 
##     Not HBCU/Not PBI/Not ANNHI/Not TRIBAL/Not AANAPII/HSI/Not NANTI 
##                                                                 367 
##     Not HBCU/Not PBI/Not ANNHI/Not TRIBAL/Not AANAPII/Not HSI/NANTI 
##                                                                  17 
## Not HBCU/Not PBI/Not ANNHI/Not TRIBAL/Not AANAPII/Not HSI/Not NANTI 
##                                                                5264 
##     Not HBCU/Not PBI/Not ANNHI/TRIBAL/Not AANAPII/Not HSI/Not NANTI 
##                                                                  29 
##             Not HBCU/PBI/Not ANNHI/Not TRIBAL/AANAPII/HSI/Not NANTI 
##                                                                   1 
##         Not HBCU/PBI/Not ANNHI/Not TRIBAL/AANAPII/Not HSI/Not NANTI 
##                                                                   1 
##         Not HBCU/PBI/Not ANNHI/Not TRIBAL/Not AANAPII/HSI/Not NANTI 
##                                                                   8 
##     Not HBCU/PBI/Not ANNHI/Not TRIBAL/Not AANAPII/Not HSI/Not NANTI 
##                                                                  84

As we can see from the various categories that were created, there are many schools that fit into multiple categories. Again, having many categories does not give us a clear picture. The solution to this problem is to collapse the categories that have multiple flags into a single one called “Multiple.”

#Clean up categories
college_scorecard_cat_demographics_complete <-  college_scorecard_cat_demographics_fused %>%
  mutate(DEMOGRAPHICS = dplyr::recode(DEMOGRAPHICS,
                            'HBCU/Not PBI/Not ANNHI/Not TRIBAL/Not AANAPII/Not HSI/Not NANTI' = 'HBCU', 
                            'Not HBCU/Not PBI/ANNHI/Not TRIBAL/Not AANAPII/HSI/NANTI' = 'ANNHI/HSI/NANTI',
                            'Not HBCU/Not PBI/ANNHI/Not TRIBAL/Not AANAPII/Not HSI/Not NANTI' = 'ANNHI',
                            'Not HBCU/Not PBI/ANNHI/TRIBAL/Not AANAPII/Not HSI/Not NANTI' = 'ANNHI/TRIBAL',
                            'Not HBCU/Not PBI/Not ANNHI/Not TRIBAL/AANAPII/Not HSI/Not NANTI' = 'AANAPII',
                            'Not HBCU/Not PBI/Not ANNHI/Not TRIBAL/Not AANAPII/HSI/Not NANTI' = 'HSI',
                            'Not HBCU/Not PBI/Not ANNHI/Not TRIBAL/Not AANAPII/Not HSI/Not NANTI' = 'None',
                            'Not HBCU/PBI/Not ANNHI/Not TRIBAL/AANAPII/HSI/Not NANTI'= 'PBI/AANAPII/HSI',
                            'Not HBCU/PBI/Not ANNHI/Not TRIBAL/Not AANAPII/HSI/Not NANTI' = 'PBI/HSI', 
                            
                            'Not HBCU/Not PBI/ANNHI/Not TRIBAL/AANAPII/Not HSI/Not NANTI' = 'ANNHI/AANAPII',
                            'Not HBCU/Not PBI/ANNHI/Not TRIBAL/Not AANAPII/Not HSI/NANTI' = 'ANNHI/NANTI',
                            'Not HBCU/Not PBI/ANNHI/TRIBAL/AANAPII/Not HSI/Not NANTI' = 'ANNHI/TRIBAL/AANAPII',
                            'Not HBCU/Not PBI/Not ANNHI/Not TRIBAL/AANAPII/HSI/Not NANTI' = 'AANAPII/HSI',
                            'Not HBCU/Not PBI/Not ANNHI/Not TRIBAL/Not AANAPII/HSI/NANTI'= 'HSI/NANTI',
                            'Not HBCU/Not PBI/Not ANNHI/Not TRIBAL/Not AANAPII/Not HSI/NANTI' = 'NANTI',
                            'Not HBCU/Not PBI/Not ANNHI/TRIBAL/Not AANAPII/Not HSI/Not NANTI' = 'TRIBAL',
                            'Not HBCU/PBI/Not ANNHI/Not TRIBAL/AANAPII/Not HSI/Not NANTI' = 'PBI/AANAPII',
                            'Not HBCU/PBI/Not ANNHI/Not TRIBAL/Not AANAPII/Not HSI/Not NANTI' = 'PBI')) %>%
  # Collapse all multi-categories into a single one
  mutate(DEMOGRAPHICS = fct_collapse(DEMOGRAPHICS, 
                                     "Multiple" = c( 'ANNHI/HSI/NANTI', 'ANNHI/TRIBAL', 'PBI/AANAPII/HSI',
                                                     'PBI/HSI', 'ANNHI/AANAPII', 'ANNHI/NANTI',
                                                     'ANNHI/TRIBAL/AANAPII', 'AANAPII/HSI', 
                                                     'HSI/NANTI', 'PBI/AANAPII'
                                                  )))

table(college_scorecard_cat_demographics_complete$DEMOGRAPHICS)

## 
##  AANAPII Multiple    ANNHI     HBCU      HSI    NANTI     None      PBI 
##       60      124        2       99      367       17     5264       84 
##   TRIBAL 
##       29

The 18 categories have been reduced to 10. We can now proceed to the next phase, which is to produce a plot that will test if there is a relationship between sex-specificity and demographics. One of the options we can use to plot a categorical variable against another categorical variable is a mosaic plot.

library(vcd)
tbl_status<- xtabs(~ DEMOGRAPHICS + SEX_SPECIFIC, college_scorecard_cat_demographics_complete)
ftable(tbl_status)

##              SEX_SPECIFIC CO-ED MEN ONLY WOMEN ONLY
## DEMOGRAPHICS                                       
## AANAPII                      59        0          1
## Multiple                    122        0          2
## ANNHI                         2        0          0
## HBCU                         96        1          2
## HSI                         365        0          2
## NANTI                        17        0          0
## None                       5176       60         28
## PBI                          84        0          0
## TRIBAL                       29        0          0

mosaic(tbl_status,
       labeling_args = list(rot_labels = c(top = 90, left = 0),
                            gp_labels=gpar(fontsize= 9), 
 offset_varnames = c(top = 2, left = 2), offset_labels = c(left = 0.5, top =2)), 
 spacing = spacing_increase(start = unit(0.45, "lines"), rate = 1), 
 margins = c(top = 0.25, bottom = 0.5),
 gp = shading_hcl,
 legend = T)

In this plot, each rectangle represents the intersection of 1 category from one variable and another category from another variable. For example, the largest rectangle represents the number of colleges that are both co-ed and have no ethnic flag. The larger the area of the rectangle, the more observations are in it. If there is no intersection between the variable categories, a circle with a line is used.

The Pearson residual colors, essentially, tell us whether we have more or fewer colleges than expected if we were to assume that demographics and sex-specificity are independent (as in there is no relationship). The presence of blue and red tells us that the variables are not independent meaning they have some kind of relationship. As we can clearly see in the plot, all of the rectangles are grey, which indicates that the two variables are independent.

Now, if we take a look at religious affiliation and plot it against sex-specificity, we might get something different.

tbl_status2<- xtabs(~ RELAFFIL +SEX_SPECIFIC, college_scorecard_cat_demographics_complete)
ftable(tbl_status2)

##                       SEX_SPECIFIC CO-ED MEN ONLY WOMEN ONLY
## RELAFFIL                                                    
## No Affiliation                      5107       28         19
## Religious Affiliation                843       33         16

mosaic(tbl_status2,
       labeling_args = list(rot_labels = c(top = 90, left = 0),
                            gp_labels=gpar(fontsize= 8), 
 offset_varnames = c(top = 2, left = 4), offset_labels = c(left = 3, top =2)), 
 spacing = spacing_increase(start = unit(0.45, "lines"), rate = 1), 
 margins = c(top = 0.25, bottom = 0.5),
 gp = shading_hcl,
 legend = T)

Now, we have some color! It looks like for men-only and women-only colleges, there are more values than expected when they are religiously affiliated. This matches some of the observations we made when we first looked at the map.

Let’s take a less confusing look by plotting the relationship between these two variables using an interactive stacked bar chart.

library(extrafont)
library(plotly)
plot_data <- college_scorecard_cat_demographics_complete %>%
  count(RELAFFIL, SEX_SPECIFIC) %>%
  add_count(SEX_SPECIFIC, wt= n) %>%
  mutate(percentage = (n/nn))




  plot_ly(x = plot_data[['SEX_SPECIFIC']], y = plot_data[['percentage']], color = plot_data[['RELAFFIL']], colors = c("#f1cb35",
"#cc8834")) %>%
  add_bars() %>%
  layout(barmode = "stack",
         font = list(family = "Georgia"), 
         yaxis = list(title = "Percentage of Colleges", tickformat = "%"),
         xaxis = list(title= "College Type" ))

From this chart, we can clearly see that around half of both men-only and women-only colleges are religious wheres 14% percent of coed colleges are religious.

To continue our examination of the relationship between sex-specificity and demographics, we will be using a logistic regression framework. Specifically, we want to answer the question “does the proportion of ethnic groups help predict sex-specificity?”

From our data set there are a number of numeric demographic variables to examine:

• UGDS_WHITE (proportion of students that are white)
• UGDS_BLACK (proportion of students that are black)
• UGDS_HISP (proportion of students that are Hispanic )
• UGDS_ASIAN (proportion of students that are Asian )
• UGDS_AIAN (proportion of students that are American Indian/ Alaskan natives)
• UGDS_NHPI (proportion of students that Native Hawaiian/Pacific Islander)
• UGDS_2MOR (proportion of students that are 2 or more ethnicities)
• UGDS_NRA (proportion of students that are non-resident aliens [ I may refer to them as international])
• UGDS_UNKN (proportion of students with unknown race)

With so many variables to consider, the easiest way to investigate their relationship with sex-specificity is to create a logistic regression summary table and see which variables have the most significant relationships. Because we need the variables to be binary, we will be using one of the data sets from earlier on in our project. There are a few clean-up steps like making sure the MENONLY and WOMENONLY columns are numeric and have no null values.

#Make character variables to be numeric 
college_scorecard_full_num_demographics<- college_scorecard_clean[, c(2:3) :=lapply(.SD, as.numeric), .SDcols= c(2:3)]

college_scorecard_dem_glance<- college_scorecard_full_num_demographics %>%
  select(UGDS, UGDS_WHITE, UGDS_BLACK, UGDS_HISP, UGDS_ASIAN,
          UGDS_AIAN, UGDS_NHPI, UGDS_2MOR, UGDS_NRA, UGDS_UNKN)

# Visualize Missing Values
library(visdat)
vis_miss(college_scorecard_dem_glance)+
  theme(axis.text.x = element_text(angle = 90, hjust = 0.65),
        text = element_text(family = "Georgia"))

Now our goal is to get rid of the missing rows so that we can run our logistic regression models.

college_scorecard_full_num_demographics2 <- college_scorecard_full_num_demographics %>%
  filter(UGDS != is.na(UGDS))


college_scorecard_dem_glance2<- college_scorecard_full_num_demographics2 %>%
  select(UGDS, UGDS_WHITE, UGDS_BLACK, UGDS_HISP, UGDS_ASIAN,
          UGDS_AIAN, UGDS_NHPI, UGDS_2MOR, UGDS_NRA, UGDS_UNKN)

vis_miss(college_scorecard_dem_glance2)+
  theme(axis.text.x = element_text(angle = 90, hjust = 0.65),
        text = element_text(family = "Georgia"))

#Exclude non-operational colleges 
college_scorecard_full_num_demographics_final<- college_scorecard_full_num_demographics2[CURROPER == "1"]

One important thing to address is the amount of colleges that were removed because of these clean up steps. We went from 6806 colleges to 5751. Because the missing values were consistent throughout these numeric demographics, the most likely reason they were left blank was because the researchers did not have access to the information.

First, we will create a table that shows the relationship between each explanatory numeric demographic variable and the MENONLY response variable.

library(broom)
#Create logistic regression model for predicting male-only colleges 
demographics_mod_male1 <- tidy(glm(MENONLY ~  UGDS_WHITE,  family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_male2 <- tidy(glm(MENONLY ~  UGDS_BLACK, family= binomial,  data = college_scorecard_full_num_demographics_final))


demographics_mod_male3 <- tidy(glm(MENONLY ~   UGDS_HISP, family= binomial,  data = college_scorecard_full_num_demographics_final))
demographics_mod_male4 <- tidy(glm(MENONLY ~   UGDS_ASIAN, family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_male5 <- tidy(glm(MENONLY ~   UGDS_AIAN,  family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_male6 <- tidy(glm(MENONLY ~  UGDS_NHPI,  family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_male7 <- tidy(glm(MENONLY ~  UGDS_2MOR,  family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_male8 <- tidy(glm(MENONLY ~  UGDS_NRA,  family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_male9 <- tidy(glm(MENONLY ~  UGDS_UNKN,  family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_male_main<-bind_rows(demographics_mod_male1, demographics_mod_male2, demographics_mod_male3,
          demographics_mod_male4, demographics_mod_male5, demographics_mod_male6,
          demographics_mod_male7, demographics_mod_male8, demographics_mod_male9)

demographics_mod_male_main

Because we cannot visualize a logistic regression that is more than 4 dimensions (3 explanatory variables and 1 response variable), we need to select the explanatory variables that are the most significant for our final model. There is also an extra incentive to get rid of some of these demographic variables because they are closely related and could create numerical instability. From this table, we can see that most of the explanatory variables independently have significant relationships with the MENONLY variable. The NRA variable and NHPI variables seem to have the least significant relationships with the MENONLY variable.

Now we will create a multivariate logistic model that takes into account all of the remaining variables to see if we can eliminate more explanatory variables.

demographics_mod_male_main2 <- glm(MENONLY ~   UGDS_WHITE +  UGDS_BLACK + UGDS_HISP + UGDS_ASIAN + 
UGDS_AIAN +  UGDS_2MOR +  UGDS_UNKN, family= binomial,  data = college_scorecard_full_num_demographics_final)

summary(demographics_mod_male_main2)

## 
## Call:
## glm(formula = MENONLY ~ UGDS_WHITE + UGDS_BLACK + UGDS_HISP + 
##     UGDS_ASIAN + UGDS_AIAN + UGDS_2MOR + UGDS_UNKN, family = binomial, 
##     data = college_scorecard_full_num_demographics_final)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6941  -0.0479  -0.0128  -0.0023   5.4707  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.4389     0.8093  -1.778 0.075415 .  
## UGDS_WHITE     0.1384     0.8362   0.166 0.868534    
## UGDS_BLACK   -10.8152     3.0340  -3.565 0.000364 ***
## UGDS_HISP     -4.9057     1.5689  -3.127 0.001767 ** 
## UGDS_ASIAN    -7.5415     5.6959  -1.324 0.185497    
## UGDS_AIAN   -201.6041   114.6773  -1.758 0.078746 .  
## UGDS_2MOR    -73.5400    20.1115  -3.657 0.000256 ***
## UGDS_UNKN    -65.3860    28.0082  -2.335 0.019568 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 648.63  on 5750  degrees of freedom
## Residual deviance: 364.70  on 5743  degrees of freedom
## AIC: 380.7
## 
## Number of Fisher Scoring iterations: 13

Here, some of these variables seem to have very large standard errors. The larger the standard error, the less precise the statistic, meaning the samples may not closely represent the population. A large standard error could be caused by having a small sample size, a large standard deviation, or a combination of both. To keep things simple, we will be removing the AIAN, 2MOR, and UNKN variables too because of their very large standard errors.

demographics_mod_male3 <- glm(MENONLY ~    UGDS_WHITE +  UGDS_BLACK + UGDS_HISP  +UGDS_ASIAN, family= binomial,  data = college_scorecard_full_num_demographics_final)

summary(demographics_mod_male3)

## 
## Call:
## glm(formula = MENONLY ~ UGDS_WHITE + UGDS_BLACK + UGDS_HISP + 
##     UGDS_ASIAN, family = binomial, data = college_scorecard_full_num_demographics_final)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5521  -0.1057  -0.0300  -0.0063   6.6810  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -6.932      2.089  -3.318 0.000907 ***
## UGDS_WHITE     5.128      2.153   2.382 0.017219 *  
## UGDS_BLACK   -16.316      5.175  -3.153 0.001618 ** 
## UGDS_HISP     -2.536      3.505  -0.724 0.469364    
## UGDS_ASIAN   -21.548     12.738  -1.692 0.090718 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 648.63  on 5750  degrees of freedom
## Residual deviance: 464.62  on 5746  degrees of freedom
## AIC: 474.62
## 
## Number of Fisher Scoring iterations: 10

Now, our table shows UGDS_ASIAN is not only statistically insignificant but also has a large standard error compared to the other three variables. Additionally, UGDS_HISP also appears to be statistically insignificant. With this in mind, our final model will only incorporate two variables (UGDS_WHITE and UGDS_BLACK).

But before we do the finalizing, we should perform a couple of steps. First, we should check if there is collinearity between our two variables so that it does not create instability in our model.

#Collinearity Check 
college_scorecard_full_num_demographics_final %>%
  plot_ly(y = ~UGDS_BLACK, x= ~UGDS_WHITE, colors = "#44668b", text = ~INSTNM) %>%
  add_markers()

#Correlation Check 
cor(college_scorecard_full_num_demographics_final$UGDS_WHITE, college_scorecard_full_num_demographics_final$UGDS_BLACK)

## [1] -0.4719334

So it looks like the variables are very correlated but not so much that they create instability.

Second, we should compare our model (with proportion of black and white students) to a model without black students so that we can see if the standard errors drastically change.

We are just going to rename the UGDS_BLACK variables so that our final plot is easier to to read.

# Rename variables 
college_scorecard_full_num_demographics_final2 <- college_scorecard_full_num_demographics_final %>%
  rename( Black_Students = UGDS_BLACK)
  
  
#  Final Male Models
demographics_mod_male_final <- glm(MENONLY ~    UGDS_WHITE + Black_Students, family= binomial,  
                                   data = college_scorecard_full_num_demographics_final2) # White and Black Students 


demographics_mod_male_final2 <- glm(MENONLY ~    UGDS_WHITE,  family= binomial,  
                                   data = college_scorecard_full_num_demographics_final2) # White Students Only
summary(demographics_mod_male_final)

## 
## Call:
## glm(formula = MENONLY ~ UGDS_WHITE + Black_Students, family = binomial, 
##     data = college_scorecard_full_num_demographics_final2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5547  -0.1080  -0.0323  -0.0089   6.8577  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -9.667      1.343  -7.199 6.05e-13 ***
## UGDS_WHITE        7.873      1.434   5.492 3.98e-08 ***
## Black_Students  -14.706      5.223  -2.816  0.00487 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 648.63  on 5750  degrees of freedom
## Residual deviance: 468.47  on 5748  degrees of freedom
## AIC: 474.47
## 
## Number of Fisher Scoring iterations: 10

summary(demographics_mod_male_final2)

## 
## Call:
## glm(formula = MENONLY ~ UGDS_WHITE, family = binomial, data = college_scorecard_full_num_demographics_final2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5465  -0.1174  -0.0404  -0.0084   5.1060  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -13.072      1.108 -11.796   <2e-16 ***
## UGDS_WHITE    11.246      1.246   9.026   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 648.63  on 5750  degrees of freedom
## Residual deviance: 478.54  on 5749  degrees of freedom
## AIC: 482.54
## 
## Number of Fisher Scoring iterations: 10

The standard errors between models don’t differ much, so we are free to use our first one (with both black and white students).

Our final model shows that there is a significant relationship between being a men-only institution and particular proportions of certain ethnicities. Specifically, we can see that there is a positive relationship between men-only institutions and white students as well as a negative relationship between men-only institutions and black students. This finding does not seem too surprising because our map above showed that many of these institutions were seminary schools, which one can safely assume are mostly white men. But for the sake of transparency let’s take a quick look at a table of these institutions:

library(knitr)
Men_only_glance <- Men_only %>% 
  select(INSTNM, UGDS_WHITE, UGDS_BLACK)

kable(Men_only_glance,  col.names = c('Name',  'Prop. White', 'Prop. Black'), align = 'cc', caption = "Table 1.1 The Proportion of White and Black Students Per Institution")

For women-only colleges, we will follow a similar procedure.

#Create logistic regression model for predicting female-only colleges 
demographics_mod_female1 <- tidy(glm(WOMENONLY ~  UGDS_WHITE,  family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_female2 <- tidy(glm(WOMENONLY ~  UGDS_BLACK, family= binomial,  data = college_scorecard_full_num_demographics_final))


demographics_mod_female3 <- tidy(glm(WOMENONLY ~   UGDS_HISP, family= binomial,  data = college_scorecard_full_num_demographics_final))
demographics_mod_female4 <- tidy(glm(WOMENONLY ~   UGDS_ASIAN, family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_female5 <- tidy(glm(WOMENONLY ~   UGDS_AIAN,  family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_female6 <- tidy(glm(WOMENONLY ~  UGDS_NHPI,  family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_female7 <- tidy(glm(WOMENONLY ~  UGDS_2MOR,  family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_female8 <- tidy(glm(WOMENONLY ~  UGDS_NRA,  family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_female9 <- tidy(glm(WOMENONLY ~  UGDS_UNKN,  family= binomial,  data = college_scorecard_full_num_demographics_final))

demographics_mod_female_main<-bind_rows(demographics_mod_female1, demographics_mod_female2, demographics_mod_female3,
          demographics_mod_female4, demographics_mod_female5, demographics_mod_female6,
          demographics_mod_female7, demographics_mod_female8, demographics_mod_female9)

demographics_mod_female_main

Based on some of the standard errors in this table, we will be getting rid of UGDS_AIAN and UGDS_NHPI. Additionally, there appears to be only two variables that are significant (UGDS_NRA, UGDS_2MOR).

demographics_mod_female_main <- glm(WOMENONLY ~  UGDS_2MOR + UGDS_NRA, family= binomial,  data = college_scorecard_full_num_demographics_final)

summary(demographics_mod_female_main)

## 
## Call:
## glm(formula = WOMENONLY ~ UGDS_2MOR + UGDS_NRA, family = binomial, 
##     data = college_scorecard_full_num_demographics_final)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5849  -0.1113  -0.1052  -0.0994   3.2663  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.3470     0.2078 -25.737   <2e-16 ***
## UGDS_2MOR     4.5186     2.4720   1.828   0.0676 .  
## UGDS_NRA      2.8175     1.2308   2.289   0.0221 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 426.91  on 5750  degrees of freedom
## Residual deviance: 421.33  on 5748  degrees of freedom
## AIC: 427.33
## 
## Number of Fisher Scoring iterations: 8

Because UGDS_2MOR did not reach our p-value threshold, we will exclude it from the final model.

demographics_mod_female_final <- glm(WOMENONLY ~ UGDS_NRA, family= binomial,  data = college_scorecard_full_num_demographics_final)

summary(demographics_mod_female_final)

## 
## Call:
## glm(formula = WOMENONLY ~ UGDS_NRA, family = binomial, data = college_scorecard_full_num_demographics_final)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.4166  -0.1084  -0.1057  -0.1057   3.2217  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.1841     0.1791  -28.94   <2e-16 ***
## UGDS_NRA      2.7833     1.2046    2.31   0.0209 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 426.91  on 5750  degrees of freedom
## Residual deviance: 423.62  on 5749  degrees of freedom
## AIC: 427.62
## 
## Number of Fisher Scoring iterations: 8

Based on this table, we can conclude that there is a significant relationship between proportion of nonresident alien students and women-only colleges. This finding was unexpected. I think a look at what is going on would be helpful:

library(knitr)
Women_only_glance <- Women_only %>% 
  select(INSTNM, UGDS_NRA)


kable(Women_only_glance, col.names = c('Name',  'Prop. Nonresident Aliens'), align = 'cc', caption = "Table 1.2 The Proportion of Nonresident Students Per Institution")

nra_avges <- college_scorecard_completely_prepped %>% 
  select(SEX_SPECIFIC, UGDS_NRA) %>% 
  group_by(SEX_SPECIFIC) %>% 
  summarize( UGDS_NRA_AVG= mean(UGDS_NRA, na.rm = TRUE))

kable(nra_avges, col.names = c('Institution Type',  'Mean Proportion Nonresident Aliens'), align = 'cc', caption = "Table 1.3 The Average Proportion of Nonresident Students Per Type of Institution")

Although an average of 4.9 % of the student body being international does not seem like much, when compared to the other types of schools, it looks like women colleges have more international students.

library(visreg)
visreg(demographics_mod_male_final, "UGDS_WHITE", by = "Black_Students",
       rug = FALSE ,
       band= FALSE,
       xlab = "Proportion of White Students", 
       ylab = "Probabiltiy of Men-Only College",
       line= list(col ="#1298d0"),
       scale="response",
       gg = TRUE)+ 
  theme_bw()+ 
  labs(title = "Logistic Regression Prediction of Being Men-Only College")+
  theme(legend.position = "none",
        text = element_text(family = "Georgia"), 
        plot.title = element_text(size = 15, margin(b = 10), hjust = 0.5, family = "Georgia"), 
        strip.background = element_rect(
        color="Black", fill="#1298d0",  linetype="solid"))

visreg(demographics_mod_female_final, "UGDS_NRA", 
       rug = FALSE ,
       band= FALSE,
       xlab = "Proportion of Nonresident Aliens", 
       ylab = "Probabiltiy of Being Women-Only College",
       line= list(col ="#cc8834"),
       scale="response",
       gg = TRUE)+ 
  theme_bw()+ 
  labs(title = "Logistic Regression Prediction of Women-Only College")+
  theme(legend.position = "none",
        text = element_text(family = "Georgia"), 
        plot.title = element_text(size = 15, margin(b = 10), hjust = 0.5, family = "Georgia") )

As we can see from the graphs, the relationships between these demographic variables and sex-specificity are still fairly small despite the statistical significance of the models. Testing the accuracy or recall of the models does not seem appropriate given how weak the relationships are.

Admission Statistics

During the previous section, we had the numeric variable be the explanatory variable and the sex-specificity (a categorical variable) be the response variable. In this section and for the rest of the sections, sex-specificity is going to become the response variable. This transition was made, firstly, because creating logistic regression models for each one of these questions is going to be too much of a pain. Secondly, because these variables under investigation seem to have, for the most part, small relationships with gender, I feel the best thing would be to use a logistic regression framework using only the most significant variables. What would probably work best for an exploratory analysis like this one would be an anova framework that allows us to get a quick glance at the data using bar charts.

Our next step is to answer, “how do admission standards differ across the sex-specific colleges?”

Here are our explanatory variables:

• ADM_RATE (Admission Rate)
• SATVRMID (SAT Verbal Median)
• SATMTMID (SAT Math Median)
• SATWRMID (SAT Writing Median)
• SAT_AVG (SAT Overall Average)
• ACTCMMID (ACT Cumulative Median)
• ACTENMID (ACT English Median)
• ACTMTMID (ACT Math Median)
• ACTWRMID (ACT Writing Median)

Because there is a strong possibility that these variables will be collinear, it would be a good idea to take a look at what the overall averages are for each sex:

college_scorecard_completely_prepped %>%
  select(SEX_SPECIFIC, 
         ADM_RATE, 
         SATVRMID, SATMTMID, SAT_AVG,
         ACTCMMID, ACTENMID, 
         ACTMTMID, ACTWRMID)%>%
  group_by(SEX_SPECIFIC)  %>%
  summarize(ADM_RATE = mean(ADM_RATE, na.rm= TRUE), 
         SATVRMID= mean(SATVRMID, na.rm= TRUE),  SATMTMID= mean(SATMTMID, na.rm= TRUE), SAT_AVG= mean(SAT_AVG, na.rm= TRUE),
         ACTCMMID= mean(ACTCMMID, na.rm= TRUE), ACTENMID=  mean(ACTENMID, na.rm= TRUE), 
         ACTMTMID= mean(ACTMTMID, na.rm= TRUE), ACTWRMID= mean(ACTWRMID, na.rm= TRUE)) %>%
  kable(align = 'cc', caption = "Table 1.4 The Average Adission Statistics Per Type of Institution")

Women-only colleges being on average more selective is not surprising. With only 35 of them and there being very well known selective institutions like Scripps College and Barnard College, these findings were more or less expected. However, what came as a bit of a surprise was the selectivity of men-only colleges. Normally, as admission rates go down scores go up. Notice that the admission rates for men-only colleges is much higher than co-ed and women-only colleges. If men-only colleges are more selective, then we should see an admission rate under 60 percent. Additionally, since none of the men-only colleges are as selective as the most selective women-only colleges, what is driving these averages up so high?

Let’s take a look:

admissions_table <- function(sex) {
  college_scorecard_completely_prepped %>%
  filter(SEX_SPECIFIC == sex)  %>%
  select(INSTNM, 
         ADM_RATE, 
         SAT_AVG, SATVRMID, SATMTMID,
         ACTCMMID, ACTENMID, 
         ACTMTMID, ACTWRMID) %>%
    kable(align = 'cc', caption = "Tables 1.4 Admission Statistics for Institution Type")
}

admissions_table("MEN ONLY")

admissions_table("WOMEN ONLY")

I knew it! Something was indeed going on here! Because most of the male colleges have NA’s, for their admission test statistics, only 5/61 contributed their numbers to those averages. This means that while most of the women-only college standards were averaged, only the more selective men-only college averages were used. This also makes sense because if most of your colleges are seminary schools, how could your average be so high?

This finding actually makes our job easier because the only valid variable to use is admission rate given the scarcity of data for the admission tests.

summary( lm( ADM_RATE ~ SEX_SPECIFIC, data= college_scorecard_completely_prepped ))

## 
## Call:
## lm(formula = ADM_RATE ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.67684 -0.12459  0.02196  0.16496  0.38945 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.676842   0.004976 136.027  < 2e-16 ***
## SEX_SPECIFICMEN ONLY    0.155375   0.033817   4.595 4.61e-06 ***
## SEX_SPECIFICWOMEN ONLY -0.067196   0.036978  -1.817   0.0693 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2168 on 1972 degrees of freedom
##   (4071 observations deleted due to missingness)
## Multiple R-squared:  0.0124, Adjusted R-squared:  0.01139 
## F-statistic: 12.38 on 2 and 1972 DF,  p-value: 4.556e-06

Here, we can clearly see that there is a significant positive relationship between admission rates and men-only institutions. Although women colleges did not cross the significance threshold, you can clearly see that there is a negative relationship between admission rate and women-only colleges.

Majors

For this section, I had two interesting questions:

Does attending a sex-specific school affect degree outcomes? Does going to a women-only college affect STEM degree outcomes?

Because the US Department of Education collects data on the percentage of students that graduate with particular majors, there are a lot of variables to work with!

However, since we are only interested in major academic divisions, we should be able to fuse most of these percentages and compare them that way. I decided to group them by three major divisions: Liberal Arts, STEM, and Professional/Vocational/Other.

Disclaimer: Many of these majors could have fit into a couple of groups, but in the interest of keeping things simple, I based my decisions mostly on this wikipedia page about academics fields.

Here is the breakdown of the variables:

college_scorecard_majors<-  college_scorecard_completely_prepped %>%
  transmute(INSTNM, 
            SEX_SPECIFIC , 
            LIB_ARTS = PCIP05 + PCIP13 + PCIP16 + PCIP23 + PCIP24 + PCIP38 + PCIP39 +PCIP42 + PCIP45 + PCIP50 + PCIP54, 
            STEM = PCIP01 + PCIP03 + PCIP11 + PCIP14 + PCIP15 + PCIP26 + PCIP27+PCIP29 + PCIP40 +  PCIP41 +  PCIP51, 
            PROF_VOC_OTHERS =  PCIP04 + PCIP09 + PCIP10 + PCIP12 + PCIP19 + PCIP22 +
              PCIP25 + PCIP30 + PCIP31 + PCIP43 + PCIP44 + PCIP46 + PCIP47 + PCIP48 + PCIP49 + PCIP52 )



kable(col.names = c("Name", 'Institution Type',  'Liberal Arts', 'STEM', "Professional & Vocation"), head(college_scorecard_majors),
      align = 'cc', caption = "Table 1.5: Sample Proportion of Degree's Conferred Per Instituion")

summary(lm(LIB_ARTS ~ SEX_SPECIFIC, data = college_scorecard_majors)) 

## 
## Call:
## lm(formula = LIB_ARTS ~ SEX_SPECIFIC, data = college_scorecard_majors)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9298 -0.2005 -0.1551  0.1407  0.7996 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            0.200534   0.003495  57.374  < 2e-16 ***
## SEX_SPECIFICMEN ONLY   0.729223   0.034722  21.002  < 2e-16 ***
## SEX_SPECIFICWOMEN ONLY 0.263843   0.044608   5.915 3.52e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2631 on 5756 degrees of freedom
##   (287 observations deleted due to missingness)
## Multiple R-squared:  0.0761, Adjusted R-squared:  0.07578 
## F-statistic: 237.1 on 2 and 5756 DF,  p-value: < 2.2e-16

summary(lm(STEM ~ SEX_SPECIFIC, data = college_scorecard_majors)) 

## 
## Call:
## lm(formula = STEM ~ SEX_SPECIFIC, data = college_scorecard_majors)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34488 -0.34488 -0.06808  0.17039  0.65512 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.344879   0.004443  77.621  < 2e-16 ***
## SEX_SPECIFICMEN ONLY   -0.327114   0.044139  -7.411 1.44e-13 ***
## SEX_SPECIFICWOMEN ONLY -0.011948   0.056706  -0.211    0.833    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3344 on 5756 degrees of freedom
##   (287 observations deleted due to missingness)
## Multiple R-squared:  0.009455,   Adjusted R-squared:  0.009111 
## F-statistic: 27.47 on 2 and 5756 DF,  p-value: 1.335e-12

summary(lm(PROF_VOC_OTHERS ~ SEX_SPECIFIC, data = college_scorecard_majors)) 

## 
## Call:
## lm(formula = PROF_VOC_OTHERS ~ SEX_SPECIFIC, data = college_scorecard_majors)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4500 -0.2691 -0.0723  0.3030  0.5501 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.449999   0.004682  96.114  < 2e-16 ***
## SEX_SPECIFICMEN ONLY   -0.432007   0.046512  -9.288  < 2e-16 ***
## SEX_SPECIFICWOMEN ONLY -0.247279   0.059754  -4.138 3.55e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3524 on 5756 degrees of freedom
##   (287 observations deleted due to missingness)
## Multiple R-squared:  0.01755,    Adjusted R-squared:  0.0172 
## F-statistic:  51.4 on 2 and 5756 DF,  p-value: < 2.2e-16

college_scorecard_majors %>% 
  group_by(SEX_SPECIFIC) %>% 
  summarize(mean_stem = mean(STEM, na.rm= TRUE), mean_lb_arts = mean(LIB_ARTS, na.rm= TRUE),mean_prof_voc = mean(PROF_VOC_OTHERS, na.rm= TRUE)) %>%
  plot_ly(x = ~SEX_SPECIFIC, y = ~mean_stem, type= 'bar', name= 'STEM Degrees', marker = list(color = "#44668b")) %>%
  add_trace(y = ~mean_lb_arts, name= 'Liberal Arts Degrees', marker = list(color = "#f1cb35"))%>%
  add_trace(y = ~mean_prof_voc, name= 'Professional and Vocational Degrees', marker = list(color = "#bdb8b8"))%>%
  layout(barmode = "group", 
         xaxis = list(title = "College Type"),
         yaxis = list(title = "Proportion of Degrees Conferred"),
         font = list(family = "Georgia"))

Earning Outcomes

Here, we want to know if earning outcomes are affected by sex-specific schools.

‘MN_EARN_WNE_P10’ (Mean earnings of students working and not enrolled 10 years after entry) ‘MD_EARN_WNE_P10’ (Median earnings of students working and not enrolled 10 years after entry)

earning_outcomes_table <- bind_rows(
tidy(lm(MN_EARN_WNE_P10 ~ SEX_SPECIFIC, data= college_scorecard_completely_prepped)),
tidy(lm(MD_EARN_WNE_P10~ SEX_SPECIFIC, data= college_scorecard_completely_prepped))
)

earning_outcomes_table

From this table, it looks like both mean and median earning outcomes are positively associated with women-only colleges, albeit slightly. Both measurements are negatively associated with men-only colleges, but they are insignificant. One thing to note is that the median outcome for men-only colleges is almost significant.

college_scorecard_completely_prepped %>% 
  group_by(SEX_SPECIFIC) %>% 
  summarize(avg_mean = mean(MN_EARN_WNE_P10, na.rm= TRUE), 
            avg_md = mean(MD_EARN_WNE_P10, na.rm= TRUE)) %>%
  plot_ly(x = ~SEX_SPECIFIC, y = ~avg_mean, type= 'bar', name= 'Mean Earning ', marker = list(color = "#44668b")) %>%
  add_trace(y = ~avg_md, name= 'Median Earnings', marker = list(color = "#f1cb35"))%>%
  layout(barmode = "group", 
         xaxis = list(title = "College Type"),
         yaxis = list(title = "Average Earnings After 10 Years"),
         font = list(family = "Georgia"))

Completion

Next, does being at a sex-specific school affect completion rates for different ethnic groups?

Variables:

• C150_4_WHITE (Completion rate for first-time, full-time students at four-year institutions (150% of expected time to completion) for white students)

• C150_4_BLACK (Completion rate for first-time, full-time students at four-year institutions (150% of expected time to completion) for black students)

• C150_4_HISP (Completion rate for first-time, full-time students at four-year institutions (150% of expected time to completion) for Hispanic students)

• C150_4_ASIAN (Completion rate for first-time, full-time students at four-year institutions (150% of expected time to completion) for Asian students)

• C150_4_AIAN (Completion rate for first-time, full-time students at four-year institutions (150% of expected time to completion) for American Indian/Alaska Native students)

• C150_4_NHPI (Completion rate for first-time, full-time students at four-year institutions (150% of expected time to completion) for Native Hawaiian/Pacific Islander students)

• C150_4_2MOR (Completion rate for first-time, full-time students at four-year institutions (150% of expected time to completion) for students of two-or-more-races)

• C150_4_NRA (Completion rate for first-time, full-time students at four-year institutions (150% of expected time to completion) for non-resident alien students)

• C150_4_UNKN (Completion rate for first-time, full-time students at four-year institutions (150% of expected time to completion) for students whose race is unknown)

summary(lm(C150_4_WHITE~ SEX_SPECIFIC, data= college_scorecard_completely_prepped))

## 
## Call:
## lm(formula = C150_4_WHITE ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.64224 -0.14239  0.01976  0.15054  0.61240 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.542386   0.004963 109.294  < 2e-16 ***
## SEX_SPECIFICMEN ONLY   -0.154786   0.031111  -4.975 7.04e-07 ***
## SEX_SPECIFICWOMEN ONLY  0.099857   0.039238   2.545    0.011 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2236 on 2113 degrees of freedom
##   (3930 observations deleted due to missingness)
## Multiple R-squared:  0.01481,    Adjusted R-squared:  0.01387 
## F-statistic: 15.88 on 2 and 2113 DF,  p-value: 1.433e-07

summary(lm(C150_4_BLACK~ SEX_SPECIFIC, data= college_scorecard_completely_prepped))

## 
## Call:
## lm(formula = C150_4_BLACK ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.52134 -0.18285 -0.02405  0.16015  0.59495 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            0.405052   0.005715  70.875  < 2e-16 ***
## SEX_SPECIFICMEN ONLY   0.184268   0.112312   1.641  0.10102    
## SEX_SPECIFICWOMEN ONLY 0.116285   0.042778   2.718  0.00662 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2508 on 1963 degrees of freedom
##   (4080 observations deleted due to missingness)
## Multiple R-squared:  0.005079,   Adjusted R-squared:  0.004065 
## F-statistic:  5.01 on 2 and 1963 DF,  p-value: 0.006753

summary(lm(C150_4_HISP~ SEX_SPECIFIC, data= college_scorecard_completely_prepped))

## 
## Call:
## lm(formula = C150_4_HISP ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.58758 -0.16029 -0.00389  0.16911  0.52788 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.476089   0.005599  85.024  < 2e-16 ***
## SEX_SPECIFICMEN ONLY   -0.003972   0.071672  -0.055  0.95581    
## SEX_SPECIFICWOMEN ONLY  0.111494   0.042211   2.641  0.00832 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2475 on 1998 degrees of freedom
##   (4045 observations deleted due to missingness)
## Multiple R-squared:  0.003483,   Adjusted R-squared:  0.002486 
## F-statistic: 3.492 on 2 and 1998 DF,  p-value: 0.03063

summary(lm(C150_4_ASIAN~ SEX_SPECIFIC, data= college_scorecard_completely_prepped))

## 
## Call:
## lm(formula = C150_4_ASIAN ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.65546 -0.17798  0.03312  0.23202  0.52488 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.566881   0.007266  78.020   <2e-16 ***
## SEX_SPECIFICMEN ONLY   -0.091765   0.122158  -0.751    0.453    
## SEX_SPECIFICWOMEN ONLY  0.088580   0.054137   1.636    0.102    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2987 on 1724 degrees of freedom
##   (4319 observations deleted due to missingness)
## Multiple R-squared:  0.001888,   Adjusted R-squared:  0.0007303 
## F-statistic: 1.631 on 2 and 1724 DF,  p-value: 0.1961

summary(lm(C150_4_AIAN~ SEX_SPECIFIC, data= college_scorecard_completely_prepped))

## 
## Call:
## lm(formula = C150_4_AIAN ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.43186 -0.43186 -0.04886  0.31814  0.64444 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.43186    0.01007  42.887   <2e-16 ***
## SEX_SPECIFICMEN ONLY    0.17529    0.26375   0.665    0.506    
## SEX_SPECIFICWOMEN ONLY -0.07630    0.09676  -0.789    0.430    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3727 on 1384 degrees of freedom
##   (4659 observations deleted due to missingness)
## Multiple R-squared:  0.0007709,  Adjusted R-squared:  -0.0006731 
## F-statistic: 0.5339 on 2 and 1384 DF,  p-value: 0.5865

Here is the summary of the summary tables:

White Student Completion Rate:

• significant negative relationship for men-only
• slightly significant positive relationship for women-only

Black Student Completion Rate:

• moderately significant positive relationship for women-only

Hispanic Student Completion Rate:

• moderately significant positive relationship for women-only

Asian Student Completion Rate:

• no significant relationship for both men-only and women-only colleges

American Indian/ Alaskan Native Student Completion Rate:

• no significant relationship for both men-only and women-only colleges

Here is what it looks like visually:

college_scorecard_completely_prepped %>% 
  group_by(SEX_SPECIFIC) %>% 
  summarize(avg_completion_w = mean(C150_4_WHITE, na.rm= TRUE), 
            avg_completion_b = mean(C150_4_BLACK, na.rm= TRUE), 
            avg_completion_h = mean(C150_4_HISP, na.rm= TRUE), 
            avg_completion_as = mean(C150_4_ASIAN, na.rm= TRUE), 
            avg_completion_ai = mean(C150_4_AIAN, na.rm= TRUE)) %>%
  plot_ly(x = ~SEX_SPECIFIC, y = ~avg_completion_w, type= 'bar', name= 'White', marker = list(color = "#44668b")) %>%
  add_trace(y = ~avg_completion_b, name= 'Black', marker = list(color = "#f1cb35"))%>%
  add_trace(y = ~avg_completion_h, name= 'Hispanic', marker = list(color = "#1298d0"))%>%
  add_trace(y = ~avg_completion_as, name= 'Asian', marker = list(color = "#cc8834"))%>%
  add_trace(y = ~avg_completion_ai, name= 'American Indian/ Alaska Native', marker = list(color = "#bdb8b8"))%>%
  layout(barmode = "group", 
         xaxis = list(title = "College Type"),
         yaxis = list(title = "Average Completion Rate"),
         font = list(family = "Georgia"))

summary(lm(C150_4_NHPI~ SEX_SPECIFIC, data= college_scorecard_completely_prepped))

## 
## Call:
## lm(formula = C150_4_NHPI ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.61904 -0.44404 -0.04404  0.55596  0.55596 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.44404    0.01384  32.076   <2e-16 ***
## SEX_SPECIFICWOMEN ONLY  0.17501    0.15521   1.128     0.26    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.409 on 878 degrees of freedom
##   (5166 observations deleted due to missingness)
## Multiple R-squared:  0.001446,   Adjusted R-squared:  0.0003085 
## F-statistic: 1.271 on 1 and 878 DF,  p-value: 0.2598

summary(lm(C150_4_2MOR~ SEX_SPECIFIC, data= college_scorecard_completely_prepped))

## 
## Call:
## lm(formula = C150_4_2MOR ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.55539 -0.20622  0.02708  0.19378  0.52708 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            0.472925   0.006913  68.414   <2e-16 ***
## SEX_SPECIFICMEN ONLY   0.227075   0.198793   1.142   0.2535    
## SEX_SPECIFICWOMEN ONLY 0.082466   0.049396   1.669   0.0952 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.281 on 1684 degrees of freedom
##   (4359 observations deleted due to missingness)
## Multiple R-squared:  0.002413,   Adjusted R-squared:  0.001228 
## F-statistic: 2.037 on 2 and 1684 DF,  p-value: 0.1308

summary(lm(C150_4_NRA~ SEX_SPECIFIC, data= college_scorecard_completely_prepped))

## 
## Call:
## lm(formula = C150_4_NRA ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.64515 -0.17095  0.03585  0.21235  0.71789 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.587646   0.007393  79.491  < 2e-16 ***
## SEX_SPECIFICMEN ONLY   -0.305537   0.062997  -4.850 1.36e-06 ***
## SEX_SPECIFICWOMEN ONLY  0.057507   0.054682   1.052    0.293    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2867 on 1550 degrees of freedom
##   (4493 observations deleted due to missingness)
## Multiple R-squared:  0.01575,    Adjusted R-squared:  0.01448 
## F-statistic:  12.4 on 2 and 1550 DF,  p-value: 4.551e-06

summary(lm(C150_4_UNKN~ SEX_SPECIFIC, data= college_scorecard_completely_prepped))

## 
## Call:
## lm(formula = C150_4_UNKN ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.59045 -0.22123  0.00957  0.21547  0.50957 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.490429   0.007379  66.459   <2e-16 ***
## SEX_SPECIFICMEN ONLY   -0.144262   0.171801  -0.840   0.4012    
## SEX_SPECIFICWOMEN ONLY  0.100025   0.058769   1.702   0.0889 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2973 on 1649 degrees of freedom
##   (4394 observations deleted due to missingness)
## Multiple R-squared:  0.002189,   Adjusted R-squared:  0.0009787 
## F-statistic: 1.809 on 2 and 1649 DF,  p-value: 0.1642

Another summary of the summaries:

Native Hawaiian/Pacific Islander Completion Rate:

• no significant relationship for both men-only and women-only colleges

2 or More Races Completion Rate:

• no significant relationship for both men-only and women-only colleges

Nonresident Aliens Completion Rate:

• very significant negative relationship for men-only

Unknown Completion Rate:

• no significant relationship for both men-only and women-only colleges

Here is the visual:

college_scorecard_completely_prepped %>% 
  group_by(SEX_SPECIFIC) %>% 
  summarize(avg_completion_nhpi = mean(C150_4_NHPI, na.rm= TRUE), 
            avg_completion_2M = mean(C150_4_2MOR, na.rm= TRUE), 
            avg_completion_nra = mean(C150_4_NRA, na.rm= TRUE), 
            avg_completion_u = mean(C150_4_UNKN, na.rm= TRUE))%>%
  plot_ly(x = ~SEX_SPECIFIC, y = ~avg_completion_nhpi, type= 'bar', name= 'Native Hawaiian/Pacific Islander', marker = list(color = "#44668b")) %>%
  add_trace(y = ~avg_completion_2M, name= '2 or More Races', marker = list(color = "#1298d0"))%>%
  add_trace(y = ~avg_completion_nra, name= 'Nonresident Aliens', marker = list(color = "#f1cb35"))%>%
  add_trace(y = ~avg_completion_u, name= 'Unknown', marker = list(color = "#cc8834"))%>%
  layout(barmode = "group", 
         xaxis = list(title = "College Type"),
         yaxis = list(title = "Average Completion Rate"),
         font = list(family = "Georgia"))

Retention

Do sex-specific schools affect retention?

Here is the variable we are looking at:

RET_FT4_POOLED (First-time, full-time student retention rate at four-year institutions)

summary(lm(RET_FT4_POOLED~ SEX_SPECIFIC, data= college_scorecard_completely_prepped))

## 
## Call:
## lm(formula = RET_FT4_POOLED ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.72061 -0.07969  0.02494  0.10271  0.27939 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            0.720612   0.003528 204.247   <2e-16 ***
## SEX_SPECIFICMEN ONLY   0.033652   0.022028   1.528    0.127    
## SEX_SPECIFICWOMEN ONLY 0.047545   0.027238   1.746    0.081 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1598 on 2137 degrees of freedom
##   (3906 observations deleted due to missingness)
## Multiple R-squared:  0.002461,   Adjusted R-squared:  0.001528 
## F-statistic: 2.636 on 2 and 2137 DF,  p-value: 0.07186

It looks like there is no significant relationship.

Faculty Employment

Next, we will be investigating the question: “Is there a difference in faculty salary and employment at different colleges?”

Variables:

• AVGFACSAL (Average faculty salary)
• PFTFAC (Proportion of faculty that is full-time)

summary(lm(AVGFACSAL~ SEX_SPECIFIC, data= college_scorecard_completely_prepped))

## 
## Call:
## lm(formula = AVGFACSAL ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6928.9 -1657.4  -289.9  1336.6 13555.1 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             6928.88      40.76 170.008  < 2e-16 ***
## SEX_SPECIFICMEN ONLY   -2217.14     325.38  -6.814 1.09e-11 ***
## SEX_SPECIFICWOMEN ONLY  1020.66     428.12   2.384   0.0172 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2521 on 3920 degrees of freedom
##   (2123 observations deleted due to missingness)
## Multiple R-squared:  0.01322,    Adjusted R-squared:  0.01271 
## F-statistic: 26.25 on 2 and 3920 DF,  p-value: 4.715e-12

summary(lm(PFTFAC~ SEX_SPECIFIC, data= college_scorecard_completely_prepped))

## 
## Call:
## lm(formula = PFTFAC ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.71468 -0.24868 -0.02378  0.28112  0.40862 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            0.591377   0.004878 121.228  < 2e-16 ***
## SEX_SPECIFICMEN ONLY   0.248299   0.050582   4.909 9.57e-07 ***
## SEX_SPECIFICWOMEN ONLY 0.012831   0.049129   0.261    0.794    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2892 on 3580 degrees of freedom
##   (2463 observations deleted due to missingness)
## Multiple R-squared:  0.006699,   Adjusted R-squared:  0.006144 
## F-statistic: 12.07 on 2 and 3580 DF,  p-value: 5.958e-06

The first summary table shows a significant negative relationship between faculty salary and being employed at a men-only institution. In contrast, there is a positive relationship (albeit slightly significant) between faculty salary and being at a women-only college. Additionally, the second table shows a very significant positive relationship between the proportion of full-time faculty and men-only institutions.

faculty_data <- college_scorecard_completely_prepped %>% 
  group_by(SEX_SPECIFIC) %>% 
  summarize(avg_fascal = mean(AVGFACSAL, na.rm= TRUE), 
            avg_pftfac = mean(PFTFAC, na.rm= TRUE)) 

faculty_data %>%
  plot_ly(x = ~SEX_SPECIFIC, y = ~avg_fascal, type= 'bar', marker = list(color = "#1298d0"))%>%
  layout(barmode = "group", 
         xaxis = list(title = "College Type"),
         yaxis = list(title = "Average Faculty Salary (Monthly)"), 
         font = list(family = "Georgia")) 

faculty_data %>%
  plot_ly(x = ~SEX_SPECIFIC, y = ~avg_pftfac, type= 'bar', marker = list(color = "#cc8834"))%>%
  layout(barmode = "group", 
         xaxis = list(title = "College Type"),
         yaxis = list(title = "Average Proportion of Full Time Faculty"),
         font = list(family = "Georgia"))

So men-only colleges pay you less on average but expect you to stick around.

Household Income

How is household income related to sex-specific schools?

Our Variable:

• MEDIAN_HH_INC (Median household income)

summary(hh_income_mod_male <- lm( MEDIAN_HH_INC ~ SEX_SPECIFIC,  data = college_scorecard_completely_prepped))

## 
## Call:
## lm(formula = MEDIAN_HH_INC ~ SEX_SPECIFIC, data = college_scorecard_completely_prepped)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -42219  -8572   -293   8398  42861 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             58009.9      200.2 289.706  < 2e-16 ***
## SEX_SPECIFICMEN ONLY    -4385.3     1807.9  -2.426   0.0153 *  
## SEX_SPECIFICWOMEN ONLY   8879.9     2240.8   3.963 7.52e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13200 on 4434 degrees of freedom
##   (1609 observations deleted due to missingness)
## Multiple R-squared:  0.004888,   Adjusted R-squared:  0.004439 
## F-statistic: 10.89 on 2 and 4434 DF,  p-value: 1.914e-05

college_scorecard_completely_prepped %>% 
  group_by(SEX_SPECIFIC) %>% 
  summarize(avg_median_hh = mean(MEDIAN_HH_INC, na.rm= TRUE)) %>%
  plot_ly(x = ~SEX_SPECIFIC, y = ~avg_median_hh, type= 'bar', marker = list(color = "#1298d0"))%>%
  layout(barmode = "group", 
         xaxis = list(title = "College Type"),
         yaxis = list(title = "Average Median Household Income"),
         font = list(family = "Georgia"))

It looks like there is a significant positive relationship between median household income and attending a women-only college. Conversely, there is a slightly negative relationship between median household income and attending a men-only college.