ANALYSIS NOTES - Primary and Secondary Effects for Educational Expectations: An International Perspective

The supplementary material contains the formulars used to construct total, primary, and secondary effects as estimated in the paper. In addition, in keeping with growing calls for reproducible research (see Mesirov, 2010 Accessible reproducible research, Science, 327, 415-416)all analysis steps to obtain the final solutions presented in the paper and links to obtain the original data are presented. The paper below provides actual input code but also associated outputs. As the exact outputs are presented any warnings produced by R are maintained. These warnings were all investigated and found to not effect the results in any way. The warning are left in tact to provide a full record for researchers who wish to replicate the results.

Formula Used to Derive Primary and Secondary Effects in Current Research

The general formula is: \( \begin{equation}y^{i} = \alpha +B y^{i} + \Gamma x^{i} + \zeta_{i}\end{equation} \)

With matricies: \[ \begin{aligned}y_{i} = \begin{bmatrix} Expect\\ Ach \end{bmatrix}\end{aligned} \]

\[ \begin{aligned}\alpha = \begin{bmatrix} \beta^{Expect}_{0}\\ \beta^{Ach}_{0} \end{bmatrix}\end{aligned} \]

\[ \begin{aligned}B= \begin{bmatrix} 0 & \beta^{Ach}_{Expect}\\ 0&0 \end{bmatrix}\end{aligned} \]

\[ \begin{aligned}\Gamma= \begin{bmatrix} \beta^{Expect}_{SES}\\ \beta^{Ach}_{SES} \end{bmatrix}\end{aligned} \]

\[ \begin{aligned} \zeta = \begin{bmatrix} e^{Expect}_{i}\\ e^{Ach}_{i} \end{bmatrix} \sim N(o, \Psi ) \text{where} \Psi = \begin{bmatrix} \sigma^{2}_{ei:Expect} & 0\\ 0&\sigma^{2}_{ei:Ach} \end{bmatrix}\end{aligned} \]

The effects of interest are then derived from the coefficents from this model as follows:

Secondary effects \( =\beta^{Expect}_{SES} \)

Primary effects \( =\beta^{Ach}_{SES} \times \beta^{Expect}_{Ach} \)

Total effects \( =\beta^{Expect}_{SES} + (\beta^{Ach}_{SES} \times \beta^{Expect}_{Ach}) \)

and where:

Expect = the expectations of obtaining a university level of education

Ach = academic achievement as measured by PISA achievement tests

SES = socioeconomic status of the child

Data Preperation

Data is taken from the PISA 2003 database which can be found at the ACER website. The data is read into R from the developed SPSS datafile (studentData.sav). In this step non-OECD countries are dropped.

University expectations variable is also created from the SISCED variable in the PISA database.

# read in spss data Export as RData

library(foreign)
studentData <- read.spss("C:/Users/30016475/Dropbox/Projects_Research/Social Forces/PISA2003/data/studentData.sav", 
    to.data.frame = TRUE, use.value.labels = FALSE)
## Warning: C:/Users/30016475/Dropbox/Projects_Research/Social
## Forces/PISA2003/data/studentData.sav: Unrecognized record type 7, subtype
## 18 encountered in system file
reducedData <- studentData[, c("CNT", "SCHOOLID", "AGE", "GRADE", "ST03Q01", 
    "HISEI", "HISCED", "IMMIG", "OECD", "SISCED", "W_FSTUWT", grep("^PV.*(MATH|READ|SCIE)$", 
        names(studentData), value = TRUE), "ESCS")]
reducedData <- reducedData[reducedData$OECD == 1, ]
# Remove non-OECD from levels
reducedData$CNT <- factor(reducedData$CNT)
# Add university expectation binary variable
reducedData$UNI <- ifelse(reducedData$SISCED == 5, 1, ifelse(reducedData$SISCED == 
    9, NA, 0))
# For pooled standadization z-score HISEI and HISCED
reducedData[, 6:7] <- scale(reducedData[, 6:7])
names(reducedData)
##  [1] "CNT"      "SCHOOLID" "AGE"      "GRADE"    "ST03Q01"  "HISEI"   
##  [7] "HISCED"   "IMMIG"    "OECD"     "SISCED"   "W_FSTUWT" "PV1MATH" 
## [13] "PV2MATH"  "PV3MATH"  "PV4MATH"  "PV5MATH"  "PV1READ"  "PV2READ" 
## [19] "PV3READ"  "PV4READ"  "PV5READ"  "PV1SCIE"  "PV2SCIE"  "PV3SCIE" 
## [25] "PV4SCIE"  "PV5SCIE"  "ESCS"     "UNI"
save(reducedData, file = "C:/Users/30016475/Dropbox/Projects_Research/Social Forces/PISA2003/data/reducedData19013013.RData")
rm(list = ls(all = TRUE))

