ASER variance over time

This notebook decomposes the variance of ASER state and district averages (and round to round changes in those averages) into a persistent component and a transitory component. We then further decompose the transitory component into sampling error and a non-sampling component.

We are unable to definitely say whether the non-sampling transitory component is due to non-sampling survey error or one-off shocks to learning outcomes, but there are several reasons to believe that survey error is more likely. In particular:

• Few policies last for only one year
• There are unlikely to be large changes in cohort ability or a big shift from rural to urban

Scroll to the end to see all the results

Decomposing variance into persistent and transitory components

Kaine and Staiger (2002) develop two methods for decomposing variance in the changes in scores into persistent and transitory components. The intuition behind these methods is that if a change in scores from one year to the next is reversed in the following year we would likely suspect that the changes in scores were not due to an actual underlying improvement in learning outcomes. Instead, we would expect that the observed temporary increase was due to measurement error or some other transitory effect (e.g. perhaps there was a heat wave which lowered scores one year).

For the first method they first assume that average test scores for a school in time t, \(y_t\), consist of a fixed component \(\alpha\), a persistent component \(v_t\) which follows a random walk, and a transitory component \(\varepsilon _t\)

\[y_t=\alpha+v_t+\varepsilon_t; v_t=v_{t-1}+u_t\]

Then \(Var(y_t-y_{t-1})= \sigma_u^2+2\sigma_{\varepsilon}^2\) and the proportion of the overall variance of the changes in y arising due to the transitory shock, the \(2\sigma_{\varepsilon}^2\) component, can be estimated as…

\[ -2*corr(\Delta y_t , \Delta y_{t-1} ) = \frac{2\sigma_{\varepsilon}^2}{(\sigma_u^2+2\sigma_{\varepsilon}^2)} \]

Similarly, we can also estimate the proportion of variance in levels (as opposed to changes) which are due to the transitory shock using the fact that…

\[ \sigma^2_{\varepsilon}=-corr(\Delta y_t , \Delta y_{t-1} )*Var(\Delta y_t ) \]

A potential downside to this method is that it relies on the assumption that the terms \(u_t\) are not serially correlated. This might be the case if, for example, a state or district implemented a program which led to sustained improvements in learning outcomes. We can test for this by looking at \(corr(\Delta y_t , \Delta y_{t-2} )\). If the \(u_t\) terms are not serially correlated, the correlation in current changes and twice lagged changes should be 0.

We find that this holds for district changes (the correlation between changes and twice lagged changes is -0.0739763) but not for the state changes (the correlation is -0.2202468). Thus, for states we use Kane and Staiger’s second method which is slightly more complicated. (We also apply the first method as a robustness check.)

The second method relies on the fact that if there is both a persistent component and a transitory component to scores, we would expect the correlation between current scores and the first lagged score to reflect both persistent and transitory shocks while the correlation between current scores and further lags would mainly reflect the transitory component. Thus, if looking at a graph of the correlation by lag number, the correlation should fall quite a bit from 1 with the first lag and then exhibit relatively steady decay after that. We do in fact see this with state scores.

Using this method, we the variance of the persistent component using the following formula:

\[ \sigma^2_{pers}=1-\sigma^2_{\varepsilon}=\frac{\sigma^2_y*\rho_1}{\rho_{pers}}\]

Where \[\rho_{pers}=1/K \sum \frac{\rho_{k+1}}{\rho_k} \]

Decomposing non-persistent variance into sampling and other sources

If we have an estimate of the sampling error for scores, we may further decompose the variance arising from transitory shocks into a component due to sampling error and a component due to other unexplained sources. That is, we assume that \(\varepsilon_t = \varepsilon_{ts} + \varepsilon_{to}\) where the first component is due to sampling and the second is due to other unexplained sources (and the two components are independent) and combine our estimates of \(\sigma^2_s\) and with estimates of \(\sigma_{\varepsilon}^2\) from above.

