Statistical Inference for Covariate-Adjusted and Interpretable Generalized Factor Model with Application to Testing Fairness
Jing Ouyang Chengyu Cui Kean Ming Tan Gongjun Xu




Statistical Inference for Covariate-Adjusted and Interpretable Generalized Factor Model with Application to Testing Fairness
    Jing Ouyang, Chengyu Cui, Kean Ming Tan, and Gongjun Xu 
 

Department of Statistics, University of Michigan
    
============================================================================================================================


    In the era of data explosion, statisticians have been developing interpretable and computationally efficient statistical methods to measure latent factors (e.g., skills, abilities, and personalities) using large-scale assessment data. In addition to understanding the latent information, the covariate effect on responses controlling for latent factors is also of great scientific interest and has wide applications, such as evaluating the fairness of educational testing, where the covariate effect reflects whether a test question is biased toward certain individual characteristics (e.g., gender and race) taking into account their latent abilities. However, the large sample size, substantial covariate dimension, and great test length pose challenges to developing efficient methods and drawing valid inferences. Moreover, to accommodate the commonly encountered discrete types of responses, nonlinear latent factor models are often assumed, bringing further complexity to the problem. To address these challenges, we consider a covariate-adjusted generalized factor model and develop novel and interpretable conditions to address the identifiability issue. Based on the identifiability conditions, we propose a joint maximum likelihood estimation method and establish estimation consistency and asymptotic normality results for the covariate effects under a practical yet challenging asymptotic regime. Furthermore, we derive estimation and inference results for latent factors and the factor loadings.  We illustrate the finite sample performance of the proposed method through extensive numerical studies and an application to an educational assessment dataset obtained from the Programme for International Student Assessment (PISA). 

   

  Keywords: Generalized factor model; Covariate adjustment; Large-scale testing fairness

    






§ INTRODUCTION



Latent factors, often referred to as hidden factors, play an increasingly important role in modern statistics to analyze large-scale complex measurement data and find wide-ranging applications across various scientific fields, including educational assessments <cit.>, macroeconomics forecasting <cit.>, and biomedical diagnosis <cit.>.
For instance, in educational testing and social sciences, latent factors are used to model unobservable traits of respondents, such as skills, personality, and attitudes <cit.>; in biology and genomics, latent factors are used to capture underlying genetic factors, gene expression patterns, or hidden biological mechanisms <cit.>.  
To uncover the latent factors and analyze large-scale complex data, various latent factor models

have been developed and extensively investigated in the existing literature <cit.>. 


In addition to measuring the latent factors, the observed covariates and the covariate effects conditional on the latent factors hold significant scientific interpretations in many applications <cit.>. 
One important application is testing fairness, which receives increasing attention in the fields of education, psychology, and social sciences <cit.>. In educational assessments, testing fairness, or measurement invariance, implies that groups from diverse backgrounds have the same probability of endorsing the test items, controlling for individual proficiency levels <cit.>. 
Testing fairness is not only of scientific interest to psychometricians and statisticians but also attracts widespread public awareness <cit.>.  
In the era of rapid technological advancements, international and large-scale educational assessments are becoming increasingly prevalent. 
One example is the Programme for International Student Assessment (PISA), which is a large-scale international assessment with substantial sample size and test length <cit.>. PISA assesses the knowledge and skills of 15-year-old students in mathematics, reading, and science domains <cit.>. In PISA 2018, over 600,000 students from 37 OECD[OECD: Organisation for Economic Co-operation and Development] countries and 42 partner countries/economies participated in the test <cit.>. 
 To assess fairness of the test designs in such large-scale assessments, it is important to develop modern and computationally efficient methodologies for interpreting the effects of observed covariates (e.g., gender and race) on the item responses, controlling for the latent factors.


However, the discrete nature of the item responses, the increasing sample size, and the large amount of test items in modern educational assessments pose great challenges for the estimation and inference for the covariate effects as well as for the latent factors.
For instance, in educational and psychological measurements, such a testing fairness issue (measurement invariance) is typically assessed by differential item functioning (DIF) analysis of item response data that aims to detect the DIF items, where a DIF item has a response distribution that depends on not only the measured latent factors but also respondents’ covariates (such as group membership). Despite many statistical methods that have been developed for DIF analysis, existing methods often require domain knowledge to pre-specify DIF-free items, namely anchor items, which may be misspecified and lead to biased estimation and inference results <cit.>. To address this limitation, researchers developed item purification methods to iteratively select anchor items through stepwise selection models <cit.>. More recently, tree-based methods <cit.>, regularized estimation methods <cit.>, item pair functioning methods <cit.>, and many other non-anchor-based methods have been proposed. However, these non-anchor-based methods do not provide valid statistical inference guarantees for testing the covariate effects. It remains an open problem to perform statistical inference on the covariate effects and the latent factors in educational assessments.

To address this open problem, we study the statistical estimation and inference for a general family of covariate-adjusted nonlinear factor models, which include the popular factor models for binary, count, continuous, and mixed-type data that commonly occur in educational assessments.
The nonlinear model setting poses great challenges for estimation and statistical inference.
Despite recent progress in the factor analysis literature, most existing studies focus on estimation and inference under linear factor models <cit.>
and covariate-adjusted linear factor models <cit.>. 
The techniques employed in linear factor model settings are not applicable here due to the nonlinearity inherent in the general models under consideration.


Recently, several researchers have also investigated the parameter estimation and inference for generalized linear factor models <cit.>. 
However, they either focus only on the overall consistency properties of the estimation or do not incorporate covariates into the models. 

In a concurrent work, motivated by applications in single-cell omics, <cit.> considered a generalized linear factor model with covariates and studied its inference theory, where the latent factors are used as surrogate variables to control for unmeasured confounding. However, they imposed relatively stringent assumptions on the sparsity of covariate effects and the dimension of covariates, and their theoretical results also rely on data-splitting. Moreover, <cit.> focused only on statistical inference on the covariate effects, while that on factors and loadings was unexplored, which is often of great interest in educational assessments. Establishing inference results for covariate effects and latent factors simultaneously under nonlinear models remains an open and challenging problem,
due to the identifiability issue from the incorporation of covariates and the nonlinearity issue in the considered general models.



To overcome these issues,  we develop a novel framework for performing statistical inference on all model parameters and latent factors under a general family of covariate-adjusted generalized factor models. 


Specifically, we propose a set of interpretable and practical identifiability conditions for identifying the model parameters, and further incorporate these conditions into the development of a computationally efficient likelihood-based estimation method. 
Under these identifiability conditions, we develop new techniques to address the aforementioned theoretical challenges and obtain estimation consistency and asymptotic normality for covariate effects under a practical yet challenging asymptotic regime. 
Furthermore, building upon these results,
we establish estimation consistency and provide valid inference results for factor loadings and latent factors that are often of scientific interest, advancing our theoretical understanding of nonlinear latent factor models. 

 