To manage memory the smaller analysis dataset is saved and then the workspace is cleared. The smaller data set is then sourced for the remaining analysis.

Creation of Principal Components

Principals components of PISA math, verbal, and science tests were constructed as follows:

# Load RData
load("C:/Users/30016475/Dropbox/Projects_Research/Social Forces/PISA2003/data/reducedData19013013.RData")
# names(reducedData) Z-score achievement
reducedData[, grep("^PV[0-9]", names(reducedData))] <- scale(reducedData[, grep("^PV[0-9]", 
    names(reducedData))])
# Principal Value analaysis
library(psych)
# Achievement scores
Ach <- list()
for (i in 1:5) {
    Ach[[i]] <- principal(reducedData[, grep(paste0("^PV", i, ".*$"), names(reducedData))], 
        1, score = TRUE)$scores
}

Ach <- data.frame(do.call(cbind, Ach))
names(Ach) <- paste0("Ach", 1:5)
# add principal components to dataframe
AchData <- data.frame(reducedData, Ach)
# Save the data with the ICCs and principal components
save(AchData, file = "C:/Users/30016475/Dropbox/Projects_Research/Social Forces/PISA2003/data/AchData19013013.RData")
rm(list = ls(all = TRUE))

Construction of Intraclass correlations (ICC)

Country ICCs created using the psych package from a set of anova models. These ICCs were created independently for all plausable values and the resulting average across the plausible values is used in further analysis.

# ICCs
load("C:/Users/30016475/Dropbox/Projects_Research/Social Forces/PISA2003/data/AchData19013013.RData")
library(multilevel)
## Warning: package 'multilevel' was built under R version 3.0.1
## Loading required package: nlme
## Loading required package: MASS
ICCdata <- split(AchData, AchData$CNT)
ICC <- list()
for (i in 1:5) {
    x <- lapply(ICCdata, function(x) ICC1(aov(x[, i + 27] ~ as.factor(SCHOOLID), 
        data = x, weights = W_FSTUWT)))
    ICC[[i]] <- data.frame(do.call(rbind, x))
}

ICCs <- do.call(cbind, ICC)
ICCs <- apply(ICCs, 1, mean)
ICCs <- data.frame(ICCs = matrix(as.numeric(ICCs), nrow = 30), CNT = names(ICCs))
names(ICCs) <- c("ICC", "CNT")

# Will need to report this in a table
ICCs[order(ICCs[, 1], decreasing = TRUE), ]
##        ICC CNT
## 22 0.57481 NLD
## 14 0.55328 HUN
## 7  0.55207 DEU
## 29 0.54890 TUR
## 17 0.53626 ITA
## 2  0.53583 AUT
## 18 0.51243 JPN
## 3  0.50623 BEL
## 11 0.47348 FRA
## 6  0.46235 CZE
## 21 0.45694 MEX
## 27 0.44607 SVK
## 19 0.40453 KOR
## 13 0.40335 GRC
## 5  0.35979 CHE
## 26 0.34198 PRT
## 20 0.31423 LUX
## 12 0.23802 GBR
## 30 0.22153 USA
## 9  0.21958 ESP
## 1  0.21646 AUS
## 4  0.17140 CAN
## 15 0.16553 IRL
## 24 0.16168 NZL
## 8  0.14245 DNK
## 25 0.14155 POL
## 28 0.09789 SWE
## 23 0.07808 NOR
## 10 0.04364 FIN
## 16 0.04143 ISL
# Aggregate the ICCs back into the data set
AchData <- merge(AchData, ICCs, by = "CNT")

Descriptive Results

Descriptive results consist of:

