arXiv:1802.08004v1 [math.ST] 22 Feb 2018

The use of sampling weights in the M-quantile
random-effects regression: an application to
PISA mathematics scores
Francesco Schirripa Spagnolo∗ , Nicola Salvati† , Antonella
D’Agostino‡ , and Ides Nicaise§
∗

Dipartimento di Studi Economici e Giuridici, Università di Napoli “Parthenope”
†

‡

Dipartimento di Economia e Management, Università di Pisa

Dipartimento di Studi Aziendali e Quantitativi, Università di Napoli “Parthenope”
§

KU Leuven

Abstract
M-quantile random-effects regression represents an interesting approach
for modelling multilevel data when the interest of researchers is focused
on the conditional quantiles. When data are based on complex survey
designs, sampling weights have to be incorporate in the analysis. A
pseudo-likelihood approach for accommodating sampling weights in the
M-quantile random-effects regression is presented. The proposed methodology is applied to the Italian sample of the “Program for International
Student Assessment 2015” survey in order to study the gender gap in
mathematics at various quantiles of the conditional distribution. Findings offer a possible explanation of the low share of females in “Science,
Technology, Engineering and Mathematics” sectors.
Keywords: Multilevel modelling; M-quantile; Pseudo-likelihood; Robust statistics; Sampling weights

1

Background

Gender differences in educational outcomes have raised concern over the last
few decades because several studies stressed that the gender gap in education
contributes to the gender segregation in the labour market and to the gender differences in wages (Brown and Corcoran, 1997; Duquet et al., 2010; World Bank,
2011).
Especially, the severe female disadvantage in mathematics is of particular
importance for its snowball effect. Although in many western countries the
share of women that participate in higher education is now large, they yet tend
to choose the STEM (Science, Technology, Engineering, and Mathematics) sectors much less frequently than men (Bradley, 2000; van Langen et al., 2006;
European Commission, 2016). Consequently, the low share of women in STEM
education is one of the reasons of the poor female presence in research and innovation decision-making positions (van Langen et al., 2006; European Commission,
1

2016). Since both research and innovation are important elements for a smarter,
more sustainable and inclusive growth, the European Commission (2012) also
stressed that the female presence among scientific decision makers is fundamental for increasing the quality of research and improving the acceptance of
innovation in the labour market.
Fig. 1 highlights the negative relationship between the female-male ratio of
people being employed in science and technology and the gap in mathematics in
secondary school measured using data of the “Program for International Student
Assessment (PISA) 2015” survey. As shown in Fig. 1, Italy is one of the
countries with the highest gender gap in mathematics and the lowest share of
female being employed in Science and Technology.

Figure 1: Plot of gender gap in mathematics and ratio of female and male human
resources in science and technology aged 25-64 in 26 OECD EU Countries (the
data about the Human Resources in Science & Technology refer to 2016. Source:
Eurostat)
Despite the gender differences in mathematics are much higher than the
OECD average (a difference of 20 points against an OECD-average difference
of 9 points) and the share of girls that does not achieve the minimum level of
competencies in mathematics is one of highest among the OECD countries, in
Italy the scores in mathematics are similar to the OECD average (Italy has an
average test score of 490, which is not significantly different from the OECD
mean) (Invalsi, 2016).
For this reason, the gender differences in mathematics have been studied
in several papers: among others, Bratti et al. (2007) used multilevel approach
for analysing the 2003 Italian wave of the PISA and they found a large gender
gap in mathematical performance; Agasisti et al. (2014) found a female penalty
on mathematics achievement stronger in Southern than Northern regions of
students attending grade 6 in the year 2011/12; Masci et al. (2016) fitted a
bivariate multilevel linear model for students attending the first year of junior
2

secondary school in the year 2012/2013 and their results highlighted that on
average males are better in mathematics, while females are better in reading;
Grilli et al. (2016) analysed the combined dataset of International Association
for the Evaluation of Educational Achievement (IEA) assessments of student
achievement in mathematics and science (TIMSS) and in reading (PIRLS) and
they found that females have a lower performance in mathematics and science,
but not in reading.
However, all these studies are also based on standard multilevel modelling
approaches and accordingly findings are valid mainly around the centre of the
distribution (Goldstein, 2011). As a consequence, this can lead to an incomplete overview of the gender gap in mathematics, as large gender gap in mathematics among high-performing students can be more serious than among lowperformers, since the best students tend to enroll more into higher education
(AlmaDiploma, 2017).
The report of Invalsi (2016) showed, based on simple raw-data, how the gender gap in mathematics really tends to be wider in the upper quantiles of the
distribution. In this perspective, what is currently lacking in the literature is
the analysis of the gender gap along the overall conditional distribution of mathematics scores, i.e. taking into account other relevant variables (e.g. students
and school-level characteristics) that can mitigate the gender effect in different
quantiles.
In other words the aim of this paper is to estimate and understand the
gender gap in math, using a distributional approach, which makes comparisons
over the entire range of mathematics scores rather than focusing on summary
measures such as the mean. Such analysis is obviously relevant as it might give
insights to policy makers for implementing policies to reduce the gap and to try
and raise the number of girls enrolled in scientific graduate programmes.
From a methodological point of view, quantile regression is the most popular approach for estimating the conditional quantiles of a response variable
(Koenker and Bassett, 1978; Koenker, 2005), but the classical implementation
of this estimation method does not allow to include into the analysis specific
random effects to take into account the hierarchical data structure. Actually, some attempts to implement quantile regression with a hierarchical data
structure have been addressed in literature. For instance, Koenker (2004) proposed a penalized fixed-effects estimation method for longitudinal data, while
Geraci and Bottai (2007, 2014) proposed a two-levels quantile regression with
mixed effects named the Linear Quantile Mixed Model (LQMM). It allows to
model data with complex dependence structures by including multiple random
effects in the linear conditional quantile functions. For a review of linear quantile
models for multilevel data see (Marino and Farcomeni, 2015). In the context of
educational research, there are not many examples of studies that applied multilevel quantile regression models. Costanzo (2015) used the LQMM to evaluate
the effect of a specific training programme (M@tabel) on the Italian sixth grade
students’ performance in mathematics at secondary schools: the author highlighted the advantages of using this model compared to the traditional linear
random effects one. Faria and Portela (2016) applied LQMM to analyse the
determinants of students’ success in Portugal using 2009 PISA survey.
A more flexible alternative to quantile regression is also the the M-quantile
(MQ) regression (Breckling and Chambers, 1988; Chambers and Tzavidis, 2006).
This approach integrates the concepts of quantile regression and expectile re3

gression (Newey and Powell, 1987) within a framework defined by a ‘quantilelike’ generalization of regression based on influence functions (M-regression).
The classical MQ regression even does not allow the analysis of multilevel data.
Recently important developments have been made in order to overcome this limitation. Tzavidis et al. (2016) extended M-quantile regression so as to include
random effects in two-level hierarchical data structure using maximum likelihood (MQRE-2L); Borgoni et al. (2016) extended the MQRE-2L to a three-level
random effects model (MQRE-3L); Alfò et al. (2017) define a finite mixture of
quantile and M-quantile regression models (FMMQ) for heterogeneous and /or
for dependent/clustered data. Although LQMM, MQRE and FMMQ allow for
modelling data with a complex dependence structure including random effects,
they do not accommodate for sampling weights in the estimation procedure.
Thus, they are not suitable for modelling survey data based on complex survey designs (e.g. when statistical units are selected with unequal probabilities
from the population) usually implemented by national and international organizations. Namely, any estimation procedure that does not take into account
sampling weights provides biased results (Cochran, 1977; Sarndal et al., 1992;
Pfeffermann, 1993).
While weighting for unequal probabilities of selection is a relatively wellestablished procedure in single level models as it can be viewed as an application of the ‘Pseudo Maximum Likelihood’ (PML) approach (Chambers, 2003;
Binder, 1983), it needs additional remarks in multilevel models because the
covariance structure of the population has to be modelled. For this reason,
the multilevel models require the knowledge of the inclusion probabilities at
each hierarchical level. Pfeffermann et al. (1998) used the PML approach and
the iterative generalized least square (IGLS) algorithm (Goldstein, 1986) in a
two-level model with a continuous response variable, obtaining the so-called
probability-weighted IGLS (PWIGLS) estimators.
Alternative approaches have been proposed by Grilli and Pratesi (2004);
Rabe-Hesketh and Skrondal (2006); Asparouhov (2006). These approaches include sampling weights in the estimation procedure. In particular, the census likelihood is estimated by weighting the sample likelihood including sampling weights in the log-likelihood function. Indeed, Rabe-Hesketh and Skrondal
(2006) proposed a pseudo-likelihood approach via adaptive quadrature for generalized linear models with dichotomous response with any number of levels;
their method also allows for stratification and PSUs that are not represented
by a random effect in the model. Grilli and Pratesi (2004) used the PML approach to develop weighted estimators in the context of ordinal and binary models. Finally, Asparouhov (2006) proposed an approximately unbiased multi-level
pseudo maximum likelihood (MPML) for general multi-level modelling and asserted that it can be used with any parametric family of distributions and any
linear or non-linear multilevel models.
The proposed multilevel approaches based on the use of the sampling weights
suffer from the same limitations as the the standard multilevel models: they
allow for modelling the centre of the distribution. This limitation can lead
to an incomplete overview of the gender gap in mathematics. Therefore, the
purpose of this paper is to extend the current modelling framework of M-quantile
random effects so that to include sampling weights and to examine the shape
of gender gap along the overall conditional distribution of mathematics scores
using PISA-OECD 2015 data.
4

