R Notebook

rm(list = ls())

library(tidyverse)

Warning: package ‘tidyverse’ was built under R version 4.1.3Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages ----------------------------------------------------------------------------------------- tidyverse 1.3.2 --v ggplot2 3.4.0      v purrr   0.3.4 
v tibble  3.1.8      v dplyr   1.0.10
v tidyr   1.2.1      v stringr 1.4.0 
v readr   2.1.3      v forcats 0.5.1 Warning: package ‘ggplot2’ was built under R version 4.1.3Warning: package ‘tibble’ was built under R version 4.1.3Warning: package ‘tidyr’ was built under R version 4.1.3Warning: package ‘readr’ was built under R version 4.1.3Warning: package ‘dplyr’ was built under R version 4.1.3-- Conflicts -------------------------------------------------------------------------------------------- tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(ggplot2)

districts = read_csv("./data/districts.csv")

Rows: 146 Columns: 27-- Column specification ----------------------------------------------------------------------------------------------------------
Delimiter: ","
chr  (2): system_name, region
dbl (25): system, alg_1, alg_2, bio, chem, ela, eng_1, eng_2, eng_3, math, science, enrollment, black, hispanic, native, el, s...
i Use `spec()` to retrieve the full column specification for this data.
i Specify the column types or set `show_col_types = FALSE` to quiet this message.

districts

There are 146 rows and 27 columns in districts.

Q2. Notice that the first row corresponds to the whole State of Tennessee. Remove this row and save the result back to districts.

districts = districts[-1, ]
head(districts)

Q3. How many districts have a proficiency rate of at least 80% for both alg_1 and eng_1?

districts%>%filter(alg_1 >= 80 & eng_1 >= 80)

count(districts%>%filter(alg_1 >= 80 & eng_1 >= 80))

There are 13 districts.

Q4. How many districts have a proficiency rate less than 50% for either alg_1 or eng_1?

districts%>%filter(alg_1 < 50  | eng_1 < 50)

count(districts%>%filter(alg_1 < 50  | eng_1 < 50))

There are 8 districts.

Q5.Which district has the lowest graduation rate?

districts%>%filter(grad == min(districts$grad, na.rm = TRUE))%>%
  select(system_name, grad)

Tennessee School for Blind has the lowest graduation rate.

#oa
districts%>%
  arrange(grad) %>% 
  head(1)

Q6. Within the Mid Cumberland region, which district has the highest ACT composite?

districts%>%filter(region == 'Mid Cumberland')%>%
  filter(act_composite == max(act_composite, na.rm = TRUE))%>%
  select(system_name, act_composite)

Q7. Create a histogram showing the distribution of graduation rates. What can you say about this distribution?

ggplot(districts, aes(grad)) +
  geom_histogram(binwidth = 5, fill = 'blue', color = 'black')+
  scale_x_continuous(breaks=seq(0, 100, 10))

NA

districts %>% 
  ggplot(aes(x=grad, fill = region)) +
  geom_histogram(bins = 50, na.rm = TRUE)

Most graduation rates fall between 80-100%. There is one outlier with a graduation rate 11%.

Q8. Create a scatter plot to compare alg_1 proficiency rates to alg_2 rates. What do you notice?

Facet this plot by region. Does anything stand out when you facet the plots?

ggplot(districts) +
  geom_point(aes(alg_1, alg_2)) +
  scale_x_continuous(breaks=seq(0, 100, 10))

The scatter plot shows there is a strong correlation between alg_1 and alg_2 proficiency rates.

ggplot(districts) + 
  geom_point(aes(alg_1, alg_2, color = factor(region))) +
  facet_wrap((.~region), nrow = 3, ncol = 3)+
  xlim(0,100)+
  ylim(0,100)

The scatter plots show Southwest/Memphis shows the strongest positive correlation.Southeast shows the weakest correlation. Upper Cumberland shows a non-linear trend.

Q9. Create a bar chart showing the total enrollment by region.

Which region has the highest total enrollment? Which has the smallest?

options(scipen = 999)
p = ggplot(districts, aes(x=region, y=enrollment, fill = region)) +
  geom_bar(stat = "identity") +
  scale_fill_brewer(palette = "Dark2")

p + coord_flip() 

NA
NA

Reordering by enrollment

library(scales)

Warning: package ‘scales’ was built under R version 4.1.3
Attaching package: ‘scales’

The following object is masked from ‘package:purrr’:

    discard