# Age of first selection
descriptives <- data.frame(n.obs = tapply(AchData$CNT, AchData$CNT, length), 
    selection = c(16, 10, 12, 16, 15, 10.5, 10, 16, 16, 16, 15, 16, 14.5, 11, 
        15, 16, 14, 15, 14, 12.5, 11.5, 12, 15.5, 16, 15, 14.5, 10.5, 16, 10, 
        16), ICC = round(tapply(AchData$ICC, AchData$CNT, mean), 3), Expect.tables = round(tapply(AchData$UNI, 
        AchData$CNT, mean, na.rm = TRUE), 3), avg.ach = round(tapply(AchData$Ach1, 
        AchData$CNT, mean), 3))
descriptives
##     n.obs selection   ICC Expect.tables avg.ach
## AUS 12551      16.0 0.216         0.613   0.305
## AUT  4597      10.0 0.536         0.251   0.084
## BEL  8796      12.0 0.506         0.350   0.256
## CAN 27953      16.0 0.171         0.637   0.223
## CHE  8420      15.0 0.360         0.169   0.103
## CZE  6320      10.5 0.462         0.437   0.342
## DEU  4660      10.0 0.552         0.197   0.110
## DNK  4218      16.0 0.142         0.247  -0.018
## ESP 10791      16.0 0.220         0.503  -0.022
## FIN  5796      16.0 0.044         0.504   0.507
## FRA  4300      15.0 0.473         0.353   0.172
## GBR  9535      16.0 0.238         0.368   0.220
## GRC  4627      14.5 0.403         0.626  -0.343
## HUN  4765      11.0 0.553         0.526  -0.038
## IRL  3880      15.0 0.166         0.534   0.169
## ISL  3350      16.0 0.041         0.359   0.063
## ITA 11639      14.0 0.536         0.477   0.098
## JPN  4707      15.0 0.512         0.501   0.334
## KOR  5444      14.0 0.405         0.774   0.448
## LUX  3923      12.5 0.314         0.408  -0.095
## MEX 29983      11.5 0.457         0.563  -0.821
## NLD  3992      12.0 0.575         0.424   0.375
## NOR  4064      15.5 0.078         0.256  -0.005
## NZL  4511      16.0 0.162         0.397   0.319
## POL  4383      15.0 0.142         0.299  -0.008
## PRT  4608      14.5 0.342         0.501  -0.258
## SVK  7346      10.5 0.446         0.447  -0.008
## SWE  4624      16.0 0.098         0.327   0.154
## TUR  4855      10.0 0.549         0.784  -0.616
## USA  5456      16.0 0.222         0.643  -0.060
cor(descriptives[, 2:5])
##               selection     ICC Expect.tables avg.ach
## selection      1.000000 -0.7873     -0.004396  0.2964
## ICC           -0.787274  1.0000      0.116549 -0.1785
## Expect.tables -0.004396  0.1165      1.000000 -0.2021
## avg.ach        0.296414 -0.1785     -0.202133  1.0000

Mplus Probit Models

First the R data file is saved for use in Mplus using the MplusAutomation package. Five data sets are created each with one plausable value for the Achievement principal component.

library(MplusAutomation)
## Warning: package 'MplusAutomation' was built under R version 3.0.1
## Loading required package: tcltk
## Loading required package: plyr
## Warning: package 'plyr' was built under R version 3.0.1
## Loading required package: boot
## Attaching package: 'boot'
## The following object is masked from 'package:psych':
## 
## logit
for (i in 1:5) {
    capture.output(prepareMplusData(AchData[, c(1:28, 34, 28 + i)], filename = paste0("C:/Users/30016475/Dropbox/Projects_Research/Social Forces/PISA2003/data/PISA2003_PV", 
        i, ".dat")))
}
## Warning: The file 'PISA2003_PV1.dat' currently exists and will be
## overwritten
## Warning: The file 'PISA2003_PV2.dat' currently exists and will be
## overwritten
## Warning: The file 'PISA2003_PV3.dat' currently exists and will be
## overwritten
## Warning: The file 'PISA2003_PV4.dat' currently exists and will be
## overwritten
## Warning: The file 'PISA2003_PV5.dat' currently exists and will be
## overwritten

All models are run using a single Mplus template file. The template file is:

TITLE: Country XXX !Outcome = Expectation of university entry
                    !Treatment mediator Interaction;
DATA: FILE = "C:/Users/30016475/Dropbox/Projects_Research/Social Forces/PISA2003/data/IMPUTATION.dat";
type=imputation;