The rest of the paper is organized as follows. In Section <ref>, we introduce the model setup of the covariate-adjusted generalized factor model. Section <ref> discusses the associated identifiability issues and further presents the proposed identifiability conditions and estimation method. Section <ref> establishes the theoretical properties for not only the covariate effects but also the latent factors and factor loadings. In Section <ref>, we perform extensive numerical studies to illustrate the performance of the proposed estimation method and the validity of the theoretical results. In Section <ref>, we analyze an educational testing dataset from Programme for International Student Assessment (PISA) and identify test items that may lead to potential bias among different test-takers.
We conclude with providing some potential future directions in Section <ref>.
  



Notation: For any integer N, let [N] = {1, …, N}. 
For any set , let # be its cardinality.
For any vector = (r_1, …, r_l)^, let _0  = # ({j: r_j ≠ 0 }),  _∞= max_j = 1, …, l |r_j|, and _q = (∑_j=1^l |r_j|^q)^1/q for q ≥ 1.
 We define 1_x^(y) to be the y-dimensional vector with x-th entry to be 1 and all other entries to be 0.
For any symmetric matrix , let λ_min() and λ_max() be the smallest and largest eigenvalues of . 
For any matrix = (a_ij)_n× l, let _∞,1 = max_j=1,…, l ∑_i = 1^n |a_ij | be the maximum absolute column sum, _1,∞ = max_i=1,…, n∑_j=1^l |a_ij| be the maximum of the absolute row sum, _max = max_i,j |a_ij| be the maximum of the absolute matrix entry,
_F = (∑_i=1^n ∑_j=1^l |a_ij|^2)^1/2 be the Frobenius norm of , and =√(λ_max(^)) be the spectral norm of .  Let ·_φ_1 be sub-exponential norm. 
 Define the notation _v = vec() ∈^nl to indicate the vectorized matrix ∈^n× l. Finally, we denote ⊗ as the Kronecker product.





§ MODEL SETUP




Consider n independent subjects with q measured responses and p_* observed covariates.  

For the ith subject, let _i ∈^q be a q-dimensional vector of responses corresponding to q measurement items and   _i^c ∈^p_* be a p_*-dimensional vector of observed covariates. Moreover, let _i be a K-dimensional vector of latent factors representing the unobservable traits such as skills and personalities, where we assume K is specified as in many educational assessments.  

 
 We assume that the q-dimensional responses _i are conditionally independent, given _i^c and _i.  
 Specifically, we model the jth response for the ith subject, Y_ij, by the following  conditional distribution:

    Y_ij∼ p_ij( y | w_ij),     where  w_ij = β_j0 + _j^⊺_i  + _jc^⊺_i^c.

 Here β_j0∈ is the intercept parameter, _jc =(β_j1,…,β_jp_*)^⊺∈^p_* are the coefficient parameters for the observed covariates, and _j = (γ_j1,…,γ_jK)^⊺∈^K are the factor loadings. 
For better presentation, we write _j= (β_j0,_jc^)^ as an assembled vector of intercept and coefficients and define _i = (1, (_i^c)^)^ with dimension p = p_* + 1, which gives 
    w_ij =  _j^⊺_i  + _j^⊺_i.

 Given w_ij, the function p_ij is some specified probability density (mass) function. 
 Here, we consider a general and flexible modeling framework by allowing different types of p_ij functions to model diverse response data in wide-ranging applications, such as binary item response data in educational and psychological assessments <cit.> and mixed types of data in educational and macroeconomic applications <cit.>; see also Remark <ref>.
 A schematic diagram of the proposed model setup is presented in Figure <ref>.







Our proposed covariate-adjusted generalized factor model in (<ref>) is motivated by applications in testing fairness.
 In the context of educational assessment, the subject's responses to questions are dependent on latent factors _i such as students' abilities and skills, and are potentially affected by observed covariates _i^c such as age, gender, and race, among others <cit.>. 
  The intercept β_j0 
  
  is often interpreted as the difficulty level of item j and referred to as the difficulty parameter in psychometrics <cit.>.
 The capability of item j to further differentiate individuals based on their latent abilities is captured by _j = (γ_j1,…,γ_jK)^⊺, which are also referred to as discrimination parameters <cit.>. 
 The effects of observed covariates _i^c on subject's response to the jth question Y_ij, conditioned on latent abilities _i, are captured by _jc =(β_j1,…,β_jp_*)^⊺, which are referred to as DIF effects in psychometrics <cit.>. This setting gives rise to the fairness problem of validating whether the response probabilities to the measurements differ across different genders, races, or countries of origin while holding their abilities and skills at the same level. 



Given the observed data from n independent subjects, we are interested in studying the relationships between _i and _i^c after adjusting for the latent factors _i in (<ref>). 
Specifically, our goal is to test the statistical hypothesis H_0: β_js = 0 versus H_a: β_js≠ 0 for s ∈ [p_*], where β_js is 
the regression coefficient for the sth covariate and the jth response, after adjusting for the latent factor _i. 
In many applications, the latent factors and factor loadings also carry important scientific interpretations such as students' abilities and test items' characteristics. This motivates us to perform statistical inference on the parameters β_j0, _j, and _i as well.  




The proposed model setup (<ref>) is general and flexible as various functions p_ij's could be used to model diverse types of response data in wide-ranging applications.
For instance, in educational assessments, logistic factor model <cit.> with 

    p_ij (y | w_ij) = exp(w_ijy)/{1+exp(w_ij)},    y ∈{0,1}

and probit factor model <cit.> with 
    p_ij  (y | w_ij) = {Φ (w_ij)}^y {1-Φ (w_ij)}^1-y,   y ∈{0,1}
 
where Φ(·) is the cumulative density function of standard normal distribution,
are widely used to model the binary responses, indicating correct or incorrect answers to the test items.

 Such types of models are often referred to as item response theory models <cit.>.
In economics and finances, linear factor models with
p_ij (y| w_ij) ∝exp{ - (y - w_ij)^2/(2σ^2)}, where y ∈ and σ^2 is the variance parameter,
are commonly used to model continuous responses, such as GDP, interest rate, and consumer index <cit.>. Moreover, 
depending on the the observed responses, different types of function
p_ij's can be used to model the response from each item 
j ∈ [q]. Therefore, mixed types of data, which are common in educational measurements <cit.> and macroeconomic applications <cit.>, can also be analyzed by our proposed model.




In addition to testing fairness, the considered model finds wide-ranging applications in the real world. For instance, in genomics, the gene expression status may depend on unmeasured confounders or latent biological factors  and also be associated with the variables of interest including medical treatment, disease status, and gender <cit.>.
The covariate-adjusted general factor model helps to investigate the effects of the variables of interest on gene expressions, controlling for the latent factors <cit.>.  This setting is also applicable to other scenarios, such as brain imaging, where the activity of a brain region may depend on measurable spatial distance from neighboring regions and latent structures due to unmodeled factors <cit.>.