The following object is masked from ‘package:readr’:

    col_factor

options(scipen = 999)
p = ggplot(districts, aes(x = reorder(region, +enrollment), y=enrollment, fill = region)) +
  geom_bar(stat = "identity", na.rm = T) +
  scale_fill_brewer(palette = "Dark2")+
  scale_y_continuous(labels = comma)

p + coord_flip() 

upper cumberland and southwest memphis do not get in order. Mid Cumberland has the highest enrollment, Northwest has the lowest enrollment.

Q10. When creating this bar chart you may have noticed that some districts have missing enrollment values.

For how many districts is this the case?

sum(is.na(districts$enrollment))

[1] 4

districts%>%filter(is.na(enrollment))%>%group_by(region, system_name)%>%tally()

Q11. What is the mean graduation rate across all districts?

What might be wrong with using just the regular mean to assess average graduation rates?

mean(districts$grad, na.rm = TRUE)

[1] 90.06562

Q12. Redo the previous question but use a weighted average (weighted.mean) graduation across all districts, weighing by enrollment.

How much does this change your answer?

Can you explain using the data the reason for the big change from using the mean?

weighted.mean(districts$grad, dplyr::coalesce(districts$enrollment,0), na.rm=TRUE)

[1] 87.33285

weighted.mean(districts$grad, coalesce(districts$enrollment,0), na.rm=TRUE)

[1] 87.33285

districts%>%
  summarize(grad_weighted_mean = weighted.mean(grad, coalesce(enrollment,0), na.rm=TRUE))

df <- districts %>% filter((!is.na(enrollment))&(!is.na(grad)))

weighted.mean(df$grad, df$enrollment)

[1] 87.33285

The unweighted graduation rate is 90%. The enrollment weighted graduation rate is 87%. Enrollment weighted – schools with larger enrollments influence the average more. The weighted average decreased–schools with larger enrollments have smaller graduation rates. Weighted-higher enrollments contribute more to the avg.

Q13. Create a boxplot showing graduation rates per region.

Does anything stand out?

ggplot(districts)+
  geom_boxplot(aes(x=region, y=grad, fill = region))+
  ggtitle('Graduation Rate by Region')+
  scale_y_continuous(breaks = seq(0,100,5))+
  coord_flip()

The boxplot shows outliers for Southwest/Memphis at ~50% and MidCumberland at ~10%

Q14. Find the weighted average of graduation rates by region using enrollment as weights.

Compare the results you get for the weighted average to what you see from the boxplots.

options(digits = 3)
districts %>% 
  group_by(region)%>%
  summarise(enr_wm = weighted.mean(grad, coalesce(enrollment,0), na.rm=T),
            mean = mean(grad, na.rm =T))

NA
NA
NA
NA

Southwest/Memphis region has a lower graduation rate of 79.5, which is consistent with its boxplot showing a small outlier and skewed left distribution. Midcumberland’s graduation rate is not that much lower than other regions even though its boxplot shows a small outlier around ~12%.

Q15. For many districts, values for alg_2 are lower than for alg_1.

districts <- districts %>% mutate(alg_diff = alg_1 - alg_2)

Create a histogram showing the distribution of differences (alg_1 - alg_2).

ggplot() + 
  geom_histogram(aes(x=alg_diff, fill="a"), alpha=.7, data=subset(districts, alg_diff <=0),  color = "blue", binwidth = 5) +
  geom_histogram(aes(x=alg_diff, fill="z"), alpha=.5, data=subset(districts, alg_diff >0),  color="black", binwidth = 5) +
  
  scale_fill_manual(name="", values=c("a" = "skyblue", "z" ="red"), labels=c("a"="Alg1 < Alg2", "z" ="Alg1 > Alg2"))+
  
  scale_x_continuous(breaks=seq(-20, 60, 10))

NA
NA

The histogram shows about 20 districts had a drop in algebra by about 15 points.

Which school had the largest drop from alg_1 to alg_2?

districts%>%filter(alg_diff == max(alg_diff, na.rm=TRUE))%>%select(system_name, alg_1, alg_2, alg_diff)

Pickett County has the largest drop from alg_1 to alg_2.

For what percentage of schools is it true that alg_2 is larger than alg_1?

districts%>%summarise(pct_sch = sum(alg_diff<0, na.rm = TRUE)/nrow(districts))

districts%>%summarise(pct_sch = sum(alg_diff<0, na.rm = TRUE)/ sum(!is.na(alg_diff)))