ASER doesn’t publish standard errors and we don’t have access to the microdata so we are unable to directly estimate the standard errors. There are two ways in which we can indirectly estimate the standard errors. First, this paper by Ramaswami and Wadhwa gives standard errors for a few states and districts. Standard errors for the proportion of class 3-5 students who can read a std 1 text by state, reported in table 1, are all roughly .01. The states for which standard errors are reported are all fairly large states and thus we should not use this value for all states, but it serves as a useful robustness check. Standard errors for districts are around .04. (To arrive at this figure, I used the margin of errors reported for std35_lang in table 2 which are roughly .12. From table 1, I assumed that the prevalence of std35_lang is around .65 and used the formula from the paper for margin of error to back out the standard errors.). So, using this approach, our estimate of \(\sigma_{s}^2\) for state averages is \(.0001=.01^2\) and for district averages is \(.0016=.04^2\). Note that since there are around 16-25 districts per state (for the states which the paper calculates standard errors) we would expect the district variance to be about 1/16th to 1/25th the state variance which is what we see.

Alternatively, we may calculate standard errors for district averages analtically using knowledge of the sampling design and an estimate of the ICC at the village level from IHDS. Using IHDS, we know that the ICC of ASER scores at the village level is around .067. (See the IHDS analysis notebook for this calculation.) Within each district, ASER samples 30 villages and interviews 20 households (or is it kids?) per village. For a variable with prevalence p = .65, the variance of the district estimates from sampling would be…

\[ DEFF = 1+(m-1)\rho \approx 1+(20-1)*.067=2.273\]

\[ \sigma_s^2=\frac{DEFF*p*(1-p)}{600}\approx 9.4708333\times 10^{-4}\] The close similarity between estimates of \(\sigma_{s}^2\) derived analytically and estimates backed out from the Ramaswami and Wadhwa paper gives us confidence.

We take as our final value of \(\sigma_{s}^2\) for districts is .0013. For states, we calculate \(\sigma_{s}^2=.0013/(N_s)\) where \(N_s\) is the number of districts in state s and then take the average across states. We note that the proportion of variance in state levels and changes due to sampling is very small and thus estimating state level sampling variance accurately is less important than for districts.

District Data

We use ASER data from the ASER district pages from years 2006 to 2014. These are the only years for which ASER has publicly released district level data (to my knowledge). There are four learning outcomes variables which all of these district pages include:

• Share class 1 and 2 students who can at least read letters
• Share class 1 and 2 students who can at least recognize numbers
• Share class 3, 4, and 5 students who can at least read level 1 text
• Share class 3, 4, and 5 students who can at least do subtraction

Note that these figures averages across all children, regardless of what type of school they attend.

We use the last two learning outcomes – those for children in class 3, 4, and 5. These learning outcomes should be, in theory, more stable. In addition, they are similar to the learning outcome used in the ASER vs. NAS analysis.

State Data

State level data has been scraped from the individual state “ASER Trends over time” reports. I first used tabula to scrape the pdf and then used an R script (“import ASER trends over time”) to clean up the tabula output. Note that I had to manually change the tabula output in several cases because it didn’t work super well.

There are a lot of different values, but I just used…

• share of class 3 kids who can read a std 1 text
• Share of class 3 kids who can do at least subtraction

These are the only variables available for class 3 kids (though they are also broken up by school type as well) in the trends over time report.

Setup

library(tidyverse)
dist_path <- "C:/Users/dougj/Documents/Data/Education/ASER District Data/Clean"
output <- "C:/Users/dougj/Dropbox/Education in India/Original research/Learning outcomes data/figures"
dists <- read_csv(file.path(dist_path, "aser_district_partial.csv"))

state_path <- "C:/Users/dougj/Documents/Data/Education/ASER trends over time"
states <- read_csv(file.path(state_path, "aser_trends.csv")) %>% select(year, State, std3_reading, std3_math)

##### IMPORTANT: TO DO ALL THE ANALYZE FOR MATH, SET THE SUBJECT BELOW TO "MATH" ### 
subject <- "math"
if (subject == "math") {
  dists$std35_std1_and_up <- dists$std35_subtraction_and_up
  states$std3_reading <- states$std3_math
}


### IMPORTANT -- THIS IS WHERE I SPECIFY DISTRICT SAMPLING VARIANCE ###
var_dist_sampling <- .0013

District analysis

Estimate \(corr(\Delta y_t , \Delta y_{t-1} )\) by using the autocorrelation of district deltas.