VARIABLE: NAMES = CNT SCHOOLID AGE GRADE ST03Q01 HISEI HISCED IMMIG OECD SISCED
     W_FSTUWT PV1MATH PV2MATH PV3MATH PV4MATH PV5MATH PV1READ PV2READ PV3READ
     PV4READ PV5READ PV1SCIE PV2SCIE PV3SCIE PV4SCIE PV5SCIE ESCS UNI ICC PVAch;
MISSING=.;

USEVARIABLES UNI HISEI PVAch int;
categorical = UNI;
useobs CNT EQ 1;
cluster=schoolid;
weight=W_FSTUWT;

DEFINE: int = HISEI*PVAch;

analysis: type=complex; estimator = MLR; LINK = PROBIT;
DEFINE: STANDARDIZE HISEI PVAch int;

MODEL:
[UNI$1] (mbeta0);
UNI on HISEI (beta2);
UNI on PVAch (beta1);
UNI on int (beta3);
[PVAch](gamma0);
PVAch on HISEI (gamma1);
PVAch(sig2);

MODEL CONSTRAINT:
new(ind dir arg11 arg10 arg01 arg00 v1 v0 
    probit11 probit10 probit01 probit00 
    tie de pie ortie orde orpie Ctotal total prop);

    dir=beta3*gamma0+beta2;
    ind=beta1*gamma1+beta3*gamma1;
    arg11=-mbeta0+beta2+(beta1+beta3)*(gamma0+gamma1);
    arg10=-mbeta0+beta2+(beta1+beta3)*gamma0;
    arg01=-mbeta0+beta1*(gamma0+gamma1);
    arg00=-mbeta0+beta1*gamma0;
    v1=(beta1+beta3)^2*sig2+1;
    v0=beta1^2*sig2+1;
    probit11=arg11/sqrt(v1);
    probit10=arg10/sqrt(v1);
    probit01=arg01/sqrt(v0);
    probit00=arg00/sqrt(v0);
    ! Phi function needed below:
    tie=phi(probit11)-phi(probit10);
    de=phi(probit10)-phi(probit00);
    pie=phi(probit01)-phi(probit00);
    ortie=(phi(probit11)/(1-phi(probit11)))/
    (phi(probit10)/(1-phi(probit10)));
    orde=(phi(probit10)/(1-phi(probit10)))/
    (phi(probit00)/(1-phi(probit00)));
    orpie=(phi(probit01)/(1-phi(probit01)))/
    (phi(probit00)/(1-phi(probit00)));
    Ctotal = de + tie;
    total= ind + dir;
    prop = (tie)/ (Ctotal+.00001);

The template file is iterative over in order to create country specific syntax files. The creation and running of files is done via the following user defined function. Before creating these files relavent directories are created in the working directory; one each for every ses.var used.

fileCreation <- function(ses.var = "HISEI") {
    temp <- readLines("C:/Users/30016475/Dropbox/Projects_Research/Social Forces/PISA2003/MPLUSforMARKDOWN/TEMPLATEFILE19012013.inp")
    temp[10:60] <- unlist(sapply(temp[10:60], function(x) gsub("HISEI", ses.var, 
        x)))
    countries <- levels(AchData$CNT)
    for (i in 1:30) {
        x <- temp
        x[1] <- gsub("XXX", countries[i], x[1])
        x[13] <- paste0("useobs CNT EQ ", i, ";")
        cat(x, sep = "\n", file = paste0("C:/Users/30016475/Dropbox/Projects_Research/Social Forces/PISA2003/MPLUSforMARKDOWN/", 
            ses.var, "/", countries[i], ".inp"))
    }
}
fileCreation(ses.var = "HISEI")  # Highest occupation prestige of parents
## Warning: incomplete final line found on
## 'C:/Users/30016475/Dropbox/Projects_Research/Social
## Forces/PISA2003/MPLUSforMARKDOWN/TEMPLATEFILE19012013.inp'
fileCreation(ses.var = "HISCED")  # Highest education level of parents
## Warning: incomplete final line found on
## 'C:/Users/30016475/Dropbox/Projects_Research/Social
## Forces/PISA2003/MPLUSforMARKDOWN/TEMPLATEFILE19012013.inp'
fileCreation(ses.var = "ESCS")  # Highest education level of parents
## Warning: incomplete final line found on
## 'C:/Users/30016475/Dropbox/Projects_Research/Social
## Forces/PISA2003/MPLUSforMARKDOWN/TEMPLATEFILE19012013.inp'