To analyze large-scale measurement data, we aim to develop a computationally efficient estimation method and to provide inference theory for quantifying uncertainty in the estimation.

Motivated by recent work in high-dimensional factor analysis, we treat the latent factors as fixed parameters and apply a joint maximum likelihood method for estimation <cit.>. 


Specifically, we let the collection of the item responses from n independent subjects be = (_1, …, _n)_n× q^⊺ and the design matrix of observed covariates to be = (_1, …, _n)_n × p^⊺. For model parameters, the discrimination parameters for all q items are denoted as 
= (_1, …, _q)_q× K^⊺, while the intercepts and the covariate effects for all q items are denoted as = (_1, …, _q)_q × p^⊺. The latent factors from all n subjects are = (_1, …, _n)_n × K^⊺.
Then, the joint log-likelihood function can be written as follows:

    L( |, , , ) = 1/nq∑_i=1^n ∑_j=1^q l_ij(β_j0 + _j^⊺_i  + _jc^⊺_i^c),

where the function l_ij (w_ij) = log p_ij (Y_ij | w_ij) is the individual log-likelihood function with w_ij = β_j0 + _j^⊺_i  + _jc^⊺_i^c. 
We aim to obtain (, , ) from maximizing the joint likelihood function L( |, , , ).

While the estimators can be computed efficiently by maximizing the joint likelihood function through an alternating maximization algorithm <cit.>,  challenges emerge for performing statistical inference on the model parameters. 

    
  * One challenge concerns the model identifiability.
Without additional constraints, the covariate effects are not identifiable due to the incorporation of covariates and their potential dependence on latent factors. The latent factors and factor loadings encounter similar identifiability issues as in traditional factor analysis <cit.>. Ensuring that the model is statistically identifiable is the fundamental prerequisite for achieving model reliability and making valid inferences <cit.>.

  * Another challenge arises from the nonlinearity of our proposed model. 

In the existing literature, most studies focus on the statistical inference for our proposed setting in the context of linear models 

<cit.>.
On the other hand, settings with general log-likelihood function l_ij (w_ij), including covariate-adjusted logistic and probit factor models, are less investigated. 
Common techniques for linear models 

are not applicable to the considered general nonlinear model setting.









Motivated by these challenges, we propose interpretable and practical identifiability conditions in Section <ref>. We then incorporate these conditions into the joint-likelihood-based estimation method in Section <ref>. Furthermore, we introduce a novel inference framework for  performing statistical inference on _j, _j, and _i in Section <ref>.





§ METHOD

	 
 
 

 §.§ Model Identifiability

 
Identifiability issues commonly occur in latent variable models <cit.>.
The proposed model in (<ref>) has two major identifiability issues. 
The first issue is that the proposed model remains unchanged after certain linear transformations of both  and , causing the covariate effects together with the intercepts, represented by , and the latent factors, denoted by , to be unidentifiable.   
 The second issue is that the model is invariant after an invertible transformation of both  and  as in the linear factor models <cit.>, causing the latent factors  and factor loadings  to be undetermined.  
 
 
Specifically, under the model setup in (<ref>), we define the joint probability distribution of responses to be P( |, , , ) = ∏_i=1^n ∏_j=1^q p_ij (Y_ij | w_ij).
The model parameters are identifiable if and only if for any response , there does not exist (, , ) ≠ (, , ) such that P( |, , , ) = P(|, , , ). 
The first issue concerning the identifiability of  and  is that for any (, , ) and any transformation matrix , there exist =, =   + ^, and =  - such that P( |, , , ) = P(|, , , ).
 This identifiability issue leads to the indeterminacy of the covariate effects and latent factors. 
The second issue is related to the identifiability of  and .

For any (, , ) and any invertible matrix , there exist = (^)^-1, =, and = such that P( |, , , ) = P(|, , , ). This causes the latent factors and factor loadings to be unidentifiable.



 

    Intuitively, the unidentifiable  =  - can be interpreted to include both direct and indirect effects of  on response . We take the intercept and covariate effect on the first item (_1) as an example and illustrate it in Figure <ref>. One part of _1 is the direct effect from  onto  (see the orange line in the left panel), whereas another part of _1 may be explained through the latent factors , as the latent factors are unobserved and there are potential correlations between latent factors and observed covariates. The latter part of _1 can be considered as the indirect effect (see the blue line in the right panel).




 
 The first identifiability issue is a new challenge introduced by the covariate adjustment in the model, whereas the second issue is common in traditional factor models <cit.>. Considering the two issues together, for any (, , ), , and , there exist transformations =  (^)^-1, =  ( + ^), and = - such that P (|, , , ) = P (|, , , ). In the rest of this subsection, we propose identifiability conditions to address these issues.

For notation convenience, throughout the rest of the paper, we define  ^* = (^*, ^*, ^*) as the true parameters.




  
Identifiability Conditions As described earlier, the correlation between the design matrix of covariates  and the latent factors ^* results in the identifiability issue of ^*.
In the psychometrics literature, the intercept β_j0^* is commonly referred to as the difficulty parameter, while _jc^* represents the effects of observed covariates, namely DIF effects, on the response to item j <cit.>.  
The different scientific interpretations motivate us to develop different identifiability conditions for β_j0^* and _jc^*, respectively.
Specifically, we propose a centering condition on ^* to ensure the identifiability of the intercept β_j0^* for all items j ∈ [q].
On the other hand, to identify the covariate effects β_jc^*, a natural idea is to impose the covariate effects _jc^* for all items j ∈ [q] to be sparse,

as shown in many regularized methods and item purification methods <cit.>. 
In <cit.>, an interpretable identifiability condition is proposed for selecting sparse covariate effects, yet this condition is specific to uni-dimensional covariates.
Motivated by <cit.>, we propose the following minimal ℓ_1 condition applicable to general cases where the covariates are multi-dimensional. To better present the identifiability conditions, we write = (_0, _1, …, _p_*) ∈^K× p and define _c = (_1, …, _p_*) ∈^K × p_* as the part applied to the covariate effects.



(i) ∑_i=1^n_i^* = 0_K.
(ii)
∑_j=1^q _jc^* _1 < ∑_j=1^q _jc^* - _c^_j^*_1 for any _c≠. 


Condition <ref>(i) assumes the latent abilities ^* are centered to ensure the identifiability of the intercepts _j0^*'s, which is commonly assumed in the item response theory literature <cit.>. 
Condition <ref>(ii) is motivated by practical applications. For instance, in educational testing, practitioners need to identify and remove biased test items, correspondingly, items with non-zero covariate effects (β_js^* ≠ 0). In practice, most of the designed items are unbiased, and therefore, it is reasonable to assume that the majority of items have no covariate effects, that is, the covariate effects _jc^*'s are sparse <cit.>.

