Comprehensive Simulation Study Research Articles

Maintaining the validity of inference from linear mixed models in stepped-wedge cluster randomized trials under misspecified random-effects structures.

Linear mixed models are commonly used in analyzing stepped-wedge cluster randomized trials. A key consideration for analyzing a stepped-wedge cluster randomized trial is accounting for the potentially complex correlation structure, which can be achieved by specifying random-effects. The simplest random effects structure is random intercept but more complex structures such as random cluster-by-period, discrete-time decay, and more recently, the random intervention structure, have been proposed. Specifying appropriate random effects in practice can be challenging: assuming more complex correlation structures may be reasonable but they are vulnerable to computational challenges. To circumvent these challenges, robust variance estimators may be applied to linear mixed models to provide consistent estimators of standard errors of fixed effect parameters in the presence of random-effects misspecification. However, there has been no empirical investigation of robust variance estimators for stepped-wedge cluster randomized trials. In this article, we review six robust variance estimators (both standard and small-sample bias-corrected robust variance estimators) that are available for linear mixed models in R, and then describe a comprehensive simulation study to examine the performance of these robust variance estimators for stepped-wedge cluster randomized trials with a continuous outcome under different data generators. For each data generator, we investigate whether the use of a robust variance estimator with either the random intercept model or the random cluster-by-period model is sufficient to provide valid statistical inference for fixed effect parameters, when these working models are subject to random-effect misspecification. Our results indicate that the random intercept and random cluster-by-period models with robust variance estimators performed adequately. The CR3 robust variance estimator (approximate jackknife) estimator, coupled with the number of clusters minus two degrees of freedom correction, consistently gave the best coverage results, but could be slightly conservative when the number of clusters was below 16. We summarize the implications of our results for the linear mixed model analysis of stepped-wedge cluster randomized trials and offer some practical recommendations on the choice of the analytic model.

Approaches to Statistical Efficiency When Comparing the Embedded Adaptive Interventions in a SMART

Sequential, multiple assignment randomized trials (SMARTs), which assist in the optimization of adaptive interventions, are growing in popularity in education and behavioral sciences. This is unsurprising, as adaptive interventions reflect the sequential, tailored nature of learning in a classroom or school. Nonetheless, as is true elsewhere in education research, observed effect sizes in education-based SMARTs are frequently small. As a consequence, statistical efficiency is of paramount importance in their analysis. The contributions of this manuscript are twofold. First, we provide an overview of adaptive interventions and SMART designs for researchers in education science. Second, we propose four techniques that have the potential to improve statistical efficiency in the analysis of SMARTs. We demonstrate the benefits of these techniques in SMART settings both through the analysis of a SMART designed to optimize an adaptive intervention for increasing cognitive behavioral therapy delivery in school settings and through a comprehensive simulation study. Each of the proposed techniques is easily implementable, either with over-the-counter statistical software or through R code provided in Supplemental Material.

A Bayesian Change Point Model for Dynamic Alternative Transcription Start Site Usage During Cellular Differentiation.

ABSTRACT An alternative transcription start site (ATSS) is a major driving force for increasing the complexity of transcripts in human tissues. As a transcriptional regulatory mechanism, ATSS has biological significance. Many studies have confirmed that ATSS plays an important role in diseases and cell development and differentiation. However, exploration of its dynamic mechanisms remains insufficient. Identifying ATSS change points during cell differentiation is critical for elucidating potential dynamic mechanisms. For relative ATSS usage as percentage data, the existing methods lack sensitivity to detect the change point for ATSS longitudinal data. In addition, some methods have strict requirements for data distribution and cannot be applied to deal with this problem. In this study, the Bayesian change point detection model was first constructed using reparameterization techniques for two parameters of a beta distribution for the percentage data type, and the posterior distributions of parameters and change points were obtained using Markov Chain Monte Carlo (MCMC) sampling. With comprehensive simulation studies, the performance of the Bayesian change point detection model is found to be consistently powerful and robust across most scenarios with different sample sizes and beta distributions. Second, differential ATSS events in the real data, whose change points were identified using our method, were clustered according to their change points. Last, for each change point, pathway and transcription factor motif analyses were performed on its differential ATSS events. The results of our analyses demonstrated the effectiveness of the Bayesian change point detection model and provided biological insights into cell differentiation.

An evaluation of computational methods for aggregate data meta-analyses of diagnostic test accuracy studies

