Science, Faculty of
Statistics, Department of
DSpace
UBCV
Lee, Gahyung (Seren)
Doctor of Philosophy - PhD
University of British Columbia
With skyrocketing improvements in the performance of modern computing machines, the area of Bayesian inference applications is much improved, and Bayesian statistical analysis is used more frequently. In Bayesian inference, most computational resources are applied to running Markov chain Monte Carlo (MCMC) algorithms to obtain samples from posterior distributions. The MCMC algorithm is the main route to implement Bayesian inference. It allows for high-dimensional and flexible sampling. However, at the same time, researchers can experience poorer computational performance when Bayesian statistical inference is performed using a specific family of models, namely partially identified models. This is because the good computational performance of the MCMC algorithm is not guaranteed. The parameters of the partially identified model are not uniquely identified, which can make off-the-shelf MCMC algorithms very inefficient. Importance sampling with transparent reparameterization (ISTP) is a good computational remedy for posterior inference with partially identified models. With the ISTP algorithm, researchers can obtain better and more stable computational performance while still obtaining samples in the original parameterization. In this thesis, we first traverse scenarios of worsening computational performance with partially identified models and compare the results of ISTP with an off-the-shelf MCMC algorithm (Chapter 2). Then, we discuss the general usability of ISTP and develop the diagnostic method for models suspected to have partial or weak identification (Chapter 3). Along with ISTP, we introduce an R package for the Bayesian inference with the partially identified model (Chapter 4). Lastly, we discuss what was completed, its limitations, and possible future improvements (Chapter 5).
Statistical models that are known or suspected to bepartially identified: Issues of parameterization,computation, and software developmentbyGahyung (Seren) LeeB.Sc, The University of British Columbia, 2020A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Statistics)The University of British Columbia(Vancouver)April 2025© Gahyung (Seren) Lee, 2025The following individuals certify that they have read, and recommend to the Fac-ulty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled:Statistical models that are known or suspected to be partially identified:Issues of parameterization, computation, and software developmentsubmitted by Gahyung (Seren) Lee in partial fulfillment of the requirements forthe degree of Doctor of Philosophy in Statistics.Examining Committee:Paul Gustafson, Professor, Department of Statistics, UBCSupervisorAlexandre Bouchard-CoËte´, Professor, Department of Statistics, UBCSupervisory Committee MemberYongjin Park, Assistant Professor, Department of Statistics, UBCSupervisory Committee MemberDaniel J. McDonald , Associate Professor, Department of Statistics, UBCUniversity ExaminerKyungchul (Kevin) Song , Professor, Vancouver School of Economics, UBCUniversity ExaminerSaman Muthukumarana, Professor, Department of Statistics, University of Mani-tobaExternal ExamineriiAbstractWith skyrocketing improvements in the performance of modern computing ma-chines, the area of Bayesian inference applications is much improved, andBayesian statistical analysis is used more frequently. In Bayesian inference,most computational resources are applied to running Markov chain Monte Carlo(MCMC) algorithms to obtain samples from posterior distributions. The MCMCalgorithm is the main route to implement Bayesian inference. It allows for high-dimensional and flexible sampling. However, at the same time, researchers canexperience poorer computational performance when Bayesian statistical inferenceis performed using a specific family of models, namely partially identified models.This is because the good computational performance of the MCMC algorithm isnot guaranteed. The parameters of the partially identified model are not uniquelyidentified, which can make off-the-shelf MCMC algorithms very inefficient. Im-portance sampling with transparent reparameterization (ISTP) is a good compu-tational remedy for posterior inference with partially identified models. With theISTP algorithm, researchers can obtain better and more stable computational per-formance while still obtaining samples in the original parameterization. In thisthesis, we first traverse scenarios of worsening computational performance withpartially identified models and compare the results of ISTP with an off-the-shelfMCMC algorithm (Chapter 2). Then, we discuss the general usability of ISTPand develop the diagnostic method for models suspected to have partial or weakidentification (Chapter 3). Along with ISTP, we introduce an R package for theBayesian inference with the partially identified model (Chapter 4). Lastly, we dis-cuss what was completed, its limitations, and possible future improvements (Chap-ter 5).iiiLay SummaryMany researchers choose Bayesian inference methods for their statistical analy-sis because of their convenience and flexibility. When researchers use Bayesianinference methods, they will usually benefit from a powerful computational tool,the Markov chain Monte Carlo (MCMC) algorithm. However, MCMC algorithmsmay not work well when used for partially identified models, despite the fact thatthe parameterization makes scientific sense. We suggest a method to alleviate thecomputational challenge of Bayesian inferences with partially identified models,explore its universal usability, and develop a corresponding package in the statisti-cal programming language R [49] to help users perform Bayesian inferences withpartially identified models.ivPrefaceThis dissertation is the original work of Seren Lee completed under the supervisionof Dr. Paul Gustafson. No Generative Artificial Intelligence tools were used in theresearch, thesis writing, or package development.A version of Chapter 2 has been preprinted [Lee S, Gustafson P. Issues of pa-rameterization and computation for posterior inference in partially identified mod-els, 2024. URL https://arxiv.org/abs/2408.10416.]. The main author, Seren Lee,wrote the main text and prepared all simulations, mock data, figures, and tables.Dr. Paul Gustafson, the research supervisor, mentored throughout the research andreviewed the final manuscript.A package pimbayes in chapter 4 is developed by the main author, SerenLee, under the guidance of Dr. Paul Gustafson. The package is published to thepublic on GitHub. Dr. Miguel Barron-Lattes and Xianda Sun provided aid forbridging this package with their packages.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Bayesian inference, partially identified model and the off-the-shelfMCMC algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Importance sampling with transparent reparameterization (ISTP) 62 Issues of parameterization and computation for posterior inferencein partially identified models . . . . . . . . . . . . . . . . . . . . . . 82.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.1.1 Partially identified model . . . . . . . . . . . . . . . . . . 8vi2.1.2 Transparent reparameterization (TP) and importance sam-pling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.1.3 MCMC diagnostic tools . . . . . . . . . . . . . . . . . . 132.2 Logistic regression on multiple variables with missing entries inthe dependent variable . . . . . . . . . . . . . . . . . . . . . . . 162.2.1 Transparent reparameterization . . . . . . . . . . . . . . . 182.2.2 Simulation methods and results . . . . . . . . . . . . . . 212.2.3 Bank data example . . . . . . . . . . . . . . . . . . . . . 302.3 Count data with missing entries . . . . . . . . . . . . . . . . . . . 362.3.1 Pseudo-TP . . . . . . . . . . . . . . . . . . . . . . . . . 372.3.2 Simulation methods and results . . . . . . . . . . . . . . 392.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 472.4.1 Performance of sampling with TP and Pseudo-TP in exam-ples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 472.4.2 Diagnostic tools for the existence of TP . . . . . . . . . . 502.4.3 Improved version of pseudo-TP . . . . . . . . . . . . . . 502.4.4 Weakly vs Partially identified model . . . . . . . . . . . . 513 Numerically diagnosing the identification of partially or weaklyidentified models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.1.1 Importance sampling with transparent reparameterization(ISTP) . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.1.2 Condition number of the matrix . . . . . . . . . . . . . . 553.1.3 Fisher information . . . . . . . . . . . . . . . . . . . . . 563.2 Representation for the identification . . . . . . . . . . . . . . . . 573.3 Exploring the existence of TP . . . . . . . . . . . . . . . . . . . . 773.3.1 A tool for diagnosing the severity of partially identifiedmodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 773.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 813.4.1 Objective threshold . . . . . . . . . . . . . . . . . . . . . 823.4.2 What can we do when warning flag is raised . . . . . . . . 82vii4 pimbayes: A comprehensive tool for Bayesian inference with par-tially identified models . . . . . . . . . . . . . . . . . . . . . . . . . . 844.1 Alternative options for partially identified models . . . . . . . . . 844.2 R package pimbayes . . . . . . . . . . . . . . . . . . . . . . . . 894.2.1 Package scope . . . . . . . . . . . . . . . . . . . . . . . 894.2.2 Technical background . . . . . . . . . . . . . . . . . . . 894.2.3 Installation and setup . . . . . . . . . . . . . . . . . . . . 904.2.4 Applying ISTP to Cross-sectional studies example . . . . 904.2.5 Applying ISTP to Saturated logistic regression model ex-ample . . . . . . . . . . . . . . . . . . . . . . . . . . . . 924.2.6 Manually performing ISTP . . . . . . . . . . . . . . . . . 934.2.7 Diagnosing the identification of the model . . . . . . . . . 944.2.8 Using non-reversible PT . . . . . . . . . . . . . . . . . . 964.3 Package usage examples with real data . . . . . . . . . . . . . . . 984.3.1 Bank data . . . . . . . . . . . . . . . . . . . . . . . . . . 984.3.2 Cross-sectional study . . . . . . . . . . . . . . . . . . . . 995 Conclusions and Discussions . . . . . . . . . . . . . . . . . . . . . . 1015.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1015.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105A Supporting Materials . . . . . . . . . . . . . . . . . . . . . . . . . . 111viiiList of TablesTable 1.1 Two-by-to table from cross-sectional study. . . . . . . . . . . . 3Table 2.1 Simulation time of methods 1 and 2 (sec) as n increases.Method 1 uses JAGS and draws a sample from Ï(Î¸ ,ymis | dn).Method 2 uses ISTP and draws a sample from Ïâ(Ï ,λ | dn). . 23Table 2.2 Simulation time of methods 1 and 2 (sec) as p increases.Method 1 uses JAGS and draws a sample from Ï(Î¸ ,ymis | dn).Method 2 uses ISTP and draws a sample from Ïâ(Ï ,λ | dn). . 26Table 2.3 Simulation time of Method 1 and 2 (sec) as Ï increases. Method1 uses JAGS and draws a sample from Ï(Î¸ ,ymis | dn). Method2 uses ISTP and draws a sample from Ïâ(Ï ,λ | dn). . . . . . . 28Table 2.4 Bank data - Simulation time of Method 1 and 2 (sec) as Ï in-creases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31Table 2.5 Simulation time of methods 1 - 3 (sec) as n increases. . . . . . 41Table 2.6 Simulation time of methods 1 - 3 (sec) as Ï increases. . . . . . 44Table 3.1 Multinomial model with four cells. . . . . . . . . . . . . . . . 58Table 3.2 Two-by-two table with two imperfect tests T1,T2. . . . . . . . . 74Table 4.1 Leptospirosis study for sheep abattoirs in New Zealand [10, 45]. 100ixList of FiguresFigure 1.1 Trace plots of Se in Example 1 . . . . . . . . . . . . . . . . . 5Figure 2.1 Two trace plots. In both plots, there are 3 chains of 15000iterations after 2000 burn-ins, respectively. . . . . . . . . . . 16Figure 2.2 Method 1 (JAGS) Scenario 1 (The number of observations) -Trace plot of β1. . . . . . . . . . . . . . . . . . . . . . . . . 24Figure 2.3 ESS and time in each method and scenario . . . . . . . . . . 25Figure 2.4 Method 1 (JAGS) Scenario 2 (The number of variables in X) -Trace plot of β1 . . . . . . . . . . . . . . . . . . . . . . . . . 27Figure 2.5 Method 1 (JAGS) Scenario 3 (Extent of MNAR) - Trace plotof β1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29Figure 2.6 Bank data - Method 1 (JAGS) Scenario 3 (Extent of MNAR) -Trace plot of β1 . . . . . . . . . . . . . . . . . . . . . . . . . 32Figure 2.7 Bank data - ESS and time in each method . . . . . . . . . . . 33Figure 2.8 Trace plots of δ12 in the simulation results with small and largeprior variances in MAR and MNAR models . . . . . . . . . . 34Figure 2.9 Posterior density plots of δ12 in the simulation results withsmall and large prior variances in MAR and MNAR models . 35Figure 2.10 Method 1 (JAGS, original) Scenario 1 (The number of obser-vations) - Trace plot of µ . . . . . . . . . . . . . . . . . . . . 42Figure 2.11 Method 2 (JAGS, Pseudo-TP) Scenario 1 (The number of ob-servations) - Trace plot of µ . . . . . . . . . . . . . . . . . . . 43Figure 2.12 ESS and time in each method and scenario . . . . . . . . . . 44xFigure 2.13 Method 1 (JAGS, original) Scenario 2 (Extent of MNAR) -Trace plot of µ . . . . . . . . . . . . . . . . . . . . . . . . . . 45Figure 2.14 Method 2 (JAGS, Pseudo-TP) Scenario 2 (Extent of MNAR) -Trace plot of µ . . . . . . . . . . . . . . . . . . . . . . . . . . 46Figure 2.15 Trace plots and f (x) graphs. . . . . . . . . . . . . . . . . . . 52Figure 3.1 Condition number of Fisher information matrix I(η) andasymptotic variance of the estimator of µ , Avar(µË), in modelM2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62Figure 3.2 MCMC trace plots of µ for each breakpoint of α1 in contourplot of M2 with fixed µ = 2,α0 = â1, and prior variance Ï =0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64Figure 3.3 MCMC trace plots of µ for each breakpoint of α1 in contourplot of M2 with fixed µ = 2,α0 =â1, and prior variance Ï = 0.1. 65Figure 3.4 MCMC trace plots of µ for each breakpoint of α1 in contourplot of M2 with fixed µ = 2,α0 =â1, and prior variance Ï = 1. 66Figure 3.5 MCMC trace plots of µ for each breakpoint of α1 in contourplot of M2 with fixed µ = 2,α0 =â1, and prior variance Ï = 10. 67Figure 3.6 The number of solutions of Equation 3.2. . . . . . . . . . . . 68Figure 3.7 The number of solutions of g(θ) = p. The vector p followsDirichlet(1,1,1,1) or Dirichlet(1,1,1,1,1). . . . . . . . . . . 69Figure 3.8 Condition number of Fisher information matrix when increas-ing k of model Mk. . . . . . . . . . . . . . . . . . . . . . . . 71Figure 3.9 Asymptotic variance of the estimator of µ , Avar(µË), when in-creasing k of model Mk. . . . . . . . . . . . . . . . . . . . . 72Figure 3.10 The number of solutions of models M2 and M3 in various θvalues. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure 3.11 The number of solutions of models M2 and M3 in various θvalues with fixed α1 = 0. . . . . . . . . . . . . . . . . . . . . 75Figure 3.12 Graphs of f (µ) under M2 of some θ values generating withα1 = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76xiFigure 4.1 Trace and density plots of µ for JAGS simulation of Example 6with various prior variances. . . . . . . . . . . . . . . . . . . 87Figure 4.2 Trace and density plots of µ for Pigeons simulation of Ex-ample 6 with various prior variances. . . . . . . . . . . . . . 88Figure 4.3 Trace and density plots of sampling with Pigeons package re-sults without extended traces option. . . . . . . . . . . . . . . 97Figure 4.4 Trace and density plots of sampling with Pigeons package re-sults with extended traces option. . . . . . . . . . . . . . . . 98xiiList of SymbolsÏâ(Ï) Convenience prior densityÏâ(Ï ,λ | dn) Convenience posterior densitylâ(dn | Ï) Convenience likelihoodhâ(Î¸) Pseudo-TP functionPl The simplex of length lxiiiList of AbbreviationsCLT central limit theoremESS effective sample sizeISTP importance sampling with transparent reparameterizationMAR missing at randomMCMC Markov chain Monte CarloMNAR missing not at randomPT parallel temperingTP transparent reparameterizationxivAcknowledgmentsFirst of all, I would like to express my genuine gratitude to my research supervi-sor, Professor Paul Gustafson, for his constant support and valuable advice. Hisresearch on Bayesian inference with partially identified models has been a signif-icant source of inspiration for my PhD work. I honestly donât think I could havemade it through this PhD journey without his kind, warm-hearted, and thoughtfulguidance. I am also deeply thankful to the late Professor Ruben Zamar, who wasthe first to see potential in me and encouraged me to pursue this path. Though he isno longer with us, his belief in me remains a meaningful source of motivation andpassion in my research. I would also like to sincerely thank the members of mydissertation committee for generously agreeing to serve on my committee and fortheir valuable time and insightful feedback throughout this process. I truly admiretheir passion for scholarship, which has greatly inspired me. Finally, I am gratefulto Pigeons team and JuliaBUGS team for their collaboration in the development ofthe R package pimbayes.I am profoundly appreciative of my great-grandmother, who lovingly raisedme. Although she passed away in 2020, her strength, wisdom, and warmth remainwith me and continue to guide me every day. I am also deeply grateful to all mygrandparents for their enduring love, life lessons, and gentle strength. Their sto-ries and values have supported and comforted me in more ways than I can express.Above all, I owe my deepest respect and endless gratitude to my parents. Theirunconditional love and unwavering belief in me have been the foundation of ev-erything I have achieved. This PhD is, in many ways, a reflection of their support.I would also like to thank my brother for his quiet strength and steady support.Though rarely spoken, I am truly grateful for his presence in my life.xvTo my dear friends and close acquaintancesâ Arim, Hyunjung, Jeong Ha &Michael, Minji & Heeseop, Yeonjoo & Minsoo, Jiyeon, and Misun, and manyothers whose names I couldnât mentionâthank you for your steady presence andencouragement. Please know that your kindness has not gone unnoticed, even ifI couldnât mention all your names here. Your smiles, empathy, and unwaveringbelief in me have been an anchor through both the brightest and most difficult mo-ments of this journey. I can only hope to be the kind of person you have been to mein your own life paths. Lastly, I would like to thank all the colleagues I met duringmy time in Vancouver since 2016. It has been a true honour to share this chapterof my academic life with such passionate, resilient, and inspiring individuals.With sincere appreciation for all the meaningful connections and support, I seethis achievement not as an endpoint, but as a foundation for the challenges andopportunities that lie ahead.xviChapter 1IntroductionThis thesis is about possible challenges in computational performance whenBayesian inference with an off-the-shelf Markov chain Monte Carlo (MCMC) al-gorithm is used with a partially identified model. We investigate the issues, suggestone potential remedy, and introduce our R package for Bayesian inference withpartially identified models in this thesis.1.1 Bayesian inference, partially identified model and theoff-the-shelf MCMC algorithmBayesian inference is one of the most famous and frequently used statistical meth-ods. Due to the faster computational speed of machines compared to the past,various strategies, techniques, and areas of application of Bayesian inference arerapidly developing. Unlike frequentist inference, Bayesian inference is based onprobabilistic uncertainty of the parameters and works well with complex, high-dimensional problems. After specifying the prior and likelihood and receivingdata, the main work is numerically approximating the posterior distribution, basedon Posterior â Priorà Likelihood. To investigate this distribution, we can havea numerical sample as a representation of it. One method for obtaining this nu-merical representation is having a Monte Carlo sample of draws from the posteriordistribution. Thanks to powerful MCMC sampling methods, we can draw MonteCarlo samples from the posterior distribution even when this posterior is not in1known, closed-form distribution family.However, while Bayesian inference could be used for statistical analysis re-gardless of the type of model the researcher specifies, good computational perfor-mance is not guaranteed for all models. There can be multiple reasons for this; oneis the nonfunctioning of the off-the-shelf MCMC algorithm, such as Gibbs sam-pling [6] with JAGS [46] or Hamiltonian Monte Carlo [43] with STAN [57], whichis used as the basis for Bayesian inference. This failure arises especially when themodel is partially identified. A partially identified model is defined where param-eters are not uniquely identified [22]. The main issue is that the failure of MCMCcould happen even if the current parameterization in the model scientifically makessense.Partial identification is also somewhat related to the term weak identification[58, 59], which is commonly used in econometrics, especially in the context ofthe generalized method of moments (GMM) for parameter estimation. Stock andWright [58] describes the impact of weak identification using the GMM momentcondition. In the identified model (point identified model), only the unique solutionθ = θ0 meets the moment condition that expectation EÏ (θ) is 0. In contrast, if themodel is nonidentified (unidentified), there is no unique solution, while the mostextreme case is that all θ has the moment condition; in other words, the expectationis zero for all θ . Weak identification lies between these two models where thereis still a unique solution θ = θ0, but there are values of θ Ì¸= θ0 for which theexpectation EÏ (θ) is almost zero. With this weakly identified model, we mayencounter issues while using GMM. Partial identification and weak identificationhave the common characteristic that they are between identified and nonidentifiedmodels and are expected to suffer issues in inference due to it.The following is a summary of an example of the partially identified modelintroduced in [45].Example 1. Data from the cross-sectional study collected by an imperfect test.The ideal data we wish we had is described in Table 1.1, where D+ and Dâindicate diseased or not and E+ and Eâ indicate exposure to the risk factor or not.In addition, variables T+,Tâ are for positive and negative results of the diagnostictest. Since the test is imperfect, positive test results do not necessarily mean dis-2eased. This also applies to negative results; there is still a chance of actually beingdiseased.D+ DâE+ x11 x12Eâ x21 x12Table 1.1: Two-by-to table from cross-sectional study.The parameters of interest are test sensitivity, Se, specificity, Sp, and probabil-ities of each cell, Ï11 = Pr(E+ â©D+),Ï12 = Pr(E+ â©Dâ),Ï21 = Pr(Eâ â©D+).The remaining cell probability, Ï22 is omitted because it is a function of the otherthree probabilities, Ï22 = 1âÏ11âÏ12âÏ21. However, the distribution of D is notdirectly learnable from data, which means the distribution of (E | D) is also notlearnable. Therefore, the actual model is specified with the distributions of (D,E)and (T | D). Then, the distribution of (E | D) is found under the assumption ofindependence between the test results and the exposure given the disease status,which is conditional independence of T and E given D. The distributions of (D,E)and (T | D) imply the distribution of (T,D), the multinomial distribution with pa-rameters ηi j instead of Ïi j. The probabilities η are the function of the originalparameters as follows:η11 = Pr(T+â©D+)= SeÏ11+(1âSp)Ï21,Î·12 = Pr(T+â©Dâ)= SeÏ12+(1âSe)Ï22,Î·21 = Pr(Tââ©D+)= SpÏ21+(1âSe)Ï11,Î·22 = Pr(Tââ©Dâ)= (1âSe)Ï12+SpÏ22.Again, the last entry η22 is the complement of the sum of other components. Theoriginal parameterization Se,Sp,Ï11,Ï22,Ï21 scientifically makes sense. However,3while the multinomial data can only describe three things, there are five unknownparameters.When we perform posterior sampling by declaring the model of (D,E) andsampling (T |D) for off-the-shelf MCMC software, on a mock data set we generateaccording to the specification the authors used in [45], the computational result ofMCMC output does not look good in Figure 1.1. These simulations are performedwith JAGS, which uses an off-the-shelf MCMC algorithm - Gibbs sampling [6],for four datasets with different total counts N. As the authors stated in [45], thecomputational performance is expected to worsen with greater N due to the poorerconvergence of the underlying Markov chain [45]. Although three chains of 15000samples after 2000 burn-in samples are used, which results in a total of 45000samples, the estimated effective sample size (ESS) is only at most 1401. ESS,which is used for diagnosing the computational performance of MCMC, is basedon the autocorrelation of the MCMC sample. Since ESS represents the amountof independent information that we could get from the posterior samples, in otherwords, the quality of the numerical representation with ESS = k is equal to that iscontained in k i.i.d. draws, larger and closer to the original sample size is ideal. Butin the present case, it is much less than that. Trace plots, which show the mixingand convergence among chains through iterations, also do not show a positive signin Figure 1.1. Chains do not mix well and have a pattern that shows a high serialautocorrelation across iteration.In addition to this example of a partially identified model, we often see modelsapplied to incomplete data, especially in biological and medical settings. One ofthe most frequently used strategies to deal with missing entries is assuming missingat random (MAR) and imputing missing entries based on that assumption. WhenMAR is assumed, the missingness depends only on the available (observed) data[35]. However, if missing not at random (MNAR) is suspected, where the missing-ness depends on the unavailable (missing) data [35], the model could be partiallyidentified. In these two different types of partially identified models, we may ex-perience issues of computational performance. However, obviously, it may not bea sufficient reason to adjust the parameters in the research. The model and parame-ters are established for scientific reasons and make sense in the context of research.Considering different computing methods is more reasonable than modifying the4Figure 1.1: Trace plots of Se in Example 1. JAGS [46, 47] is used for pos-terior sampling on the data with N = (38,380,3800,38000). The totalcount N is the sum of all entries in the table.parameters and models.Gustafson [22] introduced a transparent reparameterization (TP), one of thegood alternative strategies when partially identified models are used. If a TP isavailable, we can obtain samples from the posteriors in the parameterization weintend while getting a reasonable computational performance. Since this reparam-eterization is reversible, it is possible to seamlessly move back and forth betweenan easily interpretable parameterization (original) and a computation-friendly pa-rameterization (transparent).51.2 Importance sampling with transparentreparameterization (ISTP)Importance sampling with transparent reparameterization (ISTP) is introduced forposterior inferences with a partially identified model. Considering the definition ofthe partially identified model that the parameters cannot be uniquely identified, thealgorithm ISTP uses a different parameterization and sampling from the posterioris expressed in the new parameterization [22].A related idea is also found in [26]. As discussed in their paper, for a Markovchain, geometric ergodicity is a beneficial and preferable property, mainly dueto a Markov chain central limit theorem (CLT) and asymptotic variance calcula-tion. They introduce a variable-transforming method, change-of-variable, to ac-quire new densities with new variables that preserve geometric ergodicity, which isdesirable. After obtaining samples from the new densities, since their transformingmethod is with diffeomorphisms, transforming back to the original variables is fea-sible. As shown in their paper, for examples with conjugate priors with exponentialfamily for categorical data, such as Poisson regression or logistic regression, theposterior is are not super-exponentially light and therefore is not geometric ergodicwith a random-walk Metropolis sampler. After applying the change-of-variableformula, new densities are now super-exponentially light and MCMC chains arenow geometric ergodic.Similarly, when applying ISTP, the computational performance improves asthe new parameterization is separated according to what can and cannot be learneddirectly from the data and isolates the fully identified part of the model. However,ISTP is applicable only when it is possible to find a TP [45]. It is also discussed in[20] that not all models can have a TP, and even if a candidate TP-like parameter-ization pair exists, it may not be easy to prove that mapping between two parame-terizations satisfies the condition for a TP. Consequently, the author introduces analgorithm that uses a Monte Carlo method to find a limiting posterior distributionwhen a TP is unavailable but candidates similar to TP are available.The availability of a TP results in limited access to the algorithm, which mo-tivates us to discuss the universality of the algorithm and expand its usability. Inaddition, ultimately, we make a R package for users to easily access the algorithm6and diagnose and resolve the issues of partial identification of the model.The main motivation for this thesis stems from considering computing meth-ods that could work better than using off-the-shelf MCMC algorithms in Bayesianinference. As already discussed briefly, it is a more reasonable choice comparedto revising the model and parameters to resolve the computational issues. Then,the ultimate contribution of this thesis is to investigate whether and how the ISTPcan be made more accessible to general users and to develop a software package toconstruct a convenient platform.In this thesis, we discuss the computational issues that occurred in a specificparameterization when Bayesian inference is used and the possible remedy - ISTP- is compared to the original methods (Chapter 2). Our experiments are performedwith R [49] and julia [3] with the computing environment of a server clusterwith a 32-core CPU and 64GB RAM provided by the Digital Research Alliance ofCanada. Then, we introduce numerical methods for diagnosing the identificationof a partially identified model with a discussion of identification in the model thatis supposed to be partially or weakly identified (Chapter 3). Finally, we developa package for the suspected partially identified model that can diagnose the givenmodel and apply various methods that could relieve the worsened computationalperformance with the partially identified model (Chapter 4). Our package is de-veloped with julia packages JuliaBUGS.jl [63] and Pigeons.jl [61]. Inthe last chapter, we summarize and discuss these issues, remedies for the issues,and possible future improvements (Chapter 5).7Chapter 2Issues of parameterization andcomputation for posteriorinference in partially identifiedmodelsThis chapter explores the issues of computational performance when an off-the-shelf MCMC algorithm is used for posterior inference with a partially identifiedmodel. We also introduce an alleviation method, and discuss the comparison be-tween the simulation results of those algorithms.2.1 Background2.1.1 Partially identified modelThe term partially identified model refers to the identification of the model. Wesay that a model is well identified when distinct parameter values cannot yield thesame distribution of observable data [21, 22]. It is also called the point identifiedmodel in econometrics Lewbel [33]. Conversely, if this property does not hold, wehave a nonidentified model. This concept can be interpreted in a similar way to asolution of an equation system in mathematics: a unique solution versus infinitely8many solutions. When the equations are insufficient to yield a unique solution,there can be infinitely many solutions. Similarly, if the data are not rich enoughto learn all the parameters, the model can be nonidentified. While from a mathe-matical perspective a model is either identified or not, some nonidentified modelsmight realistically be described as being partially identified. This partially identi-fied model is also described as set identified, as in [33]. Regardless of which termis used, the main characteristic is that it is not a single point but rather a set of val-ues yielding the same distribution given the data. Lack of identification can be aninevitable consequence of what is known about, and what is not known, about thedata-generating process, but this could still lead to some level of useful inference.At least qualitatively, this describes the important feature of a partially identifiedmodel. For example, in the paper [22], the author discovers the spectrum betweentwo end points of the identification, identified model and nonidentified model, withthe dependence of priors in parameterizations.As the author showed in Section 7.1 [22], when the model is partially identified,the performance of MCMC can be poor (e.g., poor mixing). In fact, performancecan further worsen as the model deviates further from the fully identified model.A simple example of the identified model is linear regression on the data of Xand Y .Example 2. Simple linear regression:(Y | X)â¼ N(β0+β1X ,Ï2).When data consists of two variables X and Y with n i.i.d. observations, researchersare typically interested in β1. This model is fully identified. On the other hand, ifwe tweak this example to accommodate missing entries, it can result in a partiallyidentified model. This example is a variation of the example introduced in [23].Example 3. Simple linear regression with missing entries:(Y | X)â¼ N(β0+β1X ,Ï2),(R | X ,Y )â¼ Bern(expit(γ0+ γ1X + γ2Y + γ3XY )),9where R is a binary variable so that R = 1 means Y is observed, and R = 0 indi-cates that Y is missing. This data structure involves the same variables (X ,Y ) asabove, but a certain proportion of the Y entries is missing. Researchers are stillinterested in β1 but are also forced to consider other parameters, notably (γ2,γ3),which describe the relationship between the missingness of Y and Y itself. In de-tail, γi(i = 0,1,2,3) describe the missing data mechanism, with the special case ofγ2 = γ3 = 0 (missing at random (MAR), ignorable data) versus not (missing not atrandom (MNAR), nonignorable data)[35]. When we move from MAR to MNAR,we typically lose identification [35]. Therefore, introducing R and γ makes themodel partially identified, especially since Y influences R, while R explains themissingness of Y .2.1.2 Transparent reparameterization (TP) and importancesamplingTransparent reparameterizationTransparent parameterization [22], which is a mathematical tool to reshape andconvert parameters, results in a transformed set of parameters (Ï ,λ ) derived fromthe original parameter set θ by thinning out, changing, and introducing new pa-rameters. Let h be the invertible function for reparameterizing one set into another.For the function h to be a transparent reparameterizing function, two conditions areto be satisfied. Given the vector of original parameters, θ â Rk, and the result ofthe reparameterization h(θ) = (Ï ,λ ), these conditions are as follows:1. Independence between λ and the data given Ï .Let Ï(dn | Ï ,λ ) be the conditional density function of the data given param-eters. We assume that the distribution of Dn depends only on Ï . Therefore,Ï(dn | Ï ,λ ) = Ï(dn | Ï).2. Applicability of regular parametric asymptotic theory.ân-consistent estimation of Ï that is obtained through Ï(dn | Ï), i.e., anestimator ÏË constructed from Ï(dn | Ï), satisfies ÏË âÏ = Op(1/ân).10When both conditions are satisfied, we say that h(θ) = (Ï ,λ ) is a transparent repa-rameterization.The main role of a TP is to partition the parameters into two parts by identifi-cation when the model is partially identified as in Section 2.1.1. We only keep theparameters that contribute to likelihood while removing unnecessary parts to min-imize the number of elements in the first part, Ï . The second part, λ , can then bechosen arbitrarily as long as it maintains the invertibility of h and matches the totaldimension in the original parameterization. Due to the arbitrariness of λ , there willbe countless TP, if a TP exists. The choice of λ is left to the researchers. Some-times, λ can be chosen to ensure good sampling performance (e.g. i.i.d. MonteCarlo sampling by enabling a conjugacy). Or, λ can be a subset of the originalparameters that researchers can readily interpret.Finally, there is a possibility that TP does not exist. Then, a pseudo-TP can beconsidered. In this chapter, pseudo-TP is a reparameterization with a function hthat is only close to being a TP since dn and λ are weakly dependent given Ï sothat Ï(dn | Ï ,λ ) Ì¸= Ï(dn | Ï), while TP has a necessary condition that Dn â¥â¥ λ | Ï .The detailed explanation and usage of pseudo-TP are investigated in Section 2.3.Convenience priors and convenience posteriorsWhen Bayesian analysis is instantiated, there are prior distributions specified bythe researchers, ideally based on existing knowledge about the scientific contextat hand. However, those distributions are not always convenient for calculationpurposes in the Bayesian analysis. Indeed, there can be prior distributions that aremore calculation-friendly; however, they may differ from those that researchersspecify. As introduced in [22], those convenient but possibly different prior dis-tributions from those that researchers specify are called convenience prior distri-butions. Accordingly, we call the posterior distributions derived from convenienceprior distributions convenience posterior distributions. In this thesis, we place ânext to convenience density functions to differentiate them from the desired den-sity functions.11Importance samplingImportance sampling [18, 31], is a widely applicable technique in computationalstatistics that helps find a given distributionâs properties. The most noticeable fea-ture of the importance sampling is that it can be applied to distributions whichare not easy to draw a sample from. It generally can be explained in one simpleexpression:Wi âf (Xi)g(Xi), (2.1)where Xi is the i-th observation in the sample X1,X2, . . . ,Xn, f (X) is proportionalto the original density of interest, g(X) is proportional to a new, easy-to-sampledensity and Wi is weight of the observation. Rather than trying hard to draw a sam-ple from the original density f (X), we first draw an i.i.d. sample of n observationsfrom the new, easy-to-sample density g(X) and find desired properties of the orig-inal distribution by using this sample and the ratio in Equation 2.1. As weight isproportional to the gap between the shape of the actual density and the new density,weighting allows for the appropriate adjustment to the samples from the densityg(X). Consequently, additional processes other than sampling and weighting, suchas explicit normalization of densities, are not required.Specifically, importance sampling can be used to adjust for the discrepancywhen using convenience distributions rather than their desired counterparts [22].Drawing a sample from the convenience posterior obtained using the TP and thenapplying importance sampling may look like extra steps compared to simply us-ing MCMC directly in the original parameterization. However, since this methoddraws a sample from the convenience posterior density Ïâ(Ï | dn) where Ï is fullyidentified and the convenience prior density Ïâ(Î» | Ï), the computational efficiencyis expected to be better than traditional MCMC methods sampling on the partiallyidentified model that cannot uniquely determine θ .The method introduced in [22], is summarized in Algorithm 1. Given inputof the user (researcher) specified prior density function Ï0(Î¸) and the convenienceprior density function Ïâ(Ï ,λ ) = Ïâ(Î» | Ï) Â·Ïâ(Ï), this algorithm returns a samplebased approximation to the posterior of θ . Note that Ï(Ï ,λ ) is required, but sincethe researchers specify the prior Ï0(Î¸), this density will be determined from theform of the transformation using Ï(Ï ,λ ) = Ï0(hâ1(Ï ,λ ))· | âhâ1(Ï ,λ ) | [52].12Algorithm 1 ISTPInputPrior density function (user-specified) Ï(Ï ,λ )Likelihood function l(dn | Ï ,λ )Dataset dnConvenience prior distribution Ïâ(Ï ,λ ) = Ïâ(Î» | Ï) Â·Ïâ(Ï)OutputSample from desired posterior distribution Ï(Î¸ | dn)1: Calculate the marginal convenience posterior density function of Ï , Ïâ(Ï | dn).2: Find the joint convenience posterior distribution, Ïâ(Ï ,λ | dn) = Ïâ(Ï |dn)Ïâ(Î» | Ï).3: if Ïâ(Ï ,λ | dn) is composed of common distributions available in R then4: Generate m draws composed of i.i.d. Monte Carlo draws using R.5: else6: Generate m draws using an off-the-shelf MCMC package.7: end if8: Let the i-th draw be (Ïi,Î»i).9: Calculate weight Wi =Ï(Ïi,Î»i)Ïâ(Î»i|Ïi)·Ïâ(Ïi) .10: The final weight wi is normalized, wi = Wiâmj=1 Wj .11: Resample the (Ïi,Î»i) draws using weights wi to get a sample from Ï(Ï ,λ | dn)12: Return a final sample mapped back into the parameterization of θ .2.1.3 MCMC diagnostic toolsEffective Sample Size (ESS)ESS [30] refers to the adjusted sample size that reflects the amount of reliableinformation in a Monte Carlo sample. For an autocorrelated sample, the definitionisn1+2ââk=1Ïkwhere n is a sample size and Ïk is autocorrelation at lag k [17]. As seen in thedefinition, ESS is inversely proportional to autocorrelation in an MCMC sample.Theoretically, it can be larger than n, or even close to infinity with negative autocor-relation, but we will not focus on this case in this thesis. An autocorrelated samplewith ESS = n contains the same amount of information about the underlying dis-13tribution as a hypothetical i.i.d sample of size n. ESS may be much smaller thanthe actual sample size, and small ESS implies that MCMC chains are not mixingwell, which is not an ideal performance of MCMC. Finally, since we do not havepopulation-level information, such as autocorrelation, we use the estimate of ESSin practice.For importance sampling, the estimate of common ESS of all parameters isgiven by the formula [9, 18, 22]EÌSS = 1/(mâi=1w2i). (2.2)For the MCMC sampling using STAN [56] or JAGS [47], the estimation of ESS ofeach parameter can be calculated based on autocorrelation of the posterior sample,as implemented by bespoke functions in the respective packages [47, 48, 56]. Ifimportance weighting is applied after an MCMC sampling package is used, thefinal estimated ESS is calculated heuristically by multiplying losses of informationfrom each step:EÌSS =(1/(mâi=1w2i))ÃESSp/m, (2.3)where m is the actual sample size and ESSp is the estimated ESS from the package.Multivariate ESSWhen MCMC is used for posterior sampling, the ESS covered in the last sectionpertains to a single parameter. We can indeed report the ESS for the parameter ofmost interest, to diagnose MCMC performance. However, in models with multipleparameters, especially high-dimensional, this single ESS may not be representativeof overall MCMC performance. ESS estimates of some parameters can be indeedsmall, while others are closer to the actual sample size. Since the posterior distribu-tion is inherently a joint distribution of all parameters, ideally, multiple parametersshould be considered. As this issue arises from the usage of a single ESS, it can bealleviated by using an estimate of ESS for multivariate parameters. Let the matrixΣ be positive definite and defined under Markov chain central limit theorem (CLT)14[28],ân( fËnâE( f (Xi))) dâ (0,Σ)where X1,X2, . . . ,Xn is Markov chain with invari-ant probability distribution equal to the target distribution and fËn = 1n âni=1 f (Xi).Since we would like to have the Monte Carlo standard error of all parameters, f (X)is taken to be the identity function. The multivariate ESS estimator, introduced in[65], is based onESS = m · |Î|1p|Î£| 1p,where m is the actual sample size, p is the number of parameters, Î is the true co-variance matrix, and Σ is a Monte Carlo standard error. Then, a strongly consistentestimator of ESS is,EÌSS = m · |În|1p|Î£n|1p,where În is the sample covariance matrix and Σn is strongly consistent estimateof Σ. All simulations in this thesis use the R function multiESS() in mcmcsepackage [12].Trace plotA trace plot is a graph showing the sampleâs evolution versus the state (sample)number. The trace plot can show the convergence and autocorrelation. If there ispoor mixing, as in the left plot of Figure 2.1, samples are less precise representa-tions of the target distribution since they are highly autocorrelated. Also, the threechains do not agree with each other. On the contrary, the right plot in Figure 2.1does not have any prominent shape; the three chains are overlapped, and mixingwell. Even though there is no standard rule for determining the performance ofMCMC through trace plots and there is no way to consider multiple parameterssimultaneously in a single plot, it is worth using them since they are very straight-forward to interpret and compare.15Figure 2.1: Two trace plots. In both plots, there are 3 chains of 15000 itera-tions after 2000 burn-ins, respectively.2.2 Logistic regression on multiple variables withmissing entries in the dependent variableIn this section, an example of a partially identified model is introduced, a corre-sponding TP is derived, and simulation results with traditional MCMC and theISTP algorithm in multiple scenarios involving poor MCMC performance aresummarized.Example 4 is an extension and variation of the mock data example introducedin unpublished lecture notes [23]. The data are comprised of n observations ofX , a vector of independent variables with p binary entries X1,X2, . . . ,Xp, and asingle dependent binary variable Y . Researchers are interested in the associationbetween X and Y . The challenging aspect of these data is that some entries of Yare unavailable. At this point, it is not yet clear whether a model of Y versus Xresults in a partially identified model or an identified model. If Y is suspected tobe missing not at random (MNAR) i.e, the missingness of Y is influenced by bothY and X , then consequently the model is partially identified. The intuition behindthis is the same as in Example 3 - Missingness of Y is again explained by Y . After16introducing a new binary variable R, which is 1 when Y is observed and 0 when Yis missing, the model researchers have designed can be summarized as below:Example 4. Saturated logistic regression with incomplete outcome data:â¢ Dataâ There are p variables in matrix X = (X1, . . . ,Xp)Ta design matrix of X is introduced to describe a saturated model. Let avector Xdes â R2p be a generic row of the design matrix, then Xdes =(X0...0, . . . ,X1,...1) where Xi1,i2,...,ip = 1 â X1 = i1,X2 = i2, . . . ,Xp = ip.The row vector Xdes has 2p element where each element represents acombination of X1 to Xp. As these elements cover all the combinationsof X, exactly one element is 1, and the others are zero.â Y is the dependent variable, which is binary.â R is also a binary variable that describes the missingness of Y . Thevariable R = 1 if and only if Y is observed, and R = 0 if and only if Yis not available.â¢ Parametersâ α = (p0...0, . . . , p1...1)T where α belongs to the probability simplex on2p categories and pi1,i2,...,ip = Pr(X1 = i1,X2 = i2, . . . ,Xp = ip)â β â R2pâ γ â R2pâ δ â R2p .â¢ Model - Multinomial and Logistic regressionâ Xdes â¼Multinom(Î±,1)A multinomial distribution is suitable as the sum of each row is fixed at1, and the sum of the probability of each column being one is 1 whilethose columns are independent of each other.â (Y | Xdes)â¼ Bern(expit(Xdesβ ))17â (R | Xdes,Y = 0)â¼ Bern(expit(XdesÎ³))â (R | Xdes,Y = 1)â¼ Bern(expit(Xdes(expit(logit(γ+δ )))))Since both Y and R are binary, Y given X and R given (X ,Y ) are ex-plained with the logit scale of saturated linear regression in this model.In summary, the unknown parameters are α,β ,γ,δ .Then, priors are chosen forthese parameters:α â¼ Dir(12p),βi â¼Uni f (0,1) (i = 0, . . . ,2pâ1),Î³i â¼Uni f (0,1) (i = 0, . . . ,2pâ1),Î´i â¼ Normal(0,Ï2) (i = 0, . . . ,2pâ1).(2.4)While noninformative priors are chosen for α,β ,γ , the prior for δ is informativeand has a hyperparameter. Since elements of δ are log odds ratio of (R,Y ) givenX , the magnitude of δ describes the departure of missingness from MAR and,therefore, the distance from a well-identified model. Researchers set Ï2 based ontheir belief about how much deviation from MAR could exist. Therefore, in thisthesis, the prior variance of δi, Ï2, is the hyperparameter varied in the simulationto test the impact of MNAR.2.2.1 Transparent reparameterizationIf we directly apply a standard MCMC method, especially Gibbs sampling, withparameters α,β ,γ,δ , the computational performance is poor. While partial identi-fication is expected to affect performance, the dependence between elements in αleads to the issue specifically of Gibbs sampling. Since Gibbs sampling uses theconditional distribution of one parameter given other parameters, high dependencereduces computational efficiency. The elements of the simplex α are not indepen-dent, in a distributional sense, of each other. This issue is resolved by removingone element from the parameterization (e.g., the first element is the complementof the sum of the rest elements). In this thesis, whenever quantity x is constrainedto have elements summing to one, we let x(0) be the subset of x without the firstelement. Then, the parameter set, θ , consists of (α(0),β ,γ,δ ). The dimension18of θ is dim(θ) = 2pâ 1+ 3Ã 2p = 2p+2â 1. This dimension must be preservedafter applying TP. The issue with these parameters as given by θ is that some ofthem cannot be directly learned from the data. Considering Y is not available whenR = 0, while X is always available, we can directly learn the following parametersfrom the data,ε = Pr(R = 1),ζ (0) = Pr(X(0)des | R = 0),η(0) = Pr(X(0)des | R = 1),ξ = Pr(Y = 1 | Xdes,R = 1).(2.5)These parameters are taken to comprise Ï . The number of parameters in ξ isthe number of entries in Xdes, 2p, while ζ (0) and η(0) have one less. Therefore, thedimension of Ï is dim(Ï) = 1+2p+2(2pâ1) = 3Ã2pâ1. Now, we are choosingλ that makes h(θ) = (Ï ,λ ) invertible and preserves the dimension. The vector oflog odds ratios δ is one option, while infinitely many λ options are available aslong as they meet the condition. In the case of using λ = δ , the function h(θ) isgiven by the system of equations below:19ε = Pr(R = 1)= âxâXdesPr(R = 1 | x) Â·Pr(x)= âxâXdes{Pr(R = 1 | x,Y = 1) Â·Pr(Y = 1 | x) Â·Pr(x)}+ âxâXdes{(Pr(R = 1 | x,Y = 0) Â·Pr(Y = 0 | x) Â·Pr(x)}=< Î±,Î³â,Î² >+< α,γ, β¯ >,ζ = Pr(Xdes | R = 0)= Pr(Xdes,R = 0)/Pr(R = 0)= {Pr(Xdes,Y = 1,R = 0)+Pr(Xdes,Y = 0,R = 0)}/Pr(R = 0)= {β ·α · (1â Î³â)+(1âÎ² ) ·α · (1â Î³)}/(1â ε),η = Pr(Xdes | R = 1)= Pr(Xdes,R = 1)/Pr(R = 1)= {Pr(Xdes,Y = 1,R = 1)+Pr(Xdes,Y = 0,R = 1)}/Pr(R = 1)= {β ·α · Î³â+(1âÎ² ) ·α · γ}/ε,ξ = Pr(Y | Xdes,R = 1)= Pr(Y,Xdes,R = 1)/Pr(Xdes,R = 1)= (β ·α · Î³â)/(Î· · ε),(2.6)where < x,y,z > is the sum of an element-wise product of vectors, < x,y,z >=âdim(x)i=1 xi · yi · zi.20Then, the function hâ1(Ï ,λ ) is given by the system of equations below:α = Pr(Xdes) = Pr(R = 1) Â·P(Xdes | R = 1)+Pr(R = 0) Â·P(Xdes | R = 0)= ε Â·Î·+(1â ε) ·ζ ,β = Pr(Y | Xdes) = Pr(Y,Xdes)/Pr(Xdes)= (Pr(Y,Xdes,R = 0)+Pr(Y,Xdes,R = 1))/α= {expit(logit(ξ )âÎ´ ) · (1â ε) ·ζ +ξ · ε ·η}/α,γ = Pr(R = 1 | Xdes,Y = 0) = Pr(R = 1,Xdes,Y = 0)/Pr(Xdes,Y = 0)= (1âÎ¾ ) ·η · Îµ/{(1âÎ² ) ·α}.(2.7)2.2.2 Simulation methods and resultsTo compare the performance of traditional MCMC and ISTP, we set up three sce-narios that may differently burden the computation. While arguments of sampling,n, p,Ï , are varied respectively in each scenario, all other non-varied arguments arefixed to common values of n = 3000, p = 3 and Ï = 0.5. The mock data is gener-ated with the model specified in Example 4 using parameter values produced withrandom generator functions in R. To reproduce these parameter values, a seed valuefor the random number generator is fixed throughout the simulation.1. The number of observations, n. (10 levels, 3 - 10000)Having enough information is essential in estimating; however, MCMC sam-pling can take too long if the data is too large. Therefore, mock data with atiny number of observations (n = 3) to many observations (n = 10000) aretested, while other conditions remain the same.2. The number of variables in X , p. (5 levels, 1 - 5)Estimating is generally harder if we have many parameters to estimate. Sincethe model is saturated, the number of total parameters increases exponen-tially in p.3. The deviation a prior allowed from the MAR assumption, as represented bythe prior variance of δi, Ï2. (7 levels, 0.01 - 10)21This is related to how much the model is partially identified. With the widerprior, the missingness of Y is permitted to deviate further from MAR. So, weexpect that it will make the performance worse. In the simulation, a narrowprior (Ï = 0.01) to a wide prior (Ï = 10) are used.Two methods are tested with these three scenarios for generating a posteriorsample of parameters θ = (α,β ,γ,δ ). For both methods, 45000 draws are foundfor each parameter while the first method, JAGS, achieves this using three chainsof 15000 draws each after 2000 burn-ins.First, traditional MCMC is implemented with JAGS called from R [47] withmodel and prior specification declared in terms of the original parameterizationin Example 4 and Equation 2.4. Three chains are used to verify that chains aremixing well. In this method, JAGS draws a sample from the posterior distributionof parameters and missing data Y given observed data Y and X .In contrast, the second method uses ISTP and draws an i.i.d. Monte Carlosample of parameters in (Ï ,λ ) given observed data by using R base functions suchas rnorm() and rbeta(). Since Ï is newly introduced, convenience priorsare required to be set. For these parameters and data structure, conjugate priors areavailable. For ε , prior beta(1,1)will result in posterior of beta(1+nR=1,1+nR=0).For other parameters, a Dirichlet prior Dir((1, . . . ,1)) will result in posterior ofDir((1+n1, . . . ,1+nk)) where k is the number of parameters and ni is the numberof observations that belongs to condition (category) i. When finding importanceweights for resampling, the ratio between two prior densities, true and conve-nience, is evaluated. Since a sample is drawn from the convenience posterior, achange of variables [52] is required. Consequently, the determinant of the Jaco-bian of the inverse function hâ1(Ï ,λ ) is found. In all simulations in this thesis,the R function jacobian() in R package pracma [4] is used. This R pack-age and function enable the R user to access the MATLAB function jacobian()[25] for finding Jacobian matrix by using provided inputs of symbolic R functionsthat take parameters and density functions as functional arguments. Given theseinputs, jacobian() function evaluates the Jacobian matrix at given points in theparameter space, which is used to find the determinant with another R functiondet().22ResultsWhile we also consider multiESS and median ESS across all elements of θ to re-flect the performance of MCMC, we focus on the ESS and trace plot of a presumedtarget parameter β1 = logitPr(Y = 1 | X1 = 1,X2 = 0, . . . ,Xp = 0).Scenario 1 (n = 3â10000 (10 levels), p = 3,Ï = 0.5)In Figure 2.2 that shows the performance of Method 1, two things are notable- when n is larger, the trace plot of β1 has more apparent autocorrelation and theESS of β1 decreases. Considering that the posterior sample size is N = 45000, anESS of 770 is small enough to indicate that MCMC is not working well. Also,in the top left plot of Figure 2.3, the ESS/time of Method 1 is good for a smallsample size regardless of how ESS is found. However, it steeply decreases rightafter that. The main reason for this is the sampling time. ESS itself is not as steeplydecreasing as ESS/time. We also see that the sampling time increases exponentiallyfrom 4 seconds to 2 hours in Table 2.1. On the other hand, Method 2 shows smallESS/time for the small sample size but has generally stable ESS/time and is betterthan Method 1 when the sample size is large. And as in Table 2.1, time doesnâtchange much from one minute, and it is the same for ESS of 25000 after n= 1000.And most of all, the time is much smaller than that of Method 2.In conclusion, the performance of Method 1 worsens with larger n in terms ofboth ESS and time. In contrast, the performance of Method 2 is worse than Method1 when the data size is small, but it gets better when there is enough data, and theperformance does not change much in both time and ESS after that. Therefore, itshows better performance.n 3 10 30 100 300Method1 - JAGS 4 6 14 37 96Method2 - ISTP 68 67 67 67 67n 1000 2000 3000 5000 10000Method1 - JAGS 323 844 1483 2927 7292Method2 - ISTP 67 68 67 67 67Table 2.1: Simulation time of methods 1 and 2 (sec) as n increases. Method1 uses JAGS and draws a sample from Ï(Î¸ ,ymis | dn). Method 2 usesISTP and draws a sample from Ïâ(Ï ,λ | dn).23Figure 2.2: Method 1 (JAGS) Scenario 1 (The number of observations) -Trace plot of β1. The red curve denotes the true value of β1. Among 10values of n, n = 3,100,2000, and 10000 are plotted.24Figure 2.3: ESS and time in each method and scenario.The left plots are ES-S/time, and the right plots are ESS in each scenario and method. Method1 has three curves each, the red curve is with a median of ESS of all pa-rameters, the green curve is multiESS [12, 65] in Section 2.1.3, and theorange curve is ESS of β1 = logitPr(Y = 1 |X1 = 1,X2 = 0, . . . ,Xp = 0).One blue curve represents the result of Method 2.25Scenario 2 (p = 1â5 (5 levels),n = 3000,Ï = 0.5)In Figure 2.4, the trend of ESS and trace plot differs from the other two scenar-ios. ESS decreases a bit and then increases again. Similarly, the trace plot doesnâtshow much autocorrelation from the smaller to larger p. However, in Table 2.2,even p = 1 takes 7.5 minutes to complete, increasing to 100 minutes when p = 5.Therefore, even though the ESS and trace plot is reasonable, the ESS/time is al-most 0 in Figure 2.3. One noticeable thing in Figure 2.3: only multiESS shows asincreasing from p = 4 to p = 5. However, this caseâs multiESS value is unreliableas the covariance matrix used to calculate it is close to singular.On the other hand, the Method 2 performance worsens when p increases, sinceESS decreases from 30000 to a small number near 0 and time increases from 26seconds to 8 minutes. However, ESS/time is still comparable to Method 1, consid-ering the poor performance of Method 1 even in the smallest case of p = 1.While scenario 1 and scenario 2 scale up the size of the problem by increasingthe data size or parameter size, Method 1 and Method 2 show different trends inthese two scenarios. Method 1 loses performance in terms of both ESS and timein scenario 1; however, only time is the issue in scenario 2. Method 2 doesnâtlose performance in scenario 1 but worsens in scenario 2 both in ESS and time.However, this degradation is less than that of Method 1.p 1 2 3 4 5Method 1 - JAGS 459 784 1489 2901 5939Method 2 - ISTP 26 49 100 216 498Table 2.2: Simulation time of methods 1 and 2 (sec) as p increases. Method1 uses JAGS and draws a sample from Ï(Î¸ ,ymis | dn). Method 2 usesISTP and draws a sample from Ïâ(Ï ,λ | dn).Scenario 3 (Ï = 0.01â10 (7 levels),n = 3000, p = 3)Scenario 3 is about the extent of MNAR allowed for by prior variances, whichindicates how far the researcher feels the missing data mechanism could devi-ate from MAR. As the prior standard deviation, Ï , increases both Method 1 andMethod 2 show worsening performance. In Table 2.3, time doesnât change much(around 23 minutes for Method 1 and 1.7 minutes for Method 2). Therefore, theworsening of performance is mainly in ESS. Both in Figure 2.5 and Figure 2.3,26Figure 2.4: Method 1 (JAGS) Scenario 2 (The number of variables in X) -Trace plot of β1. The red curve denotes the true value of β1. Among 5values of p, p = 2,3,4 and 5 are plotted.27ESS decreases as prior SD increases. When comparing the results of Methods 1and 2, while the two methods show a trend of decreasing performance in larger Ï ,Method 2 has much better performance because it has both a shorter time of sam-pling and larger ESS. The notable thing about scenario 3 is that unlike scenarios 1and 2, Method 1 shows a complete breakdown at the extreme case Ï = 10. WhileMethod 2 also shows a steep decrease, it doesnât thoroughly crash as does Method1.Ï 0.01 0.05 0.1 0.5 1 5 10Method 1 1346 1359 1310 1325 1339 1389 1423Method 2 102 101 105 102 102 102 103Table 2.3: Simulation time of Method 1 and 2 (sec) as Ï increases. Method1 uses JAGS and draws a sample from Ï(Î¸ ,ymis | dn). Method 2 usesISTP and draws a sample from Ïâ(Ï ,λ | dn).28Figure 2.5: Method 1 (JAGS) Scenario 3 (Extent of MNAR) - Trace plot ofβ1. The red curve denotes the true value of β1. Among 7 values of thehyperparameter Ï , Ï = 0.01,0.1,1 and 10 are plotted.292.2.3 Bank data exampleBank data [40, 41] were collected from direct marketing campaigns conducted bya bank in Portugal. The campaigns were completed by telephone calls. Thereare sixteen predictors and one response variable Y that indicates whether a clientsubscribed to a term deposit.Data editingWe first choose four predictors, marital, default, housing and loan.While default, housing, loan are binary, marital is categorical withmore than two categories. Therefore, we transform it into married where it hastwo categories: married (TRUE, 1) and not -divorced, single, unknown - (FALSE,0). Then, using these four binary variables, we create the design matrix with 24 =16 columns. This design matrix is created with the rule explained in Example 4,for instance, X1111 = 1 only when all four variables are 1 (TRUE). Although thedata do not have missing entries, we revise the dataset so that it has missing entriesin the response variable. The following is a summary of the data revision andsimulation method for the first revised dataset (MNAR):â¢ Create the design matrix has 2p = 24 = 16 columns, and n = 4521 rows.â¢ Randomly generate values for the parameters γ,δ . We use priors in Equa-tion 2.4.â¢ Calculate the probabilities of Ri being 1 for each i= 1, . . . ,n. Using probabil-ities, find the values of Ri. If Ri = 0, make Yi unknown. With this process, therevised dataset is far from MAR. After applying, we find that the proportionof missingness in the dataset is about 30%.â¢ Apply methods 1 (JAGS) and 2 (ISTP).â¢ Use simulation methods similar to scenario 3. The values for n, p, n =4000, p = 4, are fixed, and Ï is varied within (0.01,0.05,0.1,0.5,1,10).In practice, we usually do not know whether we have MAR or MNAR. If we allowless deviation from MAR by using small prior variances of δ but actually havingMNAR, the point estimates are expected to be inaccurate. However, if we allowfor more deviation from MAR, the computational performance is expected to be30lower, as in Figure 2.5. In this section, we focus mainly on the computationalperformance. However, the precision and reliability of the results are also concerns.Therefore, we choose multiple Ï values for simulation to check this.ResultsIn Figure 2.6, it is clear that ESS decreases from 11645 to 299 as the prior standarddeviation Ï increases from 0.01 to 10. And the mixing in the trace plot gets worse.It is also possible to see that there is a deviation between chains in the last plot,which is for Ï = 10. Since n and p are both larger than the simulation of Scenario3, the simulation time for each method also takes longer in Table 2.4. Method 1takes about 80 minutes, and Method 2 takes about 3 minutes. It is positive that thetime decreases greatly when using Method 2, but ESS is not that good. We expectthat this is due to the value of p. We see that when p increases, the ESS of Method2 decreases substantially and is sometimes smaller than that of Method 1.Ï 0.01 0.05 0.1 0.5 1 5 10Method 1 4681 4752 4661 4737 4686 4723 4770Method 2 193 195 195 195 193 195 197Table 2.4: Bank data - Simulation time of Method 1 and 2 (sec) as Ï in-creases. Method 1 uses JAGS and draws a sample from Ï(Î¸ ,ymis | dn).Method 2 uses ISTP and draws a sample from Ïâ(Ï ,λ | dn).On the other hand, as we briefly discuss above, whether point/interval esti-mates truly represent their target parameters is also a concern for researchers. Inthe real world, we usually donât know the true value, but for our simulation, as werevise the data set, we could adjust γ and δ values. For any data created based onignorable missingness, missing completely at random (MCAR), the γ values arethe same as Pr(R = 1 | X) = Pr(R = 1), and the δ values are zero. In our secondrevised dataset, we assign Pr(R = 1) = 0.7,γ1 = γ2 = · · Â·= γ16 = 0.7. So, we havetwo datasets for the simulation: the first revised dataset (MNAR) and the second(MCAR) revised dataset. In Figure 2.8, when a small prior variance 0.01 is usedfor MAR data, the trace plot shows that the true value of δ12 is well caught in theinterval. However, this is not true for MNAR data. The computational performanceis good in both simulations, but how well the interval estimates represent the truevalues is different. When we increase the Ï value from 0.01 to 10, the compu-31Figure 2.6: Bank data - Method 1 (JAGS) Scenario 3 (Extent of MNAR) -Trace plot of β1. Among 7 hyperparameters of Ï , Ï = 0.01,0.1,1 and10 are plotted.32Figure 2.7: Bank data - ESS and time in each method. The left plots areESS/time, and the right plots are ESS in each scenario and method.Method 1 has three curves each, the red curve is with a median of ESS ofall parameters, the green curve is multiESS [12, 65] in Section 2.1.3, andthe orange curve is ESS of β1 = logitPr(Y = 1 |X1 = 1,X2 = 0, . . . ,Xp =0). One blue curve represents the result of Method 2.tational performance is much worse, but the true values are in the interval, as theposterior variances are large. We also draw density plots in Figure 2.9. The mostnotable thing is that, in results for small prior variances, prior densities and pos-terior densities are almost overlapped. When large prior variances are used, truevalues are not near the mode of priors and posteriors but in a reasonable windowfrom the mode. At the same time, all of the densities have unimodal shapes, whichmay not reflect the behaviour of partially identified models well since lack of fullidentification can lead to multimodal posterior distributions.33Figure 2.8: Trace plots of δ12 in the simulation results with small and largeprior variances in first and second revised datasets. JAGS simulation isused with the specified prior SD Ï in the plot. The red horizontal line istrue parameter value of δ12.34Figure 2.9: Density plots of δ12 in the simulation results with small and largeprior variances in MAR and MNAR data. MAR data is created withPr(R = 1) = 0.7, and MNAR data is the same data as in simulations ofthe current section. JAGS simulation is used with the specified prior SDÏ in the plot. Prior distribution is N(0,Ï2) with specified prior SD Ï .The sky blue (dotted) curve is for the prior density, and the pink curve isfor the posterior density. The black vertical line is true parameter valueof δ12.352.3 Count data with missing entriesIn the previous section, logistic regression applied to data with missing entries wasexplored. Since all variables were binary, finding a transparent reparameterizationwas not complicated. This process can be more challenging when the variables arenot binary. In this section, we deal with count outcome data, which is discrete butnot as simple as binary.Example 5. Count regression with incomplete outcome data:X â¼Multinom(Î±,1),Y | Xdes â¼ Pois(exp(Xdesβ )),R | Xdes,Y = yâ¼ Bern(expit(Xdes(γ+δy))).Also, following the same definition as in other examples, R is a binary indicatorvariable that is 1 whenever Y is available and is 0 when Y is missing. Finally, Xhas p binary elements, and a generic row of design matrix Xdes has 2p elements.The difference between Example 4 and Example 5 is whether Y is binary orcount, and therefore, we start from the assumption of the model of Poisson regres-sion for (Y | X). However, this makes a huge difference in finding an appropriateTP. There is no guarantee that Y given R = 1 belongs to the same family of dis-tributions (Poisson), as it does marginally. When Y is binary, both (Y | R = 1)and Y follow Bernoulli distribution. However, it is not for count variable Y , whichleads to non-closed form probability mass function for Y | X = 1. So, we introducepseudo-TP. To test the computational performance with a pseudo-TP, we focus onthe simpler example below:Example 6. Estimating parameters in incomplete count data:Y â¼ Pois(µ),R | Y â¼ Bern(expit(α0+α1Y )).(2.8)In fact, this example is a special case of Example 5, the model with p = 0. Itis a minimal working example that illustrates the fundamental computational chal-lenge, Y | R = 1 and Y do not belong to the same distribution. A mock dataset36is generated with fixed parameters of µ = 5, (α0,α1) = (log(2)/â5)(â5,1) us-ing R. The parameter values α0 and α1 are chosen based on the values E(Y ) =Var(Y ) = 5. In particular, the contrast with Y being one SD above/below its mean,logit(Pr(R = 1 | y+1))â logit(Pr(R = 1 | y)) = log2, induces a log odds ratio of2 for the chance of missingness. At the same time, this choice of α results in theaverage missingness rate of approximately expit(0) = 0.5.Bayesian analysis can be performed with priors as follows:µ â¼ Gamma(1,1),α0 â¼ Norm(0,10),α1 â¼ Norm(0,Ï2).This Example 6 has no X , so it is much simpler now. Nontheless, it is not clearwhether the model is fully identified or not. However, it is clear that this examplestill has nonignorable missingness [35] in Y , which may have a crucial impact oncomputational efficiency. At this point, a pseudo-TP might help, since data dependsstrongly on Ï and weakly (or almost not at all) on λ .2.3.1 Pseudo-TPThe structure of pseudo-TP is similar to the TP of Example 4, but excludes anyrelationships between Y and X . The initial parameterization is θ = (µ,α0,α1). Inthis example, the probability γ = Pr(R = 1) and the conditional expectation δ =E(Y | R= 1) can be directly learned from the data. These two parameters are takento comprise Ï = (γ,δ ). Intuitively, the data do not directly inform about α1, whichdescribes the extent of MNAR. So, we can use it as λ . As discussed in Example 5,the closed-form TP may not exist. Even if we could define h(θ) = (Ï ,λ ) whereÏ = (γ,δ ) and λ = α1, since the expression in h(θ) involves the infinite sum, wedonât have a closed-form expression for hâ1(Ï ,λ ).But with Taylor approximations for finding some necessary parameters, wemay find a reparameterization close to TP. Now we find hâ(Î¸) and hâ1â (Ï ,λ ) in-stead, as functions which define a pseudo-TP. The pseudo-TP function hâ(Î¸) is37based on the following approximations:γ = Pr(R = 1)=âây=0Pr(Y = y)Pr(R = 1 | Y = y)= E(Pr(R = 1 | Y ))= E(expit(α0+α1Y ))â expit(α0+α1E(Y ))= expit(α0+α1µ):= Î³â,Î´ = E(Y | R = 1) = âây=0 y Â·Pr(Y = y,R = 1)Pr(R = 1)=âây=0 y Â·Pr(Y = y,R = 1)âây=0 Pr(Y = y) Â·Pr(R = 1 | Y = y)=âây=0 y Â·Pr(R = 1 | Y = y) Â·Pr(Y = y)âây=0 Pr(Y = y) Â·Pr(R = 1 | Y = y)=âây=0 y · expit(α0+α1y) · eâµµyy!âây=0 expit(α0+α1y) · eâµµyy!=âây=0 y · exp(g(y)) · exp(âµ)µyy!âây=0 exp(g(y)) · exp(âµ)µyy!â âây=0 y · exp(g(µ)+gâ²(µ)(yâµ)) · exp(âµ)µyy!âây=0 exp(g(µ)+gâ²(µ)(yâµ)) · exp(âµ)µyy!=âây=0 y · exp(ygâ²(µ)) · µyy!âây=0 exp(ygâ²(µ)) · µyy!=âây=0 y · (µ exp(gâ²(µ)))yy!âây=0(µ exp(gâ² (µ)))yy!=âây=0 y · exp(âµ exp(gâ²(µ)))(µ exp(gâ² (µ)))yy!âây=0exp(âµ exp(gâ² (µ)))(µ exp(gâ² (µ)))yy!38=âây=0 y · exp(âµ exp(gâ²(µ)))(µ exp(gâ² (µ)))yy!1:= δ â,where g(y) = log(expit(α0+α1y)). By the approximation of γ , which is Î³â, above,δ â â µ(exp(α1 · (1â Î³â))). While we learn (γ,δ ) from the data, approximationsÎ³â and δ â are used for finding a closed-form inverse function.Since Taylor approximation is used multiple times, hâ(Î¸) is not a TP, but weexpect that it is close to a TP since hâ(Î¸) is a legitimate approximation that includespartitioning between parameters that intuitively are informed directly by the data,versus those that are not. Note that, the inverse function hâ1â (Ï ,λ ) is given by thesystem of equations below:µ = µâ/(exp(Î±1 · (1â Î³â))),Î±0 = logit(Î³â)âÎ±1 Â·µ.2.3.2 Simulation methods and resultsIn this section, for testing the performance of sampling with pseudo-TP, we ap-ply simulations to two scenarios introduced in Section 2.2.2, one with increasingthe number of observations and the other with further deviation from MAR. Thescenario with increasing γ is excluded since this example is inherently univariate.Scenario 1 uses data with n = 3 to n = 10000 while fixing Ï = 0.5. And in sce-nario 2, a narrow prior (Ï = 0.01) to a wide prior (Ï = 10) are used, where Ï2 isthe prior variance of α1. For all simulations in scenario 2, n = 3000.While Method 1 remains the same as Section 2.2.2 with MCMC using JAGSapplied in the original parameterization, two methods are used for testing inferenceobtained under the pseudo-TP parameterization.â¢ Method 2This method uses JAGS. But sampling is on (Ï ,λ ,Ymiss) via conveniencepriors as follows:γ â¼ Beta(1,1),39δ â¼ Gamma(1,1),α1 â¼ Norm(0,Ï2).Importance weights are required for this method since it uses conveniencepriors. The main difference between methods 1 versus 2 is that JAGS iterateson θ (Method 1) and (Ï ,λ ) with mapping back to θ (Method 2). Eventhough Method 2 could not reparameterize θ , this difference might alleviatethe computation efficiency issue arising from partial or weak identification.â¢ Method 3This method uses importance sampling with pseudo-TP. The crucial pointof this method is facilitating i.i.d. sampling from a convenience posterior,which is different from a convenience posterior for MCMC sampling inMethod 2. This is achieved by using a convenience likelihood as well asa convenience prior. The convenience likelihood is constructed under the as-sumption that (Y | R = 1) follows the Poisson distribution, while Method 2is free of the distributional assumption. Consequently, the parts of likelihoodin weight ratio canât be cancelled since the likelihood with the assumptionis no longer equivalent to the original likelihood. Hence, the weight ratio isfound with Wi =Ï(Ïi,Î»i)l(dn|Ïi,Î»i)Ïâ(Î»i|Ïi)·Ïâ(Ïi)·lâ(dn|Ïi) where lâ(dn | Ïi) is convenience likeli-hood and l(dn | Ïi,Î»i) is true likelihood. The convenience likelihood is basedon the following statistical modelRâ¼ Bern(γ),(Y | R = 1)â¼ Pois(δ ).(2.9)In the usual calculation for importance weights, unobserved data is not nec-essary. However, since the true likelihood l(dn | θ) in Equation 2.8 requiresall entries of y, imputing of Y is necessary. JAGS in Methods 1 and 2 im-putes missing entries by default, but in Method 3, we first find the conditionalexpectation E(Y | R = 0) â µ · exp(âÎ±1 · expit(α0 +α1µ)) and impute un-observed Y on the assumption that Y | R = 0 follows Poisson distribution.40ResultsFollowing the same process as in Section 2.2.2, we look into ESS, simulation time,and trace plot of target parameter, µ = E(Y ).Scenario 1In Figure 2.10, we can see that ESS of µ first steeply decreases from 11487 to303 and slightly decreases, as n increases from 3 to 10000, when using Method 1.Also, as time increases rapidly, ESS/time converges to 0 as n increases. Method 2results have larger ESS when n is larger than 2000 in Figure 2.12; however, Method2 doesnât dominate in ESS/time due to longer performance time as in Table 2.5.One noticeable thing is that in both Figure 2.10 and Figure 2.11, the bottom-rightplot (n = 10000) shows strong evidence for the hypothesis that the model is iden-tified since the posterior variance of target parameter µ is very small. It is not clearwhether the model is mathematically identified or not. And the plots are suggestiveof it being identified, albeit weakly identified. This will be discussed further in thenext section and chapter.On the other hand, Method 3 shows the lowest ESS/time and ESS among thethree methods. Even though Method 3 only took a few seconds to a few minutesto execute the algorithm, ESS is usually near 1. Its behaviour is different from theresults in Example 4, where ISTP has better performance both in terms of ESSand time.n 3 10 30 100 300Method 1 0.5 1.0 2.4 8.3 23.6Method 2 7.1 8.0 11.1 19.2 44.7Method 3 6.5 6.7 7.1 8.2 11.2n 1000 2000 3000 5000 10000Method 1 98.9 220.3 337.5 561.4 1153.9Method 2 191.5 413.3 591.5 1036.4 2026.6Method 3 21.4 36.3 51.4 83.8 163.0Table 2.5: Simulation time of methods 1 - 3 (sec) as n increases.Scenario 2The results are a bit different in scenario 2. We expect the performance willworsen when the prior distribution allows further deviation from MAR as in Ex-41Figure 2.10: Method 1 (JAGS, original) Scenario 1 (The number of observa-tions) - Trace plot of µ .ample 4. However, from Ï = 0.01 to Ï = 0.1, ESS first decreases steeply, and thetrace plot shows autocorrelation or chains are not overlapped in Figure 2.13 andFigure 2.14. Then, ESS remains similar after a small increase in both methods 1and 2. Method 2 generally shows larger ESS and larger ESS/time, but ESS/timeresults are not much larger in Method 2 due to the longer simulation time as seenin Figure 2.12 and Table 2.6. Again, Method 3 has a small ESS and ESS/time.We suspect that the reasons behind the poor performance of Method 3 are:â¢ Reparameterization is based on an approximation using the first-order Taylorexpansion.42Figure 2.11: Method 2 (JAGS, Pseudo-TP) Scenario 1 (The number of ob-servations) - Trace plot of µ .â¢ Compared to Method 2, there is a further approximation that (Y | R = 1)follows the Poisson distribution. The importance weights in Method 3 takesthis approximation into account by not cancelling the ratio of true to con-venience likelihood. However, our results show that the convenience like-lihood approximation is not close enough to the true likelihood for effec-tive importance weighting. To find the weight, ratios between two posteriordensities are used. When finding them, an imputation of Y is applied toget the likelihoods. This imputation has two approximations/assumptions.E(Y | R = 0) = µ · exp(âÎ±1 · expit(α0 +α1µ)) is found with the first orderTaylor expansion, and (Y | R = 0) is assumed to follow the Poisson distribu-43Figure 2.12: ESS and time in each method and scenarioÏ 0.01 0.05 0.1 0.5 1 5 10Method 1 351.3 349.3 345.3 343.5 345.8 339.4 342.4Method 2 659.4 654.4 658.0 664.2 656.2 642.8 656.2Method 3 50.2 50.2 50.2 51.0 51.8 57.5 60.5Table 2.6: Simulation time of methods 1 - 3 (sec) as Ï increases.tion. Most normalized weights are near 0, and only a few or one is almost1. The imputed data convenience likelihood and true likelihood are bothfor the joint distribution of (Y,R), and this problem arises from the overalldiscrepancy between two joint-likelihoods f (y,r) and f â(y,r).44Figure 2.13: Method 1 (JAGS, original) Scenario 2 (Extent of MNAR) -Trace plot of µ .45Figure 2.14: Method 2 (JAGS, Pseudo-TP) Scenario 2 (Extent of MNAR) -Trace plot of µ .462.4 DiscussionIn this chapter, we tested TP and pseudo-TP in some scenarios of a model eitherknown or suspected to be partially identified. Since we try to set a model thatmakes scientific sense, but such a model can easily and unintentionally be onlypartially identified, it would be helpful if TP and pseudo-TP could resolve the issueof poor sampling performance with traditional MCMC. When a TP is available, itis shown that there is an improvement in the computational performance by usingISTP, especially in terms of compute time. However, when TP is unavailable, andconsequently, pseudo-TP is used instead, ISTP does not work in the example wetried. The issue arises from an unbridgeable gap between actual and conveniencelikelihood. In this section, we discuss the performance of importance samplingwith TP/pseudo-TP, its challenges, and potential future improvements.2.4.1 Performance of sampling with TP and Pseudo-TP in examplesPerformance in Example 3 - Transparent reparameterization in saturatedlogistic regressionIn Section 2.2, we tested some scenarios involving a partially identified model.When an off-the-shelf MCMC algorithm is applied to draw a sample from theposterior distribution in the partially identified model, it is expected to break down.Two scenarios focus on the size of the problem. Since the MCMC method updatesthe chain of parameters based on data, the computational performance is expectedto worsen when the size of the problem is larger. Indeed, when the number ofobservations in the data, n, or the number of parameters in X , p, increase, samplingtime increases steeply. Also, in scenario 1, where n is varied, ESS is small in nottrivially small n for Method 1. The large ESS of Method 1 in small n is mainlybecause the posteriors are close to priors. When an off-the-shelf MCMC algorithmis used for sampling from a joint prior density, it usually does not result in highautocorrelation.On the other hand, as ISTP (Method 2) allows the usage of i.i.d. Monte Carlosampling in Example 4, thanks to the conjugacy, the number of observations isnot the issue for computational performance. Similarly, increasing the number of47parameters p doesnât have much impact except for the extreme case, p = 5. It wasshown that ESS is almost 0 when p = 5, where the dimension of θ is over 100.Regardless of ESS, the simulation time is generally stable around a few secondsto a minute in all scenarios. While the number of parameters in X , p, impactsthe importance weight calculation time, this will not be a big problem unless thenumber of total parameters in the problem is extremely large.The other scenario concerns the prior variance, Ï2, representing the deviationfrom MAR to MNAR. We found that both ESS/time and ESS decrease as Ï in-creases. While the results of Method 2 also show the break down of ESS in thelarge variance case of Ï2 = 100, ESS/time is larger than Method 1 for all variancevalues.There are two main concerns: (i) ESS/time and (ii) compute time. Generally,Method 2 shows better performance than Method 1 in terms of these concerns.This somewhat agrees with the comments about ISTP in [45]. When the authorsexplained this method, they mentioned that ISTP could be an effective alternativefor MCMC methods in problems with a partially identified model. However, theywere also concerned that ISTP relies on the availability of finding a TP and thismay not be easy, or even possible, especially when there are many parameters inthe problem. In this problem, even though Method 2 is not working well whenthere is a large number of parameters, it is not about finding a reparameterizationbut instead about the dimensions of parameters. In Example 4, since parametersare highly structured, finding the TP when p is large is not hard. It would still bethe same for extremely large p as long as parameters are structured. However, asthe authors commented in [45], ISTP can only be implemented when researcherscan find an appropriate TP.Performance in Example 5 - Pseudo-TP in the count dataThe best scenario is for a TP to be available for the given problem. However,as seen in Section 2.3, even a small change in the problem can make it hard orimpossible to find a TP. However, if we find a reparameterization that behavessimilarly to a TP, it is still worth considering whether it improves computationalperformance, while acknowledging the drawback arising from not being a TP.48In Example 6, we found pseudo-TP based on Taylor approximation. Then, weperformed simulations in three ways: Method 1: JAGS with original parameteri-zation, Method 2: JAGS with pseudo-TP, using a convenience prior and the actuallikelihood l(dn | hâ1(Ï ,λ )), and Method 3: importance sampling with pseudo-TP,using a convenience prior and a convenience likelihood. Method 2 performs betterthan Method 1 in terms of ESS and trace plots. However, the compute time is al-most doubled, and still, MCMC sampling with the actual likelihood l(dn | θ) usedto update the states of θ has the autocorrelation issue in parameters that explainR and X ,Y . In the actual likelihood l(dn | θ), R is still explaining the existenceof Y while Y is explaining the probability of R being 1. So, we expected that im-portance sampling with pseudo-TP can be a good option since Method 3 will usei.i.d. Monte Carlo sampling using R base functions, hence, is faster than MCMCsampling. However, to use the i.i.d. Monte Carlo sampling, another assumptionis required for the conjugate convenience priors. While q = E(Y | R = 1) is foundwith the Taylor approximation, it does not give any information regarding the dis-tribution of Y | R = 1. It is the same for Y | R = 0, which is required for evaluatingthe complete-data likelihood l(dn | θ). Therefore, simulation Method 3 is based onthe assumption of (Y | R = 1) and (Y | R = 0) following the Poisson distribution.With these assumptions, importance sampling with pseudo-TP fails, as ESS is toosmall. As briefly discussed in the previous section, the main reason is that the gapbetween true and convenience likelihood is too large. Method 2 and 3 have thetrade-off between the probably more accurate computational result by using lessapproximation and the saving of time by using conjugate convenience prior, whichenables i.i.d. Monte Carlo sampling. At least in this example, unless the issueof tiny ESS is resolved by reducing the gap between true and convenience likeli-hood, even if Method 2 yields mixed results - better ESS and trace plot but longersimulation time, Method 3 canât be an alternative for Method 2, or even Method 1.The other notable thing in Example 6 is the challenge in discriminating be-tween a partially identified versus a weakly identified model. In Section 2.1.1, wediscussed that mathematically, it is hard to distinguish those two types of mod-els. It is not clear whether this model for count data is fully identified or partiallyidentified. After performing simulations, we found that as the number of observa-tions increases, the posterior variance of the target parameter becomes really small.49This happens when α1 is learned through data, which suggests identification of themodel. However, at the same time, the ESS decreases significantly compared to theresult of small n, which is the expected behaviour of a partially identified model.This will be discussed further in the next chapter.2.4.2 Diagnostic tools for the existence of TPThere is no algorithmic resolution on whether a TP exists for a given MCMC prob-lem. Even if the investigator canât find a TP, that doesnât necessarily mean thata TP doesnât exist. From the researchersâ point of view, it is essential to decidewhether to proceed further to find a TP in terms of efficiency in the computationalprocess. For that purpose, a diagnostic tool to check the low possibility of the ex-istence of a TP could play a crucial role. In Section 2.1.2, there are two conditionsthat a transformation needs to meet to be transparent reparameterization. However,these conditions can only be checked after a putative transformation is found. So,these conditions are necessary conditions for some transformation to be a TP, notnecessary conditions for a given problem to have a TP. If we can find the necessaryconditions for the given problem to have TP, it will be easier to check whether TPis available without ad-hoc attempts to find a TP. We speculate that these necessaryconditions are related to the Jacobian matrix of TP function and Fisher informationof the parameters in the problem. In the next section, we investigate further on this.2.4.3 Improved version of pseudo-TPThe concept of a pseudo-TP has a lot of room for improvement. Since pseudo-TP is a broad term that covers any transformation that is not a TP but has someintuitive flavour of a TP, a given pseudo-TP may work well, or may fail to workbecause it relies on overly rough approximations. Therefore, the goal is to mini-mize the gap arising from approximation while still ensuring the reparameterizingtransformation is explicitly invertible.One possible way forward is to use an incremental or adaptive method. In Ex-ample 6, at least two approximations are combined: a Taylor approximation and aPoisson distributional approximation. The impact of all assumptions is not the sim-ple summation of the effects of each assumption. As seen in the simulation results50of Method 2, at least the Taylor approximation doesnât break the performance ofposterior sampling. However, simulation results show that Method 3 fails to work.The difference between methods 2 and 3 is the addition of the Poisson assumption.It is appealing to use conjugate convenience prior and, therefore, use i.i.d. MonteCarlo sampling. However, if the importance sampling with pseudo-TP can applyapproximations one by one gradually instead of applying two approximations si-multaneously, the sampling performance may get better.Another possible improvement is measuring the extent of deviation of pseudo-TP from TP. As seen in the simulation results of Example 6, it looks like ourpseudo-TP is not close enough to a TP. However, we do not have a clear ideaabout how close or far it is from a TP. This criterion, which represents the extent,could help us to expect and could give a chance to improve the performance ofpseudo-TP.2.4.4 Weakly vs Partially identified modelIn Example 6, the model is not determined to be fully identified or partially identi-fied. We patently recognize the data has the issue of nonignorable missingness, butthat doesnât necessarily indicate that the model is partially identified. Indeed, theresult with the large sample size in Figure 2.10 suggests that the model is weaklyidentified. However, this information is available only after the actual MCMCsampling is done and is not easily recognizable in the model specification stage.If we scrutinize the information that can be learned directly from the data moreclosely, we might prove the identification of the model. In this Example 6, it is pos-sible to learn the probability cy = Pr(Y = y,R = 1) = Pr(R = 1 | Y = y)Pr(Y = y)for any non-negative integer y. And, we can easily find one-to-one function be-tween (a0,a1) and (α0,α1) where ai = Pr(R= 1 |Y = y) = expit(α0+α1y). Then,it is possible to represent the probability Pr(R= 1 |Y = 2) as a function of (a0,a1).Let this function be h(a0,a1) = expit(2logit(a1)â logit(a0)). The identification ofthe model is implicated by the number of solutions (µ,a0,a1) in the system ofequations:c0 = eâµa051Figure 2.15: Trace plots and f (x) graphs. The left plots are the f (x) graphsrepresenting the number of solutions for f (µ) = 0. The right plots arethe trace plots of µ . The top two plots use the data of Example 6 withn = 100, and the bottom plots use the data with n = 10000.c1 = eâµµa1c2 =eâµµ22h(a0,a1).If we simplify this, we are interested in the number of solutions µ in the given52equation:f (µ) =eâµµ22h(c0eâµ,c1eâµµ)â c2 = 0 (2.10)To investigate the relationship between the computational performance ofMCMC and the number of solutions, we first generate the dataset using the samemethod in Example 6 and perform MCMC sampling on it. Then, by using esti-mates cËi = #(Y = i,R = 1)/n, we find a number of the solutions. In Figure 2.15,when n = 100, there is no solution for f (µ) = 0, and the trace plot shows a partialidentification. And when n = 10000, there is one solution for f (µ) = 0. In bothcases, there is poor mixing in terms of a small ESS of 300 compared to the pos-terior sample size of 45000. However, there is a sharp reduction in the posteriorvariance and consistency with the true value of µ , suggesting at least weak identi-fication. In the next chapter, we investigate further to describe the identification ofthe model numerically.53Chapter 3Numerically diagnosing theidentification of partially orweakly identified modelsIn this chapter, we discuss methods to numerically diagnose the partially identifiedmodel. The main purpose of diagnostics is to inform researchers objectively andnumerically about the complications and difficulty of the statistical problem withpartially identified models. In addition, we use this information to warn about thecomputational difficulty of finding a TP.3.1 Background3.1.1 Importance sampling with transparent reparameterization(ISTP)This section summarizes the TP and ISTP algorithm as introduced in Chapter 2.Assume that the statistical model with parameterization θ is given for the dataDn. When h(θ) = (Ï ,λ ) is invertible, which implies that the dimension of param-eterization remains the same, dim(θ) = dim(Ï)+ dim(λ ), recall that the require-ment for h(θ) being a TP consists of two factors.1. Data dn is dependent only on Ï .54Then, Ï(dn | Ï ,λ ) = Ï(dn | Ï).2. Regular parametric asymptotic theory applies.Let estimator ÏË be obtained from Ï(dn | Ï). This estimator isân-consistent,ÏË âÏ = Op(1/ân).When a TP is available, the algorithm is given in Algorithm 1.3.1.2 Condition number of the matrixWhile the term condition number has a quite general and broad meaning in Mathe-matics, we focus only on the condition number of a matrix. When there is a systemof linear equations Ax = b, then the condition number [2, 64] of a square matrix Aiscond(A) = â¥Aâ¥ · â¥Aâ1â¥, (3.1)where â¥Aâ¥ is norm of a matrix A,â¥Aâ¥= supx Ì¸=0â¥Axâ¥â¥xâ¥ .Generally, L2-norm (spectral norm) is used, but it can be any p-norm. A matrixwith a large condition number is called ill-conditioned, while a matrix with a smallcondition number is called well-conditioned. These two terms are especially rele-vant for solving linear equations. Finding a solution is much more difficult whena matrix A is ill-conditioned. In the worst case, the condition number is infinite,which means that A is a singular matrix. When A is a singular matrix, it is possiblethat no solution exists for Ax= b. Even if a solution exists, the iteration method forfinding a solution, such as the Newton-Raphson method, takes too many iterationswhen an iterative method is used [2].We suspect that this condition number is somewhat related to overall compu-tational difficulty of analysis with the given suspected partially identified models,including finding a TP and Bayesian inference using an off-the-shelf MCMC algo-rithm. This is because a procedure of analysis gives a hint of finding solutions x inAx = b where the parameters are x.553.1.3 Fisher informationIn statistical inference, the Fisher information represents the information about theparameters that we can get from a given parametric model. This thesis focuses onthe Fisher information matrix for a vector of parameters θ . The Fisher informationmatrix is given byI(θ) = E[Var{âÎ¸ log f (X ;θ)}]= E[Var{âÎ¸ log f (X ;θ)}]= E[{âÎ¸ log f (X ;Î¸)}{âÎ¸ log f (X ;θ)}T ],where f (X ;θ) is the density function of the observed data X , which is a likelihoodfunction L(θ) when we treat it is a function of θ given observed data X [37, 50].Under some regularity conditions, I(θ) is the same as the expectation of the nega-tive Hessian of the log-likelihood,I(θ) = E(ââ2Î¸ log f (X ;θ)) .Forchini and Hillier [13], Stock and Wright [58] discussed the identificationof the model and the singular and non-singular Fisher information matrix. Dur-ing the discussion, the relationship between weak identification and an almost ill-conditioned, almost singular Fisher information matrix was also covered. There-fore, in the next section, we will discover the Fisher information matrix and sus-pected partially identified models.563.2 Representation for the identificationThe last section of Chapter 2, Section 2.4.4, briefly discussed the relationship be-tween the number of solutions in the system of equations for parameters and themodelâs identification. In that example, even though the model is assumed to bepartially identified due to nonignorable missing entries, the trace plot suggests thatthe model might actually be identified in some cases. In this section, we presentmethods to represent the extent of the identification in a model that is suspected tobe only partially identified. Below is a recap of the Example 6:Y â¼ Pois(µ),R | Y â¼ Bern(expit(α0+α1Y )).Some entries of the non-negative integer variable Y are missing, with binary vari-able R indicating observation of Y . Researchers are interested in estimating theparameter µ , but there is doubt about whether Y is missing at random (MAR) ormissing not at random (MNAR).The model is presumed to be partially identified in terms of the relationshipsbetween the variables R and Y . This is because the data has missing entries and itis unclear whether the missing entries are dependent on or independent of Y. Themissingness of Y , which R indicates, is explained by Y in the model. Althoughthere are only three parameters to estimate, (θ = µ,α0,α1), the data may not pro-vide enough information, resulting in a loss of identification.There were two noteworthy things in the simulation results with this model.The off-the-shelf MCMC algorithm, JAGS with Gibbs sampling [6, 16, 46, 47],was shown to not work well with a partially identified model. Therefore, it wasworth considering a transparent reparameterization. However, since Y is a countvariable, even though the model is simple enough, there was no way to find aclosed-form expression for Pr(Y = y | R = 1) and consequently a transparent repa-rameterization. At the same time, implementing the JAGS analysis using the orig-inal model was convenient since a closed-form expression of the joint posteriordensity of (θ ,Ymiss) is available. However, simulation results showed poor MCMCperformance due to the presumed partial identification.However, interestingly, some simulation results when the number of observa-57tions is large in Figure 2.10 hinted that the model may be at least weakly identified.Although the model remains the same for all simulations, the differences arise fromthe number of observations in the data and the prior variances of α1. We start bycharacterizing properties of the model that hint at weak identification.Our strategy requires evaluating the Fisher information matrix; however, cal-culating this is not easily feasible for a model without a closed form likelihoodfunction, such as Example 6. Fortunately, it is possible to mutate this model intoanother model with an easy-to-calculate Fisher information matrix.We briefly introduced finding of the number of solutions in the equations re-lating the learnable information from the data and the original parameters θ inSection 2.4.4. When we define some probabilities pi = Pr(Y = i,R= 1), which aredirectly estimable, and ai = Pr(R = 1 | Y = i), the following system of equationsarises:p0 = eâµa0,p1 = eâµµa1,p2 =eâµµ22a2.(3.2)In this case, a2 can be expressed as a function of (a0,a1), since Pr(R = 1 | Y =i) = expit(α0+α1i). Let this function be h(a0,a1). In addition, finding (µ,a0,a1)is a sufficient condition for finding (µ,α0,α1). Then, we define a model withan aggregation of observable data R and (Y | R = 1) into four categories, suchthat g(θ) = (p0, p1, p2) parameterizes a multinomial model with four cells, as perTable 3.1.Data Counts ProbabilityCategory 0 Y = 0,R = 1 x0 p0 = Pr(Y = 0,R = 1)Category 1 Y = 1,R = 1 x1 p1 = Pr(Y = 1,R = 1)Category 2 Y = 2,R = 1 x2 p2 = Pr(Y = 2,R = 1)Category 3 Y > 2,R = 1âªR = 0 x3 p3 = 1â p0â p1â p2Table 3.1: Multinomial model with four cells.Given this transformed model, we focus on two properties: 1. Given p, thenumber of θ values that satisfy g(θ) = p, which is the number of solutions for µin Equation 2.10. 2. The properties of the Fisher information matrix I(θ).58The number of solutionsLet Pl be the simplex of length l. The question could be represented as - For agiven p = (p0, p1, p2, p3) âP4, how many θ âÎ satisfy Equation 3.2?The properties of the Fisher information matrix I(θ)We first define the transformed model for general k, where k = 2 in the currentexample. In model Mk, Y and R are aggregated into k+ 2 categories, such thatcategory membership is completely determined by the data that can be observed:(Y = 0,R = 1),(Y = 1,R = 1), . . . ,(Y = k,R = 1),(Y > k,R = 1âªR = 0).Then, a multinomial model is established with cell probabilities pâPk+2. Movingfrom the original model to Mk, the data are aggregated and some information is lost.However, although finding a closed-form Fisher information matrix for the originalmodel is hard, finding a Fisher information matrix for the multinomial likelihoodfunction is quite easy and simple. For example, for the current model M2, the like-lihood function of multinomial distribution in terms of p âP4 = (p0, p1, p2, p3)isL(p0, p1, p2, p3) =Cpx00 px11 px22 px33=Cpx00 px11 px22 px33 ,however, since p3 depends on (p0, p1, p2), considering the degrees of freedom, p3is substituted with p3 = (1â p0â p1â p2). So the Fisher information matrix ofM2, if parameterized in terms of p0, p1, p2, isI(p0, p1, p2) =ï£®ï£¯ï£°np0+ n1âp0âp1âp2nnâp0âp1âp2n1âp0âp1âp2nnâp0âp1âp2np1+ n1âp0âp1âp2n1âp0âp1âp2n1âp0âp1âp2n1âp0âp1âp2np2+ n1âp0âp1âp2ï£¹ï£ºï£» .Then, the Fisher information matrix of M2 with respect to parameterization θ isI(θ) = âÎ¸gT2 I{g2(Î¸)}âÎ¸g2, (3.3)where gk(θ) is mapping from θ = (µ,α0,α1) to (p0, p1, . . . , pk).59We could use this Fisher information matrix directly. However, we alsoconsider another intuitive parameterization that may help to understand the be-haviour of changes in identification when the model deviates further from MAR,η = (µ,α0 +α1µ,α1), by separating controls of observed data rate and extent ofdeparture from MAR, respectively. The respective parameters are the target µ ,the approximate observed data rate α0 +α1µ , and the explanation of the extentof departure from MAR α1. As shown in Section 2.3.1, the probability of ob-served, Pr(R = 1), is approximately expit(α0 +α1µ). Therefore, the parameterα0 +α1µ controls the observed data rate. Using the same rule for finding I(θ) inEquation 3.3,I(η) = âÎ·hT I{h(Î·)}âÎ·h, (3.4)where h(η) = θ .The singularity of the Fisher information matrix was shown to be related to theidentification of the models in [53, 60]. So, when the Fisher information matrix isavailable, we first check whether it is invertible. If it is invertible (nonsingular), wemay need further information about the nonsingular matrix, as we are more inter-ested in the extent, i.e., whether it is close or far from being singular. At least twonumerical properties help us to understand this; the first is the rank, discussed in[27], and the other is the condition number. In the paper [27], the authors discussedthat the identifiability of the transformation is related to the rank of Jacobian of thetransformation. In the perspective of rank and condition number both describingthe singular property of the matrix, we expect that those properties are helpful.In Section 3.1.2, we briefly discussed that a condition number equal to infinityimplies that the matrix is singular. Furthermore, if the L2 (Euclidian, Spectral)norm is used for the calculation of the condition number, the resulting value isactually the ratio between maximum and minimum singular values [2]. Therefore,a condition number based on the L2 norm of the Fisher information matrix is agood candidate for representing the identification of the model.We set up an experiment with R using Example 6 to test the relationship be-tween condition number and the extent of model identification. The experimentâsmethod is summarized below:1. Make a grid of values θ = (µ,α0,α1) with µ = 2 and 21 values of each60α0+α1µ and α1 from â5 to 5.2. Calculate the values of (p0, p1, p2) of each θ by Equation 3.2.3. Report the condition number of the Fisher information matrix in Equa-tion 3.3 with n = 1000.4. Additionally, report the asymptotic variance of an estimator of µ , Avar(µË),which is the first entry of the inverse of the Fisher information matrix.When α1 Ì¸= 0, as the absolute value of α1 increases, the missing data mechanismdeviates more from MAR. In the other direction, if the researcher assumes the dataare MAR, we can reduce the model to the two-parameter problem with α0 and µ .The model is now identified.In Figure 3.1, the most notable thing is that the condition number blows upto infinity when α1 = 0. We will return to discuss this further later in the currentsection. The value α0 +α1µ governs the missing rate. When this value is small,the condition number increases as α1, instead of |α1|, increases. We suspect thatit is because there are not many missing entries. When α0 +α1µ is large, thecondition number is similar for different values of α1, possibly because only asmall portion of the data is observed. The condition number generally increaseswith |α1| when α0 +α1µ is close to 0. The asymptotic variances Avar(µË) havepatterns quite similar to the condition numbers. Taking into account these results,the condition number and asymptotic variance are related to the identification ofthe model. However, this relationship does not seem to apply to extreme caseswhere there are very few or very many missing entries.We conduct another experiment using JAGS [46, 47] with priors:µ â¼ Gamma(1,1),α0 â¼ N(0,10),α1 â¼ N(0,Ï2),where Ï = (0.01,0.1,1,10), to test actual off-the-shelf MCMC computational per-formance. When generating the dataset, while α0 = â1 and µ = 2 are fixed forall cases, α1 increases from 0 to 5. There are 1000 observations in the data.61Figure 3.1: Condition number of Fisher information matrix I(η) and asymp-totic variance of the estimator of µ , Avar(µË), in model M2. Contourplot that shows the condition number of the Fisher information matrixI(η) and asymptotic variance of the estimator of µ , Avar(µË), of modelM2 when µ = 2. All condition numbers and variances are log-10 scaled.White points indicate infinity.62In Figure 3.2 and Figure 3.3, when (α0,α1) = (â1,0), the performance is rea-sonable. There is no obvious difference between chains. Considering the largecondition number seen in the previous experiment, this is slightly different fromwhat is expected. The possible reason is that the model with (α0,α1) = (â1,0)has a singular/almost-singular Fisher information matrix and, therefore, has fewerproblems with multiple solutions, which does not cause bad computational per-formance. However, in Figure 3.4 and Figure 3.5, when α1 increases, trace plotsshow that computational performance is not good, which somewhat agrees withthe larger condition numbers in the contour plot.On the other hand, we have not yet concluded whether M2 is identified evenwhen the Fisher information matrix of M2, I2, has a small condition number sincewe have not considered the number of solutions of the system of equations in Equa-tion 3.2. Figure 3.6 describes the result of a small experiment with fifty probabilityvectors that shows the number of solutions of g(θ) = p given p âP4. Using theprobabilities p0, p1, p2 that arise from µ = 2 and (α0,α1) = (â1,â0.5), we applysome noise to probabilities by multiplying them by Îµii.i.d.â¼ Uni f (0.95,1.05), fori = 0,1,2. Since the original values of p0, p1, p2 are smaller than 0.05, applyingthe noise will not result in an out-of-bound issue. Then, we find the number ofsolutions in a range of µ â (0,v(p0, p1,a0,a1)), by finding roots of f (µ) wheref (µ) is reduced from the following equations:a0 = eµ p0,a1 = eµ p1/µ,h(a0,a1) = eµ p2/µ2.(3.5)The range function v(p0, p1,a0,a1) ensures that all found solutions are valid, bysatisfying inequalities a0 = eµ p0 < 1 and a1 = eµ p1/µ < 1. Then, with isolatingµ , the univariate function f (µ) is,f (µ) = eâµµ2h(a0,a1)â p2.In Figure 3.6, we often have two solutions, both of which are valid. This resultprovides enough evidence to conclude that the model M2 is not identified. At the63Figure 3.2: MCMC trace plots of µ for each breakpoint of α1 in contour plotof M2 with fixed µ = 2,α0 = â1, and prior variance Ï = 0.01. Theprior variance Ï is for α1. Plot i,(i = 1,2, . . . ,21) is trace plot of thesimulation with α1 = 0.25(iâ1).64Figure 3.3: MCMC trace plots of µ for each breakpoint of α1 in contour plotof M2 with fixed µ = 2,α0 = â1, and prior variance Ï = 0.1. Theprior variance Ï is for α1. Plot i,(i = 1,2, . . . ,21) is trace plot of thesimulation with α1 = 0.25(iâ1).65Figure 3.4: MCMC trace plots of µ for each breakpoint of α1 in contour plotof M2 with fixed µ = 2,α0 = â1, and prior variance Ï = 1. The priorvariance Ï is for α1. Plot i,(i = 1,2, . . . ,21) is trace plot of the simula-tion with α1 = 0.25(iâ1).66Figure 3.5: MCMC trace plots of µ for each breakpoint of α1 in contour plotof M2 with fixed µ = 2,α0 =â1, and prior variance Ï = 10. The priorvariance Ï is for α1. Plot i,(i = 1,2, . . . ,21) is trace plot of the simula-tion with α1 = 0.25(iâ1).67Figure 3.6: The number of solutions of Equation 3.2. When there are twovalid solutions, a line is coloured with light blue and with light pinkwhen there is no valid solution.same time, there are curves that indicate that there is no valid solution. They havevalid cell probabilities (p0, p1, p2), but which are not in the image of g(), whichresults in no solution. While (p0, p1, p2) are valid cell probabilities when pi > 0and â2i=0 pi < 1, which corresponds to the probability simplex for four categories,the image of g2(θ) is only a subset ofP4.To verify that the image of the mapping is only a subset, we performed a smallexperiment for the distribution of the number of solutions in M2.1. Draw 10000 probability vectors p âP4 from Dirichlet(1,1,1,1).2. Count the number of solutions of the mapping g(θ) = p.In Figure 3.7, the left plot is for M2. In many cases, there is no valid solution, andonly a few cases have one or two solutions. So, while all are valid cell probabilities,only 7% have one solution, and 4% have two solutions. The right plot is for M3;there are no valid solutions at all. If we set the larger tolerance for finding numerical68solutions, we see a slightly larger number of valid solution cases, 65 out of 10,000,but still, it seems that the image of g3() is only a tiny subset ofP5.Figure 3.7: The number of solutions of g(θ) = p. The vector p followsDirichlet(1,1,1,1) or Dirichlet(1,1,1,1,1). The left plot is the distri-bution of the number of solutions of g2(θ) = p, pâ¼Dirichlet(1,1,1,1).The right plot is the distribution of the number of solutions of g3(θ)= p,pâ¼ Dirichlet(1,1,1,1,1).We have numerical evidence that M2 is not identified, but this does not meanthat our original model is not identified. If Mk is identified for some k, then themodels Mk+1,Mk+2, . . . are identified, which leads to the conclusion that Mâ, theoriginal count model with incomplete data, is at least weakly identified. So, weconduct the same experiment reported on in Figure 3.1, with increasing k.In Figure 3.8, as k increases, the dark purple area expands, corresponding to asmall condition number. Similarly, in Figure 3.9, the region, in terms of (µ,α0 +α1µ,α), where the asymptotic variance of µË is small, expands with larger valuesof k from 2 to 10.The right plot of Figure 3.7 shows that all cases do not have any solution. It69does not give any hint for identification, but at least it does not support the fact thatthe model is partially identified. Since those vectors p drawn from the Dirichletdistribution are not guaranteed to be in image of g2() and g3(), we tweak thisexperiment with different vectors p.1. Obtain θ values randomly sampled from µ â¼Gamma(3,1), α0 â¼N(0,2.52),Î±1 â¼ N(0,2.52) to calculate the cell probabilities p for M2 and M3. Forexample, with (µ,α0,α1) = (1.909,â1.693,1.033), cell probabilities are(p0, p1, p2, p3) = (0.023,0.096,0.160,0.138).2. Find the number of solutions g2(θ) = p. This is the number of solutions forM2.3. For each solution θ = (µ,α0,α1), check how many solutions satisfy p3 =exp(âµ)(µ)3h2 (a0,a1)/6. This is the number of solutions for M3.In this experiment, it is guaranteed that there is at least one solution, which is theoriginal θ used to calculate p. In Figure 3.10, in most cases, there is more than onesolution in g2(θ). Almost all of these cases have unique solutions in the mappingg3(θ). A few cases initially appear to have more than one solution; however, thisis due to the limitation of the iterative method to find a solution. Since finding aclosed-form solution of f (µ) is not easily feasible, we use R functions unirootand optimize. They use iterative methods to find a solution with a specifictolerance. So, although a unique solution exists, these functions can report thesolution that yields f (µ) values close enough to 0, considering the tolerance. Inaddition, the difference between the two solutions is usually small. For example,one case has two solutions (2.575757, 2.576063). It is hard to tell the differencebetween these two numbers are meaningful. Contour plots of asymptotic variancesand condition numbers and these experiment results numerically support M3, andMk with larger k being at least weakly identified. Unlike showing a counterexamplefor M2 is identified in an earlier paragraph is enough, we can only say that the resultis consistent with this example being at least weakly identified. However, we canstill conclude that the behaviours of condition numbers and asymptotic variancesalign with the identification of the model.70Figure 3.8: Condition number of Fisher information matrix when increasingk of model Mk.71Figure 3.9: Asymptotic variance of the estimator of µ , Avar(µË), when in-creasing k of model Mk.72Figure 3.10: The number of solutions of models M2 and M3 in various θvalues. The number of solutions in M2 and M3 when generatingθ = (µ,α0,α1) with µ â¼ Gamma(3,1), α0,α1 â¼ N(0,2.52).Lastly, we discuss the identification characteristic of the infinity condition num-ber when α1 = 0. A similar situation occurs in scenario 1 of the paper [19]. Belowis a summary of scenario 1:Example 7. Exposure model with imperfect tests.â¢ Variables: A binary variable E is used for an exposure where E = 1 indicatesan exposure. There are two binary variables T1,T2 that describe imperfecttests.â¢ Parameters: Sensitivity p1, p2 and specificity q1,q2 are for each test T1,T2.In addition, r = P(E = 1) is the probability of exposure. In total, there arefive parameters θ = (p1, p2,q1,q2,r).Given two-by-two table data in Table 3.2, the model is partially identified. There-fore, the author used stratification into two sub-populations to identify the model.73T1 = 1 T1 = 0T2 = 1 a bT2 = 0 c dTable 3.2: Two-by-two table with two imperfect tests T1,T2.So, in addition to (p1, p2,q1,q2), r is now divided into two probabilities r1 and r2for each sub-population with binary variable X , which is presumed to be relatedto the exposure. Then, the data is now a two-by-two-by-two table, with an addi-tion of X = 0 and X = 1 to Table 3.2. In addition, there are conditional indepen-dence assumptions between the test variables and X , T1 â¥ X | E, T2 â¥ X | E, andT1 â¥ T2 | E. With the sub-population model, an approximate root mean squarederror (ARMSE) for (pË1, pË2, qË1, qË2, rË) is measured for estimation. For the modelwith sub-populations, the author used the term âessentially identified.â As long asr1 Ì¸= r2, cell probabilities under the mapping arising from stratification could bemapped back to the original probabilities. However, when r1 = r2, the situation ischanged. Like the infinity condition number of α1 = 0 cases in Figure 3.1, ARMSEgoes to infinity when r1 = r2. As the author explained, when r1 = r2, it means thatthe sub-population arise from a completely random specification, and the informa-tion given by two sub-populations is the same, which means the Fisher informationmatrix is singular.In our Example 6, when α1 = 0, Y is missing at random. Dependence betweenR and Y is not directly estimable. A similar issue occurs in Example 7 when r1 =r2; exposure E and X are assumed to be dependent, but r1 = r2 case implies thatthey are independent, sub-population is stratified completely at random. Actually,there is not much information about α1 that the data can give. This is becausePr(R= 1 |Y ) = Pr(R= 1). Therefore, α1 = 0 makes the Fisher information matrixsingular, where mapping back is not possible, which results in the infinite conditionnumber. It is also somewhat related to the mixed results of simulations - sometimes,MCMC performance suggests a weakly identified model as in the right bottomplot of Figure 2.10, whereas usually, we see the MCMC performance results thatsuggest of a partially identified model as in the right top plot of Figure 2.10. Smallα1 generally gives the small condition number, except for α1 = 0. However, if it is740, a large condition number is given. As noted in [19], we should use a model likethis carefully, even if it shows some evidence of being identified.We also tested with the same simulation method for looking at the number ofsolutions with α1 = 0 as in Figure 3.10. In Figure 3.11, there is one solution forall cases, which is the original µ . However, at the same time, we found that thereare many cases in which a gradient optimization algorithm cannot find a solution.In Figure 3.12, we plot the corresponding f (µ) graphs for M2 for some cases.The lower right plot is for a case where the algorithm cannot find the solution.Theoretically, there must be at least one solution (the original µ), but the graphshows that f (µ) converges to zero near the solution instead of having a uniquesolution. We suspect that this kind of behaviour causes the blow-up point seen inFigure 3.1.Figure 3.11: The number of solutions of models M2 and M3 in various θ val-ues with fixed α1 = 0. The number of solutions in M2 and M3 whengenerating θ = (µ,α0,α1) with µ â¼ Gamma(3,1), α0 â¼ N(0,2.52)and α1 = 0.In this section, we explored how we can describe the identification of a given75Figure 3.12: Graphs of f (µ) of some θ values generating with α1 = 0. Eachcase is generated under θ = (µ,α0,α1) with µ â¼ Gamma(3,1), α0 â¼N(0,2.52) and α1 = 0.statistical problem involving a suspected partially identified model. We found thatthe Fisher information matrix is quite a meaningful numerical expression for iden-tification, considering both theoretical calculations of condition number and simu-lation results. At the same time, this method is quite limited to models that we cancategorize into a multinomial model framework. We need to complete this modelmapping case-by-case basis. If the Fisher information matrix of the original modelcan be easily calculated without aggregation mapping, it is possible to generalizethe identification representation.763.3 Exploring the existence of TPIn the previous section, we discovered how to explore the identification using thecondition number of the Fisher information matrix. Now in this section, we intro-duce a tool for measuring the level of computational challenge when using a givenpartially identified model issue when using an off-the-shelf MCMC algorithm. Themain goal of this tool is to diagnose the current model and give information aboutthe strength of model identification so that users can expect quality and computa-tional performance if they use an off-the-shelf MCMC algorithm.3.3.1 A tool for diagnosing the severity of partially identified modelThe obvious thing in Section 3.2 is that the condition number of the Fisher infor-mation matrix is highly related to the identification of the model. Therefore, thefirst step of a diagnostic procedure is to find the Fisher information matrix andevaluate it at the mean of the posterior samples. The tool will perform a posteriorsampling using JAGS with one chain of samples to enable this. The number ofsamples could be any desired number, with a default of 10000 samples.However, calculating a Fisher information matrix may not be easy in a partiallyidentified model. We need to differentiate the log-likelihood function twice andfind expectations, but there is no guarantee that they are in closed form. As inSection 3.2, it may be possible to transform the model to find a Fisher informationmatrix, but that is not a generalized method. So, we utilize a method to find anumerical approximation using the Monte Carlo (MC) method [55].Let the vector 1/v = [1/v1 . . .1/vn]T be the vector with reciprocals of all theelements in the vector v = [v1 . . .vn]T . Then, the method introduced in the paper issummarized in Algorithm 2. As long as the evaluatable log-likelihood function isavailable, the algorithm could approximate the Fisher information matrix withoutfinding the closed-form gradient, Hessian of the log-likelihood function, or expec-tation of Hessian.We are especially interested in the singularity and condition number of theFisher information matrix instead of the matrix itself, which eases the issue of usingapproximation instead of exact calculation. However, the singularity and scale ofa condition number of the approximated matrix are not necessarily equivalent to77Algorithm 2 Monte Carlo Resampling Method for Estimating Fisher informationmatrix [55]InputA specific point where Fisher information matrix is evaluated at, θThe number of parameters in θ , pSample size, nThe number of Hessian estimates in each pseudo dataset simulation, MThe number of pseudo dataset simulations, NLog-likelihood function, l(θ | dn)A small constant for the perturbation vector, cOutputEstimation of the Fisher information matrix1: for iâ 1 to N do2: Generate dataset din with θ .3: for k â 1 to M do4: Generate aik, aËik â Rp where all elements are i.i.d. with Rademacherdistribution.5: Perturbation vectors âik â caik and âËik â caËik.6: Gradient approximations G(θ ±δ ik | din)â(l(Î¸Â±âik+âËik|din)2)(1/âËik).7: The derivative of Gk, δGk â G(θ +δ ik | din)âG(Î¸ âÎ´ ik | din).8: Estimate of Hessian HË ik â 1/4(δGk(1/δ ik)T +(δGk(1/δ ik)T )T).9: end for10: Calculate H¯i = (âMk=1 HË ik)/M.11: end for12: Calculate FËn(Î¸) = (ââNi=1 H¯i)/N.13: Return FËn(Î¸).those of the true matrix. Therefore, considering the vulnerability of the conditionnumber to approximation error, we also use the rank of an approximated Fisherinformation matrix in this tool. While a condition number is a real number thatcan range from 1 to infinity, a rank is the number of singular values larger thantolerance (which is usually zero), which is less vulnerable to the gap induced bythe approximation. Unlike a condition number, the large rank indicates a well-conditioned matrix.In addition to the condition number C and the rank R of the Fisher informationmatrix, this tool returns the flag that indicates the expected computational diffi-78Algorithm 3 Diagnosing a partially identified modelInputThe number of parameters in θ , pSample size, nLog-likelihood function, l(θ | dn)Data, dnThreshold εOutputCondition number of a Fisher information matrix, CThe rank of a Fisher information matrix, RFlag for the high chance of absence of a TP1: Run JAGS for obtaining samples from posteriors.2: Let estimates Î¸Ë be sample mean of posteriors.3: if a closed-form Fisher information matrix is not available then4: Approximate Fisher information matrix by using. Algorithm 25: end if6: Evaluating the Fisher information matrix at Î¸Ë7: Calculating the condition number C and rank R of the Fisher information ma-trix.8: if C > ε and R < p/2 then9: Set Flag = TRUE.10: end if11: Return Flag, C, and R.culty. This flag is set to true when C is too large and R is too small compared tothe number of parameters in θ . Since there is no universal standard for the largecondition number and the small rank of an ill-conditioned matrix, we use C > 100and R < p/2.One last barrier to using this method is that it may not work well with in-complete data. This is because it is not possible to evaluate the likelihood func-tion. A possible resolution is to approximate the information matrix using the EM(Expectation-Maximization) algorithm [7, 38]. The original EM algorithm is pro-posed to find the maximum likelihood estimation (MLE) by iterations when somedata entries are missing. Although it is used in the frequentist context, it is alsoused in Bayesian inference with missing data, e.g., the Bayesian structural EMalgorithm [14]. It has been frequently discussed how to get an information ma-79trix when the EM algorithm is used to infer incomplete data [35, 36, 39]. As in[35, 36], it is possible to separate the likelihood function into observed and unob-served parts, and so on the information matrix. Similarly, in the Bayesian context,where missing data is also a parameter, the likelihood function isL(θ ,y | x) = f (x | y,θ)= f (x,y,θ)/ f (y,θ)= f (x,y | θ) f (θ)/ f (y | θ) f (θ)= L(θ | x,y)/L(θ | y).(3.6)where θ is parameters of interest, Y is unobserved, and X is observed data. UsingEquation 3.6, the final log-likelihood function is l(θ ,y) = l(θ ;x,y)â l(θ ;y) andtherefore the information matrix l(θ ,y) is I(θ ;x,y)â I(θ ;y). However, it is noteasy to find the difference between the Fisher information matrix of complete dataand unobserved data. We find that the method introduced in [51] may also work forapproximating this difference. The idea of this method is to use Monte Carlo simu-lation along with MCMC to approximate the expectation in the Fisher informationmatrix. We donât have a known solution yet, but we suspect that the following ideamight work: Whenever evaluating the log-likelihood function l(θ | dn) in Algo-rithm 2,1. Perform Monte Carlo i.i.d. sampling of unobserved data y.2. Evaluate l(θ ;x,y)â l(θ ;y) for each y.3. Average these values to approximate the likelihood of θ and the unobserveddata.803.4 DiscussionIn this chapter, we explored how to represent the extent of identification of a modelthat we suspect might only be partially identified using the Fisher information ma-trix. Since a condition number is a numerical expression that is highly related tothe invertibility of a matrix, which is a key point to solve equation Ax = b, we cal-culated the Fisher information matrix and the corresponding condition number. Tocheck the relationship between how much the model is partially identified, a condi-tion number, and computational performance, we first calculated theoretical valuesin the context of Example 6, calculated a Fisher information matrix, and performedsimulation by using off-the-shelf MCMC algorithm package JAGS [46, 47]. Sinceit is not easy to find a closed-form likelihood function of the model in Example 6,we used a transformed model with a closed-form likelihood function and an easy-to-find Fisher information matrix. Then, we introduced a diagnostic tool for check-ing the severity of a computational challenge for a given partially identified modeland whether there is empirical evidence concerning the prospect of finding a TP.We used an algorithm to approximate a Fisher information matrix [55] to seek uni-versality, i.e, an ability to apply to arbitrary models. We discussed a technique thatuses an idea of EM algorithm [7, 14, 38] to deal with a model with missing en-tries. When this tool is used, it is expected that a user can get information about theseverity of a given partially identified model and warn about whether there is a highchance of the computational burden. Following the above explorations, in this sec-tion, we discuss the challenge of this tool and possible directions for improvement.Below are two main challenges of the strategy we introduced:â¢ How to choose a good threshold?For raising a flag for computational difficulty, we set a default for the condi-tion number C > 100, and R < p/2 where p is the dimension of the parame-terization θ . While it is very obvious that C = 1 and R= p are the best casesand C = â and R = 1 are the worst cases, the middle area between the twocases is somewhat grey. At the same time, condition numbers are vulnerableto the scale of approximated values in a matrix A, unless all entries of thematrix A are scaled simultaneously, i.e. kA with some constant k. Therefore,a fixed threshold for the condition number of the approximated Fisher infor-81mation matrix has a risk of misclassification due to some large entries of thematrix.â¢ How to formally prove the result we have?We suggested using numerical expressions, such as a condition number orrank, to check the identification of a model that is presumed to be partiallyidentified. However, they are not mathematically proven to represent thisidentification. Instead, we used simulations to show they are associated.These two points challenge our strategy. However, this does not necessarily meanthat it is impossible to use the strategy in the real world. It is still a good strategyfor researchers to make a decision about their Bayesian analysis, such as whichMCMC algorithm to use. Based on these points, we discuss some further researchtopics related to the improvement strategy of the diagnostic tool.3.4.1 Objective thresholdAs discussed, the condition number is vulnerable to the scale of the matrix. Insteadof an absolute threshold, fixed for all cases, it is possible to consider a relativethreshold or standardized condition number, such as Skeelâs condition number [24,54]. The classical condition number is found with a multiplication of two norms,â¥Aâ¥ · â¥Aâ1â¥, as in Equation 3.1. On the other hand, Skeelâs condition number iscalculated by â¥|Aâ1| · |A|â¥ where |A| is element-wise absolute of the matrix. Boththese condition numbers are not affected when all elements are scaled together withthe same constant c. But, as in [54], the classical condition number is varied withrow scaling, while the Skeelâs condition number is not.3.4.2 What can we do when warning flag is raisedWhen users get a warning from the diagnosis, it does not mean that users cannotuse any off-the-shelf MCMC algorithm. Parallel tempering (PT) [62] could helpthe sampling from the posterior distribution that is not unimodal, which is highlyrelated to the partially identified model. Since the parameters are not uniquelyidentified in the partially identified model, which means that the distribution couldbe multimodal, a PT algorithm could help investigators get samples that could82adequately represent the posterior distributions. In the next chapter, we investigatethis in detail with a comparison to JAGS.83Chapter 4pimbayes: A comprehensivetool for Bayesian inference withpartially identified modelsThis chapter discusses the implementation of computational methods to diagnoseand resolve issues of the computational burden of off-the-shelf MCMC in Bayesianinference with partially identified models. Then, the R package pimbayes, whichexecutes the implementation, is introduced with examples and manuals.4.1 Alternative options for partially identified modelsIn Chapter 2, the usage of ISTP [22] is recommended to improve the computa-tional efficiency Bayesian inference with a partially identified model. In additionto this method, PT [11] may be a good option to address this issue. PT is a samplingmethod that uses N interactive chains, which is different from simply using N paral-lel chains in an off-the-shelf MCMC algorithm. While there are non-reversible andreversible PT algorithms, this chapter focuses only on non-reversible PT[62]. Apackage in julia called Pigeons.jl [61] enables the usage of non-reversiblePT. The requirement for Pigeons.jl is a model specification expressed in a lan-guage that the package can understand, for example, Turing.jl [15] or STAN[57].84This PT method is especially expected to be helpful when a TP is difficult tofind or is unavailable to researchers. This is because PT works well for multimodaland high-dimensional target distributions. When looking at the definition of a par-tially identified model, the fact that the parameters are not uniquely identified canresult in a multimodal posterior density.Below is the simulation to illustrate the idea by comparison between JAGS[46], which uses Gibbs sampling [6], and Pigeons.jl [61] results, which usesPT. Data of 1000 observations are generated under Example 6 and Equation 2.8using µ = 2,α0 = â2log(2)/â2,Î±1 = log(2)/â2. We use the same priors as inExample 6:µ â¼ Gamma(1,1),α0 â¼ Norm(0,10),α1 â¼ Norm(0,Ï2),where Ï = 0.01,0.05,0.1,0.5,1,5. Then, JAGS is run with R package runjags[8] and PT is run with julia [3] package Pigeons.jl [61] for sampling 10chains of 1024 posterior samples. Both simulations are parallelized on 32 cores.In Figure 4.1 and Figure 4.2, trace plots and density plots of µ are shown forthe simulations. In Figure 4.1, the trace plots of µ for Ï = 0.5,1 with JAGS showdeviations between chains and failure of mixing within chains, which is a badsignal regarding sampling quality. As seen in Chapter 2, the larger prior varianceshinder the computational performance. On the other hand, in Figure 4.2, the traceplots and the density plots are not very different for different prior variances. Onenotable thing is that even though both the trace plots in Figure 4.1 and Figure 4.2have 10 chains, the chains in Figure 4.1 are independent chains for the same targetand the chains in Figure 4.2 are at different temperatures, sequentially tempered[62]. In PT, the first chain iterates over the prior density, and the last chain is forthe target, which is the posterior density. Between these two chains, each chaininteracts and bridges between those two chains gradually [62]. In Figure 4.2, thelast chain is coloured with bright olive green.Differences between two figures are due to the nature of parallel temperingand Gibbs sampling. Therefore, we focus more on whether each trace plot and85ESS internally change by prior variances or not rather than comparing the twofigures. Reflecting on this result, PT is worth considering for Bayesian inferencein the partially identified model. One thing we notice during the simulation is thatmultiple hours are taken for Pigeons.jl. However, it is not an issue with thePT. Instead, the issue is with the input type we used to specify the model. We useTuring.jl [15] language for the model; however, it does not support the Markovblanket strategy [32], which will find only a subset of nodes that is informative.Since this example has many parameters where most of them are latent, executiontime can increase without a Markov blanket technique. Further investigation willbe discussed in the next section.86Figure 4.1: Trace and density plots of µ for JAGS simulation of Example 6with various prior variances.87Figure 4.2: Trace and density plots of µ for Pigeons simulation of Example 6with various prior variances. Each chain with a different colour in thetrace plot represents a chain of sequential tempered posterior densities.Trace and density curves of the last chain, which is for the target, arecoloured with bright olive green.884.2 R package pimbayes4.2.1 Package scopeR package pimbayes is designed for Bayesian inference with suspected partiallyidentified models. The functionality of the package is split into two parts:1. Diagnostic tools introduced in Section 3.3.The package helps to diagnose the statistical problem with the suspected par-tially identified model numerically. In detail, it answers two main questions- whether the model is identified or not and whether the off-the-shelf MCMCalgorithm will work well or not.2. Posterior sampling for the partially identified model.Basically, three types of posterior sampling are supported.(a) Applying ISTP in Algorithm 1 for saturated logistic regression withincomplete data (Example 4) or cross-sectional studies for imperfecttest (Example 1).(b) Applying ISTP manually when necessary elements including a TPfunction are provided.(c) Performing Non-reversible parallel tempering by Pigeons.jl [61,62] for the model provided without a transparent reparameterization(TP).In summary, if the model is a saturated logistic regression with missing entries inthe response variable or a model with cross-sectional studies for the imperfect test,users could apply all functionalities. If the model is not of these types, but usersknow a TP function, they can use them all except for 2 (a). Finally, if the model isnot of those types and TP function is unknown, the users could use 1 and 2 (c). Inthe next sections, we show the ways in which this package can be used.4.2.2 Technical backgroundFor almost all functionalities in the package pimbayes, it is necessary to parsethe userâs JAGS code. This enables the user to provide the JAGS code and data89directly without additional information, such as which part of the code correspondsto a prior distribution or a likelihood function. The package uses a parser fromthe julia package JuliaBUGS.jl [63]. In each line of the JAGS code formodel and prior specification, there are two sides: the nodes for the parametersand variables and the distributions or values for those nodes. JuliaBUGS.jlparses the JAGS code and saves this information in the form of a graph. Thepackage pimbayes uses the parsed information to identify the priors, likelihoods,and any necessary information for ISTP or diagnostic tools. At the same time,JuliaBUGS.jl provides a bridge for using Pigeons.jl. For using the non-reversible parallel tempering, an evaluation of log-likelihood and a sampling fromthe priors are required. JuliaBUGS.jl first compiles the JAGS code and dataand then provides this functionality.Lastly, in the background of all these functions, R package JuliaCall [34]supports the usage of all functionalities from julia in R.4.2.3 Installation and setupBefore installing the package, an installation of JAGS [46] is required. Then, thepackage is available via a GitHub repository.library(remotes)install_github("https://github.com/serenlee/pimbayes", dependencies = TRUE)This will install all necessary packages. Since pimbayes requires julia andits package, the users need to proceed with the initial setup. A function setupincludes the installation of the required packages from julia.setup()If a machine does not have julia,setup(installJulia = TRUE)It is also possible to specify the path of the installed julia binary.setup(JULIA_HOME = "custom_julia_path")4.2.4 Applying ISTP to Cross-sectional studies exampleTwo cases of applying ISTP are pre-coded in the package. The first is for Exam-ple 1. The transparent reparameterization is found in [45]. The parameterization Ï90is (η11,η12,η21) whereη11 = Pr(T+â©D+)= SeÏ11+(1âSp)Ï21Î·12 = Pr(T+â©Dâ)= SeÏ12+(1âSe)Ï22Î·21 = Pr(Tââ©D+)= SpÏ21+(1âSe)Ï11Î·22 = Pr(Tââ©Dâ)= (1âSe)Ï12+SpÏ22.Then, λ is (Se,Sp). The package stores a mock dataset, which is created withspecificity and sensitivity of 0.95 and Ï = (0.15,0.35,0.20,0.30). The codes forcreating the mock data are also available in the folder data-raw.data(cross)data <- list(y = cross, n = sum(cross))JAGS code with priors used in [45] for this example is as follows:jags_code <- "model{y[1:4] Ë dmulti(eta[1:4], n)eta[1:4] <- c(se*pi[1]+(1-sp)*pi[3],se*pi[2]+(1-sp)*pi[4],sp*pi[3]+(1-se)*pi[1],(1-se)*pi[2]+sp*pi[4])pi[1:4] Ë ddirch(c(1,1,1,1))se Ë dbeta(25,3)sp Ë dbeta(30,1.5)}"In this code, the names of the parameters in θ and the prior distributions are ed-itable. Then, results are returned with ESS and N samples from the posterior den-sity.results = sampling(jags_code, data, type = "cross", N = 1000)results$ess#> [1] 964.7134head(results$samples, 1)#> pi[11] pi[12] pi[21] se sp#> [1,] 0.1416467 0.3382962 0.2146765 0.9219155 0.893652791All results are returned in the type of list in R.4.2.5 Applying ISTP to Saturated logistic regression model examplePackage pimbayes also supports importance sampling with transparent reparam-eterization and convenience priors for an Example 4 of saturated logistic regressionwith incomplete outcome data introduced in Chapter 2.The data have three components: the binary response variable Y , the binaryindependent variables X = (X1, . . . ,Xp)T , and the indicator variable R, which indi-cates the missingness of Y . The binary variable R is 1 when Y is observed and 0when Y is unobserved.The following is a summary of Example 4: A design matrix of X and a rowvector Xdes â R2p of it is used. Exactly one of the 2p entries, from X0...0 to X1,...1,is one, and each entry Xi1,i2,...,ip is one if and only if the corresponding logic is true,X1 = i1,X2 = i2, . . . ,Xp = ip. Parameters of interest are α,β ,γ,δ , where the modelis specified in the model part of Example 4 and Equation 2.4.We include two mock datasets, small sat and sat. A small simulation withn = 100 and p = 2 is performed to create the data small sat.data(small_sat)head(small_sat$x)#> [,1] [,2] [,3] [,4]#> [1,] 1 0 0 0#> [2,] 0 0 1 0#> [3,] 0 0 1 0#> [4,] 0 0 1 0#> [5,] 0 0 1 0#> [6,] 0 0 0 1head(small_sat$y)#> [1] NA NA 1 1 NA 1When data and JAGS code are given with the desired posterior sample size, thefunction sampling() will return posterior samples along with ESS.n = length(small_sat$r)p = log2(dim(small_sat$x)[2])sigma = 0.5 # Hyperparameterdata <- c(small_sat, list(n = n, p2 = 2Ëp, tau = 1/sigmaË2))results = sampling(jags_code, data, type = "sat", N = 1000)results$ess#> [1] 69392head(results$samples[,1:7],3)#> px[2] px[3] px[4] py_x[1] py_x[2] py_x[3] py_x[4]#> [1,] 0.0781 0.539 0.180 0.130 0.02074 0.274 0.650#> [2,] 0.0470 0.569 0.195 0.209 0.04500 0.215 0.555#> [3,] 0.0383 0.579 0.176 0.136 0.00324 0.342 0.551head(results$samples[,8:11],3)#> pr_xy0[1] pr_xy0[2] pr_xy0[3] pr_xy0[4]#> [1,] 0.326 0.512 0.402 0.807#> [2,] 0.605 0.267 0.435 0.749#> [3,] 0.430 0.344 0.440 0.946head(results$samples[,12:15],3)#> log_or_r[1] log_or_r[2] log_or_r[3] log_or_r[4]#> [1,] -0.5529 0.1063 -0.238 -0.177#> [2,] -0.1303 0.0203 0.568 -0.809#> [3,] -0.0817 0.1216 -0.167 -0.3494.2.6 Manually performing ISTPIf required functions and variable names are provided, sampling with the ISTPwill be performed with a function istp.manual().The required arguments are as follows.1. posterior fun: A function generating samples from convenience poste-riors, which is proportional to the multiplication of convenience priors andlikelihoods.posterior_fun<-function(N){eta<- rdirichlet(N,data+1)se <- rbeta(N,25,3)sp <- rbeta(N,30,1.5)phi_lambda_matrix = cbind(se,sp,eta[,1:3])colnames(phi_lambda_matrix) = c("se","sp", paste0("eta[",1:3,"]"))return(phi_lambda_matrix)}2. h inv: hâ1(Ï ,λ ) = θ .h_cross_inv <- function(phiLambda) {eta11 = phiLambda[3]eta12 = phiLambda[4]eta21 = phiLambda[5]se = phiLambda[1]sp = phiLambda[2]eta22 = 1 - eta11 - eta12 - eta2193denom = se+sp-1pi11 = 1/denom*(sp*eta11 + (sp-1)*eta21)pi12 = 1/denom*(sp*eta12 + (sp-1)*eta22)pi21 = 1/denom*((se-1)*eta11 + se*eta21)return(c(se,sp,pi11,pi12,pi21))}3. N: Sample size (output), desired size of Monte Carlo output4. weight fun: The function that calculates importance weights.weight_fun<-function(theta_matrix, deter) {apply(theta_matrix, 1,function(theta) ifelse(valid_A(theta),(theta[1]+theta[2]-1)Ë(-2),0))}It will take theta matrix, a matrix of posterior samples of parameters inoriginal parameterization, and deter, a vector of determinants.5. var name: A vector of names of variables.When we run the code with above arguments, we observe the following outputs.results = istp.manual(posterior_fun,h_cross_inv,N = 15000,weight_fun,var_name)results$ess#> [1] 14454.38head(results$samples)#> se sp pi[11] pi[12] pi[21]#> [1,] 0.8074928 0.9417783 0.27922344 0.3398368 0.02743107#> [2,] 0.9258622 0.8667851 0.09675919 0.3384419 0.17933028#> [3,] 0.8961275 0.9853014 0.27886598 0.5304348 0.08137019#> [4,] 0.9165561 0.9081519 0.16834859 0.2533007 0.22057400#> [5,] 0.9503819 0.9579694 0.17877832 0.3971520 0.18126433#> [6,] 0.9558087 0.9830187 0.09087908 0.2064749 0.386384694.2.7 Diagnosing the identification of the modelBasically, diagnose identification takes three steps.941. Approximate the Fisher information matrix.2. Find the condition number and rank of this matrix and, if invertible, theasymptotic variance of the parameter estimators.3. If warranted, raise a flag that warns about the partial or very weak identifica-tion in a given model.The following five arguments are required.1. jags code2. data3. variable.names4. ll fun A log-likelihood function that takes data and theta as argu-ments. The type and structure of data must match the return type ofdata fun.ll_cross <- function(data, theta){se = theta[1]sp = theta[2]pi = theta[3:5]pi = c(pi, 1-sum(pi))eta = c(se*pi[1]+(1-sp)*pi[3],se*pi[2]+(1-sp)*pi[4],sp*pi[3]+(1-se)*pi[1],(1-se)*pi[2]+sp*pi[4])log(dmultinom(data, prob = eta))}5. data fun A function that creates data with a given theta (a vector of pa-rameters). data fun will take an argument args, which includes thetaand the input necessary to generate the data.data_cross <- function(args) {theta = args$thetan = args$nse = theta[1]sp = theta[2]pi = theta[3:5]pi = c(pi, 1-sum(pi))95eta = c(se*pi[1]+(1-sp)*pi[3],se*pi[2]+(1-sp)*pi[4],sp*pi[3]+(1-se)*pi[1],(1-se)*pi[2]+sp*pi[4])as.vector(rmultinom(1, size = n, prob = eta))}The above arguments are written for small cross in Section 4.2.4. We get thefollowing outputs.diagnose_identification(jags_code, data,c("se","sp","pi[1]","pi[2]","pi[3]"),ll_cross, data_cross)#> $condition_number#> [1] 266.7476#> $rank#> [1] 5#> $var#> [1] 0.09831393 -0.21314836 0.03053125 0.11506290 0.032713904.2.8 Using non-reversible PTIt is possible to obtain MCMC samples from posterior distributions using thejulia package Pigeons.jl [61], which may have better MCMC performancewith the partially identified model thanks to the non-reversible PT [62].Unlike other functions in this pimbayes package, the number of threads forthe julia session is crucial due to the multithreaded options. The defaultoption of JuliaCall [34], which bridges between julia and R, is a single-thread session. Since it is impossible to change the number of threads in the juliasession once it starts, a restart of R is required if the multithreaded option isrequired but has already been initialized without setting the number of threads.When R starts or restarts, the user can set up a fresh julia session with agiven number of threads.init_julia_session_multi(threads = 4)Next, the user can use sampling() with the desired number of threads. We usethe small number for n rounds to reduce the rendering times of this examplecode. In practice, using a larger number of n rounds is ideal to ensure a sufficientnumber of iterations, 2n.96pt = sampling(jags_code,data,"pigeons", extract_variable_names = c("mu","alpha0","alpha1"),n_chains = 4,n_rounds = 4)Using functions pigeons samples and pigeons plots, it is possible tocheck the results. The following matrix shows the first six rows of the output,which are samples from the posterior distributions.head(pigeons_samples(pt))#> mu alpha0 alpha1#> [1,] -2.0549845 3.289134 -0.43130860#> [2,] -1.3902857 6.629950 -0.01454843#> [3,] -2.1574371 1.462879 -0.82102160#> [4,] -0.6917936 -4.961777 -0.71539678#> [5,] 0.5953992 2.757618 0.14343959#> [6,] 0.7652003 1.575954 0.26909760The software also supports trace and density plots as in Figure 4.3 and Figure 4.4.pigeons_plots(pt)Figure 4.3: Trace and density plots of sampling with Pigeons package resultswithout extended traces option.There are two options for the plot. The default option of Pigeons.jl issaving the samples of target distribution only as in Figure 4.3. However, it is also97Figure 4.4: Trace and density plots of sampling with Pigeons package resultswith extended traces option.possible to save all chains with extended traces option. When the option isenabled, the plot shows each tempered chain in different colours as in Figure 4.4.4.3 Package usage examples with real dataThis section shows how to use this package with an actual data set for a saturatedlogistic regression model with incomplete data (Section 4.2.5) and data from cross-sectional studies (Section 4.2.4).4.3.1 Bank dataResultsWhile the structure of the code for the model is fixed, there are a few elements thatthe user of the package can conveniently change: the name of parameters and priordistributions for the parameters. For example, in the following code, the namesalpha, beta, gamma, delta and their priors could be changed.jags_code <- "model{for (i in 1:n) {x[i, 1:p2] Ë dmulti(alpha[1:p2],1)y[i] Ë dbern(inprod(x[i,1:p2], beta[1:p2]))98r[i] Ë dbern(pr[i])pr[i] <- (y[i]) * inprod(x[i,1:p2], pr_xy1[1:p2]) +abs(y[i]-1) * inprod(x[i,1:p2], gamma[1:p2])}for(i in 1:p2){pr_xy1[i] <- ilogit(delta[i] + logit(gamma[i]))}alpha[1:p2] Ë ddirich(Int64wrapper(ones(p2)))for(i in 1:p2) {beta[i] Ë dunif(0, 1)gamma[i] Ë dunif(0, 1)delta[i] Ë dnorm(0, tau)}}"results = pimbayes::sampling(jags_code = jags_code,data = data, "sat", N = 10000)head(results$samples[,c("alpha[2]","beta[2]","gamma[2]","delta[2]")])alpha[2] beta[2] gamma[2] delta[2][1,] 0.02223439 0.03565986 0.7827169 -0.72752236[2,] 0.02476706 0.02254182 0.8305669 -0.07641678[3,] 0.02390560 0.01900854 0.7778279 0.33015175[4,] 0.02264227 0.02434108 0.8199911 -0.87249738[5,] 0.02554273 0.02111894 0.7694983 0.06607586[6,] 0.02532275 0.03998048 0.7201092 -0.68513174The merit of using this function in the package is that the computational perfor-mance scales reasonably with the size of the dataset. The providers of the data setalso provide a full data set, which is 10 times larger. The execution times are 17secs and 46 secs for the respective datasets, which is not meaningfully differentconsidering that this includes the time to compile with JuliaBUGS in a juliasession.4.3.2 Cross-sectional studyDataset specificationThe data set in Table 4.1 has the same setting as in the mock data of Section 4.2.4but is actually gathered for studies [10, 45]. The following is the JAGS code forthis model. Similarly, the names of se, sp, pi and priors could be changed.However, the names for eta and y canât be changed since they are the names ofdata and interim variables.99D+ Dâ TotalE+ 22 25 47Eâ 82 251 333Total 104 276 380Table 4.1: Leptospirosis study for sheep abattoirs in New Zealand [10, 45].jags_code <-"model{y[1:4] Ë dmulti(eta[1:4], n)eta[1:4] <- c(se*pi[1]+(1-sp)*pi[3],se*pi[2]+(1-sp)*pi[4],sp*pi[3]+(1-se)*pi[1],(1-se)*pi[2]+sp*pi[4])pi[1:4] Ë ddirch(c(1,1,1,1))se Ë dbeta(25,3)sp Ë dbeta(30,1.5)}"We get the results and find the column means, which are posterior means.y = c(22,25,82,251)n = 380data = list(n = n, y = y)results = sampling(jags_code, data, type = "cross", N = 1000)colMeans(results$samples)# pi[11] pi[12] pi[21] se sp# 0.05868982 0.04784622 0.21837113 0.88757340 0.96090403In this section, we have discussed the functionality of the package pimbayes.The package can help users explore suspected partially identified models and per-form Bayesian inference with them. Further developments could include addingmore ISTP for various GLM models and improving computational efficiency.Also, we try to support additional platforms and inputs, such as STAN input,Turing.jl input, and python in addition to R.In the next chapter, we summarize and discuss the overall materials in this the-sis and propose possible future research topics, including pseudo-TP, mathematicalproof, and the main development points of the package.100Chapter 5Conclusions and Discussions5.1 ConclusionIn this thesis, we discuss the computation and parameterization issues whenBayesian inference is used for the models that are known to be only partially identi-fied, and models which are suspected of being either partially or weakly identified.Even though the initial parameterization scientifically makes sense, the computa-tion with it is not efficient since it is not uniquely identified. We argue for theneed for alternative methods; one of these was importance sampling with transpar-ent reparameterization (ISTP) [22]. ISTP uses TP, which reparameterizes basedon parametersâ contribution to the distribution of data. After reparameterization,the algorithm uses the convenience prior, which results in convenience posteriors.With resampling from samples from convenience posteriors, it is possible to finallyget the samples from the desired posterior while having a more efficient computa-tion. The simulation results in Chapter 2 show that ISTP generally works better,especially in terms of computational time. At the same time, the ESS is usuallylarger than that of the MCMC algorithm. It shows a good result for large-scaleproblems, such as those involving high-dimensional parameters.This algorithm is straightforward and easy to apply when a TP is available.However, the hardest part is finding a TP. It is possible that a TP does not ex-ist, and even if it exists, one cannot find a TP. So, in Chapter 3, we discuss thegeneral usability of a TP with a numerical representation of the identification of101the model. We suspect that the degree to which the partially identified model isidentified is related to the existence of a TP. So, we characterize the extent of iden-tification with a condition number and rank of the Fisher information matrix. Weset a rank and condition number threshold for the Fisher information matrix. Usingthis information, we set an algorithm for researchers to help prepare their statisti-cal analysis with Bayesian inference with the suspected partially identified model.The algorithm calculates the rank and condition number of the Fisher informationmatrix and warns about computational difficulty. Lastly, we develop a packagepimbayes that encompasses ISTP and the diagnostic algorithm. We designedthis package to be a comprehensive tool for Bayesian analysis with the partiallyidentified model. For now, the main functionality is restricted for the specific typesof the model, such as a saturated logistic regression with incomplete data, but wealso support the non-reversible parallel tempering for almost any models suspectedof being partially, or perhaps weakly, identified that can be written in a JAGS code.5.2 DiscussionIn this thesis, we discuss the issues with posterior inference for the partially identi-fied model with the off-the-shelf MCMC algorithm, how we numerically substanti-ate the issues, and how we could resolve the issues so that researchers could easilyand efficiently perform analysis with the model in its intended parameterization.The main goal of this thesis is ultimately to improve the accessibility to Bayesianinference with the partially identified model. Hence, we develop a package for re-searchers to easily access ISTP, diagnostic tools for the suspected partially iden-tified model, and non-reversible parallel tempering. In this section, we discuss thelimitations and what can be improved in the future.â¢ Performance with Pseudo-TPA difficulty of using a TP arises from the fact that there is no standard rulefor finding it, and it can be difficult to find it even if it actually exists. Asdiscussed in Chapter 2, a pseudo-TP is a good alternative in terms that itrelaxes the requirements for a TP, but it is close to having the properties of aTP. However, the performance with pseudo-TP in the simulation for Exam-ple 6 was not good enough to use it as a good alternative to posterior infer-102ence with an off-the-shelf MCMC algorithm. We only tried one example, sopseudo-TP might work better in other examples, and it is also possible thatthere is another pseudo-TP that can work better for our example. Also, theway to find pseudo-TP introduced is very specific to this example. If therewere a generic way to find a pseudo-TP that works better, it would greatlyrelieve the challenge of finding a TP.â¢ Numerical experiments versus Mathematical proofIn this thesis, especially in Chapter 3, we have multiple numerical experi-ments. Many of those results are consistent with what we expect and assume.However, while even one numerical counter-example is enough to show amodel is not identified, using numerical simulations with a big Monte Carlosample to show the consistency with the model being identified doesnât con-stitute formal proof. We could only increase the confidence by increasingthe number of numerical simulations. So, in future research, it will be goodto focus on more theoretical and mathematical proofs. A partially identifiedmodel involves many concepts of linear algebra. At the same time, MCMCalgorithms in Bayesian inference and ISTP are directly related to stochasticprocesses. Therefore, it is worth starting from those areas.â¢ Automatically find a TPIf we could find a fully automated method to find a TP, almost all concernswe have about the usability, universality, and versatility of ISTP would beresolved. The main blockades are the following:1. It is not easy to identify which parameters are directly learnable.2. TP is not unique, especially due to the λ part in h(θ) = (Ï ,λ ).3. Even when learnable parameters are available, it is not easy to find afunction h(θ) = (Ï ,λ ) without human intervention.While developing the package pimbayes, we find that the parser inJuliaBUGS is a powerful helper in removing these blocks. In the parser,information about all parameters and variables is provided. Then, all nodesare sorted according to the hierarchy. We expect that making a large matrix103of information from the graph in the parsed model could help extract therelationship between parameters.For the last blockade, we can also consider numerical approximation. Us-ing the iterative method can give a close-enough approximation to the solu-tion. Similarly, the machine may not mathematically and theoretically finda closed-form function h(θ); it may find a function that returns a numericalapproximation of the given h(θ) or hâ1(Ï ,λ ). As long as they are reliableapproximations, numerical values are enough for applying ISTP.The other possible way is using optimization method to find λ and accordingapproximated h(θ). This approach may be based on normalizing flows [44]to keep reparameterizing each step while optimizing the function that willfind λ that doesnât depend on the likelihood using the stochastic gradientdescent (SGD) method [1, 29].â¢ Using non-reversible parallel tempering with different explorer and lan-guageThe package Pigeons.jl supports many different explorers (samplers).However, many of them do not work well with integer parameters. In mostof our examples, we have integer missing entries, which are also parametersin Bayesian inference. Therefore, the default option in Pigeons.jl whenJuliaBUGS.jl model is provided with JAGS code is the slice sampler[42]. There may be more ideal samplers for problems with integer missingentries. We start with the Gibbs sampler from JuliaBUGS.jl [63]. It doesnot yet fully work for Pigeons.jl due to the incompatible type betweenthe two packages. We have a plan to continue to find the most ideal samplersfor using Pigeons.jl on these problems, including resolving the typeissue in the Gibbs sampler.The other way is using other types of input language. As previously dis-cussed, language that supports Markov blanket [32] will help resolve thecomputational issues. Blang [5] is one of the examples. Developing an in-terpreter between JAGS and Blang may also work in this case.104Bibliography[1] S.-i. Amari. Backpropagation and stochastic gradient descent method.Neurocomputing, 5(4-5):185â196, 1993. â page 104[2] A. Beck. Introduction to Nonlinear Optimization: Theory, Algorithms, andApplications with Python and MATLAB, Second Edition. SIAM,Philadelphia, PA, 2023. doi:10.1137/1.9781611977622. â pages 55, 60[3] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah. Julia: A FreshApproach to Numerical Computing. SIAM Review, 59(1):65â98, 2017.doi:10.1137/141000671. â pages 7, 85[4] H. W. Borchers. Pracma: Practical Numerical Math Functions, 2022. URLhttps://CRAN.R-project.org/package=pracma. R package version 2.4.2.Retrieved April 10, 2025. â page 22[5] A. Bouchard-CoËte´, K. Chern, D. Cubranic, S. Hosseini, J. Hume, M. Lepur,Z. Ouyang, and G. Sgarbi. Blang: Bayesian declarative modelling of generaldata structures and inference via algorithms based on distribution continua.arXiv preprint arXiv:1912.10396, 2019. â page 104[6] G. Casella and E. I. George. Explaining the Gibbs sampler. The AmericanStatistician, 46(3):167â174, 1992. â pages 2, 4, 57, 85[7] A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood fromincomplete data via the EM algorithm. Journal of the Royal StatisticalSociety. Series B (Methodological), 39(1):1â22, 1977. â pages 79, 81[8] M. J. Denwood. Runjags: An R package providing interface utilities, modeltemplates, parallel computing methods and additional distributions forMCMC models in JAGS. Journal of Statistical Software, 71(9):1â25, 2016.doi:10.18637/jss.v071.i09. â page 85105[9] A. Doucet, N. De Freitas, N. J. Gordon, et al. Sequential Monte CarloMethods in Practice, volume 1. Springer, 2001. â page 14[10] A. Dreyfus, C. Heuer, P. Wilson, J. Collins-Emerson, M. Baker, andJ. Benschop. Risk of infection and associated influenza-like disease amongabattoir workers due to two Leptospira species. Epidemiology & Infection,143(10):2095â2105, 2015. â pages ix, 99, 100[11] D. J. Earl and M. W. Deem. Parallel tempering: Theory, applications, andnew perspectives. Physical Chemistry Chemical Physics, 7(23):3910â3916,2005. â page 84[12] J. M. Flegal, J. Hughes, D. Vats, N. Dai, K. Gupta, and U. Maji. Mcmcse:Monte Carlo Standard Errors for MCMC, 2021. URLhttps://CRAN.R-project.org/package=mcmcse. R package version 1.5-0.Retrieved April 10, 2025. â pages 15, 25, 33[13] G. Forchini and G. H. Hillier. Ill-conditioned problems, Fisher information,and weak instruments. Technical report, cemmap working paper, 2005. âpage 56[14] N. Friedman. The Bayesian structural EM algorithm. arXiv preprintarXiv:1301.7373, 2013. â pages 79, 81[15] H. Ge, K. Xu, and Z. Ghahramani. Turing: a language for flexibleprobabilistic inference. In International Conference on Artificial Intelligenceand Statistics, AISTATS 2018, 9-11 April 2018, Playa Blanca, Lanzarote,Canary Islands, Spain, pages 1682â1690, 2018. â pages 84, 86[16] A. E. Gelfand and A. F. Smith. Sampling-Based Approaches to CalculatingMarginal Densities. Journal of the American Statistical Association, 85(410):398â409, 1990. â page 57[17] A. Gelman, J. Carlin, H. Stern, D. Dunson, A. Vehtari, and D. Rubin.Bayesian data analysis. Chapman and Hall/CRC, third edition, 2013. âpage 13[18] P. W. Glynn and D. L. Iglehart. Importance Sampling for StochasticSimulations. Management Science, 35(11):1367â1392, 1989. â pages12, 14[19] P. Gustafson. On Model Expansion, Model Contraction, Identifiability andPrior Information: Two Illustrative Scenarios Involving Mismeasured106Variables . Statistical Science, 20:111 â 140, 2005.doi:10.1214/088342305000000098. â pages 73, 75[20] P. Gustafson. What are the limits of posterior distributions arising fromnonidentified models, and why should we care? Journal of the AmericanStatistical Association, 104(488):1682â1695, 2009. â page 6[21] P. Gustafson. Bayesian inference for partially identified models. TheInternational Journal of Biostatistics, 6(2), 2010. â page 8[22] P. Gustafson. Bayesian Inference for Partially Identified Models: Exploringthe Limits of Limited Data. Chapman & Hall, 2015. â pages2, 5, 6, 8, 9, 10, 11, 12, 14, 84, 101[23] P. Gustafson. Summer-school. https://github.com/paulgstf/summer-school,2022. Retrieved April 10, 2025. â pages 9, 16[24] N. J. Higham. Accuracy and Stability of Numerical Algorithms. SIAM,Second edition, 2002. â page 82[25] T. M. Inc. Jacobian, 2024. URLhttps://www.mathworks.com/help/symbolic/sym.jacobian.html. RetrievedApril 10, 2025. â page 22[26] L. T. Johnson and C. J. Geyer. Variable transformation to obtain geometricergodicity in the random-walk Metropolis algorithm. The Annals ofStatistics, pages 3050â3076, 2012. â page 6[27] G. Jones, W. O. Johnson, T. E. Hanson, and R. Christensen. Identifiability ofmodels for multiple diagnostic testing in the absence of a gold standard.Biometrics, 66(3):855â863, 2010. â page 60[28] G. L. Jones. On the Markov chain central limit theorem. ProbabilitySurveys, 1:299 â 320, 2004. doi:10.1214/154957804100000051. â page 15[29] N. Ketkar. Stochastic gradient descent. In Deep learning with Python: Ahands-on introduction, pages 113â132. Springer, 2017. â page 104[30] L. Kish. Survey Sampling. John Wiley & Sons, Inc.,, 1965. â page 13[31] T. Kloek and H. K. Van Dijk. Bayesian estimates of equation systemparameters: an application of integration by Monte Carlo. Econometrica:Journal of the Econometric Society, pages 1â19, 1978. â page 12107[32] D. Koller and N. Friedman. Probabilistic Graphical Models: Principles andTechniques. MIT press, 2009. â pages 86, 104[33] A. Lewbel. The identification zoo: Meanings of identification ineconometrics. Journal of Economic Literature, 57(4):835â903, 2019. âpages 8, 9[34] C. Li. JuliaCall: an R package for seamless integration between R and Julia.The Journal of Open Source Software, 4(35):1284, 2019.doi:10.21105/joss.01284. â pages 90, 96[35] R. J. Little and D. B. Rubin. Statistical Analysis with Missing Data, volume793. John Wiley & Sons, 2019. â pages 4, 10, 37, 80[36] T. A. Louis. Finding the observed information matrix when using the EMalgorithm. Journal of the Royal Statistical Society Series B: StatisticalMethodology, 44(2):226â233, 1982. â page 80[37] A. Ly, M. Marsman, J. Verhagen, R. P. Grasman, and E.-J. Wagenmakers. Atutorial on Fisher information. Journal of Mathematical Psychology, 80:40â55, 2017. â page 56[38] G. J. McLachlan and T. Krishnan. The EM Algorithm and Extensions. JohnWiley & Sons, 2008. â pages 79, 81[39] L. Meng and J. C. Spall. Efficient computation of the Fisher informationmatrix in the EM algorithm. In 2017 51st Annual Conference on InformationSciences and Systems (CISS), pages 1â6. IEEE, 2017. â page 80[40] S. Moro, P. Cortez, and P. Rita. Bank Marketing. UCI Machine LearningRepository, 2014. doi:10.24432/C5K306. â page 30[41] S. Moro, P. Cortez, and P. Rita. A data-driven approach to predict thesuccess of bank telemarketing. Decision Support Systems, 62:22â31, 2014.â page 30[42] R. M. Neal. Slice sampling. The Annals of Statistics, 31(3):705â767, 2003.â page 104[43] R. M. Neal. MCMC using Hamiltonian dynamics. Handbook of MarkovChain Monte Carlo, 2(11):2, 2011. â page 2[44] G. Papamakarios, E. Nalisnick, D. J. Rezende, S. Mohamed, andB. Lakshminarayanan. Normalizing Flows for Probabilistic Modeling and108Inference . Journal of Machine Learning Research, 22(57):1â64, 2021. âpage 104[45] S. Pirikahu, G. Jones, and M. L. Hazelton. Bayesian credible intervals forpopulation attributable risk from caseâcontrol, cohort and cross-sectionalstudies. Australian & New Zealand Journal of Statistics, 63(4):639â657,2021. doi:10.1111/anzs.12352. â pages ix, 2, 4, 6, 48, 90, 91, 99, 100[46] M. Plummer. JAGS - Just Another Gibbs Sampler, 2017. URLhttps://sourceforge.net/projects/mcmc-jags/files/Manuals/4.x/jags user manual.pdf/. JAGS Version 4.3.0 user manual. Retrieved April 10,2025. â pages 2, 5, 57, 61, 81, 85, 90[47] M. Plummer. Rjags: Bayesian Graphical Models using MCMC, 2022. URLhttps://mcmc-jags.sourceforge.io. R package version 4-13. Retrieved April10, 2025. â pages 5, 14, 22, 57, 61, 81[48] M. Plummer, N. Best, K. Cowles, and K. Vines. CODA: ConvergenceDiagnosis and Output Analysis for MCMC. R News, 6(1):7â11, 2006. âpage 14[49] R Core Team. R: A Language and Environment for Statistical Computing. RFoundation for Statistical Computing, Vienna, Austria, 2024. URLhttps://www.R-project.org/. Retrieved April 10, 2025. â pages iv, 7[50] J. J. Rissanen. Fisher Information and Stochastic Complexity. IEEETransactions on Information Theory, 42(1):40â47, 1996. â page 56[51] M.-K. Riviere, S. Ueckert, and F. Mentre´. An MCMC method for theevaluation of the Fisher information matrix for non-linear mixed effectmodels. Biostatistics, 17(4):737â750, 2016. â page 80[52] S. Ross. A First Course in Probability. Pearson, Eighth edition, 2010. âpages 12, 22[53] T. J. Rothenberg. Identification in Parametric Models. Econometrica:Journal of the Econometric Society, pages 577â591, 1971. â page 60[54] R. D. Skeel. Scaling for Numerical Stability in Gaussian Elimination.Journal of the ACM (JACM), 26(3):494â526, 1979. â page 82[55] J. C. Spall. Monte Carlo computation of the Fisher information matrix innonstandard settings. Journal of Computational and Graphical Statistics, 14(4):889â909, 2005. â pages 77, 78, 81109[56] Stan Development Team. RStan: the R interface to Stan, 2023. URLhttps://mc-stan.org/. R package version 2.21.8. Retrieved April 10, 2025. âpage 14[57] Stan Development Team. Stan Modeling Language Users Guide andReference Manual, version 2.36, 2024. URL https://mc-stan.org. RetrievedApril 10, 2025. â pages 2, 84[58] J. H. Stock and J. H. Wright. GMM with weak identification. Econometrica,68(5):1055â1096, 2000. â pages 2, 56[59] J. H. Stock, J. H. Wright, and M. Yogo. A Survey of Weak Instruments andWeak Identification in Generalized Method of Moments. Journal ofBusiness & Economic Statistics, 20(4):518â529, 2002. â page 2[60] P. Stoica and T. L. Marzetta. Parameter estimation problems with singularinformation matrices. IEEE Transactions on Signal Processing, 49(1):87â90, 2001. â page 60[61] N. Surjanovic, M. Biron-Lattes, P. Tiede, S. Syed, T. Campbell, andA. Bouchard-CoËte´. Pigeons. jl: Distributed sampling from intractabledistributions. arXiv preprint arXiv:2308.09769, 2023. â pages7, 84, 85, 89, 96[62] S. Syed, A. Bouchard-CoËte´, G. Deligiannidis, and A. Doucet.Non-reversible parallel tempering: A scalable highly parallel mcmc scheme.Journal of the Royal Statistical Society Series B: Statistical Methodology, 84(2):321â350, dec 2021. doi:10.1111/rssb.12464. â pages 82, 84, 85, 89, 96[63] TuringLang. JuliaBUGS.jl. https://github.com/TuringLang/JuliaBUGS.jl,2021. Retrieved April 10, 2025. â pages 7, 90, 104[64] A. Van der Sluis. Condition numbers and equilibration of matrices.Numerische Mathematik, 14(1):14â23, 1969. â page 55[65] D. Vats, J. M. Flegal, and G. L. Jones. Multivariate output analysis forMarkov chain Monte Carlo. Biometrika, 106(2):321â337, 2019. â pages15, 25, 33110Appendix ASupporting MaterialsPackage pimbayes is available on GitHub.111
Thesis/Dissertation
eng
Statistics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Gustafson, Paul, 1968-
Statistical models that are known or suspected to be partially identified : issues of parameterization, computation, and software development
Text