Next, we present a sufficient and necessary condition for Condition <ref>(ii) to hold.


    Condition <ref>(ii) holds if and only if for any ∈^K∖{0_K }, 
    
    ∑_j=1^q|^_j^*|I(β_js^*= 0)>∑_j=1^q(β_js^*)^_j^* I(β_js^*≠ 0),  ∀ s∈[p_*].




Proposition <ref> implies that Condition <ref>(ii) holds when {j: β_js^* ≠ 0} is separated into {j: β_js^* > 0} and {j: β_js^* < 0} in a balanced way.
 With diversified signs of β_js^*, Proposition <ref> holds when a considerable proportion of test items have no covariate effect (β_js^* ≠ 0). 
For example, when _j^* = m 1_K^(k) with m >0, 
Condition <ref>(ii) holds if and only if ∑_j=1^q |m| { - I(β_js^*/m > 0) + I(β_js^*/m ≤  0)} > 0 and ∑_j=1^q |m| {- I(β_js^*/m ≥ 0) + I(β_js^*/m <  0)} <0. With slightly more than q/2 items correspond to β_js^* = 0, Condition <ref>(ii) holds. Moreover, if #{j: β_js^* > 0} and #{j: β_js^* < 0} are comparable, then Condition <ref>(ii) holds even when less than q/2 items correspond to β_js^* = 0 and more than q/2 items correspond to β_js^* ≠ 0. 

Though assuming a “sparse” structure, our assumption here differs from existing high-dimensional literature. 
In high-dimensional regression models,
the covariate coefficient when regressing the dependent variable on high-dimensional covariates, is often assumed to be sparse, with the proportion of the non-zero covariate coefficients asymptotically approaching zero. In our setting, Condition <ref>(ii) allows for relatively dense settings where the proportion of items with non-zero covariate effects is some positive constant.




To perform simultaneous estimation and inference on ^* and ^*, we consider the following identifiability conditions to address the second identifiability issue. 

 
 
  (i) (^*)^^* is diagonal. (ii) (^*)^⊺^* is diagonal. (iii) n^-1 (^*)^^*=q^-1(^*)^⊺^*.
 


 Condition <ref> is a set of widely used identifiability conditions in the factor analysis literature <cit.>. For practical and theoretical benefits, we impose Condition <ref> to address the identifiability issue related to . It is worth mentioning that this condition can be replaced by other identifiability conditions. 
 For true parameters satisfying any identifiability condition, we can always find a transformation such that the transformed parameters satisfy our proposed Conditions <ref>–<ref> and the proposed estimation method and theoretical results in the subsequent sections still apply, up to such a transformation. 

 

 §.§ Joint Maximum Likelihood Estimation

	


In this section, we introduce a joint-likelihood-based estimation method for the covariate effects ,  the latent factors , and factor loadings  simultaneously.


Incorporating Conditions <ref>–<ref> into the estimation procedure, 
we obtain the maximum joint-likelihood-based estimators for ^* = (^*, ^*, ^*) that satisfy the proposed identifiability conditions. 



 With Condition <ref>, we address the identifiability issue related to the transformation matrix . Specifically, for any parameters = (, , ),  there exists a matrix ^* = (_0^*, _c^*) with
 _c^* = __c ∈^K× p_*∑_j=1^q _jc - _c^_j_1 and _0^*=-n^-1∑_i=1^n(_i + _c^*_i^c ) such that the transformed matrices ^* =  +(^*)^ and
 ^* =  - ^* satisfy Condition <ref>. The transformation idea naturally leads to the following estimation methodology for ^*.
To estimate ^* and ^* that satisfy Condition <ref>,
we first obtain the maximum likelihood estimator = ( , , ) by

    = ∈Ω_ϕ -L(|, ),

where the parameter space Ω_ϕ is given as Ω_ϕ={:_max≤ C} for some large C. 

To solve (<ref>), we employ an alternating minimization algorithm. Specifically, for steps t = 0, 1,…, we compute 

    ^(t+1), ^(t+1)   = _∈^q × K, ∈^q× p - L(|, ^(t),  , ); 
    ^(t+1)   = _∈^n × K - L(|^(t+1), , ^(t+1), ),

until the quantity max{^(t+1) - ^(t)_F, ^(t+1) - ^(t)_F, ^(t+1) - ^(t)_F } is less than some pre-specified tolerance value for convergence. 
We then estimate _c by minimizing the ℓ_1-norm 
 
    _c = __c∈^K× p_*∑_j=1^q _jc - _c^_j _1.

Next, we estimate _0 = - n^-1∑_i=1^n (_i+_c_i^c) and let = (_0, _c).
 Given the estimators  , , and , we then construct  



     
    ^*  =  -    =  +^


such that Condition <ref> holds. 



Recall that Condition <ref> addresses the identifiability issue related to the invertible matrix . Specifically, 
for any parameters (, ), there exists a matrix ^* such that Condition <ref> holds for ^*=( +  (^*)^)^* and ^* =  (^*)^-. 
Let = (ϱ_1, …, ϱ_K) be a diagonal matrix that contains the K eigenvalues of (nq)^-1 (^ )^1/2 ( + ^)^ ( + ^) (^ )^1/2 and let  be a matrix that contains its corresponding eigenvectors.
We set ^* = (q^-1^)^1/2 ^-1/4. To further estimate ^* and ^*, we need to obtain an estimator for the invertible matrix ^*.
Given the maximum likelihood estimators obtained in (<ref>) and  in (<ref>), we estimate ^* via = (q^-1^)^1/2 ^-1/4 where   and  are matrices that contain the eigenvalues and eigenvectors of (nq)^-1 (^ )^1/2 (+ ^)^ ( + ^) (^ )^1/2, respectively. With   and , we now obtain the following transformed estimators that satisfy Condition <ref>:
    
    ^* =  (^)^-1  ^* = ( + ^) .


To quantify the uncertainty of the proposed estimators, we will show that the proposed estimators are asymptotically normally distributed.   Specifically, in 
Theorem <ref> of Section <ref>, we establish the asymptotic normality result for _j^*, which allows us to make inference on the covariate effects _j^*. 
Moreover, as the latent factors _i^* and factor loadings _j^* often have important interpretations in domain sciences, we are also interested in the inference on parameters _i^* and _j^*. In Theorem <ref>, we also derive the asymptotic distributions for estimators _i^* and _j^*, providing inference results for parameters _i^* and _j^*. 





§ THEORETICAL RESULTS





We propose a novel framework to establish the estimation consistency and asymptotic normality for the proposed joint-likelihood-based estimators ^* = (^*, ^*, ^*) in Section <ref>. To establish the theoretical results for ^*, we impose the following regularity assumptions.






	
There exist constants M >0, κ  >0 such that:
 



