timss.Rmd

# Trends in International Mathematics and Science Study (TIMSS) {-}

[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) <a href="https://github.com/asdfree/timss/actions"><img src="https://github.com/asdfree/timss/actions/workflows/r.yml/badge.svg" alt="Github Actions Badge"></a> 

A comparative study of student achievement in math and science across more than 50 nations.

* Grade-specific tables with one record per school, student, teacher, plus files containing student achievement, home background, student-teacher linkage, and within-country scoring reliability.

* A complex survey generalizing to fourth- and eighth-grade populations of participating countries.

* Released quadrennially since 1995.

* Funded by the [International Association for the Evaluation of Educational Achievement](https://www.iea.nl/), run at [BC](http://www.bc.edu/bc-web/schools/lsoe.html).

---

## Recommended Reading {-}

Four Example Strengths & Limitations:

✔️ [Sixty-five education systems, six benchmarking systems participated in 4th or 8th grade assessment](https://timss2023.org/encyclopedia/)

✔️ [Collects rich array of data related to teacher characteristics](https://doi.org/10.1186/s40536-024-00214-x)

❌ [Differences in sample selectivity across countries potentially undermines validity of rankings](http://doi.org/10.3386/w15949)

❌ [Low stakes examination potentially biased by non-serious test-takers](http://doi.org/10.3386/w24930)

<br>

Three Example Findings:

1. [New Zealand Year 9 student chemistry performance improved between 2019 and 2023](https://www.educationcounts.govt.nz/__data/assets/pdf_file/0005/249152/TIMSS-2023-Science-Year-9.pdf).

2. [Students who were food-insecure had lower math achievement than their peers in 2019](https://doi.org/10.1186/s40536-023-00161-z).

3. [In mathematics, U.S. 4th- and 8th-graders scored lower, on average, in 2023 than they did in 2019](https://nces.ed.gov/timss/results23/index.asp#/math/intlcompare).

<br>

Two Methodology Documents:

> [TIMSS 2019 User Guide for the International Database, 2nd Edition](https://timss2019.org/international-database/downloads/TIMSS-2019-User-Guide-for-the-International-Database-2nd-Ed.pdf)

> [Methods and Procedures: TIMSS 2019 Technical Report](https://timssandpirls.bc.edu/timss2019/methods/pdf/TIMSS-2019-MP-Technical-Report.pdf)

<br>

One Haiku:

```{r}
# brando for stella,
# gump's jenny, rock's adrian,
# students toward math test
```

---

## Function Definitions {-}

This survey uses a multiply-imputed variance estimation technique described in [Methods Chapter 14](https://timssandpirls.bc.edu/timss2019/methods/chapter-14.html). Most users do not need to study this function carefully. Define a function specific to only this dataset:

```{r eval = FALSE , results = "hide" }
timss_MIcombine <-
	function (results, variances, call = sys.call(), df.complete = Inf, ...) {
		m <- length(results)
		oldcall <- attr(results, "call")
		if (missing(variances)) {
			variances <- suppressWarnings(lapply(results, vcov))
			results <- lapply(results, coef)
		}
		vbar <- variances[[1]]
		cbar <- results[[1]]
		for (i in 2:m) {
			cbar <- cbar + results[[i]]
			vbar <- vbar + variances[[i]]
		}
		cbar <- cbar/m
		vbar <- vbar/m

		# MODIFICATION
		# evar <- var(do.call("rbind", results))
		evar <- sum( ( unlist( results ) - cbar )^2 / 4 )

		
		r <- (1 + 1/m) * evar/vbar
		df <- (m - 1) * (1 + 1/r)^2
		if (is.matrix(df)) df <- diag(df)
		if (is.finite(df.complete)) {
			dfobs <- ((df.complete + 1)/(df.complete + 3)) * df.complete *
			vbar/(vbar + evar)
			if (is.matrix(dfobs)) dfobs <- diag(dfobs)
			df <- 1/(1/dfobs + 1/df)
		}
		if (is.matrix(r)) r <- diag(r)
		rval <- list(coefficients = cbar, variance = vbar + evar *
		(m + 1)/m, call = c(oldcall, call), nimp = m, df = df,
		missinfo = (r + 2/(df + 3))/(r + 1))
		class(rval) <- "MIresult"
		rval
	}
```
---

## Download, Import, Preparation {-}

Download and unzip the 2019 fourth grade international database:

```{r eval = FALSE , results = "hide" }
library(httr)

tf <- tempfile()

this_url <- "https://timss2019.org/international-database/downloads/T19_G4_SPSS%20Data.zip"

GET( this_url , write_disk( tf ) , progress() )

unzipped_files <- unzip( tf , exdir = tempdir() )
```

Import and stack each of the student context data files for **Albania through Canada**:

```{r eval = FALSE , results = "hide" }
library(haven)

# limit unzipped files to those starting with `asg` followed by three letters followed by `m7`
asg_fns <- unzipped_files[ grepl( '^asg[a-z][a-z][a-z]m7' , basename( unzipped_files ) ) ]

# further limit asg files to the first ten countries
countries_thru_canada <- c("alb", "arm", "aus", "aut", "aze", "bhr", "bfl", "bih", "bgr", "can")

fns_thru_canada <- paste0( paste0( '^asg' , countries_thru_canada , 'm7' ) , collapse = "|" )

asg_alb_can_fns <- asg_fns[ grepl( fns_thru_canada , basename( asg_fns ) ) ]

timss_df <- NULL

for( spss_fn in asg_alb_can_fns ){

	this_tbl <- read_spss( spss_fn )
	
	this_tbl <- zap_labels( this_tbl )
	
	this_df <- data.frame( this_tbl )
	
	names( this_df ) <- tolower( names( this_df ) )
	
	timss_df <- rbind( timss_df , this_df )
	
}

# order the data.frame by unique student id
timss_df <- timss_df[ with( timss_df , order( idcntry , idstud ) ) , ]
```

### Save Locally \ {-}

Save the object at any point:

```{r eval = FALSE , results = "hide" }
# timss_fn <- file.path( path.expand( "~" ) , "TIMSS" , "this_file.rds" )
# saveRDS( timss_df , file = timss_fn , compress = FALSE )
```

Load the same object:

```{r eval = FALSE , results = "hide" }
# timss_df <- readRDS( timss_fn )
```

### Survey Design Definition {-}
Construct a multiply-imputed, complex sample survey design:

From among possibly plausible values, determine all columns that are multiply-imputed plausible values:

```{r eval = FALSE , results = "hide" }
# identify all columns ending with `01` thru `05`
ppv <- grep( "(.*)0[1-5]$" , names( timss_df ) , value = TRUE )

# remove those ending digits
ppv_prefix <- gsub( "0[1-5]$" , "" , ppv )

# identify each of the possibilities with exactly five matches (five implicates)
pv <- names( table( ppv_prefix )[ table( ppv_prefix ) == 5 ] )

# identify each of the `01` thru `05` plausible value columns
pv_columns <-
	grep( 
		paste0( "^" , pv , "0[1-5]$" , collapse = "|" ) , 
		names( timss_df ) , 
		value = TRUE 
	)
```

Extract those multiply-imputed columns into a separate data.frame, then remove them from the source:
```{r eval = FALSE , results = "hide" }
pv_wide_df <- timss_df[ c( 'idcntry' , 'idstud' , pv_columns ) ]

timss_df[ pv_columns ] <- NULL
```

Reshape these columns from one record per student to one record per student per implicate:
```{r eval = FALSE , results = "hide" }
pv_long_df <- 
	reshape( 
		pv_wide_df , 
		varying = lapply( paste0( pv , '0' ) , paste0 , 1:5 ) , 
		direction = 'long' , 
		timevar = 'implicate' , 
		idvar = c( 'idcntry' , 'idstud' ) 
	)

names( pv_long_df ) <- gsub( "01$" , "" , names( pv_long_df ) )
```

Merge the columns from the source data.frame onto the one record per student per implicate data.frame:
```{r eval = FALSE , results = "hide" }
timss_long_df <- merge( timss_df , pv_long_df )

timss_long_df <- timss_long_df[ with( timss_long_df , order( idcntry , idstud ) ) , ]

stopifnot( nrow( timss_long_df ) == nrow( pv_long_df ) )

stopifnot( nrow( timss_long_df ) / 5 == nrow( timss_df ) )
```

Divide the five plausible value implicates into a list with five data.frames based on the implicate number:
```{r eval = FALSE , results = "hide" }
timss_list <- split( timss_long_df , timss_long_df[ , 'implicate' ] )
```

Construct a replicate weights table following the estimation technique described in [Methods Chapter 14](https://timssandpirls.bc.edu/timss2019/methods/chapter-14.html):

```{r eval = FALSE , results = "hide" }
weights_df <- timss_df[ c( 'jkrep' , 'jkzone' ) ]

for( j in 1:75 ){
	for( i in 0:1 ){
		weights_df[ weights_df[ , 'jkzone' ] != j , paste0( 'rw' , i , j ) ] <- 1
		
		weights_df[ weights_df[ , 'jkzone' ] == j , paste0( 'rw' , i , j ) ] <- 
			2 * ( weights_df[ weights_df[ , 'jkzone' ] == j , 'jkrep' ] == i )
	}
}

weights_df[ c( 'jkrep' , 'jkzone' ) ] <- NULL

```

Define the design:
```{r eval = FALSE , results = "hide" }
library(survey)
library(mitools)

timss_design <- 
	svrepdesign(
		weights = ~totwgt ,
		repweights = weights_df , 
		data = imputationList( timss_list ) ,
		type = "other" ,
		scale = 0.5 ,
		rscales = rep( 1 , 150 ) ,
		combined.weights = FALSE ,
		mse = TRUE
	)
```

### Variable Recoding {-}

Add new columns to the data set:
```{r eval = FALSE , results = "hide" }
timss_design <- 
	update( 
		timss_design , 
		
		one = 1 ,
		
		countries_thru_canada = 
		
			factor( 
			
				as.numeric( idcntry ) ,
				
				levels = c(8L, 51L, 36L, 40L, 31L, 48L, 956L, 70L, 100L, 124L) ,

				labels =
					c("Albania", "Armenia", "Australia", "Austria", "Azerbaijan", "Bahrain",
					"Belgium (Flemish)", "Bosnia and Herzegovina", "Bulgaria", "Canada")
				
			) ,
		
		sex = factor( asbg01 , levels = 1:2 , labels = c( "female" , "male" ) ) ,
		
		born_in_country = ifelse( asbg07 %in% 1:2 , as.numeric( asbg07 == 1 ) , NA )

	)
```

---

## Analysis Examples with the `survey` library \ {-}

### Unweighted Counts {-}

Count the unweighted number of records in the survey sample, overall and by groups:
```{r eval = FALSE , results = "hide" }
timss_MIcombine( with( timss_design , svyby( ~ one , ~ one , unwtd.count ) ) )

timss_MIcombine( with( timss_design , svyby( ~ one , ~ sex , unwtd.count ) ) )
```

### Weighted Counts {-}
Count the weighted size of the generalizable population, overall and by groups:
```{r eval = FALSE , results = "hide" }
timss_MIcombine( with( timss_design , svytotal( ~ one ) ) )

timss_MIcombine( with( timss_design ,
	svyby( ~ one , ~ sex , svytotal )
) )
```

### Descriptive Statistics {-}

Calculate the mean (average) of a linear variable, overall and by groups:
```{r eval = FALSE , results = "hide" }
timss_MIcombine( with( timss_design , svymean( ~ asmmat , na.rm = TRUE ) ) )

timss_MIcombine( with( timss_design ,
	svyby( ~ asmmat , ~ sex , svymean , na.rm = TRUE )
) )
```

Calculate the distribution of a categorical variable, overall and by groups:
```{r eval = FALSE , results = "hide" }
timss_MIcombine( with( timss_design , svymean( ~ countries_thru_canada ) ) )

timss_MIcombine( with( timss_design ,
	svyby( ~ countries_thru_canada , ~ sex , svymean )
) )
```

Calculate the sum of a linear variable, overall and by groups:
```{r eval = FALSE , results = "hide" }
timss_MIcombine( with( timss_design , svytotal( ~ asmmat , na.rm = TRUE ) ) )

timss_MIcombine( with( timss_design ,
	svyby( ~ asmmat , ~ sex , svytotal , na.rm = TRUE )
) )
```

Calculate the weighted sum of a categorical variable, overall and by groups:
```{r eval = FALSE , results = "hide" }
timss_MIcombine( with( timss_design , svytotal( ~ countries_thru_canada ) ) )

timss_MIcombine( with( timss_design ,
	svyby( ~ countries_thru_canada , ~ sex , svytotal )
) )
```

Calculate the median (50th percentile) of a linear variable, overall and by groups:
```{r eval = FALSE , results = "hide" }
timss_MIcombine( with( timss_design ,
	svyquantile(
		~ asmmat ,
		0.5 , se = TRUE , na.rm = TRUE 
) ) )

timss_MIcombine( with( timss_design ,
	svyby(
		~ asmmat , ~ sex , svyquantile ,
		0.5 , se = TRUE ,
		ci = TRUE , na.rm = TRUE
) ) )
```

Estimate a ratio:
```{r eval = FALSE , results = "hide" }
timss_MIcombine( with( timss_design ,
	svyratio( numerator = ~ asssci , denominator = ~ asmmat )
) )
```

### Subsetting {-}

Restrict the survey design to Australia, Austria, Azerbaijan, Belgium (French):
```{r eval = FALSE , results = "hide" }
sub_timss_design <- subset( timss_design , idcntry %in% c( 36 , 40 , 31 , 956 ) )
```
Calculate the mean (average) of this subset:
```{r eval = FALSE , results = "hide" }
timss_MIcombine( with( sub_timss_design , svymean( ~ asmmat , na.rm = TRUE ) ) )
```

### Measures of Uncertainty {-}

Extract the coefficient, standard error, confidence interval, and coefficient of variation from any descriptive statistics function result, overall and by groups:
```{r eval = FALSE , results = "hide" }
this_result <-
	timss_MIcombine( with( timss_design ,
		svymean( ~ asmmat , na.rm = TRUE )
	) )

coef( this_result )
SE( this_result )
confint( this_result )
cv( this_result )

grouped_result <-
	timss_MIcombine( with( timss_design ,
		svyby( ~ asmmat , ~ sex , svymean , na.rm = TRUE )
	) )

coef( grouped_result )
SE( grouped_result )
confint( grouped_result )
cv( grouped_result )
```

Calculate the degrees of freedom of any survey design object:
```{r eval = FALSE , results = "hide" }
degf( timss_design$designs[[1]] )
```

Calculate the complex sample survey-adjusted variance of any statistic:
```{r eval = FALSE , results = "hide" }
timss_MIcombine( with( timss_design , svyvar( ~ asmmat , na.rm = TRUE ) ) )
```

Include the complex sample design effect in the result for a specific statistic:
```{r eval = FALSE , results = "hide" }
# SRS without replacement
timss_MIcombine( with( timss_design ,
	svymean( ~ asmmat , na.rm = TRUE , deff = TRUE )
) )

# SRS with replacement
timss_MIcombine( with( timss_design ,
	svymean( ~ asmmat , na.rm = TRUE , deff = "replace" )
) )
```

Compute confidence intervals for proportions using methods that may be more accurate near 0 and 1. See `?svyciprop` for alternatives:
```{r eval = FALSE , results = "hide" }
# MIsvyciprop( ~ born_in_country , timss_design ,
# 	method = "likelihood" , na.rm = TRUE )
```

### Regression Models and Tests of Association {-}

Perform a design-based t-test:
```{r eval = FALSE , results = "hide" }
# MIsvyttest( asmmat ~ born_in_country , timss_design )
```

Perform a chi-squared test of association for survey data:
```{r eval = FALSE , results = "hide" }
# MIsvychisq( ~ born_in_country + countries_thru_canada , timss_design )
```

Perform a survey-weighted generalized linear model:
```{r eval = FALSE , results = "hide" }
glm_result <- 
	timss_MIcombine( with( timss_design ,
		svyglm( asmmat ~ born_in_country + countries_thru_canada )
	) )
	
summary( glm_result )
```

---

## Replication Example {-}

This example matches the mean proficiency and standard error of the `Australia` row of the `Summary Statistics and Standard Errors for Proficiency in Overall Mathematics-Grade 4` table from the [Appendix 14A: Summary Statistics and Standard Errors for Proficiency in Grade 4 Mathematics](https://timss.bc.edu/timss2019/methods/pdf/T19_MP_Ch14-estimating-standard-errors.pdf#page=13):

```{r eval = FALSE , results = "hide" }
australia_design <- subset( timss_design , countries_thru_canada %in% "Australia" )

stopifnot( nrow( australia_design ) == 5890 )

result <- timss_MIcombine( with( australia_design , svymean( ~ asmmat ) ) )

stopifnot( round( coef( result ) , 3 ) == 515.880 )

stopifnot( round( SE( result ) , 3 ) == 2.776 )

```

This example matches the jackknife sampling, imputation, and total variances of the same row:

```{r eval = FALSE , results = "hide" }
australia_fn <- unzipped_files[ grepl( 'asgaus' , basename( unzipped_files ) ) ]
australia_tbl <- read_spss( australia_fn )
australia_tbl <- zap_labels( australia_tbl )
australia_df <- data.frame( australia_tbl )
names( australia_df ) <- tolower( names( australia_df ) )

estimate <-
	mean( c(
		with( australia_df , weighted.mean( asmmat01 , totwgt ) ) ,
		with( australia_df , weighted.mean( asmmat02 , totwgt ) ) ,
		with( australia_df , weighted.mean( asmmat03 , totwgt ) ) ,
		with( australia_df , weighted.mean( asmmat04 , totwgt ) ) ,
		with( australia_df , weighted.mean( asmmat05 , totwgt ) )
	) )

stopifnot( round( estimate , 3 ) == 515.880 )

for( k in 1:5 ){

	this_variance <- 0
	
	for( j in 1:75 ){
		for( i in 0:1 ){
			this_variance <- 
				this_variance + 
				( 
					weighted.mean( 
						australia_df[ , paste0( 'asmmat0' , k ) ] , 
						ifelse( 
							j == australia_df[ , 'jkzone' ] , 
							australia_df[ , 'totwgt' ] * 2 * ( australia_df[ , 'jkrep' ] == i ) , 
							australia_df[ , 'totwgt' ] 
						)
					) -
					weighted.mean( 
						australia_df[ , paste0( 'asmmat0' , k ) ] , 
						australia_df[ , 'totwgt' ]
					)
				)^2
		}
	}
	
	assign( paste0( 'v' , k ) , this_variance * 0.5 )

}

sampling_variance <- mean( c( v1 , v2 , v3 , v4 , v5 ) )
stopifnot( round( sampling_variance , 3 ) == 7.397 )

imputation_variance <-
	( 6 / 5 ) * 
	( 
		( ( with( australia_df , weighted.mean( asmmat01 , totwgt ) ) - estimate )^2 / 4 ) +
		( ( with( australia_df , weighted.mean( asmmat02 , totwgt ) ) - estimate )^2 / 4 ) +
		( ( with( australia_df , weighted.mean( asmmat03 , totwgt ) ) - estimate )^2 / 4 ) +
		( ( with( australia_df , weighted.mean( asmmat04 , totwgt ) ) - estimate )^2 / 4 ) +
		( ( with( australia_df , weighted.mean( asmmat05 , totwgt ) ) - estimate )^2 / 4 ) 
	)

stopifnot( round( imputation_variance , 3 ) == 0.309 )

stopifnot( round( sampling_variance + imputation_variance , 3 ) == 7.706 )

```