# Create variable for deltas and first and second lag of the delta
dist_deltas <- dists %>% 
  group_by(State, District) %>% 
  mutate(delta_reading = std35_std1_and_up - lag(std35_std1_and_up, order_by = year), delta_math = std35_subtraction_and_up - lag(std35_subtraction_and_up)) %>%
  mutate(delta_reading_lagged = lag(delta_reading, order_by = year), delta_reading_dbl_lag = lag(delta_reading, n= 2, order_by = year)) %>%
  arrange(State, District, year)


# Calculate autocorrelation of delta_reading by district and then average over all districts
# Note that since the data is already grouped by State and district, we are calculating 
# autocorrelation for each district and then averaging across all districts
corr_dist_deltas <- dist_deltas %>% 
  filter(year >= 2008) %>% 
  summarise(auto = cor(delta_reading, delta_reading_lagged)) %>%
  ungroup() %>%
  summarise(mean_auto=mean(auto)) %>% .$mean_auto


# Calculate the autocorrelation in one go for all pairs of delta and delta_lagged.
# Note that we are not doing any averaging here.
# This should spit out a very similar value and it does.
check_corr_dist_deltas <- dist_deltas %>% 
  ungroup() %>%
  filter(year >= 2008) %>% 
  summarise(auto = cor(delta_reading, delta_reading_lagged))


# Calculate the overall variance of the deltas
# Note that we have to divide by 100 here to get the actual number.
# We didn't need to do this above because we just calculated correlation and thus scale didn't matter
var_dist_deltas <- var(dist_deltas$delta_reading/100, na.rm = TRUE)

# Calculate overall variance of the district values.
# Note that we don't want to first calculate variance of dist reading scores by state and then average
# That would really underestimate the variance of the district scores.
# (To understand intuition, note districts are considered iid in the model above.)
var_dist_levels <- var(dists$std35_std1_and_up/100, na.rm = TRUE)

# Calculate variance of the transitory component
var_dist_epsilon <- -corr_dist_deltas*var_dist_deltas

Check that \(corr(\Delta y_t , \Delta y_{t-2} )=0\) and that letters and numbers are correlated.

# 
corr_dist_dbl_lag <- dist_deltas %>% filter(year >= 2009) %>%
  ungroup() %>%
  summarize(auto =cor(delta_reading,delta_reading_dbl_lag, use = "pairwise.complete.obs")) %>%
  .$auto


# Check that changes in letters and numbers is correlated -- they are
dist_deltas %>% filter(year != 2006) %>% ungroup() %>%
  summarize(cor(delta_reading,delta_math))

NA

Consolidated district level results

Based on the code above, the final values for the key parameters for the district level are…

• \(corr(\Delta y_t , \Delta y_{t-1} ) = -0.3803595\)
• \(Var(\Delta y_t )= 0.0205213\)
• \(\sigma^2_{\varepsilon} = -corr(\Delta y_t , \Delta y_{t-1} )* Var(\Delta y_t ) = 0.0078055\)
• \(\sigma^2_y = 0.0327374\)

From these figures, we can calculate:

• Share of variance in levels due to transitory shocks = 0.2384259
• Share of variance in changes due to transitory shocks = 0.7607191
• Share of the transitory variance due to sampling (for both changes and levels) = 0.1665502

State analysis - method 1

Estimate \(corr(\Delta y_t , \Delta y_{t-1} )\) at state level by average autocorrelation of state deltas.

# Create variable for year on year change for each district
state_deltas <- states %>% 
  group_by(State) %>% 
  mutate(delta_reading = std3_reading - lag(std3_reading, order_by = year), delta_math = std3_math - lag(std3_math, order_by = year)) %>%
  mutate(delta_reading_lagged = lag(delta_reading, order_by = year), delta_reading_dbl_lag = lag(delta_reading,n=2, order_by = year)) %>%
  arrange(State, year)


# Calculate average autocorrelation of delta_reading
corr_state_deltas <- state_deltas %>% 
  filter(year >= 2008) %>% 
  summarise(auto = cor(delta_reading, delta_reading_lagged, use = "pairwise.complete.obs")) %>%
  ungroup() %>%
  summarise(mean_auto=mean(auto)) %>% .$mean_auto

# Calculate correlation betwen delta letters and delta letters lagged for entire dataset (without first calculating for each state and then averaging).
# this should be about the same as when we first calculate for each state and then average --> which it is.
check_corr_state_deltas <- state_deltas %>% filter(year >= 2008) %>%
  ungroup() %>%
  summarize(auto =cor(delta_reading,delta_reading_lagged, use = "pairwise.complete.obs")) %>%
  .$auto