districts%>%summarise(pct_sch = 100*sum(na.omit(alg_diff<0))/145)

15.4% of schools have alg_2 greater than alg_1.

Is there a similar drop off for eng_2 and eng_3?

districts <- districts %>% mutate(eng23_diff = eng_2 - eng_3)

ggplot() + 
  geom_histogram(aes(x=eng23_diff, fill="a"), alpha=.9, data=subset(districts, eng23_diff <0),  color = "blue", binwidth = 5) +
  geom_histogram(aes(x=eng23_diff, fill="z"), alpha=.5, data=subset(districts, eng23_diff >0),  color="black", binwidth = 5) +
  
  scale_fill_manual(name="", values=c("a" = "skyblue", "z" ="red"), labels=c("a"="Eng2 < Eng3", "z" ="Eng2 > Eng3"))

  
  #scale_x_continuous(breaks=seq(-20, 60, 10))
  

The histogram shows almost all schools had drops in scores from Eng2 to Eng3.

districts%>%summarise(pct_sch = 100*sum(na.omit(eng23_diff>0))/145)

86% of schools showed decrease in score from eng2 to eng3.

16. You may have noticed that a lot of rows are missing values.

Which district has the largest number of missing values?

What do you notice about schools that have a lot of missing values?

districts %>% mutate(na_ct = rowSums(is.na(.)))%>%filter(na_ct == max(na_ct))

West Tenn School for Deaf has the most missing values

districts %>% mutate(na_ct = rowSums(is.na(.)))%>%arrange(desc(na_ct))

The missing values seem to come from the subject proficiency rate and act_composite columns. Maybe the schools did not participate in the state tests. They also have missing values in grad and dropout columns

###Q17. Find the correlation between graduation rate and all other variables. ### Create a horizontal bar chart showing these correlations. ### Make sure that your plot is ordered by correlation values. ### What do you notice from these correlations?

districts %>% 
  select(where(is.numeric)) %>% 
  cor(. , districts$grad, use = 'complete.obs')

                 [,1]
system        -0.0781
alg_1          0.3312
alg_2          0.3620
bio            0.3740
chem           0.2348
ela            0.2710
eng_1          0.3052
eng_2          0.3042
eng_3          0.4108
math           0.2990
science        0.3961
enrollment    -0.4164
black         -0.4199
hispanic      -0.2958
native        -0.0534
el            -0.3423
swd            0.1714
ed            -0.2130
expenditures  -0.2167
act_composite  0.2004
chronic_abs   -0.2434
suspended     -0.5017
expelled      -0.3954
grad           1.0000
dropout       -0.9218
alg_diff      -0.1934
eng23_diff    -0.2223

cor_dt <- districts %>% 
  select(where(is.numeric)) %>% 
  cor(. , districts$grad, use = 'complete.obs') %>% 
  as_tibble(cor_dt, rownames = "variable") %>% 
  rename(correlation = V1)%>%
  filter(variable != "grad")

Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
Using compatibility `.name_repair`.

cor_dt

ggplot(data = cor_dt, aes(y=reorder(variable, -correlation), x=correlation)) +
  geom_bar(stat = 'identity')+
  theme_bw()

NA
NA
NA

ggplot(data = cor_dt, aes(x=correlation, y=reorder(variable, -correlation), fill=correlation)) +
  geom_bar(stat = 'identity')+
  scale_fill_gradient(low='firebrick', high='skyblue', space='Lab')

NA
NA
NA
NA

Dropout and suspension rates are the most negatively correlated with graduation rate. Scores in English, Algebra and Sciences are the most positively correlated with graduation rates.

###Q18. Create a scatterplot for grad vs. suspended. ### Does what you see make sense given your answer from the previous part?

ggplot(districts) +
  geom_point(aes(suspended, grad),  na.rm = TRUE)+
  xlim(0,100)+
  ylim(0,100)

NA

The scatter plot shows that lower suspension rates are associated with higher graduation rates. But for lower suspension rates (<5%), there does not seem to be a strong variability in graduation rate. The negative correlation may be as strong as it is because of a few outliers.

###Q19. Create a linear regression model using lm with target variable grad and predictor variable suspended. ### What R^2 value does this model have? What is the interpretation of this number?

simple_regression_model <-lm(grad ~ suspended, data = districts)

summary(simple_regression_model)

Call:
lm(formula = grad ~ suspended, data = districts)