(i) _u^* = lim_n→∞ n^-1 (^*)^^* exists and is positive definite.  For i ∈ [n],  _i^*_2≤ M.

 (ii) _γ^* = lim_q→∞ q^-1 (^*)^^* exists and is positive definite. For j ∈ [q],  _j^*_2≤ M. 

 (iii) _x = lim_n→∞n^-1∑_i=1^n_i_i^ exists and 1/ κ^2 ≤λ_min (_x)  ≤λ_max (_x) ≤κ^ 2. For i ∈ [n],  max_i_i_∞≤ M. 
 
(iv) _ux^* = lim_n→∞ n^-1∑_i=1^n _i^* _i^ exists and _ux^*_x^-1_1,∞≤ M. The eigenvalues of (_u^*-_ux^*_x^-1(_ux^*)^)_γ^* are distinct.


	



		Assumptions <ref> is commonly used in the factor analysis literature. In particular, 
  Assumptions <ref>(i)–(ii) correspond to Assumptions A-B in <cit.> under linear factor models, ensuring the compactness of the parameter space on ^* and ^*. 
  
  Under nonlinear factor models, such conditions on compact parameter space are also commonly assumed <cit.>. Assumption <ref>(iii) is standard regularity conditions for the nonlinear setting that is needed to establish the concentration of the gradient and estimation error for the model parameters when p diverges. 
        In addition, Assumption <ref>(iv) is a crucial identification condition; similar conditions have been imposed in the existing literature such as
       Assumption G in <cit.> in the context of linear factor models and  Assumption 6 in <cit.> in the context of  nonlinear factor models without covariates.
      
  
 
	
	
	