All the models in each of the folders are run using MplusAutomation package in R

runModels("C:/Users/30016475/Dropbox/Projects_Research/Social Forces/PISA2003/MPLUSforMARKDOWN", recursive=TRUE)

Parameter Extraction and Presentation

Extraction of results for presentaton and further analysis was done using the following user defined functions.

# Run on mplus files one at a time
paths <- c("C:/Users/30016475/Dropbox/Projects_Research/Social Forces/PISA2003/MPLUSforMARKDOWN/HISEI", 
    "C:/Users/30016475/Dropbox/Projects_Research/Social Forces/PISA2003/MPLUSforMARKDOWN/HISCED", 
    "C:/Users/30016475/Dropbox/Projects_Research/Social Forces/PISA2003/MPLUSforMARKDOWN/ESCS")

parameterExtract <- function(path) {
    files <- list.files(path, full.names = TRUE)
    files <- grep(".out$", files, value = TRUE)
    # Extract relavent output
    temp <- grep("New/Additional Parameters", readLines(files[1]))
    lines <- temp + c(13, 14, 19, 21)
    output <- sapply(files, function(x) readLines(x)[lines])
    # Function for extracting relavent parameters
    f1 <- function(x) unlist(strsplit(x, "[[:space:]]+"))[c(2:4)]
    results <- apply(output, c(1, 2), f1)
    TIE <- matrix(results[, 1, 1:30], 30, 3, byrow = TRUE)
    DE <- matrix(results[, 2, 1:30], 30, 3, byrow = TRUE)
    total <- matrix(results[, 3, 1:30], 30, 3, byrow = TRUE)
    prop <- matrix(results[, 4, 1:30], 30, 3, byrow = TRUE)
    results <- cbind(total, TIE, DE, prop)
    results <- data.frame(results[, c(2:3, 5:6, 8:9, 11:12)], descriptives[, 
        3])
    names(results) <- c("total", "totalse", "ind", "indse", "dir", "dirse", 
        "prop", "propse", "ICC")
    results <- data.frame(apply(results, 2, as.numeric))
    results$CNT <- levels(AchData$CNT)
    results
}