Residuals:
   Min     1Q Median     3Q    Max 
-81.41  -2.27   1.19   3.93  12.21 

Coefficients:
            Estimate Std. Error t value            Pr(>|t|)    
(Intercept)   92.507      1.092    84.7 <0.0000000000000002 ***
suspended     -0.786      0.238    -3.3              0.0013 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.09 on 126 degrees of freedom
  (17 observations deleted due to missingness)
Multiple R-squared:  0.0795,    Adjusted R-squared:  0.0722 
F-statistic: 10.9 on 1 and 126 DF,  p-value: 0.00126

The R^2 of the model is 0.0795 or 7.95%. 7.95% of the variability in the graduation rate can be explained by suspension rate.

grad_lm <- districts %>% 
  filter(grad > 25) %>% 
  lm(grad ~ suspended, data = .)
summary(grad_lm)

Call:
lm(formula = grad ~ suspended, data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.815  -3.011   0.796   3.116  13.306 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept)   93.700      0.658  142.49 < 0.0000000000000002 ***
suspended     -0.962      0.143   -6.74        0.00000000053 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.43 on 125 degrees of freedom
Multiple R-squared:  0.266, Adjusted R-squared:  0.26 
F-statistic: 45.4 on 1 and 125 DF,  p-value: 0.000000000532

Q20. Add the regression line to your scatterplot using geom_smooth with method=‘lm’.

How do you feel about the regression line after seeing it plotted on the scatterplot?

ggplot(districts, aes(x=suspended, y=grad)) +
  geom_point()+
  geom_smooth(method = "lm")

NA
NA

The regression line does not capture most points well, there is a lot of scatter for suspended rate below 5%. This is consistent with the low R^2 value for the model.

Continued Exploration and Practice

Q21. Read in the school-level testing data for 2014, available here.

You might find the readxl library useful for this task. If you use this library, be sure to look at the na argument for the read_excel function.

library(readxl)
library(tidyverse)

testing <- read_excel("./data/data_2014_school_base.xlsx", na = c("*", "**"))

Q22. How many schools have at least 20 percent of students below bsc for Algebra I?

colnames(testing)

 [1] "year"              "system"            "system_name"       "school"            "school_name"       "subject"          
 [7] "grade"             "subgroup"          "valid_tests"       "n_below_bsc"       "n_bsc"             "n_prof"           
[13] "n_adv"             "pct_below_bsc"     "pct_bsc"           "pct_prof"          "pct_adv"           "pct_bsc_and_below"
[19] "pct_prof_adv"     

count(testing %>% filter((subject == "Algebra I")&(pct_below_bsc >= 20)))

NA

Which districts do these schools belong to?

testing %>% 
  filter((subject == "Algebra I")&(pct_below_bsc >= 20)) %>% 
  distinct(system_name)

There are 102 districts that have at least 20% below bsc for Algebra I.

count(testing %>% distinct(system_name))

There are 141 total districts in the dataset.

####Q23. How many schools have at least 20 percent of students below bsc for both Algebra I and English I?

fil_tst <- testing %>%
  filter(pct_below_bsc >= 20) %>% 
  filter(subject == "Algebra I" | subject == "English I") %>%
  select(school_name, subject, pct_below_bsc) 

fil_tst

fil_tst %>% 
  group_by(school_name) %>% 
  summarise(count = n_distinct(subject)) %>% 
  filter(count > 1)

There are 51 schools that have at least 20% below bsc in both English I and Algebra I.

###Q24. Which grade has the highest pct_adv for Algebra I?

testing %>%  
  filter(subject == "Algebra I") %>% 
  filter(pct_adv == max(pct_adv, na.rm = T)) %>% 
  group_by(grade)

8th grade has the highest Algebra I grade.

Plot the average pct_adv per grade level as a bar chart.

testing %>% 
  group_by(grade) %>%
  summarise(avg_pct_adv = mean(pct_adv, na.rm = T)) %>%
  arrange(as.numeric(grade))

Warning: NAs introduced by coercion

NA

testing %>% 
  group_by(grade) %>%
  summarise(avg_pct_adv = mean(pct_adv, na.rm = T)) %>% 
  
  ggplot(aes(x = factor(as.numeric(grade)), y=avg_pct_adv))+
  geom_bar(stat = "identity", fill = 'skyblue3', color = 'black')+
  xlab("Grade")+
  ylab("Advanced Percentage")+
  theme_minimal()

