In two previous posts I put together a data set and presented some exploratory data analysis on school achievement for national standards. After those posts I exchanged emails with a few people about the sources of data and Jeremy Greenbrook-Held pointed out Education Counts as a good source of additional variables, including number of teachers per school and proportions for different ethnic groups.

The code below call three files: Directory-Schools-Current.csv, teacher-numbers.csv and SchoolReport_data_distributable.csv, which you can download from the links.

options(stringsAsFactors = FALSE)

library(ggplot2)

# Reading School directory and dropping some address information
# for 2012 obtained from here:
# http://www.educationcounts.govt.nz/directories/list-of-nz-schools
directory <- read.csv('Directory-Schools-Current.csv')
directory <- subset(directory, select = -1*c(3:9, 11:13))

# Reading teacher numbers for 2011 obtained from here:
# http://www.educationcounts.govt.nz/statistics/schooling/teaching_staff
teachers <- read.csv('teacher-numbers.csv')

# Reading file obtained from stuff.co.nz obtained from here:
# http://schoolreport.stuff.co.nz/index.html
fairfax <- read.csv('SchoolReport_data_distributable.csv')

# Merging directory and teachers info
standards <- merge(directory, teachers, by = 'school.id', all.x = TRUE)

# Checking that the school names match
# This shows ten cases, some of them obviously the same school
# e.g. Epsom Girls Grammar School vs Epsom Girls' Grammar School
tocheck <- subset(standards, school.name.x != school.name.y)
tocheck[, c(1:2, 29)]
#school.id                          school.name.x                  school.name.y
#61          64             Epsom Girls Grammar School    Epsom Girls' Grammar School
#125        135                     Fraser High School  Hamilton's Fraser High School
#362        402                      Waiau Area School    Tuatapere Community College
#442        559                        Te Wainui a Rua           Whanganui Awa School
#920       1506 St Michael's Catholic School (Remuera)  St Michael's School (Remuera)
#929       1515                      Sunnyhills School             Sunny Hills School
#985       1573                     Willow Park School              Willowpark School
#1212      1865               Te Wharekura o Maniapoto         Te Wharekura o Oparure
#1926      2985          Sacred Heart Cathedral School Sacred Heart School (Thorndon)
#2266      3554                         Waitaha School        Waitaha Learning Centre

# Merging now with fairfax data
standards <- merge(standards, fairfax,
                   by = 'school.id', all.x = TRUE)

# Four schools have a different name
tocheck2 <- subset(standards, school.name.x != school.name)
tocheck2[, c(1:2, 32)]
#school.id                          school.name.x                   school.name
#920       1506 St Michael's Catholic School (Remuera) St Michael's School (Remuera)
#929       1515                      Sunnyhills School            Sunny Hills School
#985       1573                     Willow Park School             Willowpark School
#2266      3554                         Waitaha School       Waitaha Learning Centre

# Checking the original spreadsheets it seems that, despite
# the different name they are the same school (at least same
# city)

# Now start dropping a few observations that are of no use
# for any analysis; e.g. schools without data
standards <- subset(standards, !is.na(reading.WB))

# Dropping school 498 Te Aho o Te Kura Pounamu Wellington
# http://www.tekura.school.nz/
# It is a correspondence school without decile information
standards <- subset(standards, school.id != 498)

# Converting variables to factors
standards$decile <- factor(standards$decile)
standards$school.type <- factor(standards$school.type)
standards$education.region <- factor(standards$education.region,
                                     levels = c('Northern', 'Central North', 'Central South', 'Southern'))

# Removing 8 special schools
standards <- subset(standards, as.character(school.type) != 'Special School')
standards$school.type <- standards$school.type[drop = TRUE]

# Saving file for Eric
# write.csv(standards, 'standardsUpdated.csv',
#          row.names = FALSE, quote = FALSE)

# Create performance groups
standards$math.OK <- with(standards, math.At + math.A)
standards$reading.OK <- with(standards, reading.At + reading.A)
standards$writing.OK <- with(standards, writing.At + writing.A)