The structure of the paper is as follows: after this introduction, Section 2
provides details of PISA data, including variables used in the empirical analysis.
Section 3 provides details of the methodology proposed. In Section 4 we present
the results of applying the new methodology to PISA data for studying gender
gap in mathematics, together with details of the method that is used for model
evaluation. In Section 5 we empirically evaluate the properties of the proposed
model by using a Monte Carlo simulation study. Finally, Section 6 provides a
concluding summary and a discussion of potential areas for future research.

2

Data and sample characteristics

The empirical analyses are based on Programme for International Students Assessment (PISA) survey that is a triennial survey started in 2000 and conducted
by the OECD for measuring the extent to which the 15-year-old students have
acquired key knowledge and skills in three school subjects (mathematics, science and reading) that “are essential for full participation in modern societies”
(OECD, 2016a) in many countries and economies. In particular, our data refer
to the latest 2015 PISA wave.
PISA uses a two-stage stratified sample design. Schools having 15-year-old
students are the first-stage sampling units. A systematic Probability Proportional to Size (PPS) sampling is used to select schools. In order to improve the
precision of sample-based estimates, schools previously are divided into explicit
strata based on school characteristics. Students within sampled schools are the
second-stage sampling units. Although the two-stage sampling design used in
PISA should guarantee that the students have the same probability of selection, three factors contribute in explaining the variability of weights: i) overor under-sampling of some strata of the population; ii) lack of accuracy or no
update in defining the size of school; iii) adjustment of weights for school and
student non response. Further details on sample design and weights defined
in PISA can be found in OECD (2009, cap. 3); OECD (2017, cap. 4 and 8);
Rabe-Hesketh and Skrondal (2006).
Although the target population of PISA survey is represented by 15-years-old
students who have completed at least 6 years of formal schooling, our analysis is
based only on students enrolled in public upper secondary schools, corresponding to the 3th level of the International Standard Classification of Education
(ISCED). This choice is determined by the fact that upper secondary education
represents a sort of minimum credential for entering successfully into the labour
market and it is required to pursue their university education. For this reason,
the upper secondary level of education is considered a crucial indicator of the
output of the educational system of a Country (OECD, 2003). Our estimation
sample contains 7163 students with non-missing response variable and covariate
information. The 7163 students are clustered in 283 schools.
The response variable that we use in this paper is the students’ scores in
mathematics. In particular, we use the first plausible value generated by the
standardized procedure implemented by the PISA team (for more detail about
the methodology see OECD, 2017).
The covariates used in the models are as follows. Gender (male (baseline)
or female) whose effect is the main object of this study. Immigration status was
coded as a dummy variable (OECD, 2016a). Namely, non immigrant students

5

(students whose mother or father (or both) were born in the country where
the PISA test occurs, regardless of whether the student himself or herself was
born in that country – baseline) vs immigrant students (students whose mother
and father were both born in a country other than that where the student
performs the PISA test). The native-immigrant gap is well recognized in the
literature (see for instance Rangvid, 2007b; Ammermueller, 2007; Azzolini et al.,
2012; Hajisoteriou and Panayiotis, 2016) because usually immigrants encounter
many difficulties in achieving adequate mathematical skills, consequently they
cannot participate successfully and actively in the host society. The index of
schoolwork-related anxiety was measured by the responses of the students to the
statements regarding their worry about tests and study in general. It was standardised to have a mean of 0 and a standard deviation of 1 across OECD countries. Positive values on the index indicate that students reported higher levels
of schoolwork-related anxiety than the average student across OECD countries;
negative values indicate that students reported lower levels of anxiety than the
average student (OECD, 2016c). Grade repetition is a good indicator of the
students’ school career and it can affect negatively the academic performance
and students’ delayed entry into the labour market (OECD, 2016b, 2013; see
also Ikeda and García (2014) for a review of studies about this topic). It is
recoded into two categories: the student has never repeated a grade in any level
(baseline) and the student has repeated a grade in at least one level. Lack of
punctuality was measured by students’ answers on whether they had arrived late
for school in the two weeks before the test (no (baseline) or yes). This indicator
measures the students’ truancy and, at the same time, it may represent the lack
of interest for learning, moreover it has negative consequences on other students
because it can contribute to a disruptive learning environment. Student socioeconomic status, measured by the PISA index of economic, social and cultural
status (ESCS). It was derived from the combination of three variables related
to family background: highest parental education, highest parental occupation,
and home possessions (as a proxy of family wealth). Higher values of ESCS
indicate better socio-economic status (OECD, 2016c). This variable is strongly
associated with the students’ achievement because students with a high ESCS
have better access to educational resources provided by their family (material
resources and the educational level of their parents).
In addition to individual-level characteristics, the key explanatory variables
at school-level are as follows. The macro region that indicates where the school
is located (Northern regions (baseline) or Southern regions). Regional disparities in school achievement are indeed, generally very serious in Italy: students
attending the schools located in Northern Italy tend to have higher achievement
scores than their counterparts in Southern Italy (Invalsi, 2013, 2016). The type
of school is recoded into three categories: Lyceums (baseline), Technical and Vocational Schools and Other Schools. This classification allows to discriminate
between three groups of schools that offer different potential outcomes. Indeed,
Lyceums provide mostly a theoretical training and prepare students for tertiary
education (ISCED 5 and 6); Technical and Vocational Schools offer both a general education and a technical specialization in a specific field of study; other
schools are institutions designed to offer technical activities and allow students
to obtain a professional qualification which is immediately recognised in the
labour market.
Moreover, the school-mean of each individual-level variable has been calcu6

Table 1: Descriptive statistics of outcome variable and covariates
Variable
MATHEMATICS
ESCS
Anxiety index
Female
Immigrants
Grade repetition
Lack of punctuality
South Italy
Type of School
Lyceums
Technical and Vocational
Others
Mean gender
Mean Immigrant status
Mean Grade repetition
Mean Lack of punctuality
Mean ESCS
Number of units
Number of clusters

Min

1st Q

Median

Mean

3rd Q

Max.

140.80
-2.99
-2.51

432.14
-0.71
-0.08

498.17
-0.04
0.52

496.34
-0.05
0.48
0.51
0.07
0.13
0.35
0.44

560.72
0.67
1.05

822.64
3.56
2.55

0.15

1.07

-1.50
7163
283

-0.59

-0.27

0.42
0.54
0.04
0.50
0.09
0.16
0.36
-0.20