NA
NA

#Pct_adv only for Algebra I
testing %>%
  filter(subject=="Algebra I") %>% 
  group_by(grade) %>%
  summarise(avg_pct_adv = mean(pct_adv, na.rm = T)) %>% 
  
  ggplot(aes(x = factor(as.numeric(grade)), y=avg_pct_adv))+
  geom_bar(stat = "identity", fill = 'steelblue', color = 'black')+
  xlab("Grade")+
  ylab("Advanced Percentage for Algebra I")

NA
NA

###Q25. Find the correlation between pct_adv for Algebra I and pct_adv for Algebra II by school. #### Create a scatterplot showing Algebra II scores vs. Algebra I scores by school.

#Filter dataset for AlgI/AlgII
algebra <- testing %>% 
  filter(subject == "Algebra I" | subject == "Algebra II") %>% 
  select(school_name, subject, pct_adv)

#Pivot algebra
alg_pvt <- algebra %>%  
   
  
  pivot_wider(names_from = subject, values_from = pct_adv, values_fn = function(x) mean(x, na.rm = TRUE))

alg_pvt <- alg_pvt %>% rename('Algebra_I' = 'Algebra I',
                   'Algebra_II' = 'Algebra II')

alg_pvt

ggplot(alg_pvt)+
  geom_point(aes(Algebra_I, Algebra_II))+
  xlim(0,100)+
  ylim(0,100)

alg_pvt <- alg_pvt %>%transform(Algebra_I = as.numeric(Algebra_I), Algebra_II = as.numeric(Algebra_II))

cor(alg_pvt$Algebra_I, alg_pvt$Algebra_II, use = 'complete.obs')

[1] 0.462

Q19. Find all schools in Rutherford County that have “High School” in their name.

For these schools, create a chart (your choice) showing the differences in

pct_below_bsc, pct_bsc, pct_prof, and pct_adv for Algebra I when looking across all subgroups and grades.

testing %>% 
  filter(system_name == 'Rutherford County') %>% 
  filter(str_detect(school_name, "High School"))%>%
  filter(subject == 'Algebra I') %>% 
  
  select(school_name)%>%
  distinct(school_name)

rcs <- testing %>% 
  filter(system_name == 'Rutherford County') %>% 
  filter(str_detect(school_name, "High School"))%>%
  filter(subject == 'Algebra I') %>% 
  select(school_name, grade, subgroup, pct_below_bsc, pct_bsc, pct_prof, pct_adv, pct_bsc_and_below, pct_prof_adv)

pct_piv <- rcs %>% pivot_longer(cols = starts_with("pct"),
                     names_to = "themeasure",
                     values_to = "thenumber")

pct_piv %>% 
  filter(!themeasure %in% c("pct_bsc_and_below", "pct_prof_adv")) %>%
  group_by(school_name, themeasure) %>% 
  
  summarize(pctg = mean(thenumber, na.rm = TRUE))

`summarise()` has grouped output by 'school_name'. You can override using the `.groups` argument.

rcs_pctgs <- pct_piv %>% 
  filter(!themeasure %in% c("pct_bsc_and_below", "pct_prof_adv")) %>%
  group_by(school_name, themeasure) %>% 
  
  summarize(pctg = mean(thenumber, na.rm = TRUE))

`summarise()` has grouped output by 'school_name'. You can override using the `.groups` argument.

rcs_pctgs %>% ggplot(aes(x=themeasure, y=pctg, fill = themeasure))+
  geom_bar(position = "dodge", stat = "identity")+
  
  facet_wrap(~school_name)

rcs_pctgs$themeasure <- factor(rcs_pctgs$themeasure, levels = c("pct_adv", "pct_prof", "pct_bsc","pct_below_bsc"))

ggplot(rcs_pctgs, aes(x=pctg, y= school_name, fill=themeasure, label = pctg))+
  geom_bar(position = "stack", stat = "identity")+
  xlab("Percentage")+
  ylab("School")+
  geom_text(size=3, aes(label=paste0(sprintf("%1.0f", pctg),"%")), position = position_stack(vjust = 0.5))+
  theme(axis.text.x=element_text(angle = 45, hjust = 1))+
  coord_flip()+
  

  scale_fill_manual(values = c("pct_below_bsc" = "#d7191c",
                                "pct_bsc" = "#fdae61",
                                "pct_prof" = "#abd9e9",
                                "pct_adv" = "#2c7bb6"))