BackgroundA Generalized Linear Mixed Model (GLMM) is recommended to meta-analyze diagnostic test accuracy studies (DTAs) based on aggregate or individual participant data. Since a GLMM does not have a closed-form likelihood function or parameter solutions, computational methods are conventionally used to approximate the likelihoods and obtain parameter estimates. The most commonly used computational methods are the Iteratively Reweighted Least Squares (IRLS), the Laplace approximation (LA), and the Adaptive Gauss-Hermite quadrature (AGHQ). Despite being widely used, it has not been clear how these computational methods compare and perform in the context of an aggregate data meta-analysis (ADMA) of DTAs.MethodsWe compared and evaluated the performance of three commonly used computational methods for GLMM - the IRLS, the LA, and the AGHQ, via a comprehensive simulation study and real-life data examples, in the context of an ADMA of DTAs. By varying several parameters in our simulations, we assessed the performance of the three methods in terms of bias, root mean squared error, confidence interval (CI) width, coverage of the 95% CI, convergence rate, and computational speed.ResultsFor most of the scenarios, especially when the meta-analytic data were not sparse (i.e., there were no or negligible studies with perfect diagnosis), the three computational methods were comparable for the estimation of sensitivity and specificity. However, the LA had the largest bias and root mean squared error for pooled sensitivity and specificity when the meta-analytic data were sparse. Moreover, the AGHQ took a longer computational time to converge relative to the other two methods, although it had the best convergence rate.ConclusionsWe recommend practitioners and researchers carefully choose an appropriate computational algorithm when fitting a GLMM to an ADMA of DTAs. We do not recommend the LA for sparse meta-analytic data sets. However, either the AGHQ or the IRLS can be used regardless of the characteristics of the meta-analytic data.

CASi: A framework for cross-timepoint analysis of single-cell RNA sequencing data

Single-cell RNA sequencing (scRNA-seq) technology has been widely used to study the differences in gene expression at the single cell level, providing insights into the research of cell development, differentiation, and functional heterogeneity. Various pipelines and workflows of scRNA-seq analysis have been developed but few considered multi-timepoint data specifically. In this study, we develop CASi, a comprehensive framework for analyzing multiple timepoints’ scRNA-seq data, which provides users with: (1) cross-timepoint cell annotation, (2) detection of potentially novel cell types emerged over time, (3) visualization of cell population evolution, and (4) identification of temporal differentially expressed genes (tDEGs). Through comprehensive simulation studies and applications to a real multi-timepoint single cell dataset, we demonstrate the robust and favorable performance of the proposal versus existing methods serving similar purposes.

A Comparison of Additional Benefit Assessment Methods for Time-to-Event Endpoints Using Hazard Ratio Point Estimates or Confidence Interval Limits by Means of a Simulation Study.

For time-to-event endpoints, three additional benefit assessment methods have been developed aiming at an unbiased knowledge about the magnitude of clinical benefit of newly approved treatments. The American Society of Clinical Oncology (ASCO) defines a continuous score using the hazard ratio point estimate (HR-PE). The European Society for Medical Oncology (ESMO) and the German Institute for Quality and Efficiency in Health Care (IQWiG) developed methods with an ordinal outcome using lower and upper limits of the 95% HR confidence interval (HR-CI), respectively. We describe all three frameworks for additional benefit assessment aiming at a fair comparison across different stakeholders. Furthermore, we determine which ASCO score is consistent with which ESMO/IQWiG category. In a comprehensive simulation study with different failure time distributions and treatment effects, we compare all methods using Spearman's correlation and descriptive measures. For determination of ASCO values consistent with categories of ESMO/IQWiG, maximizing weighted Cohen's Kappa approach was used. Our research depicts a high positive relationship between ASCO/IQWiG and a low positive relationship between ASCO/ESMO. An ASCO score smaller than 17, 17 to 20, 20 to 24, and greater than 24 corresponds to ESMO categories. Using ASCO values of 21 and 38 as cutoffs represents IQWiG categories. We investigated the statistical aspects of the methods and hence implemented slightly reduced versions of all methods. IQWiG and ASCO are more conservative than ESMO, which often awards the maximal category independent of the true effect and is at risk of overcompensating with various failure time distributions. ASCO has similar characteristics as IQWiG. Delayed treatment effects and underpowered/overpowered studies influence all methods in some degree. Nevertheless, ESMO is the most liberal one. For the additional benefit assessment, the American Society of Clinical Oncology (ASCO) uses the hazard ratio point estimate (HR-PE) for their continuous score. In contrast, the European Society for Medical Oncology (ESMO) and the German Institute for Quality and Efficiency in Health Care (IQWiG) use the lower and upper 95% HR confidence interval (HR-CI) to specific thresholds, respectively. ESMO generously assigns maximal scores, while IQWiG is more conservative.This research provides the first comparison between IQWiG and ASCO and describes all three frameworks for additional benefit assessment aiming for a fair comparison across different stakeholders. Furthermore, thresholds for ASCO consistent with ESMO and IQWiG categories are determined, enabling a comparison of the methods in practice in a fair manner.IQWiG and ASCO are the more conservative methods, while ESMO awards high percentages of maximal categories, especially with various failure time distributions. ASCO has similar characteristics as IQWiG. Delayed treatment effects and under/-overpowered studies influence all methods. Nevertheless, ESMO is the most liberal one. An ASCO score smaller than 17, 17 to 20, 20 to 24, and greater than 24 correspond to the categories of ESMO. Using ASCO values of 21 and 38 as cutoffs represents categories of IQWiG.