# calculate the variance of the deltas
var_state_deltas <- var(state_deltas$delta_reading/100, na.rm = TRUE)

# calculate the variance of the levels
var_state_levels <- var(state_deltas$std3_reading/100, na.rm = TRUE)

# calculate variance of transitory component
var_state_epsilon <- -corr_state_deltas*var_state_deltas

Calculate average sampling variance. Note that since we are just calculating the overall average, I don’t have to merge on state names.

var_state_sampling <- dists %>% filter(year ==2008) %>%
  count(State) %>%
  mutate(sampling_var = var_dist_sampling/n) %>%
  summarize(mean_sampling = mean(sampling_var)) %>% .$mean_sampling

  

check that there is no correlation with double lags

# Check that changes in letters and numbers is correlated -- they are highly correlated
state_deltas %>% filter(year >= 2008) %>% ungroup() %>%
  summarize(cor(delta_reading,delta_math, use = "pairwise.complete.obs"))

corr_state_dbl_lag <- state_deltas %>% filter(year >= 2009) %>%
  ungroup() %>%
  summarize(auto =cor(delta_reading,delta_reading_dbl_lag, use = "pairwise.complete.obs")) %>%
  .$auto


temp <- state_deltas %>% filter(year >= 2009) %>%
  group_by(State) %>%
  summarize(auto =cor(delta_reading,delta_reading_dbl_lag, use = "pairwise.complete.obs"))

Consolidated state level results - method 1

Since the correlation of current levels and twice lags is not 0, we don’t use these values, but consolidating them here for reference:

• \(corr(\Delta y_t , \Delta y_{t-1} ) = -0.20926\)
• \(Var(\Delta y_t )= 0.0073836\)
• \(\sigma^2_{\varepsilon} = -corr(\Delta y_t , \Delta y_{t-1} )* Var(\Delta y_t ) = 0.0015451\)
• \(\sigma^2_y = 0.0261895\)
• \(\sigma^2_{sampling} = 1.7844227\times 10^{-4}\)

State analysis - method 2

Calculate the share of variance due to persistent effects by look at how correlation in levels decays over time.

state_lags <- states %>%
  group_by(State) %>%
  mutate(read_lag1 = lag(std3_reading, order_by = year),
         read_lag2 = lag(std3_reading, n= 2, order_by = year),
         read_lag3 = lag(std3_reading, n =3, order_by = year),
         read_lag4 = lag(std3_reading, n =4, order_by = year),
         read_lag5 = lag(std3_reading, n =5, order_by = year)) %>%
  ungroup()

rho = c()
# for each of the lags, calculate the correlation between current and the lag and store in vector rho
for (lag in seq(1,5)) {
  rho <- c(rho, cor(state_lags$std3_reading, state_lags[[paste("read_lag",as.character(lag),sep="")]], use = "pairwise.complete.obs"))
}
# Calculate variance of the transitory component by comparing rho_1 with the decay of rho after that
var_state_levels_persistent2 <- rho[1]/mean(rho/lag(rho), na.rm = TRUE)*var_state_levels

# Calculate the share of variance coming from other transitory sources.
# We will need this to calculate the variance breakup for changes
var_state_levels_other2 <- var_state_levels - var_state_levels_persistent2 - var_state_sampling

# Save a graph showing the decay of the correlation
quick_df <- tibble(lag = seq(0,5), correlation = c(1,rho))
ggplot(quick_df, aes(x= lag, y = correlation)) + 
  geom_line() +
  geom_point() +
  ylim(0,1)

ggsave(paste(subject, "- correlation_decay.png"), width = 5, height = 6 , path = output)

Consolidated state level results - method 2

With method 2, we have the same values for * \(Var(\Delta y_t )\),\(\sigma^2_y\), and \(\sigma^2_{sampling}\).

Instead of calculating \(\sigma^2_{\varepsilon}\), we calculate…

• \(\sigma^2_{pers} = 0.0252023\)

And then use these values to calculate…

• Share of variance in levels due to transitory shocks = 0.0376981
• Share of variance in changes due to transitory shocks = 0.2674285
• Share of the transitory variance due to sampling (for both changes and levels) = 0.1807383