lated for estimating the so called “compositional effect” (Sani and Grilli, 2011).
These variables are indeed proxies of the social context of the school that
might considerably influence the academic achievements because of peer interactions. Moreover, on the basis of school composition teachers can change
their instructional practices to take into account the characteristics of school
and consequently this phenomenon can affect the individual students’ performance (Rangvid, 2007a; Schneeweis and Winter-Ebmer, 2007). Therefore, how
the social context of the school affects the achievement of pupils is a relevant
element for policy makers in designing measures to control the issue of social
composition of the schools. Table 1 describes the sample by reporting sample
proportions for binary variables, means, medians and quantiles for continuous
variables.
There is an approximately even split of male and female students, only the
7% of the students’ population has a foreign background in line with the foreign
citizens who live in Italy (cfr. ISTAT, 2017), a little bit more than a third
of students reported to arrive late at school and only the 13% of them have
repeated a grade. The average score in mathematics is almost 450 and about
75% of all students exceed the level 2 of proficiency (higher than 420.07 and
less than or equal to 482.38), in line with the OECD mean (cfr. Invalsi, 2016).
The inter-quantile range is 129 a little more than OECD mean difference (125
points). More than 50% of all schools are located in Northern/Central regions
and more than half (54%) are Technical and Vocational Schools.
Finally, to support the choice of using a robust approach a simple weighted
two-level random intercept model for mathematics score has been estimated
using the student-level and the school-level weights provided in PISA database
and the covariates described above.
Fig. 2 shows the residual analysis. Clearly, several outliers and second-level
units that can be classified as influential points (filled triangle points exceed the
cut-off equal to 1) stand out.
7

(a)

(b)

Figure 2: Plot of standardised residuals (a) and Cook’s Distance (b) derived by
fitting a linear mixed-effects model for mathematics scores

3

Methodology

In this section, after a briefly reviewing of the PML in multilevel models (Section
3.1), we present the weighted two-level MQRE model (Section 3.2).

8

3.1

Pseudo-maximum likelihood approach in multilevel model

Let us have a two-level population with M level 2 units (or cluster) and Nj level
1 units within the jth cluster. Suppose that an outcome continuous variable y
is modelled using a two-level model:
yij = xTij β + zTij γ j + εij ,

i = 1, ..., Nj ,

j = 1, ..., M,

(1)

where xTij is a covariate-vector of dimension p; β is a vector of regression coefficients; zij is a vector of group indicators used to define the random part of
the model; γ j are random effects varying over clusters and we assume that they
follow a multivariate normal distribution with zero mean and covariance matrix
Ω and εij are individual random effects and we assume εij ∼ N (0, σε2 ). In the
paper we consider a two-level random-intercept model for unit i in cluster j. It
can be written as
yij = xTij β + γj + εij ,
(2)
where γj ∼ N (0, σγ2 ). When population model in Eq. 2 is evaluated at the sample, estimation of the parameters of the model can be obtained by employing
maximum likelihood method. Under the normal assumption and εij independent from γj the log-likelihood function is:
1
1
logL(β, σγ2 , σε2 ) = − log|V| − (y − Xβ)T V−1 (y − Xβ),
2
2

(3)

where y i the n × 1 response vector and n is the sample size; V = Σε + ZΣγ ZT ,
Σε = σε2 In , Σγ = σγ2 IM and Z is an n × M matrix of known positive constant.
Here Ik represents the identity matrix of dimension k. Estimates of fixed and
variance parameters can be obtained by differentiating the log-likelihood with
respect to these parameters and then solving the estimating equations defined
by setting these derivatives equal to zero. However, if the sample data are
generated by a complex sample design the different inclusion probabilities have
to take into account and the weights have to be inserted somewhere in the
estimation process.
This idea leads to use the weights within regression analysis to estimate the
model on the entire population (a census) rather than the sample. The PML
approach in multilevel framework requires to write down the census-likelihood
function and then incorporates weights at each level of the analysis in the loglikelihood equations.
The census-likelihood function can not be expressed as a simple sum of
the elementary unit contributions, but it is a sum across all levels of the data
hierarchy. In a two level model it is a function of sums across level 2 and level
1 units:



Z
Nj
M

X
X
log Lij (β, σγ2 , σε2 |γ j )  f (γ j )dγ j , (4)
log exp
log L(β, σγ2 , σε2 ) =


j=1

i=1

where f (γ j ) is the normal probability density of the level 2 random effects and
log Lij (β, σγ2 , σε2 |γ j ) is the log-likelihood contribution of the level 1 units conditioned on the level 2 random effects.
Now, suppose that the whole population is not observed, but the sample
data have been selected with the following two-stage sampling design: at the
9

first stage m clusters are selected with inclusion probabilities πj , j = 1, ...M ,
and at the sub-sequent stage elementary units are sampled with conditional
probabilities πi|j , i = 1, ...Nj . The log-likelihood function for the sample units
can be expressed as



Z

X
X
wi|j log Lij (β, σγ2 , σε2 |γ j )  f (γ j )dγ j ,
wj log exp
log Lw (β, σγ2 , σε2 ) =


i∈sj

j∈s

(5)
where s and sj indicate the sampled clusters and the units sampled in the
cluster j, respectively. As we noted, the sum for each of the levels requires the
respective conditional probabilities of selection: the log-likelihood contributions
of the level-1 units are weighted by wi|j = 1/πi|j , i.e. the inverse of the selection
probability of unit i in cluster j given that cluster j has been sampled, and the
log-likelihood contributions of the level 2 units are weighted by wj = 1/πj , i.e.
the inverse of the selection probabilities of cluster j.
In case of two-level random-intercept model the Eq. 5 can be expressed as:
log Lw (β, σγ2 , σε2 ) = −


1X 
wj log(2π|W−1 V|) + (y − Xβ)T V−1 (y − Xβ) ,
2 j∈s

(6)
where W is a diagonal matrix of the first-level sampling weights, W = diag(w1|j ,
..., wnj |j ), and V = W−1 Σε + ZΣγ ZT .
Differentiating Eq. 6 with respect to the fixed effects and to the variance
components the weighted estimation equations of the LMM are obtained.

3.2

Weighted M -quantile Random Effects Model

As we explained in Section 1, the classical multilevel regression model in Eq.
1 provides an incomplete picture of the distribution of the response variable
given the auxiliary information, because it just summarises the behaviour of
the mean of the outcome variable at each point in a set of covariates. Moreover,
the presence of outliers in the data invalidates the general assumptions of the
model; in this case the estimators of the parameters of the model under Eq. 3
could be biased and inefficient (Richardson and Welsh, 1995). Huggins (1993)
and Richardson and Welsh (1995) proposed an approach based on M-estimation
for a robust estimation of the multilevel models. Tzavidis et al. (2016) used
this idea of robust estimation in the context of multilevel model to extend Mquantile regression in order to include random effects that account for a twolevel hierarchical structure in the data and used maximum likelihood technique
to estimate the parameters of the model.
In the linear case the M-quantile regression assumes the following form:
M Qy (q|x; ψ) = xTi βψq ,

(7)

where β ψq is defined as the minimiser of
min E [|q − I(u < 0)|ρ(u)] ,
β

10

(8)

with u = (y − xT β)/σq , σq a scale parameter. In this paper we use the popular
Huber loss function (Breckling and Chambers, 1988):

2c|u| − c2
|u| > c,
ρ(u) =
(9)
u2
|u| ≤ c,
where c is a tuning constant bounded away from zero (a common choice is c
= 1.345). Assuming that ρ is (a.e.) continuously differentiable and convex, an
estimator of β ψq , β̂ ψq can be obtained as the solution of the following system
of equations
!
n
X
yi − xTi βψq
,
(10)
ψq
σ̂q
i=1

where σ̂q is a consistent estimator of σq , ψq (u) denotes an asymmetric influence
function, which is the derivative of an asymmetric loss function ρq (u) = |q −
I(u < 0)|ρ(u).
Defining the loss function as in Eq. 9, MQ regression model estimation is
obtained weighting positive residuals by q and negative residuals by 1 − q. The
Iterative Weighted Least Square (IWLS) algorithm used to fit an M-quantile
regression model guarantees convergence to a unique solution when a continuous monotone influence function (such as Huber function with c > 0) is used
(Kokic et al., 1997). An other interesting benefit is the possibility to trade efficiency for robustness by setting the value of c. For example, if c is close to
zero, the robustness of the model increases and its efficiency decreases, then
MQ moves towards quantile regression. On the other hand, setting c large, the
robustness of the model decreases and its efficiency increases, then MQ moves
towards expectile regression.
For modelling multilevel structured data Tzavidis et al. (2016) extended the
linear specification of the model in Eq. 7 to allow for the inclusion of random
effects (MQRE-2L) to account for a two-levels hierarchical structure in the data.
The M-quantile random intercept model is defined as follows:
M Qyij (q|xij , γj ; ψ) = xTij β ψq + γj ,