Marshall-olkin extended inverted kumaraswamy distribution for modeling of wind speed data

Many important factors need to be taken into account when assessing the wind energy potential of a place. These include precisely determining the distribution to be used for the analysis of wind speed data and precisely estimating the model parameters. The current study has two main goals: the first is to estimate the parameters of the Marshall Olkin extended inverted Kumaraswamy (MOEIK) distribution, and the second is to use it to model wind speed data. For parameter estimation, four renowned estimation methods are used: maximum likelihood, Anderson Darling, Cramer von Misses, and ordinary least squares. A comprehensive simulation study is used to evaluate the behavior of these derived estimators, and maximum likelihood is found to efficiently estimate the parameters. For this study, wind speed data was collected from three stations in Sanghar, Tandoghulamali, and Umerkot. The results indicated that the MOIKw distribution provides efficient fitting based on the goodness-of-fit measures among competitive distributions.

Analysis of lead free CsSnBr3 based perovskite solar cells utilizing numerical modeling

In this study, we propose several CsSnBr3-based PSC configurations using the Solar Cell Capacitance Simulator (SCAPS-1D), incorporating various efficient Electron transport layers (ETLs) such as TiO2, PCBM, WS2, SnO2, ZnO, IGZO, C60, and Hole transport layers (HTLs) like CBTS, CFTS, CuO, CuI, Spiro-OMeTAD, PEDOT:PSS, P3HT, CuSbS2, CuSCN, and Cu2O. Numerical simulation results reveal that the device structure ITO/WS2/CsSnBr3/Cu2O/Au exhibits outstanding power conversion efficiency (PCE), retaining the closest photovoltaic parameter values among 70 different configurations. In this configuration, WS2 served as the ETL, and Cu2O acted as the HTL. This device achieved an outstanding peak PCE of 20.02%. It also boasted a high open circuit voltage (Voc) of 1.23 V, a short circuit current density (Jsc) of 19.32 mA cm−2, and an impressive fill factor (FF) of 84.18%. In comparison, devices utilizing materials like TiO2, PCBM, SnO2, ZnO, IGZO, and C60 yielded PCE values of 19.72, 19.73, 19.72, 19.73, 19.72, and 15.60%, respectively. Furthermore, for the seven best-performing configurations, we investigated the effects of CsSnBr3 absorber thickness, absorber-acceptor doping density (NA), conduction band offset (CBO), ETL doping density (ND), Capacitance–Voltage (C-V), Mott–Schottky (M-S) characteristics, generation and recombination rates, series resistance (Rse), shunt resistance (Rsh), temperature, current–voltage characteristics (J-V), and quantum efficiency (QE) on performance metrics. Our findings indicate that all seven ETLs, when combined with Cu2O HTL, can serve as excellent materials for fabricating high-efficiency CsSnBr3-based PSCs with the ITO/ETL/CsSnBr3/Cu2O/Au structure. To validate our results, we compared the simulation outcomes obtained with SCAPS-1D for the best seven CsSnBr3-PSC configurations with previously published research works. This comprehensive simulation study opens a promising avenue for the cost-effective production of high-performance, lead-free CsSnBr3-based PSCs, contributing to a greener and pollution-free environment.

A comprehensive simulation study of multi-junction solar cell