resultsList <- lapply(paths, parameterExtract)
names(resultsList) <- c("HISEI", "HISCED", "ESCS")
resultsList
## $HISEI
##    total totalse   ind indse   dir dirse  prop propse   ICC CNT
## 1  0.109   0.006 0.059 0.003 0.050 0.006 0.539  0.040 0.216 AUS
## 2  0.153   0.013 0.085 0.008 0.068 0.009 0.556  0.034 0.536 AUT
## 3  0.177   0.008 0.102 0.006 0.075 0.007 0.579  0.032 0.506 BEL
## 4  0.105   0.006 0.040 0.003 0.065 0.006 0.378  0.033 0.171 CAN
## 5  0.119   0.015 0.067 0.008 0.053 0.010 0.559  0.045 0.360 CHE
## 6  0.195   0.010 0.086 0.007 0.109 0.009 0.441  0.031 0.462 CZE
## 7  0.149   0.010 0.087 0.006 0.062 0.008 0.584  0.039 0.552 DEU
## 8  0.108   0.009 0.051 0.005 0.057 0.009 0.472  0.048 0.142 DNK
## 9  0.169   0.008 0.078 0.006 0.091 0.008 0.460  0.032 0.220 ESP
## 10 0.097   0.008 0.033 0.004 0.064 0.009 0.342  0.053 0.044 FIN
## 11 0.139   0.011 0.090 0.008 0.049 0.011 0.646  0.057 0.473 FRA
## 12 0.129   0.009 0.083 0.007 0.047 0.007 0.640  0.040 0.238 GBR
## 13 0.159   0.013 0.065 0.009 0.094 0.009 0.408  0.037 0.403 GRC
## 14 0.211   0.011 0.092 0.008 0.119 0.010 0.436  0.034 0.553 HUN
## 15 0.117   0.009 0.060 0.006 0.058 0.010 0.509  0.061 0.166 IRL
## 16 0.103   0.013 0.027 0.004 0.076 0.011 0.262  0.039 0.041 ISL
## 17 0.172   0.009 0.045 0.008 0.127 0.008 0.260  0.039 0.536 ITA
## 18 0.131   0.012 0.045 0.007 0.086 0.011 0.344  0.047 0.512 JPN
## 19 0.070   0.007 0.028 0.005 0.042 0.007 0.397  0.064 0.405 KOR
## 20 0.185   0.019 0.089 0.009 0.096 0.016 0.482  0.044 0.314 LUX
## 21 0.151   0.014 0.060 0.006 0.091 0.012 0.398  0.034 0.457 MEX
## 22 0.142   0.012 0.098 0.008 0.044 0.012 0.689  0.064 0.575 NLD
## 23 0.131   0.010 0.044 0.004 0.088 0.010 0.332  0.033 0.078 NOR
## 24 0.113   0.009 0.058 0.006 0.055 0.010 0.514  0.055 0.162 NZL
## 25 0.176   0.010 0.073 0.006 0.103 0.009 0.417  0.032 0.142 POL
## 26 0.163   0.011 0.085 0.008 0.078 0.009 0.521  0.041 0.342 PRT
## 27 0.190   0.010 0.099 0.007 0.090 0.009 0.523  0.032 0.446 SVK
## 28 0.121   0.009 0.040 0.004 0.081 0.009 0.334  0.038 0.098 SWE
## 29 0.085   0.011 0.050 0.008 0.035 0.009 0.589  0.075 0.549 TUR
## 30 0.093   0.006 0.042 0.003 0.051 0.007 0.450  0.042 0.222 USA
## 
## $HISCED
##    total totalse   ind indse   dir dirse  prop propse   ICC CNT
## 1  0.121   0.007 0.042 0.004 0.079 0.006 0.348  0.029 0.216 AUS
## 2  0.154   0.012 0.055 0.007 0.099 0.009 0.356  0.034 0.536 AUT
## 3  0.176   0.010 0.079 0.006 0.097 0.009 0.450  0.033 0.506 BEL
## 4  0.118   0.006 0.027 0.002 0.091 0.007 0.227  0.022 0.171 CAN
## 5  0.158   0.017 0.054 0.009 0.104 0.012 0.343  0.036 0.360 CHE
## 6  0.202   0.010 0.089 0.007 0.113 0.009 0.440  0.032 0.462 CZE
## 7  0.158   0.012 0.086 0.008 0.072 0.010 0.546  0.042 0.552 DEU
## 8  0.133   0.015 0.053 0.005 0.079 0.013 0.402  0.040 0.142 DNK
## 9  0.159   0.009 0.064 0.006 0.095 0.007 0.402  0.029 0.220 ESP
## 10 0.098   0.007 0.027 0.003 0.071 0.008 0.278  0.038 0.044 FIN
## 11 0.138   0.011 0.076 0.007 0.062 0.009 0.554  0.047 0.473 FRA
## 12 0.150   0.010 0.066 0.006 0.084 0.009 0.439  0.036 0.238 GBR
## 13 0.140   0.012 0.055 0.007 0.085 0.008 0.396  0.032 0.403 GRC
## 14 0.224   0.011 0.096 0.008 0.127 0.010 0.430  0.030 0.553 HUN
## 15 0.119   0.009 0.050 0.005 0.070 0.009 0.416  0.047 0.166 IRL
## 16 0.172   0.013 0.033 0.004 0.139 0.013 0.194  0.025 0.041 ISL
## 17 0.153   0.009 0.045 0.007 0.108 0.008 0.296  0.039 0.536 ITA
## 18 0.152   0.014 0.051 0.009 0.101 0.012 0.334  0.052 0.512 JPN
## 19 0.092   0.007 0.031 0.005 0.061 0.006 0.333  0.046 0.405 KOR
## 20 0.169   0.021 0.072 0.009 0.097 0.016 0.428  0.037 0.314 LUX
## 21 0.164   0.015 0.057 0.006 0.107 0.011 0.345  0.025 0.457 MEX
## 22 0.104   0.012 0.066 0.007 0.038 0.010 0.634  0.071 0.575 NLD
## 23 0.136   0.013 0.030 0.004 0.106 0.012 0.218  0.026 0.078 NOR
## 24 0.143   0.011 0.053 0.005 0.090 0.010 0.369  0.035 0.162 NZL
## 25 0.180   0.010 0.073 0.006 0.107 0.009 0.404  0.029 0.142 POL
## 26 0.134   0.010 0.065 0.008 0.069 0.007 0.484  0.041 0.342 PRT
## 27 0.195   0.010 0.099 0.008 0.096 0.009 0.507  0.034 0.446 SVK
## 28 0.131   0.011 0.032 0.004 0.100 0.010 0.241  0.025 0.098 SWE
## 29 0.107   0.010 0.062 0.008 0.044 0.008 0.586  0.055 0.549 TUR
## 30 0.131   0.007 0.033 0.003 0.098 0.007 0.251  0.026 0.222 USA
## 
## $ESCS
##    total totalse   ind indse   dir dirse  prop propse   ICC CNT
## 1  0.154   0.009 0.063 0.004 0.091 0.008 0.408  0.029 0.216 AUS
## 2  0.223   0.014 0.101 0.008 0.122 0.011 0.451  0.030 0.536 AUT
## 3  0.252   0.009 0.129 0.007 0.123 0.009 0.511  0.028 0.506 BEL
## 4  0.157   0.006 0.039 0.003 0.118 0.006 0.248  0.020 0.171 CAN
## 5  0.196   0.020 0.097 0.011 0.099 0.013 0.495  0.033 0.360 CHE
## 6  0.252   0.010 0.102 0.007 0.150 0.010 0.406  0.027 0.462 CZE
## 7  0.198   0.012 0.111 0.008 0.088 0.010 0.559  0.037 0.552 DEU
## 8  0.162   0.013 0.072 0.006 0.090 0.013 0.443  0.044 0.142 DNK
## 9  0.222   0.008 0.096 0.006 0.126 0.008 0.432  0.026 0.220 ESP
## 10 0.135   0.007 0.042 0.004 0.093 0.009 0.312  0.040 0.044 FIN
## 11 0.196   0.012 0.117 0.009 0.079 0.012 0.595  0.044 0.473 FRA
## 12 0.193   0.010 0.102 0.007 0.091 0.009 0.528  0.034 0.238 GBR
## 13 0.192   0.015 0.066 0.009 0.126 0.010 0.343  0.031 0.403 GRC
## 14 0.274   0.012 0.101 0.010 0.173 0.012 0.368  0.033 0.553 HUN
## 15 0.177   0.009 0.065 0.008 0.112 0.011 0.365  0.046 0.166 IRL
## 16 0.170   0.013 0.043 0.004 0.127 0.013 0.254  0.029 0.041 ISL
## 17 0.214   0.010 0.054 0.008 0.160 0.009 0.252  0.034 0.536 ITA
## 18 0.214   0.011 0.070 0.009 0.144 0.013 0.325  0.040 0.512 JPN
## 19 0.107   0.008 0.035 0.006 0.073 0.007 0.324  0.047 0.405 KOR
## 20 0.220   0.019 0.100 0.008 0.120 0.017 0.453  0.038 0.314 LUX
## 21 0.204   0.018 0.071 0.009 0.134 0.011 0.346  0.027 0.457 MEX
## 22 0.177   0.012 0.122 0.008 0.056 0.012 0.686  0.053 0.575 NLD
## 23 0.177   0.011 0.050 0.005 0.127 0.010 0.285  0.026 0.078 NOR
## 24 0.165   0.010 0.078 0.005 0.088 0.010 0.470  0.037 0.162 NZL
## 25 0.226   0.010 0.075 0.006 0.151 0.010 0.333  0.027 0.142 POL
## 26 0.196   0.011 0.089 0.009 0.107 0.009 0.452  0.037 0.342 PRT
## 27 0.267   0.010 0.122 0.008 0.145 0.010 0.459  0.029 0.446 SVK
## 28 0.174   0.010 0.055 0.005 0.119 0.010 0.318  0.031 0.098 SWE
## 29 0.120   0.012 0.062 0.009 0.058 0.008 0.517  0.052 0.549 TUR
## 30 0.159   0.006 0.042 0.004 0.117 0.007 0.264  0.025 0.222 USA