(11)

where γj is the random effect for cluster j. For obtaining the estimation equations for the regression coefficients and the variance parameters the authors
extended the idea of asymmetric weighting of residuals by changing the estimating equation of the ML proposal II by Richardson and Welsh (1995) following
Sinha and Rao (2009).
Borgoni et al. (2016) extended the MQRE-2L to a three-level M-quantile
random effects regression (MQRE-3L). The main difference with the model
proposed by Tzavidis et al. (2016) is an additional equation for estimating the
further level and a more complex variance-covariance matrix of the outcome
variable.
The aim of this paper is to extend the idea of Tzavidis et al. (2016) and
incorporate sampling weights in the model. So, following the pseudo-likelihood
approach we propose an estimation procedure that adjusts for the effect of an
informative sampling design on estimation in MQRE-2L.
In particular our approach follows the works of Grilli and Pratesi (2004),
Asparouhov (2006) and Rabe-Hesketh and Skrondal (2006). In particular, the

11

sampling weights are inserted before the derivatives of the log-likelihood function are taken. Consequently, firstly, the weighted log-likelihood in Eq. 6 is
differentiated obtaining the weighted estimation equation for a the parameter
of the random intercept model and then, following the idea of Tzavidis et al.
(2016), we propose a robustification of the weighted estimation equations for
estimating the regression coefficients and the variance parameters obtaining in
the previous step as follows:
i
X h
−1 1/2
(12)
Uqj ψq (rqj ) = 0,
wj XTj Vqj
j∈s