This study conducts comprehensive simulation analysis of typical triple-junction solar cells using Silvaco ATLAS. Initially, modeling and simulation of the typical triple-junction solar cells under the AM1.5 solar spectrum at 300 K are performed to characterize various performance parameters of the cells under one sun illumination. Subsequently, the impact of the thickness of the top and middle cell layers on the efficiency of the cells is analyzed. Additionally, the influence of doping concentration errors in individual sub-cells on the overall power deviation of the cells is investigated. Moreover, Al0.9Ga0.1As is employed to replace AlInP as the top FSF (Front Surface Field). Furthermore, the study analyzes the variation of internal radiation efficiency of the top and middle sub-cells of the solar cells with concentration factor and bias voltage. Finally, an assessment is made on the effects of different positions and densities of interface defects, as well as the introduction of class-acceptor Gaussian distribution state levels induced by defects, on the performance of the solar cells. The results indicate that this study provides valuable insights for optimizing multi-junction solar cells and analyzing internal physical phenomena.

Confidence and Outcome Expectations in Bilateral Negotiations–A Dynamic Model

AbstractThis work proposes and studies a dynamic model of two bargaining parties exchanging offers over time, considering their confidence about the share of the “pie” they obtain, which translates into expectations regarding the outcome of the bargaining process. The model predicts the sequence of offers as well as the final agreement for given confidence parameters. A mathematical analysis of the model shows the outcome is an Asymmetric Nash Bargaining Solution with exponents determined by the bargainers’ confidence. Moreover, a compensation effect can be found between confidence and risk aversion. This work also considers that confidence levels of bargainers might change during the negotiation, and we conduct a comprehensive simulation study to analyze the effect of such changes. Through Monte-Carlo simulation, we show that a bargainer is better off if its confidence increases, but the advantage is lost if the other party’s confidence increases in a similar way. In that case, concessions are smaller and negotiations last longer. Changing confidence parameters make the outcome harder to predict, as it will depend more on the final confidence than the initial one. The simulations also show that the average size of concessions, and therefore the final agreement, depend not only on whether confidence increases or decreases, but also on the change rate, with stronger effects observed when change accelerates towards the end of the process.

Methods for online calibration of Q-matrix and item parameters for polytomous responses in cognitive diagnostic computerized adaptive testing.

The ability to rapidly provide examinees with detailed and effective diagnostic information is a critical topic in psychology. Knowing what diagnostic criteria the examinees have met enables the practitioner to seek the solution to help them in a timely manner, and this can be achieved by cognitive diagnostic computerized adaptive testing (CD-CAT). However, the pervasive challenge of replenishing items in the CD-CAT item bank limits its practical application. Online calibration is a means to address item replenishment, but in CD-CAT, most existing online calibration methods that jointly calibrate the Q-matrix and item parameters of the new items are developed only for dichotomous responses and are time-consuming. Notably, previous studies pay no attention to polytomously scored items that are frequently observed in testing, even though they can offer additional evidence for the examinees' diagnosis. To fill this gap, we propose a SCAD-based method (SCAD-EM) to calibrate the Q-matrix and item parameters of the new items with polytomous response data in order to promote the application of CD-CAT in practice. The performance of the SCAD-EM was investigated in two comprehensive simulation studies and compared against the revised single-item estimation method (SIE-BIC). Results indicated that the SCAD-EM produces a higher calibration accuracy for the category-level Q-matrix and is computationally more efficient across all conditions, but it produces a lower calibration accuracy for the item-level Q-matrix. An empirical study further demonstrated the utility of the SCAD-EM and the SIE-BIC methods in calibrating new items with a real dataset. The advantages of the proposed method, its limitations, and possible future research directions are offered at the end.

Principal simple linear regression

In this paper, we propose a new method called principal simple linear regression for predicting a continuous response variable, $Y$ using a single continuous predictor, $X$ by utilizing multiple regression lines. This method is based on the theory of principal points and offers several advantages over classical simple linear regression methods, such as the ability to predict central, dispersion, and distributional tendencies of $Y$ on $X$, and simultaneous estimation. We provide the main properties, inferences, and limiting behavior of the estimators. Additionally, we conduct a comprehensive simulation study to validate our theoretical results. The model is also applied to real datasets to demonstrate its effectiveness.

Dynamic borrowing of historical controls adjusting for covariates in vaccine efficacy clinical trials.