Plots and Correlations between Model Results and Country ICCs

# Preperation of Results Corr between total effect and ICC by country
sapply(resultsList, function(x) cor(x[, 1], x[, 9], use = "complete.obs"))
##  HISEI HISCED   ESCS 
## 0.4283 0.2562 0.4079
# Corr between proportion indirect and CNT
sapply(resultsList, function(x) cor(x[, 7], x[, 9], use = "complete.obs"))
##  HISEI HISCED   ESCS 
## 0.4051 0.6211 0.4719
# Plots for weights with country specific standadization
results <- resultsList[[1]]
plot(results[, 9], results[, 1], ylab = "Total Effect", xlab = "ICC", ylim = c(0, 
    0.3), pch = "", main = "Relationship between SES and Educational Expectations")
text(results[, 9], results[, 1], results[, 10], cex = 0.6)
abline(lm(results[, 1] ~ results[, 9]), col = "grey")
arrows(x0 = results[, 9], y0 = results[, 1] - results[, 2], x1 = results[, 9], 
    y1 = results[, 1] + results[, 2], length = 0, col = "grey")

plot of chunk unnamed-chunk-8


# Indirect Effect
plot(results[, 9], results[, 7], ylab = "Proportion of Total Effect ", xlab = "ICC", 
    ylim = c(0, 0.8), pch = "", main = "Relationship between SES and Educational Expectations \n Explained by Achievement Differentials")