i
1X h
1/2 −1
−1 1/2
−1
−
Zj ZTj Vqj
Uqj ψq (rqj ) = 0,
Zj ZTj ) − ψq (rqj )T Uqj Vqj
wj K2qj tr(Vqj
2 j∈s
−

i
1X h
1/2 −1
1/2 −1
−1
ψq (rqj ) = 0,
Wj−1 Uqj Vqj
Wj−1 ) − ψq (rqj )T Uqj Vqj
wj K2qj tr(Vqj
2 j∈s

(13)
− Xj βψq )
where Xj is the matrix of covariates in the cluster j; rqj =
is a vector of scaled residuals with components rijq ; Uqj is a diagonal matrix
with diagonal elements equal to the diagonal elements of the variance-covariance
matrix of the jth cluster Vqj = Wj−1 Σεq + Zj σγ2q ZTj ; Σεq = σε2q Inj , σγ2q and σε2q
are the M-quantile specific variance parameters; Zj is an nj × 1 vector of known
positive constants; K2qj = E[ψq (ej )ψq (ej )T ] with ej ∼ N (0, Inj ). Eq. 12 and
13 are the estimating equations of the Weighted M-quantile Random Effects
Model (henceforth Weighted-MQRE). To obtain estimators of βψq , σγ2q and σε2q
Eq. 12 and 13 are solved iteratively. For Eq. 12 a Newton-Raphson algorithm is
used and for Eq. 13 the fixed-point iterative method is implemented for getting
the estimates. Details about the algorithm can be found in the Appendix A.
Inference for the model parameters is performed using a sandwich estimator
following the proposal by Tzavidis et al. (2016) (for more details see Appendix
B).
As pointed out in Tzavidis et al. (2016) and Borgoni et al. (2016) and discussed in detail in Jones (1994), we can consider M-quantiles equivalent to
quantiles because both target the same part of the distribution. So, the estimates of β ψq can be interpreted how the effect of a unit change in x on a given
q-quantile of the distribution of y.
−1/2
Uqj (yj

3.3

Design consistency of the Weighted-MQRE regression
coefficients

In this section we prove the design consistency of β̂ψq for the regression coefficients β ψq of the Weighted-MQRE. The asymptotic design-based properties of
the regression coefficients in the M-quantile framework have been introduced by
Fabrizi et al. (2014) and the consistency of estimators in the multilevel modelling has been demonstrated by Pfeffermann et al. (1998). From a theoretical
point of view, the establishment of consistency properties in multilevel modelling
requires that the number of groups, m, and the units in each group, nj , increase;
however as pointed out in Pfeffermann et al. (1998) in practice nj are often
small and the consistency of β̂ ψq can be established when only m increases and

12

tends to M . Following the work of Fabrizi et al. (2014) and Wang and Opsomer
(2011) we define, for any q ∈ (0, 1), the population parameter Bψq as:
Bψq = inf {βψq : SM (β ψq ) ≥ 0},

(14)

PM
−1/2
ψq (Vj−1 (yj − xTj βψq )).
where SM (β ψq ) = M −1 j=1 Xj Vj
The population parameter Bψq can be estimated by
β̂ ψq = inf {βψq : ŜM (βψq ) ≥ 0},

(15)

Pm
−1/2
ψq (Vj−1 (yj − xTj β ψq )).
with ŜM (β ψq ) = M −1 j=1 wj Xj Vj
Assuming the technical conditions A1-A5 (described in Appendix C), for
each q ∈ (0, 1) and for large closed interval Θ ∈ Rp , supβψq ∈Θ |ŜM (β ψq ) −
W−MQRE

SM (β ψq )| = op (1), the estimator β̂ψq
obtained solved the Eq. 12 is design
√
W−MQRE
−1
m-consistent for Bψq , in the sense that M (β̂ ψq
−Bψq ) = Op (m−1/2 ).
Proof of this statement follows the results in Wang and Opsomer (2011) and
Fabrizi et al. (2014).

4
4.1

Application: Modelling the conditional distribution of mathematics scores
Modelling approach

In this section we apply the method proposed in the Section 3.2 to study the
effect of gender along the entire distribution of mathematics scores using PISAOECD 2015 data. In particular, taking into account the socio-economic issues
and the literature review that were described in Section 1, we are mainly interested in measuring the gender gap in mathematics, controlling for individual
and school-level characteristics. In particular, adding the cluster mean of the
individual variables we can estimate the so-called within-effect. A specific function in the statistical programming environment R (R Core Team, 2015) has
been written and it is available from the authors upon request.
The following scaled version of first-level weights named “scaling method 2”
in Pfeffermann et al. (1998) has been used in order to fit mathematics scores:
s2
wi|j

= nj wi|j

nj
X
i

wi|j

−1

,

where nj is the number of sample units in the jth cluster and wi|j are the firstlevel weights. Pfeffermann et al. (1998) have proposed
an P
alternative scaling
Pn
n
2 −1
s1
) . It is
method, named “scaling method 1”: wi|j
= wi|j ( i j wi|j )( i j wi|j
worth noting that the results using this alternative scaling method are almost
identical to those using “scaling method 2”. However, simulation results in
Pfeffermann et al. (1998) suggest that “scaling method 2” works better than
“scaling method 1” for informative weights. Moreover, the student-level weights
provided in PISA (the variable called W_FSTUWT) are student-level overall
inclusion weights (wij not wi|j ) adjusted for non-inclusion and non-participation
of students. Consequently scaling the weights is a procedure to overcome the
problem of not having the conditional sampling weights. Indeed, if student-level
13

Table 2: Results–Weighted-MQRE model for mathematics scores†
Results for the following value of q:
Variable

0.1

0.25

Intercept
Gender
−14.38∗∗∗ −14.98∗∗∗
Immigrant status
−14.58∗∗ −15.53∗∗
Grade repetition
−45.21∗∗∗ −45.21∗∗∗
Anxiety index
−5.29∗∗
−5.82∗∗∗
∗∗∗
Lack of punctuality
−12.61
−9.89∗∗∗
ESCS
2.51
2.39
Mean gender
−48.93∗∗∗ −51.64∗∗∗
Mean immigrant status
−137.76∗∗∗ −114.99∗∗
Mean grade repetition
−51.42∗
−54.71∗
Mean lack of punctuality
−41.53∗
−45.16∗∗
Mean ESCS
41.93∗∗∗
36.73∗∗∗
School type (ref. Lyceums)
Technical and Vocational
−27.54∗∗∗ −32.69∗∗∗
Others
−90.07∗∗∗ −106.40∗∗∗
South Italy
−39.62∗∗∗ −43.06∗∗∗
533.20∗∗∗

570.63∗∗∗

†Point estimates with associated p-value:

∗∗∗

0.5

0.75

608.38∗∗∗

648.20∗∗∗

−17.08∗∗∗
−14.95∗∗∗
−42.24∗∗∗
−6.74∗∗∗
−9.17∗∗∗
2.39
−53.48∗∗∗
−81.98∗∗∗
−62.36∗∗
−47.70∗∗
32.41∗∗∗

−17.90∗∗∗
−15.19∗∗∗
−38.31∗∗∗
−8.20∗∗∗
−13.82∗∗∗
3.72
−58.47∗∗∗
−78.28∗∗∗
−63.36∗∗
−47.87∗∗
28.23∗∗∗

0.9
688.40∗∗∗
−16.95∗∗∗
−18.74∗∗∗
−35.79∗∗∗
−9.66∗∗∗
−16.64∗∗∗
6.07∗∗
−66.52∗∗∗
−80.52∗∗∗
−67.51∗∗∗
−57.15∗∗∗
25.49∗∗

−35.85∗∗∗ −37.85∗∗∗ −38.18∗∗∗
−113.44∗∗∗ −106.16∗∗∗ −95.73∗∗∗
−43.21∗∗∗ −41.96∗∗∗ −42.52∗∗∗

p < 0.01;

∗∗

p < 0.05;

∗

p < 0.1

weights are rescaled, the model estimates are equivalent to those obtained when
wi|j is available; the adjustment factors do not affect the rescaled version of the
student-level weights (Rabe-Hesketh and Skrondal, 2006).

4.2

Estimation Results

The estimates from the M-quantile random effect model with sampling weights
are presented in Table 2. The estimated coefficients are significantly different
from 0 at each quantile with a different intensity for almost all the covariates;
hence they indicate the usefulness of the estimation method used for the analysis.
Findings particularly show that gender has a negative effect on mathematics score at all levels with disadvantage increasing as the mathematics score
increases, hence also controlling for other individual and school-related characteristics males still outperform females in math, and the gender gap in math is
higher at the upper tail compared with the lower tail of the distribution (the
difference between male and female students at the first decile is 14.38, while at
the quantile 0.9 is 16.95). The peak of the gender gap is at q = 0.75, where boys
score about 18 points above girls. The estimated ‘compositional effect’ of gender
(i.e. how much the percentage of females in a school affects the mathematics
scores in that school), is significantly different from zero along the entire conditional distribution and in addition it increases as the score increases. Namely,
if a student moves to a school where the share of female is 10% higher, his/her
score decreases about 5 points for q = 0.1 and about 7 points for q = 0.9.
Looking at the other covariates further interesting considerations can be
made. As expected, native students perform better than migrants. The raw
penalty for immigrant students is lower at the bottom of the distribution than
elsewhere (the estimated coefficients range from −14.58 to −18.74), indicating
a more considerable negative effect for migrant students at the top end of the
distribution. It may be the result of the difficult for immigrants to achieve
very high scores. The share of immigrants in the school shows a negative ef-

14

fect. The effect considerably decreases along the conditional distribution of the
mathematics achievement highlighting that high-performing students are less
influenced by the presence of immigrants in the school. As we expected, math
scores of students that never failed school are higher than those who did and the
effect of such covariate decreases with increasing quantiles. By contrast, being
late for school regularly has a higher effect in the upper tail of distribution,
thus the gap is more elevated between best-performing students. If we look at
the compositional effects of these two last covariates a further consideration can
be made. The negative impact of the mean of students who repeat a grade is
evident across the entire distribution; moreover this effect increases, in absolute
term, as scores increase (the estimated effect is equal to −51.42 and −67.51 for
q = 0.1 and q = 0.9, respectively). Also the effect of the percentage of students
being late at school follows a similar pattern (−41.53 and −57.15 for q = 0.1
and q = 0.9, respectively). Insofar as these two variables are good proxies for
the share of truants, findings highlight how the school environment can affect
negatively the best performing students. Namely, the presence of negligent
students can lead teachers to adopt teaching methods targeted at these more
troublesome students to the detriment of good students being less stimulated
and consequentially likely affected negatively in their performance.
Furthermore, also the index of anxiety has a greater impact on the upper quantiles, confirming that high-achieving students may be more worried
than low-achievers about getting poor grades (see among others OECD, 2016c;
Foley et al., 2017).
Looking at the individual effect of ESCS index, we found that it is significantly associated with math scores only at q = 0.90 (β0.9 = 6.07), whereas its
school compositional effect (e.g. the mean of ESCS) is always significantly different from zero. Thus, the effect of peer “socio-economic status” tends to have
a stronger impact on student performance than the corresponding individual
socio-economic status.
As a consequence, being enrolled in an ‘elite school’, i.e. a school with many
students from wealthy backgrounds, improves the math scores. Moreover, the
effect of the socio-economic context of the school tends to decrease as quantiles
increase: low performing students therefore can enhance their performance in
‘elite schools’. This last consideration leads us to claim that it is necessary
to implement targeted policy measures in order to avoid the “ghettoisation”
of schools, that is schools with a high concentration of poor students, so that
students can benefit from the school-environment and improve their scores.
Finally, both macro-region and type of school have the expected negative sign
and their effects are significantly different from zero in the overall distribution.
Students living in Southern regions have low math scores than those living in
Northern regions, geographical disparities are less evident at the lower tail of
the score distribution and then increase even if they are almost stationary as
quantiles increase. Students attending technical, vocational and other schools
perform worse than those attending Lyceums confirming findings of other studies
(see among others Bratti et al., 2007; Matteucci and Mignani, 2014).

15

5

Simulation study

A small Monte Carlo simulation study was carried out to evaluate the performance of the Weighted-MQRE regression at three quantiles, q = 0.10, 0.25, 0.50.
The focus of this simulation study is to compare the estimation of the fixed
effects of the weighted version of the MQRE and the unweighted version. Moreover, we assess the approximations of the standard errors of the fixed effects.
For both aims, finite population values yij are generated under the two-level
random intercept model:
yij = 100 + 2xij + γj + εij ,

i = 1, ..., Nj ,

j = 1, ...M.

The number of level 2 units in the population is M = 170, each with the
same number of level 1 units, Nj = 50. The auxiliary variable is uniformly
distributed in [0, 20].
The level 1 and level 2 residuals are independently generated as follows: γ ∼
0.9 N (0, 1) + 0.1 N (9, 20) and ε ∼ 0.9 N (0, 3.3) + 0.1 N (10, 75). This represents
a situation under outlier contamination in both hierarchical levels.
Once the finite population values were obtained, we adopted the following
two stage sampling design. Level 2 units are divided into three strata according
to whether γ < −1, −1 ≤ γ ≤ 1 or γ > 1 and simple random samples of
size 0.15 m, 0.65 m and 0.20 m are selected from respective strata. The overall
number of sampled cluster is m = 100. From each level 2 units, level 1 units
are partitioned into two clusters according to whether ε > 0 or ε ≤ 0 and,
also in this case, simple random samples of size 0.75 nj and 0.25 nj are selected
form respective strata. The size nj was choice to be proportional to Nj , nj =
0.3 N j. This leads to a total sample size of n = 1500. This sampling scheme is
very similar to the method that was used in Rabe-Hesketh and Skrondal (2006)
except that we sampled both level units in a different way and in particular we
have sampled more cluster and units with extreme values. With this choice, we
assure that outliers in the population are also represented in the sample and
this allows us to evaluate the performance of the models in a situation where a
robust approach is required.
According to the previous studies, by making the sampling probabilities at
level 1 and 2 dependent on the corresponding residuals, we ensure that the
sampling scheme is informative at both levels.
We replicate this scenario R = 500 times and we compare the fixed effects
of the MQRE-2L (Tzavidis et al., 2016) and the Weighted-MQRE, with tuning
constant of the Huber influence function c = 1.345. At q = 0.5 we have estimated also the Linear Mixed Model (Eq. 2). Given the informativeness of
the sampling probabilities we expect that the Weighted-MQRE performs better
than the MQRE.
For each regression parameter, performance is evaluated using the so called
“Average Relative Bias” (ARB), defined as:
ARB(θ̂) = R−1

R
X
θ̂(r) − θ
r=1

θ

× 100,

where θ̂r is the estimated parameter at quantile q for the rth replication and θ
is the corresponding ‘true’ value of this parameter.

16

Table 3 reports the simulation results for estimators of the fixed effects for
q = 0.10, 0.25, 0.50. For the Weighted-MQRE we report results with “scaling
method 2”; estimates using unscaled weights and “scaling method 1” are not reported because they are almost identical. Scaling tends to affect the estimation
of the variance components, but evaluating the estimation of the variance parameters is not a focus concerning this work (for more detail about the effect of
scaling see simulation results in Pfeffermann et al., 1998). Taking a closer look
at ARB for the fixed effects, for the slope we observe that in all models there is
almost no bias. However the use of weights have an impact on the estimation of
the intercept. Indeed, the bias in the Weighted-MQRE is lower than that of the
unweighted models. In particular, the bias tends to decrease with the quantiles.
The ARB for intercept of the Weighted-MQRE is around 0.8% for q = 0.1 and
around 0.6% for q = 0.5, while for the MQRE it is around 1.4 for the first decile
and around 1.1% for q = 0.5. The LMM appears to be the worst model, the
bias at q = 0.5 is around 2%; given the presence of the outliers in both levels,
also the weighted version of the LMM shows a severe bias (the results are not
reported here).

17

Table 3: Values of bias ARB and the average of point estimates over simulations of the fixed effects using Weighted-MQRE, MQRE and
LMM at q = 0.10, 0.25, 0.50
Method

βˆ0 ,q = 0.10

βˆ1 ,q = 0.10

βˆ0 ,q = 0.25

βˆ1 ,q = 0.25

βˆ0 ,q = 0.50

βˆ1 ,q = 0.50

βˆ0

ARB

βˆ1

ARB

βˆ0

ARB

βˆ1

ARB

βˆ0

ARB

βˆ1

Weighted-MQRE

0.797

98.452

-0.025

1.9995

0.675

99.443

-0.013

1.9997

0.626

100.627

-0.011

1.9998

MQRE

1.406

99.048

-0.025

1.9995

1.258

100.020

-0.016

1.9997

1.118

101.121

-0.013

1.9997

-

-

-

-

-

-

-

-

2.022

102.026

-0.013

1.9997

18

ARB

LMM

Table 4: Empirical and estimated standard Errors of he fixed effects for q =
0.10, 0.25, 0.50 using Weighted-MQRE with tuning constant c = 1.345

q = 0.10
q = 0.25
q = 0.50

βˆ1

βˆ0

Values of q
Empirical
standard
error

Estimated
standard
error

Empirical
standard
error

Estimated
standard
error

0.198
0.169
0.180

0.211
0.190
0.215

0.016
0.013
0.012

0.015
0.012
0.012

Having evaluated the performance of the Weighted-MQRE regression, we
want to test the performance of the sandwich estimator for the variability of
the fixed effects. Therefore, we compare the empirical standard errors and
estimated standard errors. For each estimator θ̂, at q =0.10, 0.25, 0.50, Table 4 reports
averages over simulations of the Monte Carlo standard error
q

PR
PR
(r)
(r) − θ̄)2 , where θ̄ = R−1
−1
, and the estimated
R
S(θ̂) =
r=1 θ̂
r=1 (θ̂

standard errors of the fixed effects β̂q . It can be observed that the estimated
standard error of the estimators at all values of q offers a good approximation to
the empirical standard errors. Furthermore, as we expected, inserting weights
into the estimation procedure it produces a little overestimation of the standard
errors of the intercept.
In summary, our simulation results showed that accommodating sampling
weights in the estimation procedure of the MQRE is a good strategy for dealing
with an informative sampling design and outlier contamination in the data.

6

Conclusion

In this paper, we focus on the gender effect on math outcomes at upper secondary school in Italy. Previous research has suggested that the female penalty
may help explain the low share of women in STEM education, which is in
turn one of the reasons of the poor female presence in research and innovation decision-making positions. We improve existing studies in several ways
but, most importantly, we propose a new methodology that allows to apply Mquantile regression in a multilevel framework and in complex survey design. We
also propose a simulation study in order to evaluate the performance of the fixed
effects of the Weighted-MQRE showing that the weighted estimator can reduce
the bias induced by the complex sampling design when subjects and clusters
are selected with probabilities that depend on the model’s random terms. To
do this we use data drawn from 2015 PISA survey for Italy.
In summary, the existence of a gender-related issue in the Italian education system persists not only around the mean but in the overall distribution
and, in particular, the ability of girls in school is not able to mitigate the effect
of being female. Indeed, as we pointed out in the Section 1, high-performing
students tend to continue studying and, moreover, students who obtain high
grades in mathematics tend to choose science-related subjects in tertiary edu-

19

cation (OECD, 2012). As a consequence, we can argue that the larger gender
gap among the top performers is one of the main reasons of the low share of
girls in STEM. Reducing this gap is essential to combat gender-segregation in
the tertiary education and this might in turn reduce the gender gap in labour
market outcomes. The above findings can be linked with some interesting findings of the OECD (2015) report. Student attitudes and beliefs are correlated
with educational achievement, and in particular, the lack of self-confidence in
mathematics among girls affects their mathematics’ achievement. Indeed, selfconfidence gives students the possibility to go through a trial-and-error process
that is very important in developing mathematical skills; moreover, even the
high-achieving girls have low levels of confidence in their ability and this can
explain the gender gap in mathematics. Accordingly, encouraging girls to change
their attitudes could improve gender equality in the math scores, and parents
and teachers represent the best mediators in this process, because they directly
affect their education and attitudes.
It is beyond doubt that psychological factors frame much of our tools for
reducing these gaps, but the macro-societal characteristics (i.e. societal inequality and the welfare orientation) of Italy may play an important role as well
(Marks, 2008; Cipollone et al., 2014). In this view, the welfare state of Italy,
combining elements of the continental-corporatist model and the familist model
model, is characterized by relatively low female labour market participation
anyhow (Ferrera, 2010; Katrougalos and Lazaridis, 2008; Naldini and Jurado,
2013; Lynch, 2014). Therefore, policy innovations (family-oriented policies) are
likely another important way of reducing gender gaps.
Finally, further work using other OECD countries, and comparing results
across countries would be the natural extension for applying the new methodology and analysing cross-country level effects that can add further information
for improving gender equality in Italy.

Appendix A
The steps of the estimation algorithm are as follows:
1. Let σγ2q and σε2q be known.
2. Given the variance parameter and the covariance matrix Vq , the iterative
equation for the estimates of β ψq is:
t
T
−1/2
−1
−1
βt+1
Dq (βtψq )U1/2
[XT ΓVq−1 U1/2
q Vq X]
q ψq (rq )],
ψq = β ψq + [X ΓUq

where Γ is a diagonal matrix with jth diagonal element equals to the
second level weight wj ; Dq (β tψq ) is a diagonal matrix with jth diagonal
′
′
element Dijq = ψq (rijq ) = ∂ψq (rijq )/∂rijq .
3. The estimates of β ψq are used to obtain the estimates of the variance parameters with a fixed point iterative method. As pointed out in Tzavidis et al.
(2016) this requires the change of the estimating equation in Eq. 13 to:

20

−1
T −1 1/2
Γψq (rq )T U1/2
q Vq ZZ Vq Uq ψq (rq )

 2 
σ
=0
− K2q tr ΓVq−1 ZZT Vq−1 (ZZT W−1 ) γ2q
σεq

−1
−1 −1 1/2
Vq Uq ψq (rq )
Γψq (rq )T U1/2
q Vq W

 2 
σ
−1
−1
−1
T
−1
− K2q tr ΓVq W Vq (ZZ W ) γ2q
=0
σεq

The fixed point algorithm of the estimating equation for the tth iteration
is the following:
!)−1
!
 2  (
2(t)
2(t)
σγq
σγq
σγq
= A
a
2(t)
2(t) ,
σε2q
σεq
σεq
where


 2 
K2q tr(ΓVq−1 ZZT Vq−1 ZZt ) K2q tr(ΓVq−1 ZZT Vq−1 W−1 )
σγq

A
=
σε2q
−1
−1 −1
T
−1
−1 −1
−1
K2q tr(ΓVq W Vq ZZ ) K2q tr(ΓVq W Vq W )
and


1
T 1/2 −1
T
−1 1/2
 2 
ψ
(r
)
U
V
ZZ
ΓV
U
ψ
(r
)
q
q
q
q
q
q
q
q
σ
2

a γ2q = 

σεq
1/2
1/2
1
T
−1
−1
−1
ΓVq Uq ψq (rq )
2 ψq (rq ) Uq Vq W

4. Iterate steps 2 and 3 until convergence.

Appendix B
The asymptotic covariance matrix of the estimators can be written as G−1 F(G−1 )T .
The weighted estimating equations in the Eq. 12 and Eq. 13 can be written
as:
m
X

(16)

wj Φqj (θ q ) = 0,

j=1

T

where θq = (β̂ q , σ̂ε2q , σ̂γ2q ) and Φqj (θ q ) = ΦTqjβ q , Φqjσγ2q , Φqjσε2q

T

.

Under a general response distribution D, the estimator θ̂q satisfying equations (16) is estimating a root θ q of:
m
X

ED [wj Φqj (θ q )] = 0.

j=1

21

T
where Φqj (θ q ) = ΦTqjβ q , Φqjσγ2q , Φqjσε2q
for particular choice of Φqj (θq ).
Under a general response distribution of D (i.e. not necessarily Gaussian), a
robust estimator is estimating a root θq of the expectation with respect to D of
the previous estimating equations.
Provided that
"
#
m
X
∂(wj Φqj (θ q ))
−1
ED
n
→ G,
∂θq
j=1

is the diagonal block information matrix of dimension (p + 2) × (p + 2) and
it is positive definite, and
"
#
m
X
(wj Φqj (θ q ))T
−1
ED
→ F,
n
wj Φqj (θ q )
j=1

is the the matrix of the variance of the normalised score function, a Taylor series
approximation which which holds uniformly in a neighbourhood of θq is:
θ̂q ≈ θq + G−1 n−1

d
X

wj Φqj (θq ) + op (n−1/2 ),

j=1

n → ∞.

Following Richardson and Welsh (1995) the covariance matrix can be consistently estimated by Ĝ−1 F̂(Ĝ−1 )T where the matrices Ĝ and F̂ are evaluated
at θˆq .

Appendix C
Following Fabrizi et al. (2014) we rewrite the assumption 1, 2, 8 and 9 of
Wang and Opsomer (2011).
Assumption A1. The expected number of cluster m∗ = Ed (m|y1 , y2 , ..., yM ) =
O(M δ ), where the vectors yi , i = 1, 2, ...3, M denote all the variables of interest
in each group of the M th population. The subscript d denotes the expectations
with respect to the randomization distribution induced by the sampling design.
Assumption A2. KL ≤ M πj /m∗ ≤ KU for all j, with KL and KU positive
constants.
Assumption A3. For any vector z with finite 2 + λ moments with
arbitrarPM
ily small λ > 0, we assume Vd (z̄π |y2 , ..., yM ) ≤ c1 m∗−1 (M − 1)−1 j=1 (zj −
P
z̄M )(zj − z̄M )T , for some constant c1 and z̄π = M −1 m
j=1 wj zj .
Assumption A4. The population parameter Bψq lies in a closed interval
ΘBψq on Rp .
Assumption A5. The population estimating function SM (·) and the function ψ(·) satisfy:
- the function ψ(·) is bounded;
- the population estimating function SM (βψq ) converges to S(β ψq ) uniformly on ΘBψq as M → ∞ and the equation S(β ψq ) = 0 has a unique
root in the interior of ΘBψq ;
- the limiting function S(β ψq ) is strictly increasing and absolutely continu′
ous with the finite first derivative in ΘBψq , and the derivative S (β ψq ) is
bounded away from 0 for β ψq in ΘBψq ;
22

- the population quantities:
sup
βψq ∈Θβ ψq

M α |SM (Bψq +M −α βψq )−SM (Bψq )−S(Bψq +M −α β ψq )+S(Bψq )| → 0

and
sup
β ψq ∈Θβ ψq

M −1

M
X
j=1

|[ψq (Vj−1 (yj − xTj βψq − M −α xTj β ψq ))−
−1/2

ψq (Vj−1 (yj − xTj Bψq ))]wj Xj Vj

| = Op (M −α )

where Θβψq is a large enough compact set in Rp and α ∈ (1/4, 1/2).

References
Agasisti, T., F. Ieva, and A. Paganoni (2014). Heterogeneity, school-effects
and achievement gaps across Italian regions: further evidence from statistical
modeling. MOX–Report 07/2014, MOX, Dipartimento di Matematica “F.
Brioschi”.
Alfò, M., N. Salvati, and M. G. Ranalli (2017). Finite mixtures of quantile and
M-quantile regression models. Statistics and Computing 27 (2), 547–570.
AlmaDiploma (2017). XI Indagine Esiti a distanza dei Diplomati. Rapporto
2017.
Ammermueller, A. (2007). Poor background or low returns? Why immigrant
students in Germany perform so poorly in the Programme for International
Student Assessment. Education Economics 15 (2), 215–230.
Asparouhov, T. (2006). General multi-level modeling with sampling weights.
Communications in Statistics – Theory and Methods 35 (3), 439–460.
Azzolini, D., P. Schnell, and J. R. Palmer (2012). Educational achievement gaps
between immigrant and native students in two “new” immigration countries:
Italy and Spain in comparison. The Annals of the American Academy of
Political and Social Science 643 (1), 46–77.
Binder, D. A. (1983). On the variances of asymptotically normal estimators
from complex surveys. International Statistical Review/Revue Internationale
de Statistique 51 (3), 279–292.
Borgoni, R., P. Del Bianco, N. Salvati, T. Schmid, and N. Tzavidis (2016). Modelling the distribution of health-related quality of life of advanced melanoma
patients in a longitudinal multi-centre clinical trial using M-quantile random
effects regression. Statistical methods in medical research, 1–17.
Bradley, K. (2000). The incorporation of women into higher education: Paradoxical outcomes? Sociology of education, 1–18.
Bratti, M., D. Checchi, and A. Filippin (2007). Territorial Differences in Italian
Students’ Mathematical Competencies: Evidence from PISA 2003. Discussion
paper, IZA.
23

Breckling, J. and R. Chambers (1988). M-quantiles. Biometrika 75 (4), 761–771.
Brown, C. and M. Corcoran (1997). Sex-based differences in school content and
the male-female wage gap. Journal of Labor Economics 15 (3), 431–465.
Chambers, R. and N. Tzavidis (2006). M-quantile models for small area estimation. Biometrika 93 (2), 255–268.
Chambers, R. L. (2003). Introduction to part A. In R. L. Chambers and C. J.
Skinner (Eds.), Analysis of survey data, pp. 13–27. Chichester: John Wiley
& Sons.
Cipollone, A., E. Patacchini, and G. Vallanti (2014). Female labour market
participation in Europe: novel evidence on trends and shaping factors. IZA
Journal of European Labor Studies 3 (18), 1–40.
Cochran, W. G. (1977). Sampling Techniques (3rd ed.). New York: John Wiley
& Sons.
Costanzo, A. (2015). The effect of M@tabel on italian students’ performances: A
quantile regression approach. Procedia - Social and Behavioral Sciences 197,
236 – 244.
Duquet, N., I. Glorieux, I. Laurijssen, and Y. Van Dorsselaer (2010). An unfinished job? The effect of subject choice and family formation processes on
labour market outcomes of young men and women. Journal of Education and
Work 23 (4), 319–338.
European Commission (2012). Structural change in research institutions: Enhancing excellence, gender equality and efficiency in research and innovation.
Luxemburg: Publications Office of the European Union.
European Commission (2016). She Figures 2015. Luxembourg: Publications
Office of the European Union.
Fabrizi, E., N. Salvati, M. Pratesi, and N. Tzavidis (2014). Outlier robust
model-assisted small area estimation. Biometrical Journal 56 (1), 157–175.
Faria, S. and M. C. Portela (2016). Student performance in mathematics using PISA-2009 data for Portugal. Working Paper 01/2016, Catolica Porto
Business School.
Ferrera, M. (2010). The south european countries. In F. Castel, S. Leibfried,
J. Lewis, and C. Pierson (Eds.), The Oxford Handbook of the Welfare State,
pp. 616–629. Oxford: Oxford University Press.
Foley, A. E., J. B. Herts, F. Borgonovi, S. Guerriero, S. C. Levine, and S. L.
Beilock (2017). The math anxiety-performance link: A global phenomenon.
Current Directions in Psychological Science 26 (1), 52–58.
Geraci, M. and M. Bottai (2007). Quantile regression for longitudinal data using
the asymmetric Laplace distribution. Biostatistics 8 (1), 140–154.
Geraci, M. and M. Bottai (2014). Linear quantile mixed models. Statistics and
computing 24 (3), 461–479.
24

Goldstein, H. (1986). Multilevel mixed linear model analysis using iterative
generalized least squares. Biometrika 73 (1), 43–56.
Goldstein, H. (2011). Multilevel statistical models (4th ed.). John Wiley & Sons.
Grilli, L., F. Pennoni, C. Rampichini, I. Romeo, et al. (2016). Exploiting timss
and pirls combined data: multivariate multilevel modelling of student achievement. The Annals of Applied Statistics 10 (4), 2405–2426.
Grilli, L. and M. Pratesi (2004). Weighted estimation in multilevel ordinal
and binary models in the presence of informative sampling designs. Survey
methodology 30 (1), 93–104.
Hajisoteriou, C. and A. Panayiotis (2016). The Globalisation of Intercultural
Education. Macmillan UK: Palgrave.
Huggins, R. (1993). On the robust analysis of variance components models for
pedigree data. Australian & New Zealand Journal of Statistics 35 (1), 43–57.
Ikeda, M. and E. García (2014). Grade repetition: A comparative study of
academic and non-academic consequences. OECD Journal: Economic Studies 2013/1, 269–315.
Invalsi (2013). OCSE PISA 2012. Rapporto Nazionale. A cura di Invalsi. Invalsi.
Invalsi (2016). Indagine OCSE PISA 2015. I risultati degli studenti italiani in
Science, Matematica e Lettura. Invalsi.
ISTAT (2017). Bilancio demografico nazionale. Statistiche report 13 giugno
2017, Istat.
Jones, M. C. (1994). Expectiles and M-quantiles are quantiles. Statistics &
Probability Letters 20 (2), 149–153.
Katrougalos, G. and G. Lazaridis (2008). The south european welfare states at
the dawn of the new millennium: identity and problems. Social Cohesion and
Development 3 (1), 5–25.
Koenker, R. (2004). Quantile regression for longitudinal data. Journal of Multivariate Analysis 91 (1), 74–89.
Koenker, R. (2005). Quantile Regression. Econometric Society Monographs.
Cambridge University Press.
Koenker, R. and G. Bassett (1978). Regression quantiles. Econometrica 46 (1),
33–50.
Kokic, P., R. Chambers, J. Breckling, and S. Beare (1997). A measure of
production performance. Journal of Business & Economic Statistics 15 (4),
445–451.
Lynch, J. (2014). The Italian welfare state after the financial crisis. Journal of
Modern Italian Studies 19 (4), 380–388.
Marino, M. F. and A. Farcomeni (2015). Linear quantile regression models for
longitudinal experiments: an overview. Metron 73 (2), 229–247.
25

Marks, G. N. (2008). Accounting for the gender gaps in student performance
in reading and mathematics: evidence from 31 countries. Oxford Review of
Education 34 (1), 89–109.
Masci, C., F. Ieva, T. Agasisti, and A. M. Paganoni (2016). Bivariate multilevel
models for the analysis of mathematics and reading pupils’ achievements.
Journal of Applied Statistics, 1–22.
Matteucci, M. and S. Mignani (2014). Exploring regional differences in the
reading competencies of Italian students. Evaluation review 38 (3), 251–290.
Naldini, M. and T. Jurado (2013). Family and welfare state reorientation in
Spain and inertia in Italy from a European perspective. Population Review 52 (1), 43–61.
Newey, W. K. and J. L. Powell (1987). Asymmetric least squares estimation
and testing. Econometrica, 819–847.
OECD (2003). Education at a Glance 2003: OECD Indicators. Paris: OECD
Publishing.
OECD (2009). PISA Data Analysis Manual: SPSS, Second Edition. Paris:
OECD Publishing.
OECD (2012). Closing the Gender Gap. Act Now. Paris: OECD Publishing.
OECD (2013). PISA 2012 Results: What Makes Schools Successful (Volume
IV). Paris: OECD Publishing.
OECD (2015). The ABC of Gender Equality in Education: Aptitude, Behaviour,
Confidence. Paris: OECD Publishing.
OECD (2016a). PISA 2015 Results (Volume I): Excellence and Equity in Education. Paris: OECD Publishing.
OECD (2016b). PISA 2015 Results (Volume II): Policies and practices for
succesfsul schools. Paris: OECD Publishing.
OECD (2016c). PISA 2015 Results (Volume III): Students’ Well-Being, PISA.
Paris: OECD Publishing.
OECD (2017, February). PISA 2015 Technical Report (draft format). Available
online at: http://www.oecd.org/pisa/data/2015-technical-report/.
Pfeffermann, D. (1993). The role of sampling weights when modeling survey data. International Statistical Review/Revue Internationale de Statistique 61 (2), 317–337.
Pfeffermann, D., C. J. Skinner, D. J. Holmes, H. Goldstein, and J. Rasbash
(1998). Weighting for unequal selection probabilities in multilevel models. Journal of the Royal Statistical Society, Series B (Statistical Methodology) 60 (1), 23–40.
R Core Team (2015). R: A Language and Environment for Statistical
Computing. Vienna, Austria: R Foundation for Statistical Computing.
https://www.R-project.org.
26

Rabe-Hesketh, S. and A. Skrondal (2006). Multilevel modelling of complex
survey data. Journal of the Royal Statistical Society, Series A (Statistics in
Society) 169 (4), 805–827.
Rangvid, B. S. (2007a). School composition effects in Denmark: quantile regression evidence from PISA 2000. Empirical Economics 33 (2), 359–388.
Rangvid, B. S. (2007b). Sources of immigrants’ underachievement: results from
PISA-Copenhagen. Education Economics 15 (3), 293–326.
Richardson, A. M. and A. H. Welsh (1995). Robust restricted maximum likelihood in mixed linear models. Biometrics 51 (4), 1429–1439.
Sani, C. and L. Grilli (2011). Differential variability of test scores among schools:
A multilevel analysis of the fifth-grade INVALSI test using heteroscedastic
random effects. Journal of applied quantitative methods 6 (4), 88–99.
Sarndal, C. E., B. Swensson, and B. Wretman (1992). Model Assisted Survey
Sampling. New York: Springer-Verlag.
Schneeweis, N. and R. Winter-Ebmer (2007). Peer effects in Austrian schools.
Empirical Economics 32 (2), 387–409.
Sinha, S. K. and J. Rao (2009). Robust small area estimation. Canadian Journal
of Statistics 37 (3), 381–399.
Tzavidis, N., N. Salvati, T. Schmid, E. Flouri, and E. Midouhas (2016). Longitudinal analysis of the strengths and difficulties questionnaire scores of the
millennium cohort study children in england using M-quantile random-effects
regression. Journal of the Royal Statistical Society, Series A (Statistics in
Society) 179 (2), 427–452.
van Langen, A., R. Bosker, and H. Dekkers (2006). Exploring cross-national
differences in gender gaps in education. Educational Research and Evaluation 12 (02), 155–177.
Wang, J. C. and J. D. Opsomer (2011). On asymptotic normality and variance
estimation for nondifferentiable survey estimators. Biometrika 98 (1), 91–106.
World Bank (2011). World development report 2012: gender equality and development. Washington DC: World Bank Publications.

27