Traditional vaccine efficacy trials usually use fixed designs and often require large sample sizes. Recruiting a large number of subjects can make the trial expensive, long, and difficult to conduct. A possible approach to reduce the sample size and speed up the development is to use historical controls. In this paper, we extend the robust mixture prior (RMP) approach (a well established approach for Bayesian dynamic borrowing of historical controls) to adjust for covariates. The adjustment is done using classical methods from causal inference: inverse probability of treatment weighting, G-computation and double-robust estimation. We evaluate these covariate-adjusted RMP approaches using a comprehensive simulation study and demonstrate their use by performing a retrospective analysis of a prophylactic human papillomavirus vaccine efficacy trial. Adjusting for covariates reduces the drift between current and historical controls, with a beneficial effect on bias, control of type I error and power.

Advancing unanchored simulated treatment comparisons: A novel implementation and simulation study.

Population-adjusted indirect comparisons, developed in the 2010s, enable comparisons between two treatments in different studies by balancing patient characteristics in the case where individual patient-level data (IPD) are available for only one study. Health technology assessment (HTA) bodies increasingly rely on these methods to inform funding decisions, typically using unanchored indirect comparisons (i.e., without a common comparator), due to the need to evaluate comparative efficacy and safety for single-arm trials. Unanchored matching-adjusted indirect comparison (MAIC) and unanchored simulated treatment comparison (STC) are currently the only two approaches available for population-adjusted indirect comparisons based on single-arm trials. However, there is a notable underutilisation of unanchored STC in HTA, largely due to a lack of understanding of its implementation. We therefore develop a novel way to implement unanchored STC by incorporating standardisation/marginalisation and the NORmal To Anything (NORTA) algorithm for sampling covariates. This methodology aims to derive a suitable marginal treatment effect without aggregation bias for HTA evaluations. We use a non-parametric bootstrap and propose separately calculating the standard error for the IPD study and the comparator study to ensure the appropriate quantification of the uncertainty associated with the estimated treatment effect. The performance of our proposed unanchored STC approach is evaluated through a comprehensive simulation study focused on binary outcomes. Our findings demonstrate that the proposed approach is asymptotically unbiased. We argue that unanchored STC should be considered when conducting unanchored indirect comparisons with single-arm studies, presenting a robust approach for HTA decision-making.

Recovery of High-energy Low-frequency Quasiperiodic Oscillations from Black Hole X-Ray Binary MAXI J1535–571 with a Hilbert–Huang Transform Method

We propose a method based on the Hilbert–Huang transform (HHT) to recover the high-energy waveform of low-frequency quasiperiodic oscillations (QPOs). Based on the method, we successfully obtain the modulation of the phase-folded light curve above 170 keV using the QPO phase reconstructed at lower energies in MAXI J1535–571 with Insight-HXMT observations. A comprehensive simulation study is conducted to demonstrate that such modulation indeed originates from the QPO. Thus, the highest energies turn out to significantly exceed the upper limit of ∼100 keV for QPOs reported previously using the Fourier method, marking the first opportunity to study QPO properties above 100 keV in this source. Detailed analyses of these high-energy QPO profiles reveal different QPO properties between the 30–100 and 100–200 keV energy ranges: the phase lag remains relatively stable, and the amplitude slightly increases below ∼100 keV, whereas above this threshold, soft phase lags and a decrease in amplitude are observed. Given the reports of a hard-tail detection in broad spectroscopy, we propose that the newly discovered QPO properties above 100 keV are dominated by the hard-tail component, possibly stemming from a relativistic jet. Our findings also indicate a strong correlation between the QPOs originating from the jet and corona, supporting the scenario of jet–corona coupling precession. We emphasize that our proposed HHT-based method can serve as an efficient manner in expanding the high-energy band for studying QPOs, thereby enhancing our understanding of their origin.

Spatial Functional Data analysis: Irregular spacing and Bernstein polynomials

Spatial Functional Data (SFD) analysis is an emerging statistical framework that combines Functional Data Analysis (FDA) and spatial dependency modeling. Unlike traditional statistical methods, which treat data as scalar values or vectors, SFD considers data as continuous functions, allowing for a more comprehensive understanding of their behavior and variability. This approach is well-suited for analyzing data collected over time, space, or any other continuous domain. SFD has found applications in various fields, including economics, finance, medicine, environmental science, and engineering. This study proposes new functional Gaussian models incorporating spatial dependence structures, focusing on irregularly spaced data and reflecting spatially correlated curves. The model is based on Bernstein polynomial (BP) basis functions and utilizes a Bayesian approach for estimating unknown quantities and parameters. The paper explores the advantages and limitations of the BP model in capturing complex shapes and patterns while ensuring numerical stability. The main contributions of this work include the development of an innovative model designed for SFD using BP, the presence of a random effect to address associations between irregularly spaced observations, and a comprehensive simulation study to evaluate models’ performance under various scenarios. The work also presents one real application of Temperature in Mexico City, showcasing practical illustrations of the proposed model.