# Creating a few proportions
standards$prop.euro <- with(standards, european/total.roll)
standards$prop.maori <- with(standards, maori/total.roll)
standards$prop.pacific <- with(standards, pacific.island/total.roll)
standards$prop.asian <- with(standards, asian/total.roll)

This updated data set is more comprehensive but it doesn't change the general picture presented in my previous post beyond the headlines. Now we can get some cool graphs to point out the obvious, for example the large proportion of Maori and Pacific Island students in low decile schools:

qplot(prop.maori, prop.pacific,
      data = standards, color = decile,
      xlab = 'Proportion of Maori',
      ylab = 'Proportion of Pacific Island')

Proportion of Pacific Island (vertical axis) and Maori students (horizontal axis) in schools with points colored by decile. Higher proportions for both are observed in low decile schools.

I have avoided 'proper' statistical modeling because i- there is substantial uncertainty in the data and ii- the national standards for all schools (as opposed to only 1,000 schools) will be released soon; we do't know if the published data are a random sample. In any case, a quick linear model fitting the proportion of students that meet reading standards (reading.OK) as a function of decile and weighted by total school roll—to account for the varying school sizes—will explain roughly 45% of the observed variability on reading achievement.

m1 <- lm(reading.OK ~ decile,
         data = standards,
         weights = total.roll)
summary(m1)

# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# (Intercept)  0.56164    0.01240  45.278  < 2e-16 ***
#   decile2      0.06238    0.01769   3.527 0.000439 ***
#   decile3      0.10890    0.01749   6.225 7.07e-10 ***
#   decile4      0.16434    0.01681   9.774  < 2e-16 ***
#   decile5      0.18579    0.01588  11.697  < 2e-16 ***
#   decile6      0.22849    0.01557  14.675  < 2e-16 ***
#   decile7      0.25138    0.01578  15.931  < 2e-16 ***
#   decile8      0.24849    0.01488  16.701  < 2e-16 ***
#   decile9      0.28861    0.01492  19.347  < 2e-16 ***
#   decile10     0.31074    0.01439  21.597  < 2e-16 ***
#
# Residual standard error: 1.594 on 999 degrees of freedom
# (1 observation deleted due to missingness)
# Multiple R-squared: 0.4548,  Adjusted R-squared: 0.4499
# F-statistic: 92.58 on 9 and 999 DF,  p-value: < 2.2e-16

Model fit has a few issues with distribution of residuals, we should probably use a power transformation for the response variable, but I wouldn't spend much more time before getting the full data for national standards.

par(mfrow = c(2, 2))
plot(m1)
par(mfrow = c(1, 1))

Residuals of a quick weighted linear model. The residuals show some heterogeneity of variance (top-left) and deviation from normality (top-right) with heavy tails.

Bonus plot: map of New Zealand based on school locations, colors depicting proportion of students meeting reading national standards.

qplot(longitude, latitude,
      data = standards, color = reading.OK)

New Zealand drawn using school locations; color coding is for proportion of students meeting reading national standards.

P.S. 2012-09-26 16:01. The simple model above could be fitted taking into account the order of the decile factor (using ordered()) or just fitting linear and quadratic terms for a numeric expression of decile. Anyway, that would account for 45% of the observed variability.

m1d <- lm(reading.OK ~ as.numeric(decile) + I(as.numeric(decile)^2),
          data = standards, weights = total.roll)
summary(m1d)

#Coefficients:
#                          Estimate Std. Error t value Pr(>|t|)
#(Intercept)              0.5154954  0.0136715  37.706  < 2e-16 ***
#as.numeric(decile)       0.0583659  0.0051233  11.392  < 2e-16 ***
#I(as.numeric(decile)^2) -0.0023453  0.0004251  -5.517 4.39e-08 ***
#---
#Residual standard error: 1.598 on 1006 degrees of freedom
#  (1 observation deleted due to missingness)
#Multiple R-squared: 0.4483,	Adjusted R-squared: 0.4472
#F-statistic: 408.8 on 2 and 1006 DF,  p-value: < 2.2e-16

P.S. 2012-09-26 18:17. Eric Crampton has posted preliminary analyses based on this dataset in Offsetting Behaviour.