Create bar plot of results

We now have all the inputs we need to replicate figure 4 from Kane and Staiger

# Create empty tibble to store results
# I 
df_bar <- tibble(state_or_dist = rep(c("State","District"),each =6), 
                 changes_or_levels =rep(rep(c("Changes","Levels"), each =3), times= 2),
                 bar_part = factor(rep(c("Persistent","Transitory sampling", "Transitory other"),4), levels = c("Transitory other","Transitory sampling", "Persistent")),
                 value = rep(.001,12)
                 )
# Update the values for the district changes bar
df_bar$value[(df_bar$state_or_dist == "District") & (df_bar$changes_or_levels == "Changes") & (df_bar$bar_part == "Persistent")] <- var_dist_deltas - 2*var_dist_epsilon

df_bar$value[(df_bar$state_or_dist == "District") & (df_bar$changes_or_levels == "Changes") & (df_bar$bar_part == "Transitory sampling")] <- 2*var_dist_sampling

df_bar$value[(df_bar$state_or_dist == "District") & (df_bar$changes_or_levels == "Changes") & (df_bar$bar_part == "Transitory other")] <- 2*var_dist_epsilon - 2*var_dist_sampling


# Update the values for the district levels bar
df_bar$value[(df_bar$state_or_dist == "District") & (df_bar$changes_or_levels == "Levels") & (df_bar$bar_part == "Persistent")] <- var_dist_levels - var_dist_epsilon

df_bar$value[(df_bar$state_or_dist == "District") & (df_bar$changes_or_levels == "Levels") & (df_bar$bar_part == "Transitory sampling")] <- var_dist_sampling

df_bar$value[(df_bar$state_or_dist == "District") & (df_bar$changes_or_levels == "Levels") & (df_bar$bar_part == "Transitory other")] <- var_dist_epsilon - var_dist_sampling

##### STATE DATA #####

# CHANGES
df_bar$value[(df_bar$state_or_dist == "State") & (df_bar$changes_or_levels == "Changes") & (df_bar$bar_part == "Persistent")] <- var_state_deltas-2*(var_state_sampling+var_state_levels_other2)

df_bar$value[(df_bar$state_or_dist == "State") & (df_bar$changes_or_levels == "Changes") & (df_bar$bar_part == "Transitory sampling")] <- 2*var_state_sampling

df_bar$value[(df_bar$state_or_dist == "State") & (df_bar$changes_or_levels == "Changes") & (df_bar$bar_part == "Transitory other")] <- 2*var_state_levels_other2


# LEVELS
df_bar$value[(df_bar$state_or_dist == "State") & (df_bar$changes_or_levels == "Levels") & (df_bar$bar_part == "Persistent")] <- var_state_levels_persistent2

df_bar$value[(df_bar$state_or_dist == "State") & (df_bar$changes_or_levels == "Levels") & (df_bar$bar_part == "Transitory sampling")] <- var_state_sampling

df_bar$value[(df_bar$state_or_dist == "State") & (df_bar$changes_or_levels == "Levels") & (df_bar$bar_part == "Transitory other")] <- var_state_levels_other2


# I have created this graph so that you can see it in the notebook
# The final version which combines both subjects, is created in a separate notebook
ggplot(df_bar, aes(fill = bar_part, y=value, x= changes_or_levels)) +
  geom_bar(position="stack", stat ="identity") +
  facet_grid(~ state_or_dist) + 
  scale_fill_manual(values = c("red", "orange", "blue"))+
  labs(fill = "Variance component", x = "")

ggsave(paste(subject," - variance_decomposition.png"), width = 5, height = 6 , path = output)

write_csv(df_bar, file.path(output, paste(subject, "- bar data.csv")))

Create dataframe with results

Take the output that was used to create the bar graph, reshape it, and save it as a csv

df_final <- df_bar %>%
  pivot_wider(names_from = "bar_part", values_from= "value") %>%
  mutate(total_var = `Persistent`+`Transitory sampling`+`Transitory other`) %>%
  mutate(share_pers = `Persistent`/total_var, share_samp = `Transitory sampling`/total_var, share_other = `Transitory other`/total_var)

write_csv(df_final, file.path(output, paste(subject, "- var deco.csv")))