Prospects of the Multi-channel Photometric Survey Telescope in the Cosmological Application of Type Ia Supernovae

The Multi-channel Photometric Survey Telescope (Mephisto) is a real-time, three-color photometric system designed to capture the color evolution of stars and transients accurately. This telescope system can be crucial in cosmological distance measurements of low-redshift (low-z, z ≲0.1) Type Ia supernovae (SNe Ia). To optimize the capabilities of this instrument, we perform a comprehensive simulation study before its official operation is scheduled to start. By considering the impact of atmospheric extinction, weather conditions, and the lunar phase at the observing site involving the instrumental features, we simulate light curves of SNe Ia obtained by Mephisto. The best strategy in the case of SN Ia cosmology is to take the image at an exposure time of 130 s with a cadence of 3 days. In this condition, Mephisto can obtain hundreds of high-quality SNe Ia to achieve a distance measurement better than 4.5%. Given the on-time spectral classification and monitoring of the Lijiang 2.4 m Telescope at the same observatory, Mephisto, in the whole operation, can significantly enrich the well-calibrated sample of supernovae at low-z and improve the calibration accuracy of high-z SNe Ia.

Do Recessions and Bear Markets Occur Concurrently across Countries? A Multinomial Logistic Approach

Abstract We introduce a novel multinomial logistic model for detecting and forecasting concurrent recessions and bear markets across multiple countries. Our framework leverages cross-country panel features and provides additional information for robust analysis. Through a comprehensive simulation study, we demonstrate the computational efficiency and accuracy of our model, even when handling multiple binary indicators. Applying our framework to empirical data from the United States, the UK, and Euro Area, we find that the multinomial logistic model produces superior medium-term forecasting of concurrent recession and bear market events across countries compared to multiple independent single logistic models. Additionally, our counterfactual analysis reveals that specific events, such as a recession and bear market in the United States, along with the tightening of financial conditions and a negative interest rate spread in the United States, increase the probability of concurrent and individual recession and bear market occurrences in the UK and Euro Area.

Semiparametric accelerated failure time models under unspecified random effect distributions

Accelerated failure time (AFT) models with random effects, a useful alternative to frailty models, have been widely used for analyzing clustered (or correlated) time-to-event data. In the AFT model, the distribution of the unobserved random effect is conventionally assumed to be parametric, often modeled as a normal distribution. Although it has been known that a misspecfied random-effect distribution has little effect on regression parameter estimates, in some cases, the impact caused by such misspecification is not negligible. Particularly when our focus extends to quantities associated with random effects, the problem could become worse. In this paper, we propose a semi-parametric maximum likelihood approach in which the random-effect distribution under the AFT models is left unspecified. We provide a feasible algorithm to estimate the random-effect distribution as well as model parameters. Through comprehensive simulation studies, our results demonstrate the effectiveness of this proposed method across a range of random-effect distribution types (discrete or continuous) and under conditions of heavy censoring. The efficacy of the approach is further illustrated through simulation studies and real-world data examples.

High-dimensional covariate-augmented overdispersed poisson factor model.

The current Poisson factor models often assume that the factors are unknown, which overlooks the explanatory potential of certain observable covariates. This study focuses on high dimensional settings, where the number of the count response variables and/or covariates can diverge as the sample size increases. A covariate-augmented overdispersed Poisson factor model is proposed to jointly perform a high-dimensional Poisson factor analysis and estimate a large coefficient matrix for overdispersed count data. A group of identifiability conditions is provided to theoretically guarantee computational identifiability. We incorporate the interdependence of both response variables and covariates by imposing a low-rank constraint on the large coefficient matrix. To address the computation challenges posed by nonlinearity, two high-dimensional latent matrices, and the low-rank constraint, we propose a novel variational estimation scheme that combines Laplace and Taylor approximations. We also develop a criterion based on a singular value ratio to determine the number of factors and the rank of the coefficient matrix. Comprehensive simulation studies demonstrate that the proposed method outperforms the state-of-the-art methods in estimation accuracy and computational efficiency. The practical merit of our method is demonstrated by an application to the CITE-seq dataset. A flexible implementation of our proposed method is available in the R package COAP.