text(results[, 9], results[, 7], results[, 10], cex = 0.6)
abline(lm(results[, 7] ~ results[, 9]), col = "grey")
arrows(x0 = results[, 9], y0 = results[, 7] - results[, 8], x1 = results[, 9], 
    y1 = results[, 7] + results[, 8], length = 0, col = "grey")

plot of chunk unnamed-chunk-8

Plots with highest parent education rather than highest parent occupational prestige

results <- resultsList[[2]]
plot(results[, 9], results[, 1], ylab = "Total Effect", xlab = "ICC", ylim = c(0, 
    0.3), pch = "", main = "Relationship between SES and Educational Expectations")
text(results[, 9], results[, 1], results[, 10], cex = 0.6)
abline(lm(results[, 1] ~ results[, 9]), col = "grey")
arrows(x0 = results[, 9], y0 = results[, 1] - results[, 2], x1 = results[, 9], 
    y1 = results[, 1] + results[, 2], length = 0, col = "grey")

plot of chunk unnamed-chunk-9


# Indirect Effect
plot(results[, 9], results[, 7], ylab = "Proportion of Total Effect ", xlab = "ICC", 
    ylim = c(0, 0.8), pch = "", main = "Relationship between SES and Educational Expectations \n Explained by Achievement Differentials")
text(results[, 9], results[, 7], results[, 10], cex = 0.6)
abline(lm(results[, 7] ~ results[, 9]), col = "grey")
arrows(x0 = results[, 9], y0 = results[, 7] - results[, 8], x1 = results[, 9], 
    y1 = results[, 7] + results[, 8], length = 0, col = "grey")

plot of chunk unnamed-chunk-9

Plots with PISA INDEX OF ECONOMIC, SOCIAL AND CULTURAL STATUS (ESCS)

results <- resultsList[[3]]
plot(results[, 9], results[, 1], ylab = "Total Effect", xlab = "ICC", ylim = c(0, 
    0.3), pch = "", main = "Relationship between SES and Educational Expectations")
text(results[, 9], results[, 1], results[, 10], cex = 0.6)
abline(lm(results[, 1] ~ results[, 9]), col = "grey")
arrows(x0 = results[, 9], y0 = results[, 1] - results[, 2], x1 = results[, 9], 
    y1 = results[, 1] + results[, 2], length = 0, col = "grey")

plot of chunk unnamed-chunk-10


# Indirect Effect
plot(results[, 9], results[, 7], ylab = "Proportion of Total Effect ", xlab = "ICC", 
    ylim = c(0, 0.8), pch = "", main = "Relationship between SES and Educational Expectations \n Explained by Achievement Differentials")
text(results[, 9], results[, 7], results[, 10], cex = 0.6)
abline(lm(results[, 7] ~ results[, 9]), col = "grey")
arrows(x0 = results[, 9], y0 = results[, 7] - results[, 8], x1 = results[, 9], 
    y1 = results[, 7] + results[, 8], length = 0, col = "grey")

plot of chunk unnamed-chunk-10

Relationships Excluding Outliers

In the models it was clear that Italy was an outliers. While we report on the sample of countries as a whole here, we also provide analysis without Italy here .

# Deletion of Italy outlier. Corr between proportion indirect and CNT
sapply(resultsList, function(x) cor(x[!c(x$CNT == "ITA"), 7], x[!c(x$CNT == 
    "ITA"), 9]))
##  HISEI HISCED   ESCS 
## 0.5282 0.6802 0.5632