45-piaac.Rmd

# Programme for the International Assessment of Adult Competencies (PIAAC) {-}

[![Build Status](https://travis-ci.org/asdfree/piaac.svg?branch=master)](https://travis-ci.org/asdfree/piaac) [![Build status](https://ci.appveyor.com/api/projects/status/github/asdfree/piaac?svg=TRUE)](https://ci.appveyor.com/project/ajdamico/piaac)

The Programme for the International Assessment of Adult Competencies (PIAAC) offers cross-national comparisons for the serious study of advanced-nation labor markets.

* One row per sampled adult.

* A multiply-imputed, complex sample survey designed to generalize to the population aged 16 to 65 across thirty three OECD nations.

* No expected release timeline.

* Administered by the [Organisation for Economic Co-operation and Development](http://www.oecd.org/).

## Simplified Download and Importation {-}

The R `lodown` package easily downloads and imports all available PIAAC microdata by simply specifying `"piaac"` with an `output_dir =` parameter in the `lodown()` function. Depending on your internet connection and computer processing speed, you might prefer to run this step overnight.

```{r eval = FALSE }
library(lodown)
lodown( "piaac" , output_dir = file.path( path.expand( "~" ) , "PIAAC" ) )
```

`lodown` also provides a catalog of available microdata extracts with the `get_catalog()` function. After requesting the PIAAC catalog, you could pass a subsetted catalog through the `lodown()` function in order to download and import specific extracts (rather than all available extracts).

```{r eval = FALSE , results = "hide" }
library(lodown)
# examine all available PIAAC microdata files
piaac_cat <-
	get_catalog( "piaac" ,
		output_dir = file.path( path.expand( "~" ) , "PIAAC" ) )

# download the microdata to your local computer
piaac_cat <- lodown( "piaac" , piaac_cat )
```

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

Construct a multiply-imputed, complex sample survey design:

```{r eval = FALSE }

```

```{r eval = FALSE }
library(survey)
library(mitools)

piaac_design <- readRDS( file.path( path.expand( "~" ) , "PIAAC" , "prgusap1 design.rds" ) )
```

### Variable Recoding {-}

Add new columns to the data set:
```{r eval = FALSE }
piaac_design <-
	update(
		piaac_design ,
		
		one = 1 ,
		
		sex = factor( gender_r , labels = c( "male" , "female" ) ) ,

		age_categories = 
			factor( 
				ageg10lfs , 
				levels = 1:5 , 
				labels = c( "24 or less" , "25-34" , "35-44" , "45-54" , "55 plus" ) 
			) ,
		
		working_at_paid_job_last_week = as.numeric( c_q01a == 1 )
		
	)
```

### Unweighted Counts {-}

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

MIcombine( with( piaac_design , svyby( ~ one , ~ age_categories , unwtd.count ) ) )
```

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

MIcombine( with( piaac_design ,
	svyby( ~ one , ~ age_categories , svytotal )
) )
```

### Descriptive Statistics {-}

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

MIcombine( with( piaac_design ,
	svyby( ~ pvnum , ~ age_categories , svymean , na.rm = TRUE )
) )
```

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

MIcombine( with( piaac_design ,
	svyby( ~ sex , ~ age_categories , svymean )
) )
```

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

MIcombine( with( piaac_design ,
	svyby( ~ pvnum , ~ age_categories , svytotal , na.rm = TRUE )
) )
```

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

MIcombine( with( piaac_design ,
	svyby( ~ sex , ~ age_categories , svytotal )
) )
```

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

MIcombine( with( piaac_design ,
	svyby(
		~ pvnum , ~ age_categories , svyquantile ,
		0.5 , se = TRUE ,
		keep.var = TRUE , ci = TRUE , na.rm = TRUE
) ) )
```

Estimate a ratio:
```{r eval = FALSE , results = "hide" }
MIcombine( with( piaac_design ,
	svyratio( numerator = ~ pvnum , denominator = ~ pvlit , na.rm = TRUE )
) )
```

### Subsetting {-}

Restrict the survey design to self-reported fair or poor health:
```{r eval = FALSE , results = "hide" }
sub_piaac_design <- subset( piaac_design , i_q08 %in% 4:5 )
```
Calculate the mean (average) of this subset:
```{r eval = FALSE , results = "hide" }
MIcombine( with( sub_piaac_design , svymean( ~ pvnum , 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 <-
	MIcombine( with( piaac_design ,
		svymean( ~ pvnum , na.rm = TRUE )
	) )

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

grouped_result <-
	MIcombine( with( piaac_design ,
		svyby( ~ pvnum , ~ age_categories , 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( piaac_design$designs[[1]] )
```

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

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

# SRS with replacement
MIcombine( with( piaac_design ,
	svymean( ~ pvnum , 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( ~ working_at_paid_job_last_week , piaac_design ,
	method = "likelihood" )
```

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

Perform a design-based t-test:
```{r eval = FALSE , results = "hide" }
MIsvyttest( pvnum ~ working_at_paid_job_last_week , piaac_design )
```

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

Perform a survey-weighted generalized linear model:
```{r eval = FALSE , results = "hide" }
glm_result <- 
	MIcombine( with( piaac_design ,
		svyglm( pvnum ~ working_at_paid_job_last_week + sex )
	) )
	
summary( glm_result )
```

---

## Replication Example {-}

The [OECD's Technical Report Table 18.9 on PDF page 455](https://www.oecd.org/skills/piaac/_Technical Report_17OCT13.pdf#page=455) includes statistics and standard errors for the three PIAAC domains. This code precisely replicates the Austria row shown in that official table.

```{r eval = FALSE , results = "hide" }

austria_design <-
	readRDS( file.path( path.expand( "~" ) , "PIAAC" , "prgautp1 design.rds" ) )

austria_pvlit <-
	MIcombine( with( austria_design , svymean( ~ pvlit , na.rm = TRUE ) ) )
	
austria_pvnum <-
	MIcombine( with( austria_design , svymean( ~ pvnum , na.rm = TRUE ) ) )

austria_pvpsl <-
	MIcombine( with( austria_design , svymean( ~ pvpsl , na.rm = TRUE ) ) )
	
# confirm each estimate and standard error matches the published statistics
stopifnot( round( coef( austria_pvlit ) ) == 269 )
stopifnot( round( SE( austria_pvlit ) , 1 ) == 0.7 )
stopifnot( round( coef( austria_pvnum ) ) == 275 )
stopifnot( round( SE( austria_pvnum ) , 1 ) == 0.9 )
stopifnot( round( coef( austria_pvpsl ) ) == 284 )
stopifnot( round( SE( austria_pvpsl ) , 1 ) == 0.7 )

```