For any i ∈ [n] and j ∈ [q], assume that l_ij (·) is three times differentiable, and we denote the first, second, and third order derivatives of l_ij(w_ij) with respect to w_ij as l_ij^' (w_ij), l_ij^'' (w_ij), and l_ij^''' (w_ij), respectively. There exist M > 0 and ξ≥ 4 such that  (|l_ij^' (w_ij) |^ξ) ≤ M and 
 |l_ij^' (w_ij)| is sub-exponential with l_ij^' (w_ij)_φ_1≤ M. Furthermore, we assume {l_ij^' (w_ij^*)} = 0. 
 Within a compact space of w_ij, we have b_L ≤ -l_ij^'' (w_ij) ≤ b_U and |l_ij^'''(w_ij)| ≤ b_U for b_U > b_L > 0. 

There exist M > 0 and ξ≥ 4 such that  (|l_ij^' (w_ij) |^ξ) ≤ M and 
 |l_ij^' (w_ij)| is sub-exponential with l_ij^' (w_ij)_φ_1≤ M. we assume {l_ij^' (w_ij^*)} = 0 and {X_ij l_ij^' (w_ij^*)} = 0.
 
				 
	(ii)  There exist M > 0 and ξ≥ 4 such that  (|l_ij^' (w_ij) |^ξ) ≤ M and 
 |l_ij^' (w_ij)| is sub-exponential with l_ij^' (w_ij)_φ_1≤ M. We assume {l_ij^' (w_ij^*)} = 0 and {X_ij l_ij^' (w_ij^*)} = 0 for all i ∈ [n] and j ∈ [q]. 
			 
	(iii) Within a compact space of w_ij, b_L ≤ -l_ij^'' (w_ij) ≤ b_U and |l_ij^'''(w_ij)| ≤ b_U for b_U > b_L > 0. 
		
  
  There exist constants b_U > b_L > 0 such that b_L ≤ -l_ij^'' (w_ij) ≤ b_U within a compact space of w_ij. 
  
	
	(iv) 
		 |l_ij^'''(w_ij)| ≤ b_U  within a compact space of w_ij.	



Assumption <ref> assumes smoothness on the log-likelihood function l_ij(w_ij). 
In particular, it assumes sub-exponential distributions and finite fourth-moments of the first order derivatives l^'_ij (w_ij). 
For commonly used linear or nonlinear factor models, the assumption is not restrictive and can be satisfied with a large ξ. 
For instance, consider the logistic
model with l_ij^'(w_ij) = Y_ij - exp(w_ij)/{1+exp(w_ij)}, we have |l_ij^'(w_ij)|≤ 1 and ξ can be taken as ∞.
The boundedness conditions for l_ij^''(w_ij) and l_ij^'''(w_ij) 
are necessary to guarantee the convexity of the joint likelihood function. 
In a special case of linear factor models, 

l_ij^''(w_ij) is a constant and the boundedness conditions naturally hold.

For popular nonlinear models such as logistic factor models, probit factor models, and Poisson factor models, the boundedness of l_ij^''(w_ij) and l_ij^'''(w_ij) can also be easily verified. 


    

    For ξ  specified in Assumption 2 and a sufficiently small ϵ >0, 
we assume   as n,q,p →∞, 
    
    p/√(n ∧(pq))   (nq)^ϵ+3/ξ→ 0.




Assumption <ref> is needed to ensure that the derivative of the likelihood function equals zero at the maximum likelihood estimator with high probability, a key property in the theoretical analysis. 
In particular, we need the estimation errors of all model parameters to converge to 0 uniformly with high probability. Such uniform convergence results involve delicate analysis of the convexity of the objective function, for which technically we need Assumption <ref>. For most of the popularly used generalized factor models, ξ can be taken as any large value as discussed above,  thus

(nq)^ϵ + 3/ξ is of a smaller order of √(n ∧ (pq)), given small ϵ.

Specifically, Assumption <ref> implies 
p=o(n^1/2∧ q) up to a small order term, an asymptotic regime that is reasonable for many educational assessments. 



Next, we impose additional assumptions crucial to establishing the theoretical properties of the proposed estimators. 

One challenge for theoretical analysis is to handle the dependence between the latent factors ^* and the design matrix . To address this challenge, we employ the following transformed ^0 that are orthogonal with , which plays an important role in establishing the theoretical results (see Supplementary Materials for details). In particular, for i ∈ [n], we let _i^0 = (^)^(_i^* -  ^_i). Here ^ = (q^-1(^*)^^*)^1/2 ^* (^*)^-1/4 and ^ = (^*)^(^)^-1, where ^* = (ϱ_1^*, …, ϱ_K^*) with diagonal elements being the K eigenvalues of 
(nq)^-1 ((^*)^^* )^1/2 (^*)^ (_n - _x) ^*  ((^*)^^* )^1/2
with _x =  (^)^-1^ and ^* containing the matrix of corresponding eigenvectors. 
Under this transformation for _i^0, we further define _j^0 = (^)^-1_j^* and _j^0 = _j^* + (^)^_j^* for j ∈ [q], and write _i^0 = ((_i^0)^ _i^)^ and w_ij^0 = (_j^0)^_i^0 + (_j^0)^_i. 
These transformed parameters _j^0's, _i^0's, and _j^0's give the same joint likelihood value as that of the true parameters _j^*'s, _i^*'s and _j^*'s, which facilitate our theoretical understanding of the joint-likelihood-based estimators. 





(i) For any j ∈ [q], -n^-1∑_i=1^n l_ij^''(w_ij^0) _i^0 (_i^0)^p→_jz^0 for some positive definite matrix _jz^0 and n^-1/2∑_i=1^n l_ij^'(w_ij^0) _i^0 d→(0, _jz^0).

(ii) For any i ∈ [n], -q^-1∑_j=1^q l_ij^''(w_ij^0) _j^0 (_j^0)^p→_iγ^0 for some positive definite matrix _iγ^0 and q^-1/2∑_j=1^q l_ij^'(w_ij^0) _j^0d→(0, _iγ^0).




    



Assumption <ref> is a generalization of Assumption F(3)-(4) in <cit.> for linear models to the nonlinear setting. 
  

Specifically, we need Assumption <ref>(i)
to derive the asymptotic distributions of the estimators _j^* and _j^*, and Assumption <ref>(ii) is used for establishing the asymptotic distribution of _i^*.
Note that these assumptions are imposed on the log-likelihood derivative functions evaluated at the true parameters w_ij^0, _i^0, and _j^0.

In general, for the popular generalized factor models, such assumptions hold with mild conditions. For example, under linear models, l_ij^'(w_ij) is the random error and l_ij^''(w_ij) is a constant. Then _jz^0 and _iγ^0 naturally exist and are positive definite followed by Assumption <ref>. The limiting distributions of n^-1/2∑_i=1^n l_ij^'(w_ij^0) _i^0 and q^-1/2∑_j=1^q l_ij^'(w_ij^0) _j^0 can be derived by the central limit theorem under standard regularity conditions. Under logistic and probit models, l_ij^'(w_ij) and l_ij^''(w_ij) are both finite inside a compact parameters space and similar arguments can be applied to show the validity of Assumption <ref>.


We present the following assumption to establish the theoretical properties of the transformed matrix  as defined in (<ref>). 
In particular, we define ^0=(^)^^ and write ^0 = (_0^0, …, _p_*^0)^. Note that the estimation problem of (<ref>) is related to the median regression problem with measurement errors. To understand the properties of this estimator, following existing M-estimation literature <cit.>, we define ψ_js^0()=_j^0{β_js^0+(_j^0)^(-_s^0)} and
    χ_s()=∑_j=1^qψ_js^0() for j ∈ [q] and s ∈ [p_*]. We further define a perturbed version of ψ_js^0(), denoted as ψ_js(,_js), as follows:
     
    ψ_js(,_js)=(_j^0+[_js/√(n)]_[1:K]) {β_js^0+[_js/√(n)]_K+1-(_j^0+[_js/√(n)]_[1:K])^(-_s^0)}, s∈[p_*]

    where the perturbation

    _js   =[        _K         0;         0 (_s^(p))^ ](-∑_i=1^nl_ij^''(w_ij^0)_i^0(_i^0)^)^-1(√(n)∑_i=1^nl_ij^'(w_ij^0)_i^0),
  
    is asymptotically normally distributed by Assumption <ref>. 
    We define χ̂_s()=∑_j=1^qψ_js(,_js).
    




      
    For χ_s(), we assume that there exists some constant c>0 such that min_≠|q^-1χ_s()|> c holds for all s∈[p_*]. Assume there exists _s0 for each s∈[p_*] such that χ̂_s(_s0) = 0 with p√(n)_s0→0. In a neighbourhood of _s0, χ̂_s() has a nonsingular derivative such that {q^-1∇_χ̂_s(_s0)}^-1=O(1) and q^-1|∇_χ̂_s()-∇_χ̂_s(_s0)|≤ k|-_s0|. We assume ι_nq,p:=max{_s0,q^-1∑_j=1^qψ_js(_s0,_js)}=o ((p√(n))^-1). 



Assumption <ref> is crucial in addressing the theoretical difficulties of establishing the consistent estimation for ^0, a challenging problem related to median regression with weakly dependent measurement errors. 

In Assumption <ref>, we treat the minimizer of |∑_j=1^qψ(,_js)| as an M-estimator and adopt the Bahadur representation results in <cit.> for the theoretical analysis. 

For an ideal case where _js are independent and normally distributed with finite variances, which corresponds to the setting in median regression with measurement errors <cit.>, these assumptions can be easily verified.
Assumption <ref> discusses beyond such an ideal case and covers general settings. 
In addition to independent and Gaussian measurement errors, this condition also accommodates the case when _js are asymptotically normal and weakly dependent with finite variances, as implied by Assumption <ref> and the conditional independence of Y_ij. 


We want to emphasize that Assumption <ref> allows for both sparse and dense settings of the covariate effects. Consider an example of K=p=1 and γ_j=1 for j∈ [q]. Suppose β_js^* is zero for all j ∈ [q_1] and nonzero otherwise. 
Then this condition is satisfied as long as #{j: β_js^* > 0} and #{j: β_js^* < 0} are comparable, even when the sparsity level q_1 is small.
 
 

Under the proposed assumptions, we next present our main theoretical results.


    
	Suppose the true parameters ^* = (^*, ^*, ^*) satisfy identifiability conditions <ref>–<ref>. Under Assumptions <ref>–<ref>, we have
 
 
    q^-1^* - ^* _F^2     = O_p( p^2log qp/n + plog n/q);

 if we further assume  p^3/2(nq)^ϵ+3/ξ(p^1/2 n^-1/2 +q^-1/2)=o(1), then we have
 
    n^-1^* - ^*_F^2     =  O_p( plog qp/n + log n/q); 
    
              q^-1^* - ^* _F^2    =O_p( plog qp/n + log n/q).



Theorem <ref> presents the average convergence rates of ^*. 
Consider an oracle case with ^* and ^* known, the estimation of ^* reduces to an M-estimation problem. For M-estimators under general parametric models, 
it can be shown that the optimal convergence rates in squared ℓ_2-norm is O_p(p/n) under p (log p)^3 /n → 0 <cit.>. 
In terms of our average convergence rate on ^*, the first term in (<ref>), n^-1 p^2log (qp), approximately matches the convergence rate O_p(p/n) up to a relatively small order term of plog (qp). The second term in (<ref>), q^-1plog n, is mainly due to the estimation error for the latent factor ^*. 
In educational applications, it is common to assume the number of subjects n is much larger than the number of items q.
Under such a practical setting with n ≫ q and p relatively small, the term q^-1log n in (<ref>) dominates in the derived convergence rate of ^*, which matches with the optimal convergence rate O_p(q^-1) for factor models without covariates <cit.>  up to a small order term. 



  The additional condition  p^3/2(nq)^ϵ+3/ξ(p^1/2 n^-1/2 +q^-1/2)=o(1) in Theorem <ref> is used to handle the challenges related to the invertible matrix  that affects the theoretical properties of ^* and ^*. It is needed for establishing the estimation consistency of ^* and ^* but not for that of ^*. With sufficiently large ξ and small ϵ, this assumption is approximately p = o(n^1/4∧ q^1/3) up to a small order term.
  



 
One challenge in establishing the estimation consistency for ^* arises from the unrestricted dependence structure between ^* and . If we consider the ideal case where the columns of ^* and  are orthogonal, i.e., (^*)^ = 0_K× p, then we can achieve comparable or superior convergence rates with less stringent assumptions. Specifically, with Assumptions <ref>–<ref> only, we can obtain the same convergence rates for  ^* and ^* as in (<ref>) and (<ref>), respectively. Moreover, with Assumptions <ref>–<ref>, the average convergence rate for the consistent estimator of ^* is O_p( n^-1plog qp + q^-1log n), which is tighter than (<ref>) by a factor of p.



 With estimation consistency results established, we next derive the asymptotic normal distributions for the estimators, which enable us to perform statistical inference on the true parameters.



Suppose the true parameters ^* = (^*, ^*, ^*) satisfy identifiability conditions <ref>–<ref>. Under Assumptions <ref>–<ref>, 
we have the asymptotic distributions  
as follows. Denote ζ_nq,p^-2 = n^-1plog qp+q^-1log n. 
If p^3/2√(n)(nq)^3/ξζ_nq,p^-2→ 0, for any j ∈ [q] and ∈^p with _2 = 1, 
  
    √(n)^ (_β,j^*)^-1/2(_j^*-_j^*)    d→(0,1),
 
      where _β, j^* =  (-(^0)^, _p) (_jz^0)^-1_jz^0(_jz^0)^-1(-(^0)^, _p)^, 
      and  for any j ∈ [q], 
        
    √(n)(_γ, j^*)^-1/2(_j^*-_j^*)      d→(, _K),

where _γ, j^* = ^(_K,  0)  (_jz^0)^-1_jz^0 (_jz^0)^-1 (_K,  0)^ (^)^. 
Furthermore, for any i ∈ [n], if q= O(n) and p^3/2√(q)(nq)^3/ξζ_nq,p^-2→ 0,
 
    √(q)(_u, i^*)^-1/2(_i^*-_i^*)     d→(, _K ),

  where _u,i^* = (^)^- (_iγ^0)^-1_iγ^0(_iγ^0)^-1 (^)^-1. 
  
    

    The asymptotic covariance matrices in Theorem <ref>  can be consistently estimated. Due to the space limitations, we defer the construction of the consistent estimators _β, j^*,  _γ, j^*, and _u, i^* to Supplementary Materials. 
Theorem <ref> provides the asymptotic distributions for all individual estimators. 

In particular, with the asymptotic distributions and the consistent estimators _β, j^* for the asymptotic covariance matrices, we can perform hypothesis testing on β_js^* for j ∈ [q] and s ∈ [p_*]. We reject the null hypothesis β_js^* = 0 at significance level α if | √(n) (σ̂_β,js^*)^-1β̂_js^* | >Φ^-1 ( 1- α/2), where (σ̂_β,js^*)^2 is the (s+1)-th diagonal entry in _β, j^*. 


For the asymptotic normality of _j^*, the condition p^3/2√(n)(nq)^3/ξ(n^-1plog qp+q^-1log n )→ 0 together with Assumption <ref> gives p = o{n^1/5∧ (q^2/n)^1/3} up to a small order term, and further implies n ≪ q^2, which is consistent with established conditions in the existing factor analysis literature <cit.>.

For the asymptotic normality of _i^*, the additional condition that q = O(n) is a reasonable assumption in educational applications where the number of items q is much fewer than the number of subjects n.
In this case, the scaling conditions imply p = o{ q^1/3∧ (n^2/q)^1/5} up to a small order term. 
Similarly for the asymptotic normality of _j^*, 
the proposed conditions give p = o{n^1/5∧ (q^2/n)^1/3} up to a small order term. 

 


 Similar to the discussion in Remark <ref>, the challenges arising from the unrestricted dependence between ^* and  also affect the derivation of the asymptotic distributions for the proposed estimators. If we consider the ideal case with (^*)^ = 0_K× p,  we can establish the asymptotic normality for all individual estimators under Assumptions <ref>–<ref> only and weaker scaling conditions. 
 Specifically, when (^*)^ = 0_K× p, 
 the scaling condition becomes p√(n)(nq)^3/ξ(n^-1plog qp+q^-1log n )→ 0 for deriving asymptotic normality of _j^* and _j^*, which is milder than that for (<ref>) and (<ref>).
 







§ SIMULATION STUDY


In this section, we study the finite-sample performance of the proposed joint-likelihood-based estimator. 
We focus on the logistic latent factor model in (<ref>) with p_ij (y | w_ij) = exp(w_ijy)/{1+exp(w_ij)}, where w_ij = (_j^*)^_i^* + (_j^*)^_i.  
The logistic latent factor model is commonly used in the context of educational assessment and is also referred to as the item response theory model <cit.>. 

We apply the proposed method to estimate ^* and 

perform statistical inference on testing the null hypothesis β_js^* = 0.  





We start with presenting the data generating process. 
We set the number of subjects n = {300, 500,  1000, 1500, 2000}, the number of items q = {100, 300, 500}, the covariate dimension p = {5, 10, 30}, and the factor dimension K = 2, respectively. 
We jointly generate _i^c and _i^* from (0, ) where _ij = τ^|i - j| with τ∈{0, 0.2, 0.5, 0.7}. 

In addition, we set the loading matrix _[,k]^* = 1_k^(K)⊗_k, where  ⊗ is the Kronecker product and _k is a (q/K)-dimensional vector with each entry generated independently and identically from Unif[0.5, 1.5].
For the covariate effects ^*, we set the intercept terms to equal β_j0^* = 0. 
For the remaining entries in ^*, we consider the following two settings: (1) sparse setting: β_js^* = ρ for s = 1, …, p and j = 5s-4, …, 5s and other β_js^* are set to zero; (2) dense setting: β_js^* = ρ for s = 1, …, p and j =  R_s q/5 + 1, …, (R_s + 1) q/5 with  R_s = s - 5 ⌊s/5⌋, and other β_js^* are set to zero. 
Here, the signal strength is set as ρ∈{0.3, 0.5}. Intuitively, in the sparse setting, we set 5 items to be biased for each covariate whereas in the dense setting, 20% of items are biased items for each covariate. 





For better empirical stability, after reaching convergence in the proposed  alternating maximization algorithm and transforming the obtained MLEs into ones that satisfy Conditions <ref>–<ref>, we repeat another round of maximization and transformation. 
We take the significance level at 5% and calculate the averaged type I error based on all the entries β_js^* = 0 and the averaged power for all non-zero entries, over 100 replications.
The averaged hypothesis testing results are presented in Figures <ref>–<ref> for p = 5 and p = 30, across different settings. 
Additional numerical results for p=10 are presented in the Supplementary Materials.
 























From Figures <ref>–<ref>, we observe that the type I errors are well controlled at the significance level 5%, which is consistent with the asymptotic properties of ^* in Theorem <ref>. 
Moreover, the power increases to one as the sample size n increases across all of the settings we consider. 
Comparing the left panel (ρ = 0.3) to the right panel (ρ = 0.5) in Figures 3–6, we see that the power increases as we increase the signal strength ρ. 
Comparing the plots in Figures <ref>–<ref> to the corresponding plots in Figures <ref>–<ref>, we see that the powers under the sparse setting (Figures <ref>–<ref>) are generally higher than that of the dense setting (Figures <ref>–<ref>).  
Nonetheless, our proposed method is generally stable under both sparse and dense settings. 
In addition, we observe similar results when we increase the
covariate dimension p from p=5 (Figures <ref> and <ref>) to p = 30 (Figures <ref> and <ref>).


We refer the reader to the Supplementary Materials for additional numerical results for p=10.  
Moreover, we observe similar results when we increase the test length q from q=100 (top row) to q=500 (bottom row) in 


Figures <ref>–<ref>. 

In terms of the correlation between  and ^*, we observe that while the power converges to one as we increase the sample size, the power decreases as the correlation τ increases. 







§ DATA APPLICATION


We apply our proposed method to analyze the Programme for International Student Assessment (PISA) 2018 data[The data can be downloaded from: https://www.oecd.org/pisa/data/2018database/].  PISA is a worldwide testing program that compares the academic performances of 15-year-old students across many countries <cit.>. More than 600,000 students from 79 countries/economies, representing a population of 31 million 15-year-olds, participated in this program.
The PISA 2018 used the computer-based assessment mode and the assessment lasted two hours for each student, with test items mainly evaluating students' proficiency in mathematics, reading, and science domains. 
A total of 930 minutes of test items were used and each student took different combinations of the test items. In addition to the assessment questions, background questionnaires were provided to collect students' information.



In this study, we focus on PISA 2018 data from Taipei. The observed responses are binary, indicating whether students' responses to the test items are correct, and we use the popular item response theory model with the logit link <cit.>. Due to the block design nature of the large-scale assessment, each student was only assigned to a subset of the test items, and for the Taipei data, 86% response matrix is unobserved. Note that this missingness can be considered as conditionally independent of the responses given the students' characteristics. 
Our proposed method and inference results naturally accommodate such missing data and can be directly applied. 
Specifically, to accommodate the incomplete responses, we can modify the 

joint log-likelihood function in (<ref>) into L^obs(|, , ,) = ∑_i=1^n ∑_j ∈_i l_ij(_j^⊺_i + _j^⊺_i),
where _i defines the set of questions to which the responses from student i are observed. 

In this study, we include gender and 8 variables for school strata as covariates (p_*=9). These variables record whether the school is public, in an urban place, etc. After data preprocessing, we have n = 6063 students and q = 194 questions. 
Following the existing literature <cit.>, we take K = 3 to interpret the three latent abilities measured by the math, reading, and science questions. 


We apply the proposed method to estimate the effects of gender and school strata variables on students' responses. We obtain the estimators of the gender effect for each PISA question and construct the corresponding 95% confidence intervals.
The constructed 95% confidence intervals for the gender coefficients are presented in Figure <ref>. There are 10 questions highlighted in red as their estimated gender effect is statistically significant after the Bonferroni correction.

Among the reading items, there is only one significant item and the corresponding confidence interval is below zero, indicating that this question is biased towards female test-takers, conditioning on the students' latent abilities.
Most of the confidence intervals corresponding to the biased items in the math and science sections are above zero, indicating that
these questions are biased towards male test-takers.

In social science research, it is documented that female students typically score better than male students during reading tests, while male students often outperform female students during math and science tests <cit.>. Our results indicate that there may exist potential measurement biases resulting in such an observed  gender gap in educational testing. Our proposed method offers a useful tool to identify such biased test items, thereby contributing to enhancing testing fairness by providing practitioners with valuable information for item calibration.




 











 

 To further illustrate the estimation results, Table <ref> lists the p-values for testing the gender effect for each of the identified 10 significant questions, along with
 the proportions of female and male test-takers who answered each question correctly.
We can see that the signs of the estimated gender effect by our proposed method align with the disparities in the reported proportions between females and males.

For example, 
the estimated gender effect corresponding to the item “CM496Q01S Cash Withdrawal” is positive with a p-value of 2.77× 10^-7, implying that this question is statistically significantly biased towards male test-takers. This is consistent with the observation that in Table <ref>, 58.44% of male students correctly answered this question, which exceeds the proportion of females,  51.29%.


 
Besides gender effects, we estimate the effects of school strata on the students' response and present the
 point and interval estimation results in the left panel of Figure <ref>. All the detected biased questions are from math and science sections, with 6 questions for significant effects of whether attending public school and 5 questions for whether residing in rural areas.

  To further investigate the importance of controlling for the latent ability factors, we compare results from our proposed method with the latent factors, to the results from directly regressing responses on covariates without latent factors.
From the right panel of Figure <ref>, we can see that without conditioning on the latent factors, there are excessive items detected for the covariate of whether the school is public or private. On the other hand, there are no biased items detected if we only apply generalized linear regression to estimate the effect of the covariate of whether the school is in rural areas. 









§ DISCUSSION



In this work, we study the covariate-adjusted generalized factor model that has wide interdisciplinary applications such as educational assessments and psychological measurements. 
In particular, new identifiability issues arise due to the incorporation of covariates in the model setup. To address the issues and identify the model parameters, we propose novel and interpretable conditions,  crucial for developing the estimation approach and inference results. With model identifiability guaranteed, we propose a computationally efficient joint-likelihood-based estimation method for model parameters.
Theoretically, we obtain the estimation consistency and asymptotic normality for not only the covariate effects but also latent factors and factor loadings. 


There are several future directions motivated by the proposed method. 
In this manuscript, we focus on the case in which p grows at a slower rate than the number of subjects n and the number of items q, a common setting in educational assessments.  
It is interesting to further develop estimation and inference results under the high-dimensional setting in which p is larger than n and q. 


Moreover, in this manuscript, we assume that the dimension of the latent factors K is fixed and known. 
One possible generalization is to allow K to grow with n and q.  
Intuitively, an increasing latent dimension K makes the identifiability and inference issues more challenging due to the increasing degree of freedom of the transformation matrix. 
With the theoretical results in this work, another interesting related problem is to further develop simultaneous inference on group-wise covariate coefficients, which we leave for future investigation.






agsm




