Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Linear clustering with application to single nucleotide polymorphism genotyping Yan, Guohua 2008

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
ubc_2008_fall_yan_guohua.pdf [ 3.13MB ]
Metadata
JSON: 1.0066454.json
JSON-LD: 1.0066454+ld.json
RDF/XML (Pretty): 1.0066454.xml
RDF/JSON: 1.0066454+rdf.json
Turtle: 1.0066454+rdf-turtle.txt
N-Triples: 1.0066454+rdf-ntriples.txt
Original Record: 1.0066454 +original-record.json
Full Text
1.0066454.txt
Citation
1.0066454.ris

Full Text

Linear Clustering with Application to Single Nucleotide Polymorphism Genotyping by Guohua Yan B.Sc., Liaocheng University, 1992 M.Sc., Beijing Normal University, 1995 M.Sc., The University of Windsor, 2003 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Statistics)  The University of British Columbia (Vancouver) June, 2008 c Guohua Yan 2008  Abstract Single nucleotide polymorphisms (SNPs) have been increasingly popular for a wide range of genetic studies. A high-throughput genotyping technologies usually involves a statistical genotype calling algorithm. Most calling algorithms in the literature, using methods such as k-means and mixturemodels, rely on elliptical structures of the genotyping data; they may fail when the minor allele homozygous cluster is small or absent, or when the data have extreme tails or linear patterns. We propose an automatic genotype calling algorithm by further developing a linear grouping algorithm (Van Aelst et al., 2006). The proposed algorithm clusters unnormalized data points around lines as against around centroids. In addition, we associate a quality value, silhouette width, with each DNA sample and a whole plate as well. This algorithm shows promise for genotyping data generated from TaqMan technology (Applied Biosystems). A key feature of the proposed algorithm is that it applies to unnormalized fluorescent signals when the TaqMan SNP assay is used. The algorithm could also be potentially adapted to other fluorescence-based SNP genotyping technologies such as Invader Assay. Motivated by the SNP genotyping problem, we propose a partial likelihood approach to linear clustering which explores potential linear clusters in a data set. Instead of fully modelling the data, we assume only the signed orthogonal distance from each data point to a hyperplane is normally distributed. Its relationships with several existing clustering methods are discussed. Some existing methods to determine the number of components in a data set are adapted to this linear clustering setting. Several simulated and real data sets are analyzed for comparison and illustration purpose. We also investigate some asymptotic properties of the partial likelihood approach. A Bayesian version of this methodology is helpful if some clusters are sparse but there is strong prior information about their approximate locations or properties. We propose a Bayesian hierarchical approach which is particularly appropriate for identifying sparse linear clusters. We show that the sparse cluster in SNP genotyping datasets can be successfully identified after a careful specification of the prior distributions. ii  Table of Contents Abstract  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  List of Tables  vi  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . .  x  Dedication  xi  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  Statement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 SNP genotyping . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 SNPs and their applications . . . . . . . . . . . . 1.1.2 TaqMan SNP genotyping technology . . . . . . . 1.2 Review of several clustering algorithms . . . . . . . . . . 1.2.1 MCLUST: normal mixture model-based clustering 1.2.2 MIXREG: mixture of linear regressions . . . . . . 1.2.3 LGA: linear grouping algorithm . . . . . . . . . . 1.3 Outline of the thesis . . . . . . . . . . . . . . . . . . . . .  . . . . . . . . .  1 1 1 2 3 3 5 5 6  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  8  2 Automatic SNP genotype calling . . . . . . . . . . 2.1 Background . . . . . . . . . . . . . . . . . . . . . 2.1.1 TaqMan SNP genotyping . . . . . . . . . . 2.1.2 Review of some genotype calling algorithms 2.1.3 ROX-normalized versus unnormalized data 2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Conclusions . . . . . . . . . . . . . . . . . . . . . .  . . . . . . . . . . . . . . . . . . . . . . . . . . .  . . . . . . . . .  . . . . . . .  . . . . . . .  10 10 10 12 13 17 17 iii  Table of Contents 2.4  . . . . . .  18 18 18 19 20 20  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  22  3 A partial likelihood approach to linear clustering . . . . . 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Partial likelihood-based objective function for linear clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Partial likelihood for orthogonal regression . . . . . . 3.2.2 Partial likelihood for linear clustering . . . . . . . . . 3.3 The EM algorithm . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Asymptotic properties . . . . . . . . . . . . . . . . . . . . . . 3.5 Relationships with other clustering methods . . . . . . . . . 3.5.1 With LGA . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 With normal mixture models . . . . . . . . . . . . . . 3.5.3 With mixture of ordinary regressions . . . . . . . . . 3.6 Choosing the number of linear clusters . . . . . . . . . . . . 3.6.1 Bootstrapping the partial likelihood ratio . . . . . . . 3.6.2 Information criteria . . . . . . . . . . . . . . . . . . . 3.7 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7.1 Simulated data I . . . . . . . . . . . . . . . . . . . . . 3.7.2 Simulated data II . . . . . . . . . . . . . . . . . . . . 3.7.3 Australia rock crab data . . . . . . . . . . . . . . . . 3.7.4 Taqman single nucleotide polymorphism genotyping data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . .  24 24 25 25 26 28 30 32 32 33 33 33 33 34 35 35 37 37  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  47  4 Bayesian linear clustering . . . . . . . . . 4.1 Introduction . . . . . . . . . . . . . . . 4.2 Model specification . . . . . . . . . . . 4.2.1 Model for one cluster: orthogonal 4.2.2 Model for linear clustering . . . 4.3 Sampling algorithms . . . . . . . . . . .  51 51 55 55 55 57  2.5  Methods . . . . . . . . . . . . . . . . . . . . . 2.4.1 Data preprocessing . . . . . . . . . . . 2.4.2 Fitting lines using LGA . . . . . . . . . 2.4.3 Genotype assignment . . . . . . . . . . 2.4.4 Quality assessment of DNA samples and Discussion . . . . . . . . . . . . . . . . . . . .  . . . . . . . . . . . . . . . . plates . . . .  . . . . . . . . . . . . . . . . . . regression . . . . . . . . . . . .  . . . . . .  . . . . . .  . . . . . . . . . . . . . . . . . . . . . . .  . . . . . .  . . . . . .  . . . . . .  42 45  iv  Table of Contents . . . . .  57 60 62 64 76  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  78  5 Consistency and asymptotic normality . . 5.1 Introduction . . . . . . . . . . . . . . . . . 5.2 Population version of the objective function 5.3 Consistency . . . . . . . . . . . . . . . . . 5.4 Asymptotical normality . . . . . . . . . . .  . . . . .  80 80 82 85 94  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  96  6 Future Work . . . . . . . . . 6.1 Robustness consideration 6.2 Asymptotics . . . . . . . 6.3 Model Extension . . . . . 6.4 Variable/model selection 6.5 Bayesian computation . .  . . . . . .  97 97 97 98 98 98  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  99  4.4 4.5 4.6  4.3.1 Gibbs sampling . . . . . . . . . . . . . . . . . . . 4.3.2 Metropolis-Hastings sampling . . . . . . . . . . . Australia rock crab data . . . . . . . . . . . . . . . . . . Taqman single nucleotide polymorphism genotyping data Discussion . . . . . . . . . . . . . . . . . . . . . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . .  . . . . . .  . . . . .  . . . . . .  . . . . .  . . . . . .  . . . . .  . . . . . .  . . . . .  . . . . . .  . . . . .  . . . . . .  . . . . .  . . . . . .  . . . . .  . . . . . .  . . . . .  . . . . .  . . . . . .  Appendix A Glossary of some genetic terms . . . . . . . . . . . . . . . . . 100  v  List of Tables 3.1  3.2 3.3  3.4 3.5  Criteria for choosing the numbers of clusters, K, when MLC is applied to simulated data I. “Boot p-value” is the p-value using the bootstrapping partial likelihood ratio test for H0 : K = K0 vs H1 : K = K0 + 1 where 99 bootstrapping samples are used . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Misclassification matrices of simulated data II from MLC, MCLUST, LGA and MIXREG. . . . . . . . . . . . . . . . . . Criteria for choosing the number of clusters, K, when MLC is applied to the blue crab data. “Boot p-value” is the pvalue using the bootstrapping partial likelihood ratio test for H0 : K = K0 vs H1 : K = K0 + 1 where 99 bootstrapping samples are used. . . . . . . . . . . . . . . . . . . . . . . . . . Misclassification matrices of the blue crab data (using log scales) from MLC, MCLUST, LGA and MIXREG. . . . . . . Largest membership probabilities of the 7 points labelled in Figure 3.5 by MCLUST and MLC. . . . . . . . . . . . . . . .  37 38  39 40 43  vi  List of Figures 2.1  2.2 2.3 2.4  2.5  2.6  3.1  3.2 3.3 3.4  3.5  Scatterplots of ROX-normalized data for (a): a good plate; (b) & (c): messy plates; (d): a plate with only a few points in the YY cluster; (e): a plate with no points in the YY cluster; (f): a plate with only one cluster. . . . . . . . . . . . . . . . . Scatterplots of unnormalized and ROX-normalized data for a plate. Both points and letters represent samples. . . . . . . . Probability of calling a genotype versus ROX . . . . . . . . . Silhouette width versus ROX value. Left panel: k-means results using ROX-normalized data; right panel: LGA result using unnormalized data. . . . . . . . . . . . . . . . . . . . . Control points in scatterplots of ROX-normalized and unnormalized data for a plate. Letter “P” represents a positive control point and ”N” a negative control point. . . . . . . . Left panel: the calling result of SDS software displayed in the scatterplot of unnormalized data. Right panel: the scoring results of LGA using unnormalized data. The symbol ◦ represents noncalls. . . . . . . . . . . . . . . . . . . . . . . . Comparison of clustering results of simulated data I. Upperleft panel displays the true classification; upper-right, lowerleft and lower-right are that of MLC, LGA and MCLUST. . . Pairwise scatterplots of simulated data II. . . . . . . . . . . . Pairwise scatterplots of the blue crab data. . . . . . . . . . . Clustering results of the blue crab data (using log scales) from MLC, MCLUST, LGA and MIXREG. CL is used as the response in MIXREG. . . . . . . . . . . . . . . . . . . . . . . . Clustering results from four methods applied to a Taqman data set plate. ◦, + and denote the labels assigned for the three linear clusters (genotypes) by each method, × denotes points assigned as negative controls and denotes points assigned to the background cluster. . . . . . . . . . . . . . . . .  11 14 15  16  16  17  36 38 40  41  44  vii  List of Figures 4.1 4.2  4.3 4.4 4.5  4.6 4.7  4.8  4.9  4.10  4.11  4.12 4.13 4.14  Scatterplot of one plate in which the variant allele homozygous cluster has only five points (upper-left). . . . . . . . . . Clustering results of several clustering algorithms. For kmeans and MCLUST, the number of clusters are set to 4; for the remaining algorithms, the number of lines are set to 3. . . Pairwise scatter plots of the blue crab data. . . . . . . . . . . Evolution of sampled values for θ for 11, 000 iterations using Algorithm 1 for blue crab data. . . . . . . . . . . . . . . . . . Evolution of log unnormalized posterior density of sampled values for (θ, z1 , . . ., zn ) for 11, 000 iterations using Algorithm 1 for blue crab data. . . . . . . . . . . . . . . . . . . . . . . . Autocorrelations of sampled values for θ for 11, 000 iterations using Algorithm 1 for blue crab data. . . . . . . . . . . . . . . Density curves of sampled values for θ for 10, 000 iterations using Algorithm 1, after a burn-in of 1000 iterations for blue crab data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Unnormalized marginal posterior density values of θ along ˆ1 + θ ˆ 2 connecting the two estimated a line segment (1 − t)θ ˆ ˆ modes θ 1 and θ 2 . One mode is estimated by maximum a posteriori estimation from 10, 000 iterations using Algorithm 1, after a burn-in of 1000 iterations and the other is obtained by permutation. The blue crab data is used. . . . . . . . . . . Estimated posterior probability of being classified into Class 1 of the 100 crabs from 10, 000 iterations using Algorithm 1, after a burn-in of 1000 iterations. The solid vertical line separates these crabs by their true classes and the horizontal dashed line corresponds to a probability of 0.5. . . . . . . . . Evolution of log unnormalized posterior density of sampled values for η for 11, 000 iterations using Algorithm 2 for the SNP genotyping data. . . . . . . . . . . . . . . . . . . . . . . Evolution of sampled values for θ for 11, 000 iterations using Algorithm 2 for the SNP genotyping data. The first three panels are for slopes. . . . . . . . . . . . . . . . . . . . . . . . Autocorrelations of sampled values for θ for 11, 000 iterations using Algorithm 2 for the SNP genotyping data. . . . . . . . Density curves of sampled values for θ for 11, 000 iterations using Algorithm 2 for the SNP genotyping data. . . . . . . . Clustering results of the SNP genotyping data in which points are classified into the cluster with the largest posterior membership probability. . . . . . . . . . . . . . . . . . . . . . . . .  53  54 63 65  66 67  68  69  69  72  73 74 75  76 viii  List of Figures 4.15 Clustering results of the Bayesian approach for a plate without a sparse cluster. . . . . . . . . . . . . . . . . . . . . . . .  77  ix  Acknowledgements Foremost I would like to thank my supervisors Dr. William J. Welch and Dr. Ruben H. Zamar, for their excellent guidance, continuous encouragement, great patience and generous support. I am indebted to them for many things I have learned which lay the foundation of my future career. I am very grateful to other members of my supervisory committee, Dr. Paul Gustafson and Dr. Arnaud Doucet for their valuable comments and suggestions. I am also very grateful to Dr. Harry Joe and Dr. Lang Wu for their great help and valuable advice both in statistics and in my job application process. Also I would like to thank them for teaching me playing table tennis although they must have been frustrated that I am a hopeless learner. Many thanks to Dr. Mat´ıas Salibi´an-Barrera, Dr. Raphael Gottardo, Dr. Nancy Heckman and Dr. Jiahua Chen for their help and valuable suggestions. I also thank Dr. Jason Loeppky for his help on my proposal. Special thanks to Ms. Loubna Akhabir and Ms. Treena McDonald for providing me the data and explaining to me the genomic background. I would also like to thank Ms. Elaine Salameh, Ms. Christine Graham, Ms. Rhoda Morgan, Ms. Peggy Ng and Ms. Viena Tran for their help with administrative matters. Furthermore, I thank my fellow students for their help throughout my program. I choose not to list their names here but their help and friendship will be remembered forever. Finally, I would like to acknowledge the support of the University of British Columbia through a University Graduate Fellowship.  x  To my parents, my wife and my daughter  xi  Statement of Co-Authorship This thesis is completed under the supervision of Dr. William J. Welch and Dr. Ruben H. Zamar. Chapter 2 is co-authored with Dr. William J. Welch, Dr. Ruben H. Zamar, Ms. Loubna Akhabir and Ms. Treena McDonald. I conducted the data analysis and prepared a draft of the manuscript. Chapter 3 is co-authored with Dr. William J. Welch and Dr. Ruben H. Zamar. My main contribution is the derivation of the partial likelihood. I conducted the data analysis and prepared a draft of the manuscript. Chapter 4 is co-authored with Dr. William J. Welch and Dr. Ruben H. Zamar. My main contributions are the formulation of the model and the specification of the prior distribution. I conducted the data analysis and prepared the manuscript. Chapter 5 is co-authored with Dr. William J. Welch and Dr. Ruben H. Zamar. I derived the asymptotic properties and prepared the manuscript.  xii  Chapter 1  Introduction This thesis work investigates methods to detect linearly shaped clusters in a dataset. It is motivated by a clustering problem in SNP (single nucleotide polymorphism) genotyping and much effort of this thesis has been devoted to automatic genotype calling algorithms in TaqMan SNP genotyping technology. In this chapter, we introduce the SNP genotyping problem, review several clustering algorithms used in the rest of the thesis and give an outline of the thesis.  1.1 1.1.1  SNP genotyping SNPs and their applications  A single nucleotide polymorphism (SNP, pronounced as “snip”) is a singlebase variation in a genome. The genetic code is specified by the four nucleotide “letters”: A (adenine), C (cytosine), T (thymine) and G (guanine). There are two complimentary DNA strands. It is sufficient to consider only one. SNP variation occurs when a single nucleotide, such as an A, is replaced by one of the other three letters C, G or T at a particular base of the target strand. An example of a SNP is the alteration of the DNA segment AAGGTTA to ATGGTTA, where the second A in the first snippet is replaced with a T (www.ncbi.nlm.nih.gov). A SNP variation usually involves only two nucleotides; the two possible nucleotides are called two SNP alleles. Throughout we shall use “X” to generically represent the common wild-type allele and “Y” the variant allele. The DNA sequences of any two individuals are mostly identical and SNPs are found about every 250 − 350 base pairs in the human genome (Niu et al., 2002). Approximately 3 to 5 percent of a person’s DNA sequence codes for the production of proteins, most SNPs are found outside of these “coding sequences”. SNPs found within a coding sequence, cSNPs, are of particular interest to researchers because they are more likely to alter the biological function of a protein (www.ncbi.nlm.nih.gov). Moreover, the abundance of SNPs makes them useful markers for genetic association studies that work 1  Chapter 1. Introduction to localize the genes involved in complex diseases or adverse drug reactions. The popularity of SNPs is also due to their usual biallelic property which makes them amenable to automated genotyping. (Ranade et al., 2001).  1.1.2  TaqMan SNP genotyping technology  Determination of the alleles at a single nucleotide polymorphism site is called genotyping. The nature of large-scale association studies requires rapid, reliable and cost-effective SNP genotyping. Various SNP genotyping technologies have been developed. The choice of a technology depends on whether a few SNPs are to be typed in many individuals or many different SNPs are to be examined in a few individuals (Callegaro et al., 2006). The TaqMan SNP Genotyping Assay (Applied Biosystems) is a widely used fluorescence-based high-throughput genotyping technology suitable for the former case. In this method, the region flanking the SNP site of interest is amplified in the presence of two probes each specific for one or the other allele. Probes haves a fluor, called “reporter” at one end but do not fluoresce when free in solution because they have a “quencher” at the other end that absorbs fluorescence from the reporter. During the amplification in which many copies of the same sequence are produced, the probe specifically basepaired with the target is unwound, its reporter liberated from the quencher and the fluorescence is increased. The presence of two probes, each labelled with a different fluor, allows the detection of both alleles in a single tube (De La Vega et al., 2005; Ranade et al., 2001). In the TaqMan SNP genotyping technology, DNA samples of many individuals are arranged in a 96- or 384-well plate and they are amplified simultaneously. For each individual, two quantitative fluorescent signals are measured at the end of amplification for the two alternative SNP alleles, indicating their presence or absence. The pair of signals for an individual forms a point in a scatterplot. Ideally there are four distinct clusters in a scatterplot. The NTC (no template control) cluster lies in the lower-left corner, close to the origin, containing negative control cases (a few cases that do not have DNA samples) and samples that fail to amplify. In the lower-right, upper-left, and upper-right corners are three clusters presumably containing samples of wild-type homozygotes, variant homozygotes, and heterozygotes, respectively. Genotype calls for individual samples are made by a clustering algorithm in the propriety Sequence Detection Software (Applied Biosystems). However, considerable manual intervention of an expert operator is required to assess the data quality, to set fluorescent signal thresholds and to decide the 2  Chapter 1. Introduction genotypes, especially when the variant allele Y is rare. The accuracy of SNP genotype calls is critical to later studies. Even the slightest amount of genotyping errors can lead to serious consequences on haplotype analysis, linkage analysis, genetics distance estimation and background linkage disequilibrium estimation (Kang et al., 2004). For a brief review of clustering algorithms in fluorescence-based SNP genotyping, see Chapter 2. One aim of this thesis to develop a reliable automated SNP genotype calling algorithm in the TaqMan SNP genotyping technology that involves minimal manual intervention and should also work well even the variant allele homozygous cluster YY is sparse. We conclude this subsection by noting that any genotyping calling algorithm proposed in the TaqMan SNP genotyping technology should be applicable to other fluorescence-based genotyping method, such as Invader Assay (Mein et al., 2000). In next section, we review several competing clustering algorithms used in the rest of the thesis.  1.2 1.2.1  Review of several clustering algorithms MCLUST: normal mixture model-based clustering  Finite mixture model has long been proposed for clustering problems. See for example, Banfield and Raftery (1993), Fraley and Raftery (1998), Fraley and Raftery (2002), MacLachlan and Peel (2000). In this approach, data are assumed arising from a mixture of probability distributions; each component of the mixture is regarded as a cluster. This model-based method has some advantage over a heuristic approach. First of all, it is somehow data dependent via a variety of parameter restrictions coupled with a model selection mechanism, for example in the packages MCLUST (Fraley and AE, 2006) and EMMIX (McLachlan et al., 1999). Second, some early proposed heuristic clustering criteria, including the most widely used kmeans method, were later formulated in model frameworks. A statistical model helps people understand for what data sets a particular clustering algorithm is likely to work well. Third, within a model framework, many statistical procedures are readily applicable. For example, maximum likelihood is naturally used as an optimizing criterion; the estimated probabilities of a data point conforming to the components is an appealing measure of uncertainty of the resulting classification; a component in a mixture model, which has a clear meaning in the assumed model, can be interpreted as a cluster; the determination of clustering methods and the number of clusters is then a model selection problem. 3  Chapter 1. Introduction MCLUST is a contributed R package for normal mixture modeling and model-based clustering. It provides a number of functions for normal mixture model-bases clustering, classification likelihood-based hierarchical clustering, normal mixture model-based density estimation and discriminant analysis etc. In this subsection, we review only the normal mixture modelbased clustering which is relevant to the thesis. Given observations x1 , . . ., xn , MCLUST assumes a normal mixture model K  p(xi |θ) =  k=1  pk φk (xi |µk , Σk ), i = 1, . . . , n,  where K is the number of components, pk is the mixing proportion, pk > 0, K k=1 pk = 1, φ(·|µk , Σk ) is the multivariate normal density with mean vector µk and covariance matrix Σk and θ is the collection of all parameters. Banfield and Raftery (1993) proposed a general framework for geometric cross-cluster constraint by parameterizing covariance matrices through eigenvalue decomposition Σk = λk Dk Ak DkT , where Dk is the orthogonal matrix of eigenvectors, Ak is a diagonal matrix whose elements are proportional to the eigenvalues, and λk is an associated constant of proportionality. The orientation of principal components of Σk is determined by Dk , while Ak determines the shape of the density contours; λk specifies the volume of the corresponding ellipsoid. Characteristics (orientation, shape and volume) of distributions are usually estimated from the data and can be allowed to vary between clusters or constrained to be the same for all clusters. This parameterization is very flexible; it includes but is not restricted to earlier proposals such as equal-volume spherical variance (Σk = λI) which has a close relationship with the k-means method, constant variance (Σk = Σ) and the most general unconstrained variance. The EM algorithm is used for maximum likelihood estimation; details are omitted. All models corresponding to the above parameterization possibilities are usually tried. Selection of the number K of components/clusters as well as models is through the Bayesian Information Criterion (BIC), which relates to the Bayes factor. The BIC has the form BIC = 2l − log(n)v, where l is the log-likelihood of the model and v is the number of independent parameters in the model. 4  Chapter 1. Introduction For noisy data, MCLUST adds a first order Poisson process to the normal mixture model, K p0 + p(xi |θ) = pk φk (xi |µk , Σk ), V k=1  where V is the hyper-volume of the data region, p0 > 0 and  1.2.2  K k=0 pk  = 1.  MIXREG: mixture of linear regressions  Mixtures of linear regressions may be regarded as a special case of mixture models. The response of interest is assumed univariate. Let (y1 , x1 ), . . ., (yn , xn ) be the observations. The mixture of regression models is K  p(yi |θ, xi ) =  k=1  pk φ(yi |xTi β k , σk2 ), i = 1, . . . , n,  where K is the number of components, pk is a mixing proportion, pk > 0, K 2 T k=1 pk = 1, β k is the vector of regression coefficients, φ(·|xi β k , σk ) is the T 2 univariate normal density with mean xi β k and variance σk . For ease of notation, we assume that an intercept is already included in β k if necessary. The EM algorithm can be used for maximum likelihood estimation. MIXREG is also a contributed R package (Turner, 2006); the computational details are in Turner (2000). The package allows for the selection of equal variances (σk2 = σ 2 ) or unequal variances. Model selection can be performed using BIC outside the package MIXREG.  1.2.3  LGA: linear grouping algorithm  Van Aelst et al. (2006) proposed a linear grouping algorithm to detect linear patterns in a dataset. It combines ideas of k-means, orthogonal regression and resampling and can uncover linear patterns when most traditional algorithms do not work well. LGA is the corresponding contributed R package (Harrington, 2007). Suppose that we are to uncover k linear groups in a data set with n points in d dimensions. The linear grouping algorithm works as follows: 1. Initialization. Starting values are generated by randomly selecting k mutually exclusive subsets of d points (d-subsets). For each of these d-subsets, a hyperplane through the d points is computed. 5  Chapter 1. Introduction 2. Updating group assignment. Each point is assigned to its closest hyperplane in the sense of orthogonal distance; each hyperplane is recomputed from the updated grouping using orthogonal regression. 3. Iterative refinement. The objective function is set to be the aggregated sum of the squared orthogonal distances of the data points from their closest hyperplane. Repeat step 2 until no significant improvement is attained. A moderate number of iterations (say, 10) is usually sufficient. 4. Resampling. Repeat steps 1 − 3 a number of times (say, 100) and select the solution with the lowest aggregated sum of squared orthogonal distances. This solution is refined further as in step 3 until no improvement is achieved. The idea of silhouette width (Rousseeuw, 1987) is adapted to the linear grouping scenario to measure the strength of the group assignment. Denote by s1 (i) and s2 (i) the orthogonal distances of point i from its closest and second closest hyperplanes, respectively. The silhouette width for point i is defined as s1 (i) . w(i) = 1 − s2 (i) The GAP statistic (Tibshirani et al., 2001) is used to determine the number of linear groups in a data set.  1.3  Outline of the thesis  In Chapter 2, we propose a genotype calling algorithm by further adapting the linear grouping algorithm (Van Aelst et al., 2006). A key feature of this algorithm is that it applies to unnormalized fluorescent signals. Whereas there is no explicit statistical model in Chapter 2, a partial likelihood approach to linear clustering is proposed in Chapter 3. It borrows ideas from normal mixture models and mixtures of regressions and provides an extension to the heuristic linear grouping algorithm. The SNP genotyping problem is revisited in this model framework. The method in Chapter 3 is more flexible than that in Chapter 2, but still cannot handle sparse (homozygous variant) clusters well. In Chapter 4, we move to a Bayesian hierarchical approach to the SNP genotyping problem. With careful specification of the prior structures, the Bayesian approach is able to handle a sparse variant allele homozygous cluster. 6  Chapter 1. Introduction Chapter 5 includes the technical details of the asymptotic results for the partial likelihood approach in Chapter 3. Chapter 6 summarizes a few possible future research directions.  7  Bibliography Banfield, J. and A. Raftery (1993). Model-Based Gaussian and NonGaussian Clustering. Biometrics 49 (3), 803–821. Callegaro, A., R. Spinelli, L. Beltrame, S. Bicciato, L. Caristina, S. Censuales, G. De Bellis, and C. Battaglia (2006). Algorithm for automatic genotype calling of single nucleotide polymorphisms using the full course of TaqMan real-time data. Nucleic Acids Research 34 (7), e56. De La Vega, F., K. Lazaruk, M. Rhodes, and M. Wenz (2005). Assessment of two flexible and compatible SNP genotyping platforms: TaqMan snp genotyping assays and the SNPlex genotyping system. Mutation research. Fundamental and molecular mechanisms of mutagenesis 573 (1-2), 111–135. Fraley, C. and R. AE (2006). mclust Version 3 for R: Normal Mixture Modeling and Model-based Clustering. Technical report, Technical Report 504, University of Washington, Department of Statistics. Fraley, C. and A. Raftery (1998). How Many Clusters? Which Clustering Method? Answers Via Model-Based Cluster Analysis. The Computer Journal 41 (8), 578. Fraley, C. and A. Raftery (2002). Model-Based Clustering, Discriminant Analysis, and Density Estimation. Journal of the American Statistical Association 97 (458), 611–632. Harrington, J. (2007). lga: Tools for linear grouping analysis (LGA). R package version 1.0-0. Kang, H., Z. Qin, T. Niu, and J. Liu (2004). Incorporating Genotyping Uncertainty in Haplotype Inference for Single-Nucleotide Polymorphisms. The American Journal of Human Genetics 74 (3), 495–510. MacLachlan, G. and D. Peel (2000). Finite mixture models. J. Wiley.  8  Bibliography McLachlan, G., D. Peel, K. Basford, and P. Adams (1999). The EMMIX software for the fitting of mixtures of normal and t-components. Journal of Statistical Software 4 (2). Mein, C., B. Barratt, M. Dunn, T. Siegmund, A. Smith, L. Esposito, S. Nutland, H. Stevens, A. Wilson, M. Phillips, et al. (2000). Evaluation of Single Nucleotide Polymorphism Typing with Invader on PCR Amplicons and Its Automation. Genome Research 10 (3), 330–343. Niu, T., Z. Qin, X. Xu, and J. Liu (2002). Bayesian Haplotype Inference for Multiple Linked Single-Nucleotide Polymorphisms. Am. J. Hum. Genet 70, 157–169. Ranade, K., M. Chang, C. Ting, D. Pei, C. Hsiao, M. Olivier, R. Pesich, J. Hebert, Y. Chen, V. Dzau, et al. (2001). High-Throughput Genotyping with Single Nucleotide Polymorphisms. Genome Research 11 (7), 1262– 1268. Rousseeuw, P. (1987). Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics 20, 53–65. Tibshirani, R., G. Walther, and T. Hastie (2001). Estimating the number of clusters in a data set via the gap statistic. Journal of the Royal Statistical Society. Series B (Statistical Methodology) 63 (2), 411–423. Turner, T. (2000). Estimating the propagation rate of a viral infection of potato plants via mixtures of regressions. Journal of the Royal Statistical Society Series C(Applied Statistics) 49 (3), 371–384. Turner, T. (2006). mixreg: Functions to fit mixtures of regressions. R package version 0.0-2. Van Aelst, S., X. Wang, R. Zamar, and R. Zhu (2006). Linear grouping using orthogonal regression. Computational Statistics and Data Analysis 50 (5), 1287–1312.  9  Chapter 2  Automatic genotype calling of single nucleotide polymorphisms using a linear grouping algorithm 2.1 2.1.1  Background TaqMan SNP genotyping  Single nucleotide polymorphisms (SNPs) make up about 90% of all human genetic variations. They occur approximately every 100 to 300 bases along the 3.2-billion-base human genome. They have been investigated as popular biological markers for a wide range of genetic studies. Such studies require reliable and cost-effective high-throughput genotyping technologies. Various SNP genotyping technologies have been developed and commercialized. Shi (2001) provides a review of some SNP genotyping technologies. De La Vega et al. (2005) give an elaborate assessment of two popular SNP genotyping technologies, TaqMan SNP Genotyping Assay and the SNPlex Genotyping System (Applied Biosystems). Kang et al. (2004) give a brief overview of three widely used high-throughput technologies, the TaqMan SNP Genotyping Assay, the OLA (Oligonucleotide Ligation Assay) and the MassARRAY system. Although they employ different mechanisms, the most popular high-throughput technologies share the same conceptual framework: for each DNA sample, two quantitative signal intensities are measured after amplification for two alternative SNP alleles, indicating their presence or absence; these pairs of signal intensities are used to call genotypes by an exA version of this chapter will be submitted for publication. Authors: Guohua Yan, William J. Welch, Ruben H. Zamar, Loubna Akhabir and Treena McDonald.  10  Chapter 2. Automatic SNP genotype calling pert manually or by a clustering algorithm. The two alleles are generically called X and Y hereafter. 3.5  (a)  (b)  (c)  Allele.Y.Rn 2.0 2.5  5  4  3.0  6  5  (YY)  0.5  Allele.Y.Rn 4 3  1  2  1.0  2  1.5  Allele.Y.Rn 3  (XY)  (XX)  (NTC) 0.0  0.5  1.0 1.5 Allele.X.Rn  2.0  2.5  0.5  1.0  2.5  3.0  3.5  3.0  3.5  1.4 Allele.Y.Rn 0.8 1.0  2.5 Allele.Y.Rn 2.0  0.6  1.5  0.4  1.0 2.5 3.0 Allele.X.Rn  2.0 2.5 Allele.X.Rn  1.2  4.0 3.5 Allele.Y.Rn 2.5 3.0 2.0  2.0  1.5  (f)  1.5 1.0  1.5  1.0  (e) 3.0  (d)  1.5 2.0 Allele.X.Rn  2.0  2.5 3.0 Allele.X.Rn  3.5  0.5  1.0  1.5 2.0 Allele.X.Rn  2.5  3.0  Figure 2.1: Scatterplots of ROX-normalized data for (a): a good plate; (b) & (c): messy plates; (d): a plate with only a few points in the YY cluster; (e): a plate with no points in the YY cluster; (f): a plate with only one cluster. Figure 2.1 illustrates several scenarios that are typical in genotyping data from the TaqMan SNP Genotyping Assays (Applied Biosystems) at the James Hogg iCAPTURE Centre. In this technology, two alleles are labelled by fluorescent dyes VIC and FAM; in addition, a third dye ROX, which is assumed unchanged during the amplification, is used to normalize the data. The proprietary ABI PRISM r 7900HT Sequence Detection System Plate Utility Software (referred to as SDS system hereafter) actually uses VIC/ROX and FAM/ROX (referred to as ROX-normalized data hereafter) to make genotype calls. In Figure 2.1, ROX-normalized signals are shown. Ideally, there are four distinct clusters in a scatterplot as in Figure 2.1 (a). The NTC (no template control) cluster lies in the lower-left corner, close to the origin, containing negative control cases (a few cases that do not have DNA samples) and samples that fail to amplify. In the lowerright, upper-left and upper-right corners are three clusters labelled XX, YY 11  Chapter 2. Automatic SNP genotype calling and XY, presumably containing samples of wild-type homozygotes, variant homozygotes, and heterozygotes, respectively. Owing to various artifacts, however, segregation can be poor with points lying between clusters, which make the calls less trustworthy (Figure 2.1 (b) and (c)). Furthermore, in some plates there are only a few points or even no points in the variant allele homozygous cluster YY, which often causes classical clustering algorithms to fail or to make genotype calls incorrectly (Figure 2.1 (d) and (e)). To our knowledge, the proprietary clustering algorithm incorporated in the SDS system cannot make genotyping calls in this situation and one has to call manually. In extreme cases, there is only one cluster present (Figure 2.1 (f)), which usually indicates that something has gone wrong in the amplification procedure.  2.1.2  Review of some genotype calling algorithms  For a small project, it is possible to make genotype calls manually. In most cases, it is not hard for an expert to perform this job, and the “eyeballing” procedure usually gives reasonable results due to its sophisticated incorporation of prior information. For large-scale studies, however, manual scoring can become a daunting challenge. Furthermore, humans are likely to make errors due to fatigue or oversight when a job becomes routine and different readers may have different views (Kang et al., 2004; van den Oord et al., 2003). van den Oord et al. (2003) conducted a study on the error rates of manual scoring and three statistical procedures and report that these statistical procedures uniformly outperform the manual procedure in error rates. Perhaps the most widely used clustering method in the SNP genotyping literature is the classical k-means, possibly with minor modification (Akula et al., 2002; Olivier et al., 2002; Ranade et al., 2001). Olivier et al. (2002) notice that k-means algorithms often split one genotype group of the data, especially when the variant allele homozygous cluster has only a few data points or when there is no sharp distinction between one of the genotype clusters and the NTC cluster. They develop a new algorithm, the CA algorithm, that initially examines all data points and determines the approximate location of the centroid for each cluster by dividing the space into sections for each expected cluster. They report that the heuristic algorithm works better than classical k-means. Ranade et al. (2001) and Akula et al. (2002) assign a “quality score” to each sample. After the genotype assignments are made, they assume that each genotype cluster is bivariate normally distributed. For each sample, 12  Chapter 2. Automatic SNP genotype calling they calculate the probability density value of the normal distribution for its genotype cluster as the quality score that this sample is called. Akula et al. (2002) also assume equal covariance matrices across the clusters, arguing that this assumption gives conservative but better results than otherwise. Lovmar et al. (2005) propose using “silhouette scores” to assess the quality of SNP genotype clusters, the idea of which originates from Rousseeuw (1987). They report that the measure is satisfactory and empirically the genotypes can be unequivocally assigned without manual inspection when the silhouette score is greater than 0.65. van den Oord et al. (2003) propose a mixture model approach, in which each cluster is assumed to follow a normal distribution. The number of component clusters and the initial values of the parameters are set upon inspecting the scatterplot. This is a more delicate approach than the kmeans approach since it considers the elliptical structures in the assignments of genotypes and the k-means algorithm can be regarded as a special case of the mixture models approach. Kang et al. (2004) go one step further in this direction. They assume a bivariate t-mixture model, in which a t distribution with small degrees of freedom is assumed for each cluster. They argue that clusters from SNP genotyping data usually have heavier tails than the normal distribution and report that the t-mixture model is less sensitive to outlying points and has advantages over the normal mixture model or the k-means approach. Fujisawa et al. (2004) also use a model-based approach in which only angle data are used.  2.1.3  ROX-normalized versus unnormalized data  The datasets illustrated in Figure 2.1 are after ROX-normalization. An important feature of our approach is that it uses data without ROX-normalization. We now discuss four reasons why we prefer to work with the unnormalized data. First, the motivation for ROX-normalization is presumably to form spherical clusters for which classical clustering algorithms such as k-means can be used. Inspecting the ROX-normalized scatterplots, there are quite a few plates for which the normalization does not produce reasonable spherical structures as expected, e.g., Figure 2.1 (c). Second, the incorporation of ROX intensities is aiming to correct systematic biases in the chemistry such as plate to plate variation. However, there are scenarios where some points may be pulled away from their home clusters merely by this correction and fail to be called a genotype; by the same token, points which fail to amplify may be assigned a genotype incorrectly 13  Chapter 2. Automatic SNP genotype calling due to this correction. In the left panel of Figure 2.2, points a, b, c are away Unnormalized data  4  20000  ROX−normalized data  d  c  b d  10000  2  Allele.Y.Rn 3  Allele.Y 15000  c  e  e  a  1  5000  a b  0.5  1.0  1.5 Allele.X.Rn  2.0  2.5  2000  4000  6000 8000 Allele.X  12000  Figure 2.2: Scatterplots of unnormalized and ROX-normalized data for a plate. Both points and letters represent samples. from the heterozygous cluster when the data are ROX-normalized, but are within a linear arranged cluster in the right panel when the unnormalized data are used. Similarly, point d is away from the variant allele homozygous cluster in the left panel, but is much closer to the variant allele homozygous cluster in the right panel. The proprietary software assigns very low quality values for these points. Intuitively, these points “should” be classified into one of the genotype clusters, or at least with higher “confidence” than their “quality” values indicate. Point e is outlying in both panels; in this case, normalization will not help any way. Third, the undisclosed, proprietary algorithm used by the SDS system for clustering, and hence calling genotypes, seems to have difficulty with samples with extreme (small or large) ROX values, as evidenced in Figure 2.3. The curve in Figure 2.3 shows that the empirical chance of calling a genotype by the proprietary algorithm decreases when the ROX value is too low or too high. Unnormalized data, like in Figure 2.2, often show well-separated clusters, but they are along lines. The linear grouping algorithm (LGA) (Van Aelst et al., 2006) can identify such clusters. When we apply LGA to the unnormalized data, the problem of low-call rate for extreme ROX values is much alleviated. Figure 2.4 contrasts the silhouette width (a measure of calling quality) versus ROX value when the LGA algorithm is applied to unnormal14  0.0  Empirical chance of calling 0.2 0.4 0.6 0.8  1.0  Chapter 2. Automatic SNP genotype calling  0  5000  10000  15000  ROX  Figure 2.3: Probability of calling a genotype versus ROX ized data and when the k-means method is used for ROX-normalized data. For k-means, the definition of silhouette width in Rousseeuw (1987) is used; for the LGA, the definition of silhouette width in Van Aelst et al. (2006) is used (see also equation (2.1)). Last, the positions of positive control points in scatterplots of ROXnormalized data suggest the use of unnormalized data. In a TaqMan SNP genotyping assay, individual samples are usually arranged in a 96- or 384well plate and are amplified simultaneously. For quality assessment purpose, some wells have DNA samples with known genotypes (in our example 12 in a 384-well plate), called positive controls, while some other wells have no DNA samples (typically 8 in a 384-well plate), called negative controls or no template controls (NTC). The positions of control points in one plate are shown in Figure 2.5. The positive control points are represented by “P”. Note that quite a few control points are outside of their corresponding genotype clusters in the left panel where ROX-normalized data are used. In the right panel, however, these points are well within their corresponding linear clusters. In conclusion, we shall use unnormalized data hereafter in our method of making genotype calls.  15  0.0  0.0  0.2  0.2  LGA Silhouette Width 0.4 0.6  0.8  kmeans silhouette Width 0.4 0.6 0.8  1.0  1.0  Chapter 2. Automatic SNP genotype calling  0  1750 3250 4750 6250 7750 9500 ROX  15500  0  1750 3250 4750 6250 7750 9500 ROX  15500  Figure 2.4: Silhouette width versus ROX value. Left panel: k-means results using ROX-normalized data; right panel: LGA result using unnormalized data.  Unnormalized data 20000  4  ROX−normalized data  P P P P  P  Allele.Y 15000  Allele.Y.Rn 3  P  P P P  2  10000  P  P P  P  PP P  5000  P  1  P P N N N  0.5  1.0  1.5 Allele.X.Rn  2.0  2.5  PPP  P P  N N N N N N  2000  4000  6000 8000 Allele.X  12000  Figure 2.5: Control points in scatterplots of ROX-normalized and unnormalized data for a plate. Letter “P” represents a positive control point and ”N” a negative control point.  16  Chapter 2. Automatic SNP genotype calling  2.2  Results  As an illustration, we apply the LGA approach to the unnormalized data in Figure 2.2. The silhouette threshold is set to be 0.75, which means a point is called if its distance from the nearest line is one quarter of its distance from the second nearest line. We set the signal threshold empirically. Let mx and my be the median values for X signals (“Allele.X”) and Y signals (“Allele.Y”) respectively. Any points with X signal less than 0.5mx and Y signal less than 0.5my are not called. The results are in Figure 2.6. More points are called compared with the SDS software.  5000  Allele.Y 10000 15000  20000  25000  Unnormalized data  0  0  5000  Allele.Y 10000 15000  20000  25000  Unnormalized data  0  5000  10000 Allele.X  15000  0  5000  10000  15000  Allele.X  Figure 2.6: Left panel: the calling result of SDS software displayed in the scatterplot of unnormalized data. Right panel: the scoring results of LGA using unnormalized data. The symbol ◦ represents noncalls.  2.3  Conclusions  We have proposed an automatic genotype calling algorithm by taking advantage of a linear grouping algorithm (Van Aelst et al., 2006). The proposed method uses unnormalized signals and clusters points around lines as against centroids. In addition, we associate a quality value, silhouette width (Rousseeuw, 1987; Van Aelst et al., 2006), with each DNA sample and a whole plate as well. This algorithm shows promise for genotyping data from TaqMan technology (Applied Biosystems). The algorithm is reliable. 17  Chapter 2. Automatic SNP genotype calling It could be potentially adapted to other fluorescent-based SNP genotyping technologies as well, such as the Invader Assay.  2.4  Methods  In the proposed calling algorithm, unnormalized data are used and each genotype cluster is represented by a straight line.  2.4.1  Data preprocessing  In a TaqMan SNP genotyping assay, the quality of a plate is monitored by making sure that all positive controls go to the correct genotype clusters and negative controls stay in the NTC clusters. Without well to well contamination, the negative controls should have low signals for both VIC and FAM dyes (horizontal and vertical coordinates in a scatterplot, respectively) since there are no samples in these wells to be amplified. It is advantageous that we empirically discard negative control points and a portion (5% or 10%, say) of points with very low signals (close to the origin) prior to applying LGA to locate the lines. Negative control points should not contribute to the lines which decide the orientations of the genotype clusters. (Some points with low signals are assigned to lines afterwards.)  2.4.2  Fitting lines using LGA  In SNP genotyping, there are at most three clusters, ignoring the NTC cluster. Two clusters are possible when there are very few or no points in the upper-left, YY cluster. We shall fit three lines and two lines to the data separately using LGA. LGA minimizes aggregated sum of squared orthogonal distances of points to their closest hyperplanes. We made some modification specific to the genotyping setting. Grid based initialization. LGA usually depends on multiple starts to have large probability of attaining the global minimum. In the genotyping setting, it is computationally affordable to try a sequence of initializations such that the global minimum is very likely attained. We set the intercepts to 0 for all lines and choose the slopes from a set of values {s1 , · · · , sm }, where atan(si ) are an equally spaced grid from 0 to π/2, since the slopes should always be positive in the SNP genotyping setting. Given a set of intercepts and slopes, points are assigned to their nearest lines; these lines are recalculated. These two steps are iterated until convergence. Denote by s1 (i) and s2 (i) the orthogonal distances of point i from its 18  Chapter 2. Automatic SNP genotype calling closest and second closest lines, respectively. The silhouette width for point i is defined as s1 (i) . (2.1) w(i) = 1 − s2 (i) From the multiple starts, we choose the solution that has the largest average silhouette width. (This criterion respects small clusters and penalizes two lines which are too close.) In addition, we shall restrict the slope for the XX cluster be smaller than that of the XY cluster and the slope of the XY cluster be smaller than that of the YY cluster. Another restriction is added such that there are few points in YY than XX and XY. For each LGA, we associate the lines with smallest, second-smallest and largest slopes to XX, XY, and YY, respectively. We simply discard solutions in which any slope is negative or there are more points in the cluster with the largest slope than one of the other clusters. Number of lines. We also choose between the best two line solution and the best three line solution according to this average silhouette width. If the three line solution has large average silhouette width, each line represents a genotype cluster; if the two line solution is chosen, they correspond to XX and XY.  2.4.3  Genotype assignment  For genotype assignment, we could empirically set a signal threshold beforehand, below which a point is not called, labelled as “Undetermined”, and assigned to the “NTC” cluster. The threshold can be set based on previous genotyping practice in a laboratory. In this case, a point beyond the threshold is assigned to its closest line and its genotype is called accordingly; a silhouette width value is calculated for each point to measure adequacy of the fitted line structures. A cutoff point for the silhouette widths is also predetermined such that outlying points with low silhouettes are also labelled as “Undetermined”. Our code adapts LGA to deal with thresholding of points as follows. Let mx = median(xi ), and my = median(yi ), where xi and yi are the X signal and Y signal of sample i respectively. Points with xi ≤ 0.5mx and yi ≤ 0.5my are not called. 19  Chapter 2. Automatic SNP genotype calling As an alternative to an empirical signal threshold, we could penalize points closer to the origin by modifying the definition of silhouette width. For example, s1 (i) + c , (2.2) w(i) = 1 − s2 (i) + c where the tuning parameter c is also set empirically such that points too close to the origin are not called. Another version we tried is as follows. Let ri =  x2i + yi2 ,  be the distance of point i from the origin. Let m = min{ri }, and M = median{ri }. Let ci =  min(ri , M ) − M−m 2  m 2  1 2  .  The adjusted value for silhouette width si is s∗i = ci si .  (2.3)  This downweights the silhouette score for points with weak xi and yi signals.  2.4.4  Quality assessment of DNA samples and plates  The silhouette width for a DNA sample serves as a quality value of the genotype call. In addition, for each plate, the average silhouette width can be computed. An expert may decide upon a threshold such that plates with lower average silhouette width are regarded as unreliable. Meanwhile, the average silhouette width of positive controls in a plate and the average silhouette width of negative controls are also indications of the quality of genotype calls for a specific plate.  2.5  Discussion  The existing methods assume spherical or elliptical clusters. The unnormalized signals tend to fall on lines, however, forcing normalization to try to 20  Chapter 2. Automatic SNP genotype calling produce clusters with elliptical shapes, but the ROX normalization signal may be noisy. We take a quite different approach, identifying the linear clusters associated with the unnormalized signals. As seen in Figure 2.6, fewer uncalled samples results. We call only samples with at least one strong signal. The definition of “strong” will vary from plate to plate. We normalize internally by not calling points having weak signals relative to median signals, rather than externally with respect to ROX.  21  Bibliography Akula, N., Y. Chen, K. Hennessy, T. Schulze, G. Singh, and F. McMahon (2002). Utility and accuracy of template-directed dye-terminator incorporation with fluorescence-polarization detection for genotyping single nucleotide polymorphisms. Biotechniques 32 (5), 1072–1078. De La Vega, F., K. Lazaruk, M. Rhodes, and M. Wenz (2005). Assessment of two flexible and compatible SNP genotyping platforms: TaqMan snp genotyping assays and the SNPlex genotyping system. Mutation research. Fundamental and molecular mechanisms of mutagenesis 573 (1-2), 111–135. Fujisawa, H., S. Eguchi, M. Ushijima, S. Miyata, Y. Miki, T. Muto, and M. Matsuura (2004). Genotyping of single nucleotide polymorphism using model-based clustering. Bioinformatics 20, 5. Kang, H., Z. Qin, T. Niu, and J. Liu (2004). Incorporating Genotyping Uncertainty in Haplotype Inference for Single-Nucleotide Polymorphisms. The American Journal of Human Genetics 74 (3), 495–510. Lovmar, L., A. Ahlford, M. Jonsson, and A. Syv¨anen (2005). Silhouette scores for assessment of SNP genotype clusters. feedback . Olivier, M., L. Chuang, M. Chang, Y. Chen, D. Pei, K. Ranade, A. de Witte, J. Allen, N. Tran, D. Curb, et al. (2002). High-throughput genotyping of single nucleotide polymorphisms using new biplex invader technology. Nucleic Acids Research 30 (12), e53. Ranade, K., M. Chang, C. Ting, D. Pei, C. Hsiao, M. Olivier, R. Pesich, J. Hebert, Y. Chen, V. Dzau, et al. (2001). High-Throughput Genotyping with Single Nucleotide Polymorphisms. Genome Research 11 (7), 1262– 1268. Rousseeuw, P. (1987). Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics 20, 53–65. 22  Bibliography Shi, M. (2001). Enabling large-scale pharmacogenetic studies by highthroughput mutation detection and genotyping technologies. Clin Chem 47 (2), 164–172. Van Aelst, S., X. Wang, R. Zamar, and R. Zhu (2006). Linear grouping using orthogonal regression. Computational Statistics and Data Analysis 50 (5), 1287–1312. van den Oord, E., Y. Jiang, B. Riley, K. Kendler, and X. Chen (2003). FD-TDI SNP Scoring by Manual and Statistical Procedures: A Study of Error Rates and Types. BioTechniques 34 (3), 610–624.  23  Chapter 3  A partial likelihood approach to linear clustering 3.1  Introduction  In this paper we are concerned with detecting linearly shaped clusters in a data set. This work is motivated by a clustering problem in a SNP (single nucleotide polymorphism) genotyping setting. Data of this type are bivariate, consisting of two signals, and conform mostly to linear clusters. More detailed discussion of these genotyping data appears in Section 3.7. Van Aelst et al. (2006) proposed a method called the Linear Grouping Algorithm (LGA) for linear clustering problems. For earlier approaches to linear clustering, see the references therein. There is a need for algorithms specialized in linear clustering. First, the usual mixture of normal or t distributions (see for example, Banfield and Raftery, 1993; Fraley and Raftery, 1998, 2002; MacLachlan and Peel, 2000, among others) aims at elliptical structures in a data set; when a data set displays highly linear patterns, especially when these patterns are not well separated from each other, we need a different criterion for detecting linear patterns. Second, most earlier linear clustering approaches assume that a response variable is available supervising the search for linear clusters, see for example, Spath (1982), DeSarbo and Cron (1988) and Turner (2000). Van Aelst et al. (2006) have a simulated data set illustrating that the selection of a response variable is not trivial. In that data set, two linear clusters are obvious in the first two dimensions; the third dimension is merely noise. If the second variable is selected as response, the mixture of regressions method by Spath (1982) successfully recovers the linear structures; it fails when the third variable is used as the response. Chen et al. (2001) proposed a method to detect linear A version of this chapter has been submitted for publication. Authors: Guohua Yan, William J. Welch and Ruben H. Zamar.  24  Chapter 3. A partial likelihood approach to linear clustering structures one by one without assuming a response variable in the computer vision context. We propose in this paper a partial likelihood approach and compare it with existing algorithms such as LGA (Harrington, 2007), MCLUST (Fraley and Raftery, 2007) and MIXREG (Turner, 2006). Our approach borrows ideas from finite mixture model-based clustering and from mixture of regressions. Each linear cluster is characterized by a hyperplane; the orthogonal deviation of a data point from that hyperplane is assumed to be normally distributed but the position of the data on the hyperplane is not modelled. With this strategy, all variables are treated symmetrically, as compared with the mixture of regressions approach. The rest of the paper is organized as follows. In Section 3.2, we describe a partial likelihood-based objective function, leading to a model-based linear clustering algorithm, MLC. In Section 3.3, an EM algorithm to maximize the partial likelihood function is presented. In Section 3.4, we discuss the asymptotic properties of MLC. In Section 3.5, we briefly discuss its relationships with existing algorithms such as MCLUST, LGA and MIXREG. In Section 3.6, we propose several methods for determining the number of linear clusters. In Section 3.7, several real and simulated datasets are used to compare MLC with LGA, MCLUST and MIXREG. Some closing remarks are given in Section 3.8.  3.2 3.2.1  Partial likelihood-based objective function for linear clustering Partial likelihood for orthogonal regression  First, we formulate a partial likelihood representation for orthogonal regression (see for example, Fuller, 1987). To model a single linear cluster, assume that the vector-valued data x1 , . . . , xn are independently drawn from a random mechanism represented by a random vector X in a d-dimensional space. We do not model the distribution of X except for the deviation from an also unknown hyperplane. Let {x : a x − b = 0} be a hyperplane. For identifiability purpose, we assume that a a = 1 and that the first nonzero element of a is positive. We assume that the signed orthogonal deviation a X − b from the hyperplane {x : a x − b = 0}, is normally distributed, i.e., a X − b ∼ N (0, σ 2 ).  (3.1)  25  Chapter 3. A partial likelihood approach to linear clustering Since a and b are unknown, a X − b is also unknown. Nevertheless, we can still estimate a, b and σ 2 through the following partial likelihood function n i=1  N (a xi − b; 0, σ 2 ),  (3.2)  where N (· ; 0, σ 2 ) is the density function of the normal distribution with mean 0 and variance σ 2 . ¯ and S be the sample mean and sample covariance matrix of Let x x1 , . . . , xn respectively. Then the maximum partial likelihood estimate σ ˆ 2 of σ 2 is the smallest eigenvalue of S; the maximum partial likelihood estimate ˆ of a is the (standardized) eigenvector associated with σ a ˆ 2 , which is not necessarily unique. Finally, the maximum partial likelihood estimate ˆb of b ¯ . See for example Fuller (1987). ˆx is ˆb = a Note that model (3.1) and the partial likelihood (3.2) treat the components of x symmetrically.  3.2.2  Partial likelihood for linear clustering  Now we assume that the data x1 , . . . , xn are independently drawn from a more complicated random mechanism, still represented by a random vector X in a d-dimensional space which we do not model fully. The data now lie around K hyperplanes {x : ak x = bk }, k = 1, . . . , K. Let Z = (Z1 , . . . , ZK ) be a random vector indicating these hyperplanes, where Zk = 1 with probability pk for k = 1, . . . , K. Let p = (p1 , . . . , pK ) . We assume that, conditional on Zk = 1, ak X − bk ∼ N (0, σk2 ), k = 1, . . . , K. Let z1 , . . . , zn be the corresponding unobservable indicators for the data x1 , . . . , xn . Let κ be the collection of component parameters, κ = (a1 , b1 , 2 ) , and θ = (κ , p ) . σ12 , . . ., aK , bK , σK When the indicators z1 , . . . , zn are regarded as unknown parameters, the partial likelihood function for parameters (θ , z1 , . . . , zn ) is n  L(θ, z1 , . . . , zn |x1 , . . . , xn ) =  K  [pk N (ak xi − bk ; 0, σk2 )]zik .  (3.3)  i=1 k=1  This is a so-called classification likelihood. See for example Scott and Symons (1971) and Banfield and Raftery (1993) in the context of mixture 26  Chapter 3. A partial likelihood approach to linear clustering models for elliptical clusters. In this approach, the parameters θ and indicators z1 , . . . , zn are chosen to maximize (3.3) and data points are classified according to these indicators. In another view (see for example, Fraley and Raftery, 1998, 2002), the indicators z1 , . . ., zn are regarded as realizations of random vectors Z1 , . . . , Zn , which in turn are an independent and identically distributed sample from Z. After integrating out the indicators z1 , . . . , zn , the partial likelihood function for θ is n  K  L(θ|x1 , . . . , xn ) = i=1 k=1  pk N (ak xi − bk ; 0, σk2 ).  (3.4)  The latter is a mixture likelihood approach. Celeux and Govaert (1993) did simulation studies of classification likelihood approaches and mixture likelihood approaches in normal mixture models for elliptical clusters and reported that no likelihood method uniformly outperforms the others. It is noted that the classification approach cannot consistently estimate parameters due to the “all-or-nothing” classification bias (Bryant and Williamson, 1978; Ganesalingam, 1989; Marriott, 1975). In our experience with linear clustering, both approaches give practically very similar clustering results. We shall pursue in this paper the mixture likelihood (3.4), as many existing model selection criteria like BIC can be used and allows for a more feasible asymptotic theory (see Section 3.4). As in the usual normal mixture model for elliptical clusters, the partial likelihood function (3.4) is unbounded: when a cluster consists of only points lying on a hyperplane, the contribution of each of these points to the partial likelihood tends to infinity as the variance tends to zero. The infinity occurs on the boundary of the parameter space. Hathaway (1985) proposed a constrained formulation of the maximum likelihood estimation in the univariate normal mixture model. Specifically, the author added a constraint on the standard deviations, min (σi /σj ) ≥ c > 0,  1≤i=j≤K  (3.5)  where c is a known constant determined a priori. We shall incorporate the constraint (3.5) into the partial likelihood function (3.4) (See Step 2 in the EM algorithm in Section 3.3). The solution may depend on the choice of the constant c. A usual strategy is to decrease c gradually and monitor the resulting solutions (Hathaway, 1986). [Chen and Kalbfleisch (1996); Ciuperca et al. (2003) proposed penalized approaches to this unboundedness problem which are also readily applicable to linear clustering although not 27  Chapter 3. A partial likelihood approach to linear clustering adopted in this paper; on the other hand, the Hathaway’s constraint can also be regarded as a form of penalization on the likelihood.] The partial likelihood function (3.4) naturally brings the clustering problem into a finite mixture model framework. Standard EM algorithms can be ˆ adapted to maximize (3.4); once the maximum partial likelihood estimate θ is obtained, data point xi can be assigned to the component with the largest posterior probability. The probabilities are given by w ˆik =  pˆk N (ˆ ak xi − ˆbk ; 0, σ ˆk2 ) , i = 1, . . . , n; k = 1, . . . , K, K pˆk N (ˆ a xi − ˆbk ; 0, σ ˆ2) k=1  k  (3.6)  k  which also serve as a measure of uncertainty of classifying data point xi .  3.3  The EM algorithm  Now we describe an EM algorithm for maximizing the partial likelihood (3.4). It is straightforward and works well in our experience. The completed log partial likelihood is K  n  l(θ|x1 , . . . , xn , z1 , . . . , zn ) = i=1 k=1  zik {log(pk ) + log(N (ak xi − bk ; 0, σk2 ))}. (3.7)  In the E-step, with θ  (t−1)  from the previous iteration, we have  (t)  wik ≡ E(Zik |θ (t−1) , x1 , . . . , xn ) (t−1)  =  pk  (t−1)  N (a k  (t−1)  xi − bk  (t−1) (t−1) K N (a k xi k=1 pk  −  ; 0, σk2  (t−1)  )  (t−1) (t−1) bk ; 0, σk2 )  .  (3.8)  In the M-step, we have (t)  pk = Let (t) ¯k x  and (t) Σk  n  = i=1  (t) n i=1 wik  n  .  =  (t) n i=1 wik xi , (t) n i=1 wik  (t)  (t)  (3.9)  (t)  ¯ (k) )(xi − x ¯ (k) ) . wik (xi − x  28  Chapter 3. A partial likelihood approach to linear clustering Then (t)  (t)  ak is the eigenvector associated with the smallest eigenvalue of Σk , (3.10) (t) (t) (t) ¯ (k) , (3.11) bk = a k x and σk2  (t)  (t) 2 (t) (t) n i=1 wik (a k xi − bk ) . (t) n w i=1 ik  =  (3.12)  If all σk2 in equation (3.4) are assumed to be equal, then the common value is estimated by σ  2 (t)  =  n i=1  (t) (t) K k=1 wik (a k xi  n  (t)  − bk )2  .  (3.13)  The MLC algorithm is described as follows. 1. Initialize with θ (0) . To initialize the EM algorithm, we adopt the strategy of Van Aelst et al. (2006). We randomly select K mutually exclusive subsets of d + 1 observations from x1 , . . . , xn . For k = 1, . . . , K, (0) (0) we compute the maximum partial likelihood estimates ak , bk and (σk2 )(0) from the kth subset of d+1 observations (see Subsection 3.2.1). (0) The proportions pk are all set as 1/K. The initial values are then (0) (0) 2 )(0) , p(0) , . . . , p(0) ) . θ (0) = ((a1 )(0) , b1 , (σ12 )(0) , . . . , (aK )(0) , bK , (σK 1 K 2. If constraint (3.5) is not satisfied, go back to step 1; otherwise go to step 3. 3. Update θ (0) for a predefined number of iterations or until the improvement in the partial likelihood (3.7) is less than a predefined threshold. At iteration t, (t)  • E-step. Update wik by equation (3.8). (t)  (t)  (t)  • M-step. Update pk , ak , bk and (σk2 )(t) by equations (3.9), (3.10), (3.11) and (3.12), respectively. If c = 1, which corresponds to equal variances across clusters, use (3.13) in place of (3.12). 4. Repeat steps 1 – 3 for a predefined number of times. The final clusˆ are the solution that ˆ1 , . . . , z ˆn and parameter estimates θ ter labels z has the largest completed log partial likelihood (3.4) and satisfies constraint (3.5). 29  Chapter 3. A partial likelihood approach to linear clustering For the classification likelihood (3.3), we can use a so-called CEM (classification EM) algorithm. A C-step is inserted between the E-step and M-step (t) of step 3. In the C-step, xi is assigned to the cluster with the largest wik , (t) (t) i.e., the maximum wik is replaced by the value 1 and all other wik take the (t) value 0. If the maximum value of wik is not unique, choose one of the tying clusters at random.  3.4  Asymptotic properties  Let  K  f (x; θ) = k=1  pk N (ak x − bk ; 0, σk2 ),  which is the contribution of one observation to the partial likelihood (3.4). The constrained parameter space is   2 2   θ = (a1 , b1 , σ1 , . . . , aK , bK , σK , p1 , . . . , pK ) : ak ak = 1, −∞ < bk < ∞, σk > 0, k = 1, . . . , K. . Θc =   p = 1. mini,j σi /σj ≥ c > 0, 0 < pk < 1, K k=1 k  Assuming that the data arise from a probability distribution measure P , we have the following results.  Theorem 3.1 Let P be an absolutely continuous probability measure with finite second moments. Let g(θ) ≡  log f (x; θ)dP (x).  Then the supremum of g over Θc is finite and attainable. Since the function f (·, θ) is not a density function, this theorem and the following ones cannot be proved by verifying some regularity conditions. The proofs find their roots in Wald (1949), Redner (1981), Hathaway (1983, 1985) and Garc´ıa-Escudero et al. (2007) and are available on request. Theorem 3.2 Let P be an absolutely continuous probability measure with finite second moments. Let X1 , . . . , Xn be a random sample from P . Then a global maximizer of the partial likelihood function L(θ |X1 , . . ., Xn ) in (3.4) over Θc exists almost surely if n ≥ Kd + 1. 30  Chapter 3. A partial likelihood approach to linear clustering Note that points in Θc are not identifiable for f (x; · ). The function f (x; · ) remains the same if we permute the labels 1, . . ., K; pk1 and pk2 are not identifiable if (ak1 , bk1 , σk21 ) = (ak2 , bk2 , σk22 ). Thus the consistency result is in a quotient topology space. Let ∼ be an equivalent relation on Θc such that θ 1 ∼ θ 2 if and only if f (x; θ 1 ) = f (x; θ 2 ) almost surely in P . Denote by Θqc the quotient topological space consisting of all equivalent classes of ∼. For a point θ 0 that maximizes g(θ) = log f (x; θ)dP (x), its equivalent class in Θqc is denoted by θ q0 . Theorem 3.3 (Consistency). Let X1 , . . . , Xn be a sample from an absolutely continuous probability measure P with finite second moments. Let ˆ (n) be a global maximizer of the partial likelihood function L(θ|X1 , . . ., θ ˆ (n) → θ q almost surely in the topological space Xn ) in (3.4) over Θc . Then θ 0  Θqc .  Let v1 (θ) = E  ∂ log f (X; θ) ∂θ  and v2 (θ) = E  ∂ log f (X; θ) ∂θ  ∂ 2 log f (X; θ) ∂θ∂θ  ,  ,  where the expectations are taken with respect to P . Theorem 3.4 (Asymptotic normality). Let X1 , . . . , Xn be a sample from an absolutely continuous probability measure P with finite sixth moments. ˆ (n) be a subsequence of global maximizers of the partial likelihood funcLet θ tion L(θ|X1 , . . . , Xn ) in (3.4) over Θc , which tends to an interior point θ 0 . Then √ L ˆ (n) − θ 0 ) → n(θ N (0, v(θ 0 )),  where v(θ 0 ) = [v2 (θ 0 )]+ v1 (θ 0 )[v2 (θ 0 )]+ and A+ is the Moore-Penrose inverse of A.  ˆ Denote w ˆ In equation (3.6), w ˆik is a function of xi and θ. ˆik = hk (xi , θ), k = 1, . . . , K. By the Delta method, we have Corollary 3.1 Let X1 , . . . , Xn be a sample from an absolutely continuous ˆ (n) be a subsequence probability measure P with finite sixth moments. Let θ of global maximizers of the partial likelihood function L(θ|X1 , . . . , Xn ) in  31  Chapter 3. A partial likelihood approach to linear clustering (3.4) over Θc , which tends to an interior point θ 0 . Let x be a data point. Then for k = 1, . . . , K, √  L (0) (0) ˆ (n) ) − hk (x, θ 0 )) → n(hk (x, θ N (0, [hk (x, θ 0 )] v(θ 0 )[hk (x, θ 0 )]),  where (0)  hk (x, θ 0 ) =  ∂hk (x, θ) ∂θ  θ = θ0  .  Using this corollary, we can build approximate confidence intervals for ˆ and hence evaluate the clustering of a data point. wik by replacing θ 0 with θ  3.5 3.5.1  Relationships with other clustering methods With LGA  In the LGA of Van Aelst et al. (2006), the objective function is the aggregated sum of squared orthogonal distances of the data points to their closest hyperplanes. Using our notation, LGA minimizes n  K  d(κ, z1:n ) = i=1 k=1  zik (ak xi − bk )2 .  (3.14)  In the classification likelihood (3.3), if we assume that the variances σk2 are equal across all components (or equivalently, c = 1 in (3.5)) and the mixing proportions pk are equal, the log completed partial likelihood is n  K  l(κ, z1:n ) = i=1 k=1  zik log[N (ak xi − bk ; 0, σ 2 )].  (3.15)  It is straightforward to check that the minimization of (3.14) and the maximization of (3.15) are equivalent. Hence the iterative procedure in LGA coincides with the CEM algorithm, which is discussed in Section 3.3, for maximizing (3.15). In this sense, the proposed partial likelihood approach, or the classification partial likelihood (3.3), is an extension and a model for LGA. The proposed approach permits one dispersion parameter σk2 for each cluster while LGA implicitly uses one dispersion parameter for all clusters. Furthermore, in our model framework, we are able to use the membership probability wik (see (3.6)) to assess the uncertainty of clustering a data point while an ad hoc measure, silhouette width, is used in LGA. 32  Chapter 3. A partial likelihood approach to linear clustering  3.5.2  With normal mixture models  In the M-step in Section 3.3, we can see that the proposed approach is closely related to normal mixture models. The difference resides in the Estep: the proposed approach weighs distances to hyperplanes while the latter compares distances to cluster centres. In the use of normal mixture models one is forced to estimate the means and also the covariance matrix of each component. Banfield and Raftery (1993) reparameterized the covariance matrices and lead to a variety of parsimonious covariance matrix structures. The unrestricted covariance matrix requires K(d + d(d + 1)/2) parameters and the most parsimonious model, equal spherical covariance matrices, requires only Kd + 1 parameters. The proposed approach can be regarded as a parsimonious simplification of another type which is appropriate for highly linear data sets. It requires d + 1 parameters per cluster, i.e., d first order parameters for the location of a hyperplane and one second order parameter for the dispersion around the hyperplane.  3.5.3  With mixture of ordinary regressions  As mentioned in the introduction, a mixture of ordinary regressions approach requires a response variable for guidance. A response variable suitable for clustering purpose may not exist. Furthermore, due to the “regression to the mean” phenomenon, the regression hyperplanes are different when different variables are selected as the response. The proposed approach treats each variable symmetrically, and therefore is more suitable in a clustering setting.  3.6  Choosing the number of linear clusters  In some problems, the number of linear clusters is obvious from subjectmatter knowledge; in other problems, it has to be estimated from the data. Too many components exploit the randomness of the data and make it hard to interpret the clustering; too few clusters lose information and may be misleading.  3.6.1  Bootstrapping the partial likelihood ratio  A natural approach is to apply the partial likelihood ratio test to a sequences of hypotheses K = K0 against K = K0 + 1, starting from K0 = 1. However, 33  Chapter 3. A partial likelihood approach to linear clustering −2 ln λ, where λ is the likelihood ratio, does not have a usual asymptotic χ2 distribution under the null hypothesis. The reason is that the regularity conditions are violated as the parameter space under the null hypothesis is on the edge of the parameter space under the alternative hypothesis when the former is embedded into the latter. Much theoretic research has been done on likelihood ratios for mixture problems (Chen and Chen, 2001a,c,b; Chen et al., 2001; Chen, 1998; Chen and Kalbfleisch, 2005; Everitt and Hand, 1981; Thode Jr et al., 1988, among others). We adapt the bootstrap approach in McLachlan (1987) to our setting. The log likelihood ratio statistic −2 log λ for the test of the null hypothesis of K = K1 versus the alternative K = K2 can be bootstrapped as follows. ˆ be the maximum partial likelihood estimate of all parameters from the Let θ original sample under the null hypothesis. For i = 1, . . . , n, sample zi from the multinomial distribution with probabilities (ˆ p1 , . . ., pˆK1 ); if zik = 1, sample orthogonal distance ei from N (0, σ ˆk2 ). The ith data point in the bootstrap sample is xi + (ei − ˆbk )ˆ ak . Suppose B copies of n samples are generated. The value of −2 log λi is computed after fitting mixture models for K = K1 and K = K2 , i = 1, . . . , B. The significance level of the test is obtained by comparing −2 log λ with the empirical distribution of −2 log λ1 , . . . , −2 log λB .  3.6.2  Information criteria  The partial likelihood approach enables us to use an approximate Bayes factor to determine the number of linear clusters. Banfield and Raftery (1993) used a heuristically derived approximation to twice the log Bayes factor, called AWE (approximate weight of evidence) to determine the number of clusters. Fraley and Raftery (2002) used BIC (Bayesian information criterion) as a more reliable approximation to twice the log Bayes factor, BIC(K) = 2L(K) − vK log(n),  (3.16)  where L(K) is the log partial likelihood of a mixture of K components and vK is the number of independent parameters to be estimated in a K component mixture. As stated in Fraley and Raftery (2002), although the regularity conditions that underly the approximation are not satisfied in finite mixture models, there are some theoretical and practical support of its use for determining the number of clusters. Keribin (2000) showed that the estimator of the number of clusters based on BIC is consistent; Fraley and Raftery (1998, 2002) included a range of applications in which estimating the number of clusters based on BIC has given good results. 34  Chapter 3. A partial likelihood approach to linear clustering A large number of criteria have been proposed in the literature, including NEC (Biernacki et al., 1999; Celeux and Soromenho, 1996), ICL (Biernacki et al., 2000), and cross validation likelihood (Smyth, 2000). The performance of some of these criteria were compared by Biernacki et al. (1999). As we are adopting a partial likelihood approach, in principle almost all methods proposed in a model framework are approximately applicable although we do not fully model the data. As pointed out by Everitt et al. (2001), “it is advisable not to depend on a single rule for selecting the number of groups, but to synthesize the results of several techniques”. In our experience, though, all the methods mentioned above work well if there do exist obvious linear patterns in a data set.  3.7 3.7.1  Examples Simulated data I  One data set of size 300 in two dimensions is simulated to illustrate the different behaviours of MCLUST, LGA and the proposed MLC. The xi1 are uniformly sampled from the interval (0, 15); for the first 250 points, xi2 = xi1 + 3 i and for the last 50 points, xi2 = 8 + 1.5xi1 + i , where i are independent and identically standard normal distributed. The data are shown in the top-left panel of Figure 3.1. For this data set, MCLUST chooses a model of three clusters with different volumes, different shapes and different orientations by the BIC criterion. LGA favours one cluster using the gap statistic. MLC selects two clusters by BIC, ICL, NEC and the bootstrapping likelihood ratio test (see Table 3.1). We set the threshold c = 0.05. The result of MIXREG, which is omitted, is similar to that of MLC as we generate the data in favour of it. Figure 3.1 displays the clustering results from MCLUST, LGA and MLC when the numbers of clusters are all set as 2. This simulated data helps us understand how these algorithms work. MCLUST tries to find elliptical structures, but the two linear clusters have data distributions along the lines not well modelled by normal distributions. For LGA, when two clusters are roughly parallel and touching, it imposes a structure to the data set and tends to partition the data set into bands with equal width.  35  Chapter 3. A partial likelihood approach to linear clustering  20 x2 10 0  0  10  x2  20  30  MLC  30  True classification  5  10  15  0  5 x1  LGA  MCLUST  10  15  10  15  20 x2 10 0  0  10  x2  20  30  x1  30  0  0  5  10 x1  15  0  5 x1  Figure 3.1: Comparison of clustering results of simulated data I. Upper-left panel displays the true classification; upper-right, lower-left and lower-right are that of MLC, LGA and MCLUST.  36  Chapter 3. A partial likelihood approach to linear clustering  Table 3.1: Criteria for choosing the numbers of clusters, K, when MLC is applied to simulated data I. “Boot p-value” is the p-value using the bootstrapping partial likelihood ratio test for H0 : K = K0 vs H1 : K = K0 + 1 where 99 bootstrapping samples are used Criterion K BIC ICL NEC Boot p-value 1 -1541.61 – – 0 2 -1461.63 -1477.28 0.14 1 3 -1468.25 -1712.08 0.69 –  3.7.2  Simulated data II  Another dataset is simulated to demonstrate the difference between MIXREG and the other three algorithms, MLC, MCLUST and LGA. This dataset has 200 data points in four dimensional space. For the first 100 data points, the first two components, xi1 and xi2 , are linearly related while the third and the fourth components, xi3 and xi4 , are randomly scattered: xi3 , xi4 ∼ N (0, 52 ), for i = 1, . . . , 100. (3.17) Here ij ∼ N (0, 1) and ti ∼ N (0, 52 ). For the remaining 100 data points, the first two components are randomly scattered while the last two components are linearly related: xi1 = ti +  i1 ,  xi2 = 0.5ti +  xi1 , xi2 ∼ N (0, 52 ), xi3 = ti +  i2 ,  for i = 101, . . . , 200. (3.18) 2 Again, ij ∼ N (0, 1) and ti ∼ N (0, 5 ). The pairwise scatterplots of the simulated data are in Figure 3.2. From the way the dataset is simulated, it is not possible to find a natural response variable to discriminate the two linear clusters. As a result, MIXREG does not work well for this dataset and the two clusters are not well separated. The misclassification matrices of the four algorithms are in Table 3.2. MLC, MCLUST and LGA have the same good performance here. Note that (3.17) and (3.18) generate (elongated) multivariate normal clusters in x1 and x2 or in x3 and x4 , respectively. The assumptions of MCLUST are therefore satisfied, but MLC nonetheless has the same performance.  3.7.3  i3 ,  xi4 = 0.5ti +  i4 ,  Australia rock crab data  Now we consider the rock crab data set of Campbell and Mahon (1974) on the genus Leptograpsus. We focus on the 100 blue crabs, 50 of which are 37  Chapter 3. A partial likelihood approach to linear clustering  2  2  10  2 2  22 2 2 22 2 2 2 2 1 22 2 222 1 1 2 22 1 1 1 22 222 2 22 1 11211 11 21 2 2 221 1211111 12111 1111 11 1 1 22 2 112 11112 1121111211121111111 2 2 12111 1 11 11 1 11 222 2 11 2 1 2 2 12 2 1 122 12 121 1 2 2 111 212 1 2 12 2 2 222 2 2 1 1 12 2 2 2 2 2 2 2 22 2 2 2 1 2 2 2  2  2  −10  −5  0  5  2  2  2  2  20  2  1 11 1 2 2  x2  1  2  2 2 2  1  2 1 1 1 2 22 2 1 1 121 12 122 1 21 11 2 2 2 11 11 2 2 1 2 2 22 111 2112 11 1 1 1 1 1 22222 1 2 2 1 22 1 11 22 21 1 1 21111 12 22 1 2 2 2 2 2 21 11 112111111 21 1 2 1 11 1 2 1 122 12 111 21 211 2 2 1 2 2 22 1 2 2 21 2 2 2 1 222 112 2 21 1 11 2 1 2 2 221 1 1 2 112 2 1 1 2 1 2 21 2 1 22 11 1 1 1  10 5  1  2  2 2 2 2 2 22 2 22 2 2 2 1 22 122 2 22 222 12 1 2 1 1 12 2 1 11 1 1 1 1 22 2 12 2 2 1 2 111 11 2 11111 1 11 1 1 11 2 2 122 1 2 2121111 11112 1111 1 1 1 1 11 1 2 2 2 1 2 1 121 1 21 2121 1 1 12 1 1 2 2 2 21 1 2 1 11 11 1 1 21 12 1 21 1 1 1 1 2 1 12 1 2 2 2 2 1 22 2 12 12 2 2 2 2 2 22 2 2 2 22 1 2 2 2 2 2 2 2 2  1  2 2 2 2 22 2 2212 2 2 2 22 212222 2 22 2 1 111 2 2 2 1 1 1 1 1 2222 2 1 11 1 1 1212211 211 2 111 1 1 1 11 1 1 1 11 1 1 1 12 12112211222111111111 21 2 11 1 11 1 1 112 21112 1212111 2122 1 1111 1 1 1 2122 2 1 1 2 22 1 1 1 22 22 11112 1 2 2 22 2 121 2 2 1 22 2 2 1 2 2 22222 22 2 2 2 2 2  1  2 1 1 1 2 22 1 1 1 12 222 1 2 22 1 1 1 1 1 1 1 2221 1 1 1 2 2 11 11 111122122211222212 1 2 22 12 11 1 1 111 1 21222122 2 21 1 1 111222122 221 1 1 1 1 122 1 1 1 1 1 11 1 1212222222 11 1 11 122 222222222 1 1 1 1 1 2212222122 21 1 1 1 11 1 222 1 1 22 1 2 2 11 21 2 1 1 22 1 2  1  1  2 2 2 2  x3  2  1  2  2  20  2 2  1  1  1  0  5  10  15  1  15  2  2  2  2  2  1 2 1 21 2 122 1 11 1 2122 12 2 1 2 12 1 1 1 2 1 21 1 1 2 1 2 1 1 22 111122 22 1 2212 2 1 1 12 1 1 2 2 2 21 1 22 11 1 112 1 2 1 2 2 2 2 12 1 1 2 2 1 1212 2221 1111212 1 1 1 211 1 2 212 1 112 1 211 1 1 2 2 2 22 222 2 1 12211 222121 1 2 2 12 1 222 1 21 2 1 1 2 2 2 1 12 2 2 1 1 12 1 2 1 1 2 2 2 1 2 1 21 2 1 1 2 1 2 2 2 2  2  1  2  1 21  1 2 2 22 2 1 12 2 11 2 1 21 1 1 1211 11 1 1 11122212 11 2 2 2 1 22 2 1 11 1 1 1 2 2 12222 1221 11 1 1 1211 2121221 1 12 11122 2 1 1 2 21 221 221 1 1111 1 1 2 2 2221 2111222 2 1 2 1 12 2122 21 1 1 11 1 11 11221222 12 1 1 1 2 2 1 2 1 1 2 1 2222 2 21 21 1 1 2 2 2 2222 2 12 1 2 222 1 12 2 1  2  2  10  0  x1  2 2 2 1 2 11 1 22 1 11 2 1 1 1 112 11211 2 1 1 2 1 1 1 1 2 2 2 2 1 1 12 22 2 221 12 12 2 11 21 1 2 12 2 21 1 2 1 111111 1 1 12 2212 2 1 1 1 2 1 2 1 1 2 1 2 1 1 222 111222221 1 22 12 1 2112 21 2 2 2 1 2 1 2 2 2 121 1211 11 1 12 2 1 1 2 1 2 121 2 2221 2 22 21 1 1222 1 2 2 21211 2 2 2 2 21 2 2 1 221 2 1 2 2 1 2  1  2  5  2  5  2  1 2  1  0 1  −5  1 2 2 1 2 21 2 2 11 11 1 1 1 11 2 1 22 11111 1121 1 2 11111111 2 2 2 2 2 2 2 111111 2 1 2 2 22 1 1 2 2 11 12 1 111 2 2 2 22 2 21 111111121 121 2 111 2 1 1 2 1 2 2 2 2 2 2 2 22 22 2 122 11 111 2 1111 1 2 2 2 2 2 2 21 12 111211112 2 2 1 12 1 2 22 222 2 11 2 2 2 2 2 2 12 2 2 2 2 1122 2 1 1 2 2  2  2  −15 −10 −5 2  1  −15  10 1  10  5  0  0 2  −5  −5  −10  −10  −15  −5  2  1 1 1 1 111 1 1 11 1 1 1 12 112 1 2 1 1 1 222 2 2 2 2 12 2221 1211 1 1 12 2 2 22 1212211 221212112212 1221222 2 2 2 2 2 2 2 1 22 21 1 1 2 12 2 12 1 2 1 2 21 2 2 2 1 2 2222 1212112222 2 1 22221221112 1 1 1 1 22 1 22 21 12 2 12 1 2 21 2121 2 1112 1 2 1 1112 21 112 1 21 1 1 1 1 1 1 1 1 1 1 1 1 1  −15  −10  −5  22  2 2  1 1 1 1 1 111 1 1 11 11 1 1 1 21 22 12 1 1 2 2 2 222 2 2 222111 111 1 2 2 22 11 1 2112 111 11 2 2 22 2222 222 22 2 222 22221 21111 1 222 22 21 2 2 2 222 1222 2 22111 22211211 1112122 212 11 22 212 2 2 2 2 2 1 112 22 12111122 1211 1 1 1 2 1 211 12 11 2 2 1 1 1 1 11 1 1 1 1 1 1 1  1 11 1 11 1 1 1 1 2 1 11 2 1 22 11 1 2 2222 2 1 1 11 111 1 122112 212122222222 12222 22 1 2222 21121 2212212221122221 2 1 1 1 1 2 2 12 1 2 22222222222122222121222112 21 112 1 1 1 1 2 2 222 2 221212 222 1 1 11 11 1 21 1 1 1 1 1 2 1 1 1 1 1 1 22 2 2 1 1 1 1 11 1 1 1 1 11 1 11 11  1  0  5  10  x4  1  −10  −5  0  5  10  Figure 3.2: Pairwise scatterplots of simulated data II.  Table 3.2: Misclassification matrices of simulated data II from MLC, MCLUST, LGA and MIXREG.  True class 1 2  MLC, MCLUST, LGA 1 2 88 12 9 91  MIXREG Resp= x1 or x2 Resp= x3 1 2 1 2 26 74 5 95 25 75 10 90  Resp= x4 1 2 63 37 75 25  38  Chapter 3. A partial likelihood approach to linear clustering  Table 3.3: Criteria for choosing the number of clusters, K, when MLC is applied to the blue crab data. “Boot p-value” is the p-value using the bootstrapping partial likelihood ratio test for H0 : K = K0 vs H1 : K = K0 + 1 where 99 bootstrapping samples are used. Criterion K BIC ICL NEC Boot p-value 1 476.13 – – 0 2 518.88 500.77 0.24 0.63 3 506.35 446.83 0.51 – males and 50 are females. Each crab has five measurements, FL, the width of frontal lip, RW, the rear width, CL, the length along the midline, CW, the maximum width of the carapace, and BD, the body depth in mm. As in Peel and McLachlan (2000), we are interested in the classification of male and female crabs. Peel and McLachlan (2000) had 18 crabs misclassified applying a mixture of two t distributions using all the five measurements. Using a mixture of two normal distributions resulted in one more crab being misclassified. Inspecting all pairwise scatter plots of the data set, RW and CL display two fairly well separated lines as in Figure 3.3. We apply MLC, MCLUST, LGA and MIXREG to these two measurements. For this data set, MCLUST chooses two clusters with different volumes, different shapes and different orientations by the BIC criterion. LGA chooses two clusters by the gap statistic. MLC also chooses two clusters by the various criteria in Table 3.3. MCLUST has 13 cases misclassified. LGA and MLC both have 7 cases misclassified. Inspecting the posterior probabilities of classification, MLC assigns roughly equal probabilities to the two clusters for two of the misclassified cases. In other words, there is extreme uncertainty about the assignment of these two cases, and their misclassifications are not surprising. There are no such diagnostic probabilities for LGA. MIXREG has eight cases misclassified when RW is regarded as the response; it has seven cases misclassified if CL is taken as the response. If we take logarithms of RW and CL, only five crabs are misclassified by LGA and MLC. MIXREG has five cases misclassified, if CL is taken as the response; it has 15 cases if RW is taken as the response. The results are summarized in Table 3.4. The clustering results are displayed in Figure 3.4.  39  Chapter 3. A partial likelihood approach to linear clustering  12  16  20  30  40  50  16  20  8  12  16  8  12  FL  35  45  8  RW  40  50  15  25  CL  14  18  20  30  CW  6  10  BD 8  12  16  20  15  25  35  45  6  10  14  18  Figure 3.3: Pairwise scatterplots of the blue crab data.  Table 3.4: Misclassification matrices of the blue crab data (using log scales) from MLC, MCLUST, LGA and MIXREG. MLC MCLUST LGA MIXREG MIXREG True class 1 2  1 50 5  2 0 45  1 50 14  2 0 36  1 50 5  2 0 45  (Response: CL)  (Response: RW)  1 50 5  1 38 3  2 0 45  2 12 47  40  Chapter 3. A partial likelihood approach to linear clustering  3.6 3.4 CL 3.2 3.0 2.8  2.8  3.0  3.2  CL  3.4  3.6  3.8  MCLUST  3.8  MLC  2.0  2.2  2.4 RW  2.6  2.8  2.0  2.2  2.6  2.8  2.6  2.8  3.6 3.4 CL 3.2 3.0 2.8  2.8  3.0  3.2  CL  3.4  3.6  3.8  MIXREG  3.8  LGA  2.4 RW  2.0  2.2  2.4 RW  2.6  2.8  2.0  2.2  2.4 RW  Figure 3.4: Clustering results of the blue crab data (using log scales) from MLC, MCLUST, LGA and MIXREG. CL is used as the response in MIXREG.  41  Chapter 3. A partial likelihood approach to linear clustering  3.7.4  Taqman single nucleotide polymorphism genotyping data  A single nucleotide polymorphism (SNP) is a single-base variation in a genome. The genetic code is specified by the four nucleotide “letters”: A (adenine), C (cytosine), T (thymine) and G (guanine). A SNP involves usually only two possibilities of these four letters, referred to as allele X and allele Y . An individual has two copies of genetic information passed from both parents. So the genotypes at this specific genetic position has three possible types: homozygote for allele X, homozygote for allele Y and heterozygote for both alleles. SNP genotyping, determination of genotypes of individuals, plays an increasing role in genetics studies. Taqman assay is a fluorescence-based high-throughput genotyping technology. Two fluorescent signals, x1 and x2 , called Allele.X and Allele.Y below, are used to detect the presence or absence of the two alleles, and their intensities determine the genotypes. (There is a third fluorescent signal presumably for normalizing the signal intensities. However, its effectiveness is suspicious in some situations. We choose to use signals without this normalization.) Blood samples of many individuals are assayed in wells of a plate simultaneously. For quality controls purpose, there are usually some blood samples with known genotypes, called positive controls, and there are some wells without genetic material, called negative controls. In the scatterplot of a SNP data set, there are typically four clusters: three clusters for the three possible genotypes and one for negative controls. In addition, there may be some failed samples. A clustering algorithm is usually employed to classify the samples. Algorithms commonly used include k-means algorithms (Olivier et al., 2002; Ranade et al., 2001) and modelbased algorithms (Fujisawa et al., 2004; Kang et al., 2004). Finding an algorithm that works for all plates/SNPs is a challenge. We fit a mixture of five components to the SNP genotyping data, for the three genotype clusters, the negative controls, and failed samples, respectively. The three genotype components are modelled linearly with the mixture likelihood (3.4). The fourth component is modelled with a bivariate normal distribution, and the fifth is a uniform component for possible outlying points. In the model fitting, we use the known labels for the negative controls. All the negative controls are assigned to the fourth component a priori, as are all points with both signals less than the corresponding maximum intensities of negative controls. Denote the set of these points by N . We 42  Chapter 3. A partial likelihood approach to linear clustering  Table 3.5: Largest membership probabilities of the Figure 3.5 by MCLUST and MLC. Points 1 2 3 4 5 MLC 0.99 1.00 1.00 0.99 0.98 MCLUST 0.96 0.93 1.00 1.00 1.00  7 points labelled in  6 0.86 0.97  7 0.99 0.50  do not use the known labels of the positive controls; these samples are used only for evaluating the clustering results. Hence, the modified mixture likelihood for the problem is L(κ, p1 , . . . , p5 |x1 , . . . , xn ) 3  = i∈N  ×  k=1  pk N (ak xi − bk ; 0, σk2 ) + p4 N (xi ; µ, Σ) + p5 U(xi ; R)  N (xi ; µ, Σ),  (3.19)  i∈N  where U is the uniform distribution on the rectangular region R given by (min(xi1 ), max(xi1 )) × (min(xi2 ), max(xi2 )). This mix of linear, elliptical and uniform components demonstrates the flexibility of our modelling approach. Figure 3.5 displays the clustering results of one plate using four different algorithms: MLC, MCLUST, LGA and MIXREG. Note that the majority of the points lie on three lines, hence the need for linear clustering. The results of modified MLC and MCLUST are displayed in the upper-left and the upper-right panels of Figure 3.5. All the points assigned to one of the three genotype clusters by MCLUST are assigned exactly the same way by MLC. Points 1 to 7, however, are assigned by MCLUST as background noise while MLC classifies all these points to one of the genotype clusters. Indeed, we know that these points should belong to the genotype clusters, and furthermore, points 1, 3, 6 and 7 are positive controls, all classified correctly by MLC. Table 3.5 displays the largest membership probabilities of these seven points by MCLUST and MLC, respectively. For the first six of the seven points, the two methods are very confident of their different assignments. We speculate here that, MCLUST does not place these points in the genotype clusters because they appear to be outliers with respect to the normality assumptions. In contrast, MLC classifies a point only using the proximity to a line and does not model the position along a line. 43  Chapter 3. A partial likelihood approach to linear clustering  MCLUST  3 4  1  5 6 7  2  3 4 6 7  0.0  0.5  0.5  Allele.Y  12  5  1.0  Allele.Y 1.0  1.5  1.5  2.0  2.0  MLC  0.0  0.5  1.0 Allele.X  1.5  0.5  1.5  1.5 Allele.Y 1.0 0.5  0.5  1.0  Allele.Y  1.5  2.0  MIXREG  2.0  LGA  1.0 Allele.X  0.5  1.0 Allele.X  1.5  0.5  1.0 Allele.X  1.5  Figure 3.5: Clustering results from four methods applied to a Taqman data set plate. ◦, + and denote the labels assigned for the three linear clusters (genotypes) by each method, × denotes points assigned as negative controls and denotes points assigned to the background cluster.  44  Chapter 3. A partial likelihood approach to linear clustering The lower-left panel of Figure 3.5 displays the results from LGA. Since LGA is formulated to deal with only linear clusters, we set the number of clusters to three. (If four or five clusters are chosen, the genotype clusters are mixed up.) The restriction to linear clusters means that LGA has difficulty separating the three genotype clusters from the negative controls and the failed samples. In the lower-right panel, we choose the signal for Allele Y as the response to implement MIXREG and set the number of regression lines to 3. MIXREG fails completely here. It finds only one of the three linear clusters. We obtained similar results if the signal for Allele X is chosen as the response.  3.8  Discussion  We have proposed a flexible model-based approach to linear clustering. It is flexible in the sense that only the deviation from a hyperplane is modelled parametrically; the position on the hyperplane is not modelled. The advantage of this approach is illustrated in the genotyping example, where the distribution along the line is complex and difficult to model. Furthermore, as was also illustrated in this example, we can incorporate elliptical clusters as necessary and a background cluster, borrowing from standard model-based clustering. Robustness to outliers is desirable, as the assumption of normal deviations around hyperplanes is sensitive to large deviations in the orthogonal direction. In addition to the inclusion of a uniform background cluster (Banfield and Raftery, 1993), one option would be to use a heavier tailed distribution, for example, Student’s t distribution with small degrees of freedom or with degrees of freedom depending on the data. This would adapt Peel and McLachlan (2000)’s EM algorithm for t mixture models from the elliptical context to the linear. The adaptation is straightforward but computationally more expensive. Further ideas include estimating the component covariance matrices in the M-step in a robust way, for example, trimming off some points. With a x = b, we are specifying a hyperplane in d − 1 dimensions. With little effort, this could be generalized to a mixture of partial likelihoods, each of which specifies a hyperplane of dimension q < d, n  K  l(κ, p|x1:n ) = i=1 k=1  pk N (Ak xi − bk ; 0, Σk ),  (3.20)  where A is of dimension d × (d − q), b is a vector of dimension d − q, 45  Chapter 3. A partial likelihood approach to linear clustering and Σk is a (d − q) × (d − q) covariance matrix for the deviation from the hyperplane. In the extreme case of a 0-dimension hyperplane, which is a point, we have the usual mixture of multivariate normal distributions. A mixture of components with various dimensions could be considered. A Bayesian version of this methodology would be helpful if some clusters are sparse but there is strong prior information about their approximate locations or properties (e.g., the parameters defining lines).  46  Bibliography Banfield, J. and A. Raftery (1993). Model-Based Gaussian and NonGaussian Clustering. Biometrics 49 (3), 803–821. Biernacki, C., G. Celeux, and G. Govaert (1999). An improvement of the NEC criterion for assessing the number of clusters in a mixture model. Pattern Recognition Letters 20 (3), 267–272. Biernacki, C., G. Celeux, and G. Govaert (2000). Assessing a mixture model for clustering with the integrated completed likelihood. IEEE Transactions on Pattern Analysis and Machine Intelligence 22 (7), 719–725. Bryant, P. and J. Williamson (1978). Asymptotic behaviour of classification maximum likelihood estimates. Biometrika 65 (2), 273–281. Campbell, N. and R. Mahon (1974). A multivariate study of variation in two species of rock crab of genus Leptograpsus. Australian Journal of Zoology 22 (3), 417–425. Celeux, G. and G. Govaert (1993). Comparison of the mixture and the classification maximum likelihood in cluster analysis. Journal of statistical computation and simulation(Print) 47 (3-4), 127–146. Celeux, G. and G. Soromenho (1996). An entropy criterion for assessing the number of clusters in a mixture model. Journal of Classification 13 (2), 195–212. Chen, H. and J. Chen (2001a). Large sample distribution of the likelihood ratio test for normal mixtures. Statistics and Probability Letters 52 (2), 125–133. Chen, H. and J. Chen (2001c). The Likelihood Ratio Test for Homogeneity in Finite Mixture Models. The Canadian Journal of Statistics/La Revue Canadienne de Statistique 29 (2), 201–215. Chen, H. and J. Chen (2001b). The likelihood ratio test for homogeneity in the finite mixture models. Canadian Journal of Statistics 29 (2), 201–215. 47  Bibliography Chen, H., J. Chen, and J. Kalbfleisch (2001). A modified likelihood ratio test for homogeneity in finite mixture models. Journal of the Royal Statistical Society: Series B (Methodological) 63 (1), 19–29. Chen, H., P. Meer, and D. Tyler (2001). Robust regression for data with multiple structures. 2001 IEEE Conference on Computer Vision and Pattern Recognition 1, 1069–1075. Chen, J. (1998). Penalized likelihood-ratio test for finite mixture models with multinomial observations. The Canadian Journal of Statistics 26 (4), 583–599. Chen, J. and J. Kalbfleisch (1996). Penalized Minimum-Distance Estimates in Finite Mixture Models. The Canadian Journal of Statistics/La Revue Canadienne de Statistique 24 (2), 167–175. Chen, J. and J. Kalbfleisch (2005). Modified likelihood ratio test in finite mixture models with a structural parameter. Journal of Statistical Planning and Inference 129 (1-2), 93–107. Ciuperca, G., A. Ridolfi, and J. Idier (2003). Penalized Maximum Likelihood Estimator for Normal Mixtures. Scandinavian Journal of Statistics 30 (1), 45–59. DeSarbo, W. and W. Cron (1988). A maximum likelihood methodology for clusterwise linear regression. Journal of Classification 5 (2), 249–282. Everitt, B. and D. Hand (1981). Finite mixture distributions. Monographs on Applied Probability and Statistics, London: Chapman and Hall . Everitt, B., S. Landau, and M. Leese (2001). Cluster Analysis. Hodder Arnold. Fraley, C. and A. Raftery (1998). How Many Clusters? Which Clustering Method? Answers Via Model-Based Cluster Analysis. The Computer Journal 41 (8), 578. Fraley, C. and A. Raftery (2002). Model-Based Clustering, Discriminant Analysis, and Density Estimation. Journal of the American Statistical Association 97 (458), 611–632. Fraley, C. and A. Raftery (2007). mclust: Model-Based Clustering / Normal Mixture Modeling. R package version 3.1-1. 48  Bibliography Fujisawa, H., S. Eguchi, M. Ushijima, S. Miyata, Y. Miki, T. Muto, and M. Matsuura (2004). Genotyping of single nucleotide polymorphism using model-based clustering. Bioinformatics 20, 5. Fuller, W. (1987). Measurement error models. Wiley Series in Probability and Mathematical Statistics, New York: Wiley, 1987 . Ganesalingam, S. (1989). Classification and Mixture Approaches to Clustering via Maximum Likelihood. Applied Statistics 38 (3), 455–466. Garc´ıa-Escudero, L., A. Gordaliza, R. San Mart´ın, S. van Aelst, and R. Zamar (2007). Robust linear clustering. Preprint. Harrington, J. (2007). lga: Tools for linear grouping analysis (LGA). R package version 1.0-0. Hathaway, R. (1983). Constrained maximum-likelihood estimation for a mixture of m univariate normal distributions. Ph. D. thesis, Rice University. Hathaway, R. (1985). A Constrained Formulation of Maximum-Likelihood Estimation for Normal Mixture Distributions. The Annals of Statistics 13 (2), 795–800. Hathaway, R. (1986). A constrained EM algorithm for univariate normal mixtures. Journal of Statistical Computation and Simulation 23 (3), 211– 230. Kang, H., Z. Qin, T. Niu, and J. Liu (2004). Incorporating Genotyping Uncertainty in Haplotype Inference for Single-Nucleotide Polymorphisms. The American Journal of Human Genetics 74 (3), 495–510. Keribin, C. (2000). Sankhya: The Indian Journal of Statistics 62, 49–66. MacLachlan, G. and D. Peel (2000). Finite mixture models. J. Wiley. Marriott, F. (1975). Separating mixtures of normal distributions. Biometrics 31 (3), 767–769. McLachlan, G. (1987). On Bootstrapping the Likelihood Ratio Test Stastistic for the Number of Components in a Normal Mixture. Applied Statistics 36 (3), 318–324.  49  Bibliography Olivier, M., L. Chuang, M. Chang, Y. Chen, D. Pei, K. Ranade, A. de Witte, J. Allen, N. Tran, D. Curb, et al. (2002). High-throughput genotyping of single nucleotide polymorphisms using new biplex invader technology. Nucleic Acids Research 30 (12), e53. Peel, D. and G. McLachlan (2000). Robust mixture modelling using the t distribution. Statistics and Computing 10 (4), 339–348. Ranade, K., M. Chang, C. Ting, D. Pei, C. Hsiao, M. Olivier, R. Pesich, J. Hebert, Y. Chen, V. Dzau, et al. (2001). High-Throughput Genotyping with Single Nucleotide Polymorphisms. Genome Research 11 (7), 1262– 1268. Redner, R. (1981). Note on the Consistency of the Maximum Likelihood Estimate for Nonidentifiable Distributions. The Annals of Statistics 9 (1), 225–228. Scott, A. and M. Symons (1971). Clustering Methods Based on Likelihood Ratio Criteria. Biometrics 27 (2), 387–397. Smyth, P. (2000). Model selection for probabilistic clustering using crossvalidated likelihood. Statistics and Computing 10 (1), 63–72. Spath, H. (1982). Fast algorithm for clusterwise linear regression. COMPUT. 29 (2), 175–181. Thode Jr, H., S. Finch, and N. Mendell (1988). Simulated Percentage Points for the Null Distribution of the Likelihood Ratio Test for a Mixture of Two Normals. Biometrics 44 (4), 1195–1201. Turner, T. (2000). Estimating the propagation rate of a viral infection of potato plants via mixtures of regressions. Journal of the Royal Statistical Society Series C(Applied Statistics) 49 (3), 371–384. Turner, T. (2006). mixreg: Functions to fit mixtures of regressions. R package version 0.0-2. Van Aelst, S., X. Wang, R. Zamar, and R. Zhu (2006). Linear grouping using orthogonal regression. Computational Statistics and Data Analysis 50 (5), 1287–1312. Wald, A. (1949). Note on the Consistency of the Maximum Likelihood Estimate. The Annals of Mathematical Statistics 20 (4), 595–601.  50  Chapter 4  Bayesian linear clustering with application to single nucleotide polymorphism genotyping 4.1  Introduction  This paper was motivated by a clustering problem in single nucleotide polymorphism (SNP) genotyping. A single nucleotide polymorphism (SNP, pronounced as “snip”) is a single-base variation in a genome. The genetic code of life is specified by the four nucleotide “letters”: A (adenine), C (cytosine), G (guanine) and T (thymine). There are two complementary DNA strands. It is sufficient to consider only one. SNP variation occurs when a single nucleotide, such as an A, is replaced by one of the other three letters C, G or T. One SNP usually involves only two letters, referred to generically throughout this paper as allele X and allele Y . An individual has two copies of genetic information passed from both parents. So the genotypes at a specific genetic position have three possibilities: homozygote for allele X, homozygote for allele Y and heterozygote for both alleles. SNP genotyping, determination of genotypes of individuals, plays an increasing role in genetics studies. For a small project, it is possible to make genotype calls manually. In most cases, it is not hard for an expert to perform this job, and the “eyeballing” procedure usually gives reasonable results due to its sophisticated incorporation of prior information. For large-scale studies, however, manual A version of this chapter will be submitted for publication. Authors: Guohua Yan, William J. Welch and Ruben H. Zamar.  51  Chapter 4. Bayesian linear clustering scoring can become a daunting challenge. Typical SNP genotyping applications involve thousands of patients and hundreds of SNPs. Hence, reliable automated genotyping methods are highly needed. TaqMan SNP Genotyping Assay (Applied Biosystems) is a fluorescencebased high-throughput genotyping technology. Blood samples of many individuals are arranged in a 96- or 384-well plate and are assayed simultaneously. Two fluorescent signals, x1 and x2 , also called Allele.X and Allele.Y below, are used to detect the presence or absence of the two alleles. There is a third fluorescent signal presumably for normalizing the signal intensities. However, its effectiveness is suspicious in some situations. See Chapter 2. For quality controls purpose, some blood samples with known genotypes are also included in a plate; these samples are called positive controls. As well, some wells do not have genetic material; these are called negative controls. In the scatterplot of a SNP data set, there are typically four clusters, as in Figure 4.1. In the lower-right, upper-left, and upper-right corners are three clusters, presumably containing samples of wild-type homozygotes, variant homozygotes, and heterozygotes, respectively. In the lower-left corner, the cluster may contain negative controls and/or some failed samples. A clustering algorithm is usually employed to call the SNP genotypes. Algorithms commonly used in the literature include k-means (Olivier et al., 2002; Ranade et al., 2001) and model-based algorithms (Fujisawa et al., 2004; Kang et al., 2004). The proprietary software Sequence Detection System is included in a thermal cycler. Usually human intervention of an expert operator is usually needed to review the genotype calling results. Finding an algorithm that works well for all plates/SNPs is a challenge, especially when the variant allele homozygous genotype cluster is sparse, in which standard clustering algorithms often fail. Several clustering algorithms k-means, MCLUST (Fraley and AE, 2006), LGA (Harrington, 2007), MIXREG (Turner, 2006) and MLC (Yan et al., 2008) are applied to the genotyping data in Figure 4.1. Figure 4.2 displays their corresponding clustering results. We can see that all of these algorithms fail to identify the five points in the upper-left corner as variant allele homozygotes while an expert operator would do with confidence. These algorithms fail because the variant allele homozygous genotype cluster is very sparse. This motivates us to adopt a Bayesian approach to incorporate available prior information. In addition, we found that it is convenient to model the genotype clusters as linear structures. Most clustering methods and algorithms cluster data points around “centers”. However, some data sets, as in the SNP genotype setting, form groups 52  0.5  Allele.Y 1.0  1.5  2.0  Chapter 4. Bayesian linear clustering  0.2  0.4  0.6 0.8 Allele.X  1.0  1.2  1.4  Figure 4.1: Scatterplot of one plate in which the variant allele homozygous cluster has only five points (upper-left). through linear relationships and standard clustering techniques are usually not able to find these linear patterns. By linear clustering, we mean detecting these linear clusters in a data set. Van Aelst et al. (2006) proposed a method called Linear Grouping Algorithm (LGA) for linear clustering problems without assuming a response variable. Yan et al. (2008) proposed a partial likelihood approach which models only the signed orthogonal distance from each data point to a hypothesized hyperplane and treats all variables symmetrically. In this paper we introduce a hierarchical modeling approach, which is particularly appropriate for identifying sparse linear clusters. We show that the sparse cluster in our SNP genotyping dataset can be successfully identified after a careful specification of the prior distributions. The rest of this paper is organized as follows. In Section 2, we describe a hierarchical model framework for linear clustering. In Section 3, we discuss in details sampling issues. Our approach to linear clustering is illustrated with a relatively simple dataset (crab dataset) in Section 4. In Section 5, we revisit the SNP genotyping motivating dataset and show that a careful specification of the prior distributions is critical for the success of the clustering algorithm. A brief discussion follows in Section 6.  53  Chapter 4. Bayesian linear clustering  k−means 2.0  3  3  4 4 44 44 444 4 4 4 4 44  0.2  1.5  3  0.4  0.6 0.8 Allele.X  3  1  1 111 1 11111111 1111 111 11 1 1111 1 11111 1 11 11 1111 1 11 11 1111 11 1 1 1 1 1 11 1 1 1 1 1 1 1 11 1 2 2 2 22 2 2 222222 1 22 222 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2222 22 2 2 22222222 222222222 2 222222 222 2  4  3  3  3 3 33 3 33 3 3333 3 33 3 33 33 3 3333 3 33333333 33 33 3 3333 3 33 3333 33 3 3 33 33 333 3 3 3 33 3 33 4 2  Allele.Y 1.0  2.0 Allele.Y 1.0  1.5  3  1.0  1.2  0.5  3  3  0.5  MCLUST  3  3 4 44 44 444 4 4 4 4 44  1.4  3 3 33 3 33 3 3333 3 33 3 33 33 3 3333 3 33333333 33 33 3 3333 3 33 3333 33 3 3 33 33 333 3 3 3 33 3 33 3 3  3  0.2  0.4  LGA 2.0  2 2 22 22222 2 2 2 2 22  0.2  0.4  0.6 0.8 Allele.X  1.5 3  1  1 1 11 11111 11111111 1111 1111111 3 111 111 11 11 1 111 11 311113 11 1 11 13 11 1111 1 1 1 3 3 3 3 33 1 1 1 1 3 1 1 1 1 1 111 11 11 11 1 111 1 1 11 1 11 1111 11 111 111 1 1111 111 111 111 111111111 1 111111 111 1  2  3 3  Allele.Y 1.0  2.0 1.5 Allele.Y 1.0 0.5  2  1.0  1.2  1.4  3  3  2 2 22 2 22 2 2222 2 22 2 22 22 2 2222 2 22222222 22 22 2 2222 2 32 2222 33 2 2 33 2222 22 3 2 22 2 22 2 3  1.0  1.2  1.4  0.5  3 3  0.6 0.8 Allele.X  MIXREG/MLC  3  3  3  3 2 22 2222 2222 2 222222 22 22 21122 1 2 11121 21 11 22 111 1111 11 2111 11 1 12 1 1 1 11 1 1 1 1 1 1 1 11 1 1 1 1 11 1 1 11111 1 11 111 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1113 1 11 1 1 11111111 111111111 1 111111 111 1  2 2 22 22222 2 2 2 2 22  2 2 22 2 22 2 2222 2 22 3 22 22 2 2222 2 22222222 22 22 2 2222 2 22 2222 22 2 2 22 2222 22 3 2 22 2 22 2 3  2  0.2  3  3 1 11 11111 11111111 1111 1111111 3 11131 111 11 11 1 111 11 11111 1 1 1 1111 111 11 1 1 1 1 1 1 11 1 1 1 1 1 1 1 1 1 1111 111 11 11 11 1 111 1 1 11 1 11 11 111 111 3 1113 111 111 111 111111111 1 111111 111 1  0.4  0.6 0.8 Allele.X  1.0  1.2  1.4  Figure 4.2: Clustering results of several clustering algorithms. For k-means and MCLUST, the number of clusters are set to 4; for the remaining algorithms, the number of lines are set to 3.  54  Chapter 4. Bayesian linear clustering  4.2 4.2.1  Model specification Model for one cluster: orthogonal regression  We first introduce a Bayesian approach for orthogonal regression (see e.g., Fuller, 1987). To model a single cluster, assume that the vector-valued data x1 , . . . , xn are drawn by a random mechanism represented by a random vector X in a d-dimensional space. We do not model the distribution of X except for its deviation from an unknown hyperplane. Let {x : a x − b = 0} be a hyperplane. For identifiability purpose, we assume that a a = 1 and that the first nonzero element of a is positive. Our prior information about this hyperplane is summarized by a prior distribution π(a, b). Given a and b, the signed orthogonal deviation a X−b from the hyperplane {x : a x − b = 0} has a density function p(· |σ). Let π(σ) be the prior distribution for σ. Denote θ = (a , b, σ) . The posterior distribution of θ is n  π(θ|x1 , . . . , xn ) ∝  i=1  p(a xi − b|σ) π(a, b)π(σ).  (4.1)  ¯ and S be the sample mean and sample covariance matrix of Let x x1 , . . . , xn respectively. If we take p(·|σ) to be the normal density function N (· |0, σ 2 ) and set π(a, b) ∝ 1 and π(σ) ∝ 1, the maximum a posteriori ˆ is the estimator of θ is (ˆ a , ˆb, σ ˆ ) , where σ ˆ 2 is the smallest eigenvalue of S, a ¯ . This is a conˆx (standardized) eigenvector associated with σ ˆ 2 and ˆb = a nection with orthogonal regression, see for example, Fuller (1987); for the ˆ when X is normal, distributions of the eigenvalue σ ˆ 2 and the eigenvector a see Anderson (2003).  4.2.2  Model for linear clustering  Now we assume that the data x1 , . . . , xn are drawn by a more complicated random mechanism, still represented by random vector X in a d-dimensional space which we do not model fully. The data now lie around K hyperplanes {x : ak x = bk }, k = 1, . . . , K, for which our prior knowledge is summarized in π(a1 , b1 , . . . , aK , bK ). Let Z = (Z1 , . . . , ZK ) be a random vector indicating these hyperplanes, Zk = 1 with probability pk for k = 1, . . . , K, Zk = 0 or 1 and K k=1 Zk = 1. Let p = (p1 , . . . , pK ) . We denote by π(p) the prior distribution for p. We assume that, conditional on Zk = 1, the signed orthogonal deviation ak X − bk has a density function p(·|σk ) for k = 1, . . . , K. Let π(σ1 , . . . , σK ) be the prior distribution for (σ1 , . . . , σK ) . Let z1 , . . . , zn be the corresponding unobservable indicators for the data x1 , . . . , xn , where 55  Chapter 4. Bayesian linear clustering zi = (zi1 , . . . , ziK ) for i = 1, . . . , n. Let θ be the collection of all the parameters, θ = (a1 , b1 , σ1 , . . . , aK , bK , σK , p ) . Formally, the posterior distribution of the unobservable cluster indicators z1 , . . ., zn and θ is π(z1 , . . . , zn , θ|x1 , . . . , xn ) K  n  [pk p(ak xi − bk |σk )]zik  ∝  π(a1 , b1 , . . . , aK , bK )π(σ1 , . . . , σK )π(p).  i=1 k=1  (4.2) This formula is the basis for our linear clustering, in which we are primarily interested in the cluster indicators z1 , . . . , zn and the elements of θ are nuisance parameters. A Gibbs sampling scheme based on this formula is detailed in next section. In addition, if we sum out the indicators z1 , . . . , zn , the marginal distribution for θ is π(θ|x1 , . . . , xn ) n  ∝  K  i=1 k=1  pk p(ak xi − bk |σk ) π(a1 , b1 , . . . , aK , bK )π(σ1 , . . . , σK )π(p). (4.3)  This marginal distribution can be used if we would like to set up other sampling schemes such as Metropolis-Hastings, tempering MCMC or SMC methods. Once a MCMC sample for θ is available, it is straightforward to sample cluster labels at each iteration. At iteration t, with θ (t) , sample zi from multinomial distribution (t)  (t)  (t)  (t)  Pr(zi = k) ∝ pk p((a )k xi − bk |σk ). (1)  (4.4) (T )  After a sequence of iterations, we have a sample (zi , . . . , zi ) from posterior marginal distribution π(zi |x1 , . . . , xn ), for i = 1, . . . , n. The point xi can be classified into the cluster with the largest frequency in the posterior sample; and the largest frequency provides a good measure of the cluster membership.  56  Chapter 4. Bayesian linear clustering  4.3  Sampling algorithms  4.3.1  Gibbs sampling  The model framework in (4.2) is conceptually straightforward; the posterior sampling, however, can be difficult, as in any Bayesian mixture modelling (Celeux et al., 2000). We present here a sampling algorithm for normal orthogonal deviations and, to make the algorithm more efficient, we integrate some parameters by using conjugate priors. For ease of presentation, we assume the data are in a 2-dimensional space. The density functions for normal orthogonal deviations are p(ak x − bk |σk ) = N (ak x − bk ; 0, σk2 ),  (4.5)  for k = 1, . . . , K, where N (· ; 0, σ 2 ) is the normal density function with mean 2 ) follow independent inverse 0 and variance σk2 . The priors for (σ12 , . . . , σK Gamma distributions, K 2 π(σ12 , . . . , σK )  IG(σk2 ; δ1k , δ2k ),  =  (4.6)  k=1  where IG(· ; δk1 , δk2 ) is the density function of inverse Gamma with shape parameter δk1 and rate parameter δk2 . In the case of equal variances, σ12 = 2 = σ 2 , we assume . . . = σK π(σ 2 ) = IG(σ 2 ; δ1 , δ2 ),  (4.7)  For b1 , . . . , bK , we set K 2 π(b1 , . . . , bK |σ12 , . . . , σK )=  N (bk ; κ1k , σk2 /κ2k ).  (4.8)  k=1  With these specifications, we consider the following special case of the posterior distribution (4.2), π(z1 , . . . , zn , θ|x1 , . . . , xn ) n  ∝ ×  K  i=1 k=1 K  [pk N (ak xi − bk ; 0, σk2 )]zik  N (bk ; κ1k , σk2 /κ2k )IG(σk2 ; δ1k , δ2k ) π(a1 , . . . , aK )π(p).  (4.9)  k=1  57  Chapter 4. Bayesian linear clustering Conditional on θ and the data x1 , . . . , xn , z1 , . . . , zn are independently distributed, 2 zi |(θ, xi ) ∼ M (1, (p1 N (a1 xi − b1 ; 0, σ12 ), . . . , pK N (aK xi − bK ; 0, σK )) ), (4.10) for i = 1, . . . , n, where M (n, p) is the probability mass function for multinomial distribution with n trials and probability vector p. Let n  nk = i=1  and  1 ¯k = zik , x nk  n  zik xi ,  (4.11)  ¯ k )(xi − x ¯k) . zik (xi − x  (4.12)  i=1  n  Ak = i=1  The full conditional distribution of bk is bk |e.e. = bk |(z1 , . . . , zn , ak , σk2 , x1 , . . . , xn ) ∼ N(  ¯k σk2 κ1k κ2k + nk ak x , ), nk + κ2k nk + κ2k  (4.13)  for k = 1, . . . , K, where “e.e.” stands for “everything else”. We integrate out b1 , . . . , bK from (4.9) and get 2 π(z1 , . . . , zn , a1 , σ12 , . . . , aK , σK , p|x1 , . . . , xn ) K  ∝ ×  pnk k (σk2 )−nk /2 (nk + κ2k )−1/2 exp(− k=1 K  ∗ δ2k ) 2σk2  IG(σk2 ; δ1k , δ2k ) π(a1 , . . . , aK )π(p),  (4.14)  k=1  where ∗ ¯ k )2 + κ21k κ2k − δ2k = ak Ak ak + nk (ak x  ¯ k )2 (κ1k κ2k + nk ak x . nk + κ2k  (4.15)  The distribution of σk2 conditioning on everything else except bk is σk2 |(z1 , . . . , zn , ak , x1 , . . . , xn ) ∼ IG(δ1k +  δ∗ nk , δ2k + 2k ). 2 2  (4.16)  58  Chapter 4. Bayesian linear clustering 2 = σ 2 , then the distribution of σ 2 condiIf we assume that σ12 = . . . = σK tioning on z1 , . . . , zn , a1 , . . . , aK , p is  n σ |(z1 , . . . , zn , a1 , . . . , aK , x1 , . . . , xn ) ∼ IG(δ1 + , δ2 + 2  K  2  k=1  ∗ δ2k ). 2  (4.17)  2 from (4.14) and get We further integrate out σ12 , . . . , σK  π(z1 , . . . , zn , a1 , . . . , aK , p|x1 , . . . , xn ) K  pnk k (nk + κ2k )−1/2 Γ(δ1k +  ∝  k=1  nk δ∗ nk )(δ2k + 2k )−(δ1k + 2 ) 2 2  × π(a1 , . . . , aK )π(p).  (4.18)  For proportions p, we use a Dirichlet prior, π(p) = D(p; α),  (4.19)  where D(· ; α) is the Dirichlet density function with parameter α = (α1 , . . ., αK ) . The full conditional distribution of p is p|e.e. = p|(z1 , . . . , zn ) ∼ D(α∗ ),  (4.20)  ∗ ) with α∗ = α + n for k = 1, . . . , K. Therefore, where α∗ = (α1∗ , . . . , αK k k  π(z1 , . . . , zn , a1 , . . . , aK |x1 , . . . , xn ) ∝  1 Γ(  K  K ∗ i=1 αk ) k=1  Γ(αk∗ )(nk + κ2k )−1/2 Γ(δ1k +  nk δ∗ nk )(δ2k + 2k )−(δ1k + 2 ) 2 2  × π(a1 , . . . , aK ).  (4.21)  Note that ak ak = 1 for k = 1, . . . , K. We re-parameterize ak as ak = (  1 1 + a2k  ,  ak 1 + a2k  ),  for k = 1, . . . , K and use normal priors for (a1 , . . . , aK ) , K 2 N (ak ; ν1k , ν2k ).  π(a1 , . . . , aK ) =  (4.22)  k=1  59  Chapter 4. Bayesian linear clustering The conditional distribution of ak on cluster indicators is nk δ∗ nk 2 ). )(δ2k + 2k )−(δ1k + 2 ) N (ak ; ν1k , ν2k 2 2 (4.23) From the above, discussion, a possible sampling algorithm is as follows.  π(ak |z1 , . . . , zn , x1 , . . . , xn ) ∝ Γ(δ1k +  Algorithm 1: 1. Initialize the algorithm with θ (0) . 2. At iteration t, (t)  (t)  (a) Sample cluster indicators z1 , . . . , zn from multinomial distribution (4.10), where θ is replaced with θ (t−1) . (b)  i. Sample p(t) from Dirichlet distribution (4.20), where the clus(t) (t) ter indicators are replaced with z1 , . . ., zn . (t)  (t)  ii. Sample a1 , . . . , aK using a random walk Metropolis-Hastings scheme from (4.23), where cluster indicators are re(t) (t) placed by z1 , . . ., zn , respectively. 2 )(t) from inverse Gamma distribution iii. Sample (σ12 )(t) , . . . , (σK (t) (4.16), where z1 , . . . , zn , a1 , . . . , aK are replaced by z1 , . . ., (t) (t) (t) zn , a1 , . . ., aK , respectively; in the case of equal variances 2 = σ 2 , sample (σ 2 )(t) from (4.17). σ12 = . . . = σK iv. Sample b1 , . . . , bK from normal distribution (4.13), with all parameters and cluster indicators replaced with those in the current iteration.  4.3.2  Metropolis-Hastings sampling  Alternatively, we can sum out the indicators z1 , . . . , zn and sample θ from the marginal posterior distribution π(θ|x1 , . . . , xn ) in (4.3). Standard random walk Metropolis-Hastings algorithms can be used directly to sample from the marginal posterior distribution π(θ|x1 , . . . , xn ). However, we found that it is advantageous to reparameterize θ so that the constraints such as K k=1 pk = 1 and a a = 1 are implicitly incorporated (see Section 4.5). Since the proposal distribution and the posterior distribution have the same support, the sampling schemes on the transformed parameters are more efficient. Suppose we reparameterize θ as θ = φ(η), 60  Chapter 4. Bayesian linear clustering such that φ is a one-to-one map and that the domain of η is the dimension of η. Let J be the Jacobian J=  q,  where q is  ∂θ , ∂η  with the understanding that only q free components of θ are used. Then the marginal posterior distribution of η is π(η|x1 , . . . , xn ) = π(θ|x1 , . . . , xn )|J|.  (4.24)  To sample η, we first run a Gibbs sampler, within which a random walk Metropolis step is used if a full conditional distribution is not of an explicit form. Algorithm 2.0: 1. Initialize the algorithm with η (0) . (t−1)  2. At iteration t, update η1 ponent,  (t−1)  , . . . , ηq  (t−1)  (a) Sample ν from N (ηj  , s2j ). (t)  (b) Compute r = min{1,  sequentially. For the jth com-  (t)  (t−1)  (t−1)  π(η1 ,...,ηj−1 ,ν,ηj+1 ,...,ηq  |x1 ,...,xn ) (t) (t) (t−1) (t−1) (t−1) π(η1 ,...,ηj−1 ,ηj ,ηj+1 ,...,ηq |x1 ,...,xn )  }.  (t)  (c) Update ηj = ν with probability r. ˆ of the covariance After a burn-in period, we can get a crude estimate Σ matrix of η. In the second stage of the algorithm, we update η using a Metropolis-Hastings sampler in one block. Algorithm 2: 1. Initialize the algorithm with η (0) . 2. At iteration t, ˆ (a) Sample ν from N (η (t−1) , s2 Σ). 1 ,...,xn ) (b) Compute r = min{1, π(ηπ(ν|x (t−1) |x ,...,x ) }. 1  n  (c) Update η (t) = ν with probability r.  Here we use an adaptive strategy; see page 307, Gelman et al. (2004). 61  Chapter 4. Bayesian linear clustering  4.4  Australia rock crab data  Now we consider the rock crab data set of Campbell and Mahon (1974) on the genus Leptograpsus. We focus on the 100 blue crabs, 50 of which are males and 50 are females. Each crab has five measurements, FL, the width of frontal lip, RW, the rear width, CL, the length along the midline, CW, the maximum width of the carapace, and BD, the body depth in mm. As in Peel and McLachlan (2000), we are interested in the classification of male and female crabs. Peel and McLachlan (2000) had 18 crabs misclassified applying a mixture of two t distributions using all the five measurements. Using a mixture of two normal distributions resulted in one more crab being misclassified. Inspecting all pairwise scatter plots of the data set, RW and CL display two fairly well separated lines as in Figure 4.3. We shall classify these crabs using only RW and CL. The number of clusters is set to 2. We run the Gibbs sampler “Algorithm 1” based on (4.2). To initialize the sampler, we adopt the strategy of Van Aelst et al. (2006). We randomly select K = 2 mutually exclusive subsets of d + 1 = 3 observations from x1 , . . . , xn . For k = 1, . . . , K, we compute the maxi(0) (0) mum partial likelihood estimates ak , bk and (σk2 )(0) from the kth subset of d + 1 observations (see Subsection 2 of Yan et al. (2008)). The (0) proportions pk are all set as 1/K. The initial values for θ are then (0) (0) 2 )(0) , p(0) , . . . , p(0) ) . We then θ (0) = ((a1 )(0) , b1 , (σ12 )(0) , . . . , (aK )(0) , bK , (σK 1 K monitor the evolution of cluster indicators z1 , . . . , zn . If the number of points in one cluster has fallen below d + 1, the sampler will get stuck and never recover from it; we then re-initialize the sampler with the above strategy. With the above initialization strategy, we actually impose a constraint on the distribution of cluster labels such that any allocation leading to less than d + 1 points for a cluster is not allowed, i.e., has a probability of zero. As a result, we can then specify the prior distributions of component parameters with vague priors. Specifically, we set the priors as follows. ak ∼ N (0, 1000002 ), π(bk ) ∝ 1, σk2 ∼ IG(0.0001, 0.0001), pk ∝ 1, for k = 1, 2. To alleviate the effect of autocorrelations, we keep only one simulated point for every ten iterations. Hereafter the iteration numbers are based on the “thinned” samples. Figure 4.4 displays the evolution of sampled values of θ for 11, 000 iterations. Figure 4.5 displays the evolution of log unnormalized posterior density values of sampled values for (θ, z1 , . . ., zn ) for 11, 000 iterations. 62  Chapter 4. Bayesian linear clustering  12  16  20  30  40  50  16  20  8  12  16  8  12  FL  35  45  8  RW  40  50  15  25  CL  14  18  20  30  CW  6  10  BD 8  12  16  20  15  25  35  45  6  10  14  18  Figure 4.3: Pairwise scatter plots of the blue crab data.  63  Chapter 4. Bayesian linear clustering From these two figures we found no reason to suspect that the sampler has reached a stationary distribution. Figure 4.6 displays the autocorrelations of sampled values for θ for 11, 000 iterations. The autocorrelations decay fast, especially for variances σ12 and σ22 and proportions p1 and p2 . Further checking the density curves in Figure 4.7 of sampled values for θ for 10, 000 iterations after a burn-in of 1000 iterations, we observe that all the estimated density curves are roughly unimodal. There is no obvious reason for the existence of minor modes. In Figure 4.8, we demonstrate the isolation of the two modes of the marginal posterior distribution of θ. One mode is estimated by maximum a posteriori estimation; the other one is obtained by permutation of cluster indices. The horizontal axis t ranges from 0 to 1, the vertical axis is the unnormalized marginal posterior density ˆ 1 + tθ ˆ 2 , where θ ˆ 1 and θ ˆ 2 are the two estimated modes. From value of (1 − t)θ this plot, we can see that the two modes are very peaky and isolated. It is then understandable that label-switching does not occur. While Gibbs sampling is generally criticized for not traversing the whole support of the posterior distribution for θ, it is reasonable to think that it has fully explored the support around one isolated mode in this data set. Figure 4.9 shows the posterior probabilities of being classified into Class 1 of the 100 crabs estimated from 10, 000 iterations using the Gibbs sampler Algorithm 1, after a burn-in of 1000 iterations. If we were to classify all the crabs, i.e., classifying a crab to the class with larger posterior probability, five crabs are misclassified, as indicated in the upper-left corner of Figure 4.9. The number of misclassification is the same as that of the partial likelihood approach in Chapter 3; this result is expected as we used essentially noninformative priors in this example. When informative priors are necessary, as in our next example, the Bayesian approach does have an advantage. We have also run the adaptive Metropolis-Hastings algorithm for the simulation of marginal posterior distribution for θ and get very similar results. It runs much slower than the Gibbs sampler though.  4.5  Taqman single nucleotide polymorphism genotyping data  Now we analyze the SNP genotyping data in Figure 4.1 using the linear clustering strategy. We fit a mixture of five components to the SNP genotyping data, to represent the three genotype clusters, the negative controls, and the failed samples, respectively. The three genotype components are modelled 64  Chapter 4. Bayesian linear clustering  a2  −1.50  −1.20  a2 −1.10  a1 −1.35  −1.00  −1.20  a1  0  2000  4000  6000 Iteration  8000  10000  0  2000  4000  8000  10000  8000  10000  8000  10000  8000  10000  b2  −0.15  −0.25  b2 −0.05  b1 −0.15  b1  6000 Iteration  0  2000  4000  6000 Iteration  8000  10000  0  2000  4000  σ22  0.00005  1e−04  σ2 2  σ2 1 0.00015  3e−04  σ21  6000 Iteration  0  2000  4000  6000 Iteration  8000  10000  0  2000  4000  p2  p2 0.4 0.2  0.3  p1 0.5  0.6  0.7  p1  6000 Iteration  0  2000  4000  6000 Iteration  8000  10000  0  2000  4000  6000 Iteration  Figure 4.4: Evolution of sampled values for θ for 11, 000 iterations using Algorithm 1 for blue crab data. 65  170  Log unnormalized posterior density 180 190 200 210  Chapter 4. Bayesian linear clustering  0  2000  4000  6000 Iteration  8000  10000  Figure 4.5: Evolution of log unnormalized posterior density of sampled values for (θ, z1 , . . ., zn ) for 11, 000 iterations using Algorithm 1 for blue crab data. linearly. The fourth component is modelled with a bivariate distribution with diagonal covariance matrix, and the fifth is a uniform component for possible outlying points. In the model fitting, we use the known labels for the negative controls. All the negative controls are assigned to the fourth component a priori, as are all points with both signals less than the corresponding maximum intensities of negative controls. Denote the set of these points by N . Hence, the modified mixture likelihood for the problem is L(θ|x1 , . . . , xn )    3  1 xi1 + ak xi2 bk pk t 2 ( = − ) + p4 t2 (xi ; µ, Σ) + p5 U(xi ; R)   σk σ 1 + a2 σk i∈N  ×  k=1  t2 (xi ; µ, Σ),  k  k  (4.25)  i∈N  where t2 is the density for Student’s t distribution with 2 degrees of freedom, σ11 0 and U is the uniform distribution on the rectangular reΣ= 0 σ22 gion R given by (min(xi1 ), max(xi1 )) × (min(xi2 ), max(xi2 )). For robustness consideration, we use Student’s t distribution for orthogonal deviation of a point from a line and also for the elliptical distribution for negative control 66  Chapter 4. Bayesian linear clustering  a2  ACF 0.4 0.0  0.0  ACF 0.4  0.8  0.8  a1  0  10  20 Lag  30  40  0  10  30  40  30  40  30  40  30  40  b2  ACF 0.4 0.0  0.0  ACF 0.4  0.8  0.8  b1  20 Lag  0  10  20 Lag  30  40  0  10  σ22  ACF 0.4  0.0  0.0  ACF 0.4  0.8  0.8  σ21  20 Lag  0  10  20 Lag  30  40  0  10  ACF 0.4 0.0  0.0  ACF 0.4  0.8  p2  0.8  p1  20 Lag  0  10  20 Lag  30  40  0  10  20 Lag  Figure 4.6: Autocorrelations of sampled values for θ for 11, 000 iterations using Algorithm 1 for blue crab data. 67  Chapter 4. Bayesian linear clustering  a2  0  0  2  2  Density 4 6 8  Density 4 6 8  10  10  a1  −1.5  −1.4  −1.3  −1.2  −1.25 −1.20 −1.15 −1.10 −1.05 −1.00 −0.95  b2  0  0  5  Density 5 10  Density 10 15  15  20  b1  −0.25  −0.20  −0.15  −0.10  −0.15  −0.10  0.00  σ22  Density 5000 0  0  Density 10000  15000  20000  σ21  −0.05  0.00005  0.00010  0.00015  0.00020  0.00005  0.00015  0.00035  6 Density 3 4 5 2 1 0  0  1  2  Density 3 4 5  6  7  p2  7  p1  0.00025  0.3  0.4  0.5  0.6  0.7  0.8  0.2  0.3  0.4  0.5  0.6  0.7  Figure 4.7: Density curves of sampled values for θ for 10, 000 iterations using Algorithm 1, after a burn-in of 1000 iterations for blue crab data. 68  0.0e+00  Unnormalized density 5.0e+110 1.0e+111  1.5e+111  Chapter 4. Bayesian linear clustering  0.0  0.2  0.4  0.6  0.8  1.0  t  Estimated posterior probability of being in Class 1 0.0 0.2 0.4 0.6 0.8 1.0  Figure 4.8: Unnormalized marginal posterior density values of θ along a ˆ1 + θ ˆ 2 connecting the two estimated modes θ ˆ 1 and θ ˆ2. line segment (1 − t)θ One mode is estimated by maximum a posteriori estimation from 10, 000 iterations using Algorithm 1, after a burn-in of 1000 iterations and the other is obtained by permutation. The blue crab data is used.  0  20  40 60 Index of Observations  80  100  Figure 4.9: Estimated posterior probability of being classified into Class 1 of the 100 crabs from 10, 000 iterations using Algorithm 1, after a burn-in of 1000 iterations. The solid vertical line separates these crabs by their true classes and the horizontal dashed line corresponds to a probability of 0.5.  69  Chapter 4. Bayesian linear clustering points. This mix of linear, elliptical and uniform components demonstrates the flexibility of our modelling approach. Based on the subject matter knowledge, we add various constraints in the specification of priors. For the slopes −1/ak of lines, 0<−  1 1 1 1 1 1 1 <− <− , − + >− + , a1 a2 a3 a3 a2 a2 a1  (4.26)  which maintains the order of positive slopes of the three genotype clusters and requires that the “gap” between the two lines of the heterozygotic clusters and the variant allele homozygotic clusters cannot be too small relative to the “gap” between the lines of the wild type allele homozygotic clusters and the heterozygotic clusters. These constraints on slopes are natural for SNP genotyping data. To implement Algorithm 2, we reparameterize (a1 , a2 , a3 ) into (η1 , η2 , η3 ),   a1 = −1/ exp(η1 ), a2 = −1/(exp(η1 ) + exp(η2 )), (4.27)  a3 = −1/(exp(η1 ) + 2 exp(η2 ) + exp(η3 )). It is obvious that the constraint (4.26) is satisfied. In the case of few points in the variant allele homozygotic cluster, the orthogonal variance σ32 may not be estimated reliably. Hence we require σ32 = σ12 ,  (4.28)  and then apply a log transformation to (σ12 , σ22 ), σ12 = exp(η7 ),  σ22 = exp(η8 ),  (4.29)  For variances (σ11 , σ22 ), a log transformation is applied, σ11 = exp(η11 ),  σ22 = exp(η12 ).  (4.30)  For the proportions, we add the constraint p 3 < p1 , p3 < p2 ,  (4.31)  which requires the proportion of variant allele homozygotic points is less than that of the other two genotype clusters. We reparameterize the proportions p first with q for constraint (4.31),  p1 = q1 + q3 /3,      p2 = q2 + q3 /3, p3 = q3 /3, (4.32)   p = q ,  4 4   p5 = q 5 . 70  Chapter 4. Bayesian linear clustering and then apply a logit transformation to q,   q1 = exp(η13 )/(1 + (exp(η13 ) + . . . + exp(η16 ))), ...,  q4 = exp(η16 )/(1 + (exp(η13 ) + . . . + exp(η16 ))).  (4.33)  With the constraints (4.26), (4.28) and (4.31), we can fully identify the five components of the mixture. Label-switching is effectively prevented by these informative constraints. In addition, conforming to these constraints, noninformative priors cause no problem to guarantee a proper posterior distribution. From (4.27), (4.29), (4.30) and (4.33), we obtain a one-to-one map from θ = (a1 , a2 , a3 , b1 , b2 , b3 , σ12 , σ22 , σ32 , µ1 , µ2 , σ11 , σ22 , p1 , . . . , p5 )  to η = (η1 , η2 , η3 , b1 , b2 , b3 , η7 , η8 , µ1 , µ2 , η11 , η12 , η13 , η14 , η15 , η16 ), in which the constraints (4.26), (4.28) and (4.31) are satisfied. The absolute value of the Jacobian of the transformation is |J| =  ∂θ ∂η  =|a23 (a1 − a2 )(a2 − 2a1 )/a1 |σ12 σ22 σ11 σ22 (p1 − p3 )(p2 − p3 )p3 p4 p5 .  (4.34)  We set the priors as follows. ηk ∼ N (0, 1010 ), bk ∼ N (0, 10−6 ), σk2 ∼ IG(0.0001, 0.0001), for k = 1, 2, 3, and π(µ, σ11 , σ22 ) ∝  1 , π(q) ∝ 1. σ11 σ22  The priors for b1 , b2 , b3 is strong with the implication that all three lines should roughly pass through the origin; all priors for others parameters are vague except the constraints specified above. The five-component mixture model is applied to the plate shown in Figure 4.1. For an algorithm to work in this relatively high dimensional problem, it is essential that the initial state be close to the substantial support of the posterior distribution. We first run the “optim” function in the R language multiple times, which is initialized with standard normal N (0, I16 ). Then we initialize the algorithm with the solution with the largest posterior density, say η (0) . 71  740  Log unnormalized posterior density 745 750 755 760 765  Chapter 4. Bayesian linear clustering  0  2000  4000  6000 Iteration  8000  10000  Figure 4.10: Evolution of log unnormalized posterior density of sampled values for η for 11, 000 iterations using Algorithm 2 for the SNP genotyping data. We first run Algorithm 2.0, which is relatively slow, to get a rough estimate of the covariance matrix of η. Next we run the faster Algorithm 2. To alleviate the effect of autocorrelations, we keep only one simulated point for every ten iterations. Hereafter the iteration numbers are based on the “thinned” samples. Figure 4.10 displays the evolution of log unnormalized posterior density values of sampled values for η for 11, 000 iterations. Figure 4.11 displays the evolution of sampled values of θ for 11, 000 iterations. The first three panels are for slopes of the three lines −1/a1 , −1/a2 and −1/a3 . From these two figures we conclude empirically that the sampler has reached a stationary distribution. Figure 4.12 displays the autocorrelations of sampled values for θ for 11, 000 iterations. There are high autocorrelations in most parameters indicating that the sampler is not very efficient; more efficient transformations/samplers shall be investigated in future research. Further checking the density curves in Figure 4.13 of sampled values for θ for 11, 000 iterations, we observe that all the estimated density curves are roughly unimodal. The clustering results of this plate is displayed in Figure 4.14 where each point is classified into the cluster with the largest posterior membership probability. We note that the Bayesian method does identify several points into the variant allele homozygous genotype clusters (represented by “4” in Figure 4.14). Quite a few points are classified as background noise 72  Chapter 4. Bayesian linear clustering  8000  10000  − 1 a3 20 25  4000  6000 Iteration  8000  10000  b2  2000  4000  6000 Iteration  8000  10000  0  2000  4000  6000 Iteration  8000  10000  4000  6000 Iteration  1e−03  8000  10000  0  2000  4000  6000 Iteration  8000  10000  8000  10000  2000  4000  6000 Iteration  8000  10000  8000  10000  4000  6000 Iteration  8000  4000  10000  6000 Iteration  8000  10000  2000  4000  6000 Iteration  8000  10000  0  2000  4000  6000 Iteration  8000  10000  8000  10000  p5  p5 0.04  0.020  8000  0.02 0.00  0.000  6000 Iteration  0  10000  p4 0.010  0.04 p3 0.02  4000  10000  p2 0.20  2000  p4  0.00  2000  8000  0.15  0  p3  0  2000  0.06  6000 Iteration  6000 Iteration  p2  0.70  4000  0  p1  p1 0.75 0.80  σ22 0.0025 0.0005  2000  4000  σ11 6e−04  0 0.85  σ22  0  2000  0.25  6000 Iteration  0  2e−04  0.23  4000  10000  σ11  µ2 0.25 0.27  µ1 0.110 0.100  2000  8000  b3  µ2  0.120  µ1  0  6000 Iteration  σ23  2e−04  2000  4000  σ22  σ2 2 6e−04  σ2 1 0.00012  0  2000  2 σ3 0.00012  0  0  b3 0.00  0.02 −0.02  b2 0.00  b1 −0.020 −0.035  2000  b1  σ21  0.00006  0  0.02  6000 Iteration  −0.02  4000  0.00006  2000  10  1.80  0  −0.005  0.42  15  − 1 a2 1.90  − 1 a1 0.44  −1 a3  30  −1 a2  2.00  −1 a1  0  2000  4000  6000 Iteration  8000  10000  0  2000  4000  6000 Iteration  Figure 4.11: Evolution of sampled values for θ for 11, 000 iterations using Algorithm 2 for the SNP genotyping data. The first three panels are for slopes. 73  Chapter 4. Bayesian linear clustering − 1 a3  10  20 Lag  30  40  ACF 0.6 0.8 0.4  10  1.0  30  40  10  20 Lag  30  40  10  20 Lag  30  40  40  20 Lag  30  40  40  20 Lag  30  40  20 Lag  30  40  40  0  10  20 Lag  30  40  30  40  0.8 ACF 0.4 0.0  ACF 0.4  40  30  p5  0.0  30  20 Lag  0.8  10  0.8  0.8 ACF 0.4  20 Lag  10  p4  0.0  10  0  ACF 0.4  0  p3  0  40  0.0  ACF 0.4  40  30  p2  0.0  30  20 Lag  0.8  10  0.8  0.8 ACF 0.4  20 Lag  10  p1  0.0  10  0  ACF 0.4  0  σ22  0  40  0.0  ACF 0.4 0.0  30  30  σ11  0.8  0.8 ACF 0.4  20 Lag  20 Lag  ACF 0.4  10  µ2  0.0  10  10  0.0  0  µ1  0  40  0.8  0.8 0.0  30  30  σ23  ACF 0.4  0.8 ACF 0.4  20 Lag  0  σ22  0.0  10  20 Lag  ACF 0.85  0  σ21  0  10  0.75  0.80  0.4  0  0  b3  ACF 0.90  ACF 0.6 0.8  20 Lag  b2  1.00  b1  0  0.95  0  0.65  0.4  0.75  ACF 0.85  ACF 0.6 0.8  0.95  1.0  − 1 a2  1.0  − 1 a1  0  10  20 Lag  30  40  0  10  20 Lag  Figure 4.12: Autocorrelations of sampled values for θ for 11, 000 iterations using Algorithm 2 for the SNP genotyping data. 74  Chapter 4. Bayesian linear clustering − 1 a3  − 1 a2  0.00  0  0  20  Density 0.10  Density 5 10  Density 40 60  80  15  − 1 a1  0.42  0.43  0.44  0.45  0.46  1.80  1.85  1.90  2.00  10  15  20  25  30  b3  b2  −0.03  −0.02  −0.01  0.00  Density 30 10 0  0  0  10  20  Density 30  Density 40 60  80  50  50  b1  1.95  −0.02 −0.01  0.00  0.02  0.03  −0.02  0.00010  0.00014  0.120  0.125  Density 10000 25000 6000  0.24  0.25  0.27  0.28  0.29  2e−04  4e−04  0.004  Density 10 15  0.70  0.75  0.80  0.85  0.15  0.20  p4  p5  Density 20 40  250  0.00  0.01  0.02  0.03  0.04  0.05  0  0  0  50  Density 150  60  0.25  60  p3  Density 20 40  8e−04  0  0  0.003  6e−04  5  5  Density 10  Density 400 800  0.002  0e+00  p2  15  1200  0.26  p1  0  0.001  0.00014  0  0.23  σ22  0.000  0.00010  Density 2000 4000  Density 30 10 0  0.115  0.02  σ11  50  120 Density 40 80  0.110  0.00006  µ2  0  0.105  0.01  0  0.0002 0.0004 0.0006 0.0008 0.0010 0.0012  µ1  0.100  0.00  σ23  0  Density 10000 25000 0  0.00006  −0.01  σ22 Density 1000 3000  σ21  0.01  0.000  0.005  0.010  0.015  0.020  0.00  0.01  0.02  0.03  0.04  0.05  0.06  Figure 4.13: Density curves of sampled values for θ for 11, 000 iterations using Algorithm 2 for the SNP genotyping data. 75  Chapter 4. Bayesian linear clustering  2.0  3  5  3  0.5  Allele.Y 1.0  1.5  3  3  4 44 4 44 4 4 4 44 44  2 222 222 2 222 2 222222222 2 22 2222 22 2 2 222 2 222 222 222 222 222 22 222 5 2 2 22 2 2 5 2  4  0.2  5  1 111 111111111 111 1111111 1111111 1111111 1 11 11 11 1 1 1 1 1 1 11 1 1 1 1 1 1 1 1 11 1 1 1 1 1 111111 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 11 1 1 1 1111111111 11 111 11111111111111 1111  0.4  0.6 0.8 Allele.X  1.0  1.2  1.4  Figure 4.14: Clustering results of the SNP genotyping data in which points are classified into the cluster with the largest posterior membership probability. (represented by “5”); this is conservative because the orthogonal variation in the wild-type allele homozygous genotype cluster XX is very small and the orthogonal variation of cluster YY is linked to that of cluster XX. By depleting points in the YY cluster, we have observed that the Bayesian method also works if only one point is present in the YY cluster, which is due to the informative priors and constraints. Some other plates with sparse clusters are also analyzed using the Bayesian approach and the clustering results are satisfactory; the clustering results are omitted. For plates without sparse clusters, the effect of the above prior specification is minimal, we usually obtain similar clustering results to that of the partial likelihood approach in Chapter 3. Figure 4.15 shows the clustering results of the Bayesian approach for the plate analyzed in Chapter 3.  4.6  Discussion  We have proposed a Bayesian approach to linear clustering. It is flexible in the sense that only the deviation from a hyperplane is modelled parametrically; the position on the hyperplane is not modelled. The advantage of this approach is illustrated in the genotyping example, where the distribution along the line is complex and difficult to model. Furthermore, as was also illustrated in the SNP genotyping, we can incorporate elliptical clusters as 76  3 33 33 3 3 3 333 3 3  22 22  3  0.5  2 2 2 2 22 2 22 222 22 222222 2 22 222222222 2  5  3 3  1.0  Allele.Y  1.5  2.0  Chapter 4. Bayesian linear clustering  4 4 44 4 4 4 4 44 4 4 4 4 4 444 44 44 4 0.5  2  1 1  1  1 11 1 11 1111 11111 1 1111 11 1 1 1 1 1 1 1111111 1 1111111 111111111 1 1 111 1 11 1 1  1.0 Allele.X  1.5  Figure 4.15: Clustering results of the Bayesian approach for a plate without a sparse cluster. necessary and a background cluster, borrowing from standard model-based clustering. In our examples, label-switching is prevented either by a Gibbs sampler applied to a posterior distribution with isolated modes or by informative priors. In more general situations, we may need the ideas of tempering MCMC or Sequential Monte Carlo to explore the whole support of the posterior distribution and to deal with the label-switching problem. In this paper, the number of linear clusters are assumed known. In the situation of unknown number of clusters, our first thought is to investigate the feasibility of the Reversible Jump MCMC of Richardson and Green (1997). This may imply heavy computational burden. A related problem is the scalability of the Bayesian approach to large datasets and high dimensions. We leave these problems for further research.  77  Bibliography Anderson, T. (2003). An introduction to multivariate statistical analysis. Wiley Series in Probability and Statistics. Campbell, N. and R. Mahon (1974). A multivariate study of variation in two species of rock crab of genus Leptograpsus. Australian Journal of Zoology 22 (3), 417–425. Celeux, G., M. Hurn, and C. Robert (2000). Computational and Inferential Difficulties with Mixture Posterior Distributions. Journal of the American Statistical Association 95 (451). Fraley, C. and R. AE (2006). mclust Version 3 for R: Normal Mixture Modeling and Model-based Clustering. Technical report, Technical Report 504, University of Washington, Department of Statistics. Fujisawa, H., S. Eguchi, M. Ushijima, S. Miyata, Y. Miki, T. Muto, and M. Matsuura (2004). Genotyping of single nucleotide polymorphism using model-based clustering. Bioinformatics 20, 5. Fuller, W. (1987). Measurement error models. Wiley Series in Probability and Mathematical Statistics, New York: Wiley, 1987 . Gelman, A., J. Carlin, H. Stern, and D. Rubin (2004). Bayesian data analysis. Boca Raton, FL: Chapman and Hall/CRC. Harrington, J. (2007). lga: Tools for linear grouping analysis (LGA). R package version 1.0-0. Kang, H., Z. Qin, T. Niu, and J. Liu (2004). Incorporating Genotyping Uncertainty in Haplotype Inference for Single-Nucleotide Polymorphisms. The American Journal of Human Genetics 74 (3), 495–510. Olivier, M., L. Chuang, M. Chang, Y. Chen, D. Pei, K. Ranade, A. de Witte, J. Allen, N. Tran, D. Curb, et al. (2002). High-throughput genotyping of single nucleotide polymorphisms using new biplex invader technology. Nucleic Acids Research 30 (12), e53. 78  Bibliography Peel, D. and G. McLachlan (2000). Robust mixture modelling using the t distribution. Statistics and Computing 10 (4), 339–348. Ranade, K., M. Chang, C. Ting, D. Pei, C. Hsiao, M. Olivier, R. Pesich, J. Hebert, Y. Chen, V. Dzau, et al. (2001). High-Throughput Genotyping with Single Nucleotide Polymorphisms. Genome Research 11 (7), 1262– 1268. Richardson, S. and P. Green (1997). On Bayesian Analysis of Mixtures with an Unknown Number of Components. Journal of the Royal Statistical Society. Series B (Methodological) 59 (4), 731–792. Turner, T. (2006). mixreg: Functions to fit mixtures of regressions. R package version 0.0-2. Van Aelst, S., X. Wang, R. Zamar, and R. Zhu (2006). Linear grouping using orthogonal regression. Computational Statistics and Data Analysis 50 (5), 1287–1312. Yan, G., W. Welch, and R. Zamar (2008). A likelihood approach to linear clustering. Preprint.  79  Chapter 5  Consistency and asymptotic normality in a partial likelihood approach to linear clustering 5.1  Introduction  By linear clustering, we mean detecting linearly shaped clusters in a data set. We proposed in Yan et al. (2008) a parsimonious partial likelihood approach to linear clustering. Assume that the data x1 , . . . , xn are drawn by a random mechanism, represented by random vector X in a d-dimensional space which we do not model fully. Assume that The data lie around K hyperplanes {x : ak x = bk }, k = 1, . . . , K. Let Z = (Z1 , . . . , ZK ) be a random vector indicating these hyperplanes and Zk = 1 with probability pk for k = 1, . . . , K. Let p = (p1 , . . . , pK ) . We assume that, conditional on Zk = 1, ak X − bk ∼ N (0, σk2 ), k = 1, . . . , K. Let z1 , . . . , zn be the corresponding unobservable indicators for the data x1 , . . . , xn . Let κ be the collection of component parameters, κ = (a1 , b1 , 2 ) , and θ = (κ , p ) . The indicators z , . . . , z can be σ12 , . . ., aK , bK , σK 1 n regarded as realizations of random vectors Z1 , . . . , Zn , which in turn are an independent and identically distributed sample from Z. After integrating out the indicators z1 , . . . , zn , the partial likelihood function for parameters A version of this chapter will be submitted for publication. Authors: Guohua Yan, William J. Welch and Ruben H. Zamar.  80  Chapter 5. Consistency and asymptotic normality θ is  n  K  L(θ|x1 , . . . , xn ) = i=1 k=1  pk N (ak xi − bk ; 0, σk2 ).  (5.1)  As in the usual normal mixture model, the partial likelihood function (5.1) is unbounded: when a cluster consists only of points lying on a hyperplane, the contribution of each of these points to the partial likelihood tends to infinity. The infinity occurs on the boundary of the parameter space. We adopt the constraint of Hathaway (1985). Specifically, the following constraint is imposed on the standard deviations, min (σi /σj ) ≥ c > 0,  1≤i=j≤K  (5.2)  where c is a known constant determined a priori. In this article, constraint (5.2) is assumed whenever we refer to the partial likelihood function (5.1). The partial likelihood function (5.1) naturally brings the clustering problem into a finite mixture model framework. An EM algorithm can be used to ˆ is obtained, maximize (5.1); once an maximum partial likelihood estimate θ data point xi can be assigned to the component with the largest posterior probability. The probabilities are given by w ˆik =  pˆk N (ˆ ak xi − ˆbk ; 0, σ ˆk2 ) , i = 1, . . . , n; k = 1, . . . , K, K pˆk N (ˆ a xi − ˆbk ; 0, σ ˆ2) k=1  k  (5.3)  k  which also serve as a measure of uncertainty of classifying data point xi . We shall investigate the asymptotic properties of solutions to the partial likelihood function (5.1). As (5.1), regarded as a function of x, is not a density function for X, classical results of asymptotics on maximum likelihood estimators cannot be used directly. We borrow ideas in the formulation and in the proofs of results from Garc´ıa-Escudero et al. (2007), Wald (1949), Redner (1981) and Hathaway (1983, 1985). The rest of this article is organized as follows. Section 5.2 discusses the population counterpart of the partial likelihood function (5.1) and establishes the existence of its maximum. Section 5.3 proves the consistency of a maximum of the partial likelihood function (5.1) to the maximum of the population counterpart. The asymptotic normality of a solution of the partial likelihood function (5.1) is investigated in section 5.4.  81  Chapter 5. Consistency and asymptotic normality  5.2  Population version of the objective function  To motivate the population version of the objective function (5.1), we first review the connection between maximum likelihood estimation and the Kullback-Leibler divergence. See for example Kullback and Leibler (1951) and Eguchi and Copas (2006). Suppose that probability distributions P and Q are absolutely continuous with respect to the Lebesgue measure λ. Let p(x) and q(x) be the density functions (Radon-Nikodym derivatives) of P and Q respectively with respect to λ. Then the Kullback-Leibler divergence from P to Q is defined as DKL (P Q) =  log  p(x) dP (x). q(x)  Given a random sample x1 , x2 , . . . , xn from the underlying distribution P , let Pn be the empirical distribution. Now let Qθ be a statistical model having density f (x; θ) with respect to λ, where θ is a collection of unknown parameters. The Kullback-Leibler divergence from P to Qθ is DKL (P Qθ ) =  [log p(x) − log f (x; θ)]dP (x).  The empirical version of DKL (P ||Qθ ) is DKL (Pn  1 Qθ ) = n  n i=1  [log(1/n) − log f (xi ; θ)]  1 = log(1/n) − n  n  log f (xi ; θ). i=1  Apart from the factor 1/n, the second term is just the log likelihood function. Hence, maximizing the likelihood function is equivalent to minimizing the Kullback-Leibler divergence DKL (Pn Qθ ). If P = Qθ0 for some θ0 , then a maximum likelihood estimator is strongly consistent under some regular conditions, i.e., it converges almost surely to arg minθ DKL (P Qθ ). Back to the linear clustering setting, let K  f (x; θ) =  pk fk (x; θ),  (5.4)  k=1  where fk (x; θ) = √  (a x − bk )2 1 exp − k 2 2σk 2πσk  , 82  Chapter 5. Consistency and asymptotic normality and x is a generic point in d . We still use P to denote the underlying distribution of the data and use p(x) to denote the density of P with respect to the Lebesgue measure λd . Although f (x; θ) cannot be a density as a function of x, we can nevertheless define DKL (P f (·; θ)) =  [log p(x) − log f (x; θ)]dP (x),  or equivalently define g(θ) =  log f (x; θ)dP (x),  (5.5)  and show that the maximum partial likelihood solution of (5.1) converges to the set {θ 0 : g(θ 0 ) = max g(θ)}. θ  With constraint (5.2), The parameter space is   2 2   θ = (a1 , b1 , σ1 , . . . , aK , bK , σK , p1 , . . . , pK ) : 0 < p < 1, a a = 1, −∞ < b < ∞, σ > 0, k = 1, . . . , K. . Θc = k k k k k   K p = 1, min σ /σ ≥ c > 0. i,j i j k k=1  We show in this section that the supremum of g is finite and attainable. First we prove the following lemma. Lemma 5.1 Let P be an absolutely continuous probability measure with finite second moments. Let 0 < c < 1, if c ≤ σk ≤ 1/c, for all k = 1, . . . , K. Let (a x − bk )2 s(κ) ≡ min k 2 dP (x). k σk Then s(κ) attains its infimum s0 > 0 under constraint (5.2) at some κ0 . Proof It is obvious that 0 < s(κ) < ∞. First we show that it is possible to bound bk , k = 1, . . . , K and the infimum is not missed. Suppose that there is a positive number r < |bk |, for k = 1, . . . , K. Let r0 > 0 such that Prob(|X| ≤ r0 ) > 0. Then when r > r0 , s(θ) ≥ c2 (r − r0 )2 Prob(|X| ≤ r0 ) → ∞,  as r → ∞.  83  Chapter 5. Consistency and asymptotic normality Now without loss of generality, we assume that |bk | ≤ v for some v and for k = 1, . . . , K − 1. Let r = c2 (|bK | − v). Then when |bK | is getting large, s(κ)  min  ≥  |x|≤r 1≤k≤K  (ak x − bk )2 dP (x) σk2  (ak x − bk )2 dP (x) σk2 |x|≤r 1≤k≤K−1  =  min  (ak x − bk )2 dP (x), 1≤k≤K−1 σk2  →  min  ≥  min  1≤k≤K  as |bK | → ∞  (ak x − bk )2 dP (x), σk2  where the |bK | in the last term is arbitrary number bounded from below by v. Since (a x − bk )2 2(|x|2 + v 2 ) ≤ , min k 2 k c2 σk by Lebesgue’s dominated convergence theorem, s is continuous. Now it is constrained on a compact set, the infimum of s is attainable for some κ0 as s0 = s(κ0 ) > 0. Theorem 5.1 Let P be an absolutely continuous probability measure with finite second moments. Then the supremum of g, which is defined in (5.5), over Θc is finite and attainable. Proof Let σk2 = τk2 σ 2 for k = 1, . . . , K. For every x and θ, we have log f (x, θ) ≤ log max fk (x, θ) = max log fk (x, θ). k  k  Therefore, g(θ) ≤ =  max log fk (x, θ)dP (x) k  − min k  (a x − bk )2 1 log(2πτk2 σ 2 ) + k 2 2 2 2τk σ  1 ≤ − log(2πc2 ) − log σ − 2 1 ≤ − log(2πc2 ) − log σ − 2  1 2σ 2 s0 , 2σ 2  min k  dP (x)  (ak x − bk )2 dP (x) τk2  84  Chapter 5. Consistency and asymptotic normality where Lemma 5.1 is used. Clearly, the last term attains its maximum at σ 2 = s0 and it goes to −∞ as σ → 0 or ∞. Take an arbitrary θ 0 , there exist positive numbers s1 and s2 such that g(θ) < g(θ 0 ) uniformly for θ such that σ < s1 /c or σ > cs2 . Now if we bound s1 /c ≤ σ ≤ cs2 , we have log f (x, θ) ≥ log min fk (x, θ) k  = min log fk (x, θ) (a x − bk )2 1 = min(− log(2πτk2 σ 2 ) − k 2 2 ) 2 2τk σ 2(|x|2 + maxk |bk |2 ) 1 ≥ − log(2πs22 ) − , 2 2s21 and  1 log f (x, θ) ≤ log max fk (x, θ) ≤ − log(2πs21 ). k 2 By Lebesgue’s dominated convergence theorem, g is continuous. If bk → ±∞ for some k, then by Fatou’s lemma, we have lim sup g(θ) ≤ bk →±∞  lim sup log f (x, θ)dP (x) =  log  bk →±∞  pk fk (x, θ)dP (x). k =k  Therefore, there exist real numbers t1 and t2 such that we would not miss the supremum of g, if we bound θ in the compact set {θ : s1 ≤ σk2 ≤ s2 , t1 ≤ bk ≤ t2 , for all k = 1, . . . , K}. Since g is continuous, its supremum is attainable. Remark In the proof it is not ruled out that some pk may be 0 for g to attain its maximum.  5.3  Consistency  First we prove that there exists a global maximizer of the partial likelihood function (5.1) over Θc . Theorem 5.2 Let {x1 , . . . , xn } be a set of points in d that cannot fit in K hyperplanes. That is, there do not exist {(ak , bk ) : k = 1, . . . , K} such that, for every xi , ak xi = bk for some k (which may depend on i). Let n  l(θ) =  log f (xi , θ). i=1  85  Chapter 5. Consistency and asymptotic normality Then a constrained global maximizer of l(θ) over Θc exists. Proof Let B(0, r) = {x : |x| ≤ r} be a ball that contains all the points x1 , . . . , xn . For any θ, if there is some bk such that |bk | > r, then (ak xi − bk )2 ≥ (ak xi − rsgn(bk ))2 . So l(θ) is not decreased if bk is replaced with rsgn(bk ) whenever |bk | ≥ r. So we can bound θ such that |bk | ≤ r, for all k = 1, . . . , K. Let s0 =  min  {(a1 ,b1 ,...,aK ,bK ):ak ak =1,|bk |≤r}  max min(ak xi − bk )2 . i  k  Since the points x1 , . . . , xn cannot fit in K hyperplanes, we have s0 > 0. There exists a point xj such that (ak xj − bk )2 ≥ s0 for all k = 1, . . . , K. Note also that 1/σk ≤ 1/(cσ1 ) and that 1/σk ≥ c/σ1 , we have l(θ) ≤  log( i=j  k  pk √  1 ) + log( 2πσk  k  pk √  1 ) + log( 2πcσ 1 k √ cs0 = −n log( 2πcσ1 ) − 2 . 2σ1  ≤ (n − 1) log(  pk √  k  1 s0 exp(− 2 )) 2σk 2πσk pk √  1 cs0 exp(− 2 )) 2σ 2πcσ1 1  The last term tends to −∞ as σ12 tends to zero. If some σk2 tends to zero, so do all other σk2 s. Thus, there exists a positive number s1 such that l(θ) ≤ l(θ 0 ) whenever a σk2 < s1 . On the other hand, √ 1 1 l(θ) ≤ log( pk √ ) ≤ n log(max √ ) = −n min log( 2πσk ). k k 2πσk 2πσk i k If some σk2 tends to ∞, so do all other σk2 ’s and then l(θ) decreases to −∞. Similarly, there is a positive number s2 such that l(θ) ≤ l(θ 0 ) whenever a σk2 > s2 . Therefore, we can bound θ in a compact set {θ : |bk | ≤ r, s1 ≤ σk2 ≤ s2 , for all k = 1, . . . , K}. As l is continuous, a global maximizer over Θc exists. Corollary 5.1 Let P be an absolutely continuous probability measure with finite second moments. Let X1 , . . . , Xn be a sample from P . Then a constrained global maximizer of l(θ|X1 , . . . , Xn ) exists almost surely if n ≥ Kd + 1. 86  Chapter 5. Consistency and asymptotic normality Proof As P is absolute continuous, X1 , . . . , Xn are distinct with probability 1. We need Kd points to determine K hyperplanes. Each of the remaining points lies on one of these hyperplanes with probability 0, since the set of points in K hyperplanes has zero probability. To prove the consistency result, we need some preliminary results. Lemma 5.2 If the probability measure P has finite second moments, then | log f (x, θ)|dP (x) < ∞. Proof | log f (x, θ)| dP (x) [log f (x, θ)]+ dP (x) +  = ≤  [log max fk (x, θ)]+ dP (x) k  {x:f (x,θ)≥1}  [log min fk (x, θ)]− dP (x)  +  k  {x:f (x,θ)<1}  ≤  [log f (x, θ)]− dP (x)  [log fk (x, θ)]+ dP (x) {x:f (x,θ)≥1} k  [log fk (x, θ)]− dP (x)  + {x:f (x,θ)<1} k  = k  ≤  | log fk (x, θ)|dP (x) √ [| log( 2πσk )| +  k  (ak x − bk )2 dP (x)] 2σk2  <∞.  Lemma 5.3 Let w(x, θ, ρ) =  sup  f (x, θ ).  {θ :|θ −θ|≤ρ}  Then [log w(x, θ, ρ)]+ dP (x) < ∞. 87  Chapter 5. Consistency and asymptotic normality Proof In fact, w(x, θ, ρ) ≤ ≤  fk (x, θ)  sup {θ :|θ −θ|≤ρ} k  sup {θ :|θ −θ|≤ρ} k  √  1 2πσk  ≡ M (θ, ρ). Therefore, [log w(x, θ, ρ)]+ dP (x) ≤  [log M (θ, ρ)]+ dP (x) < ∞.  Lemma 5.4 lim  ρ→0  log w(x, θ, ρ)dP (x) =  log f (x, θ)dP (x).  Proof Since log w(x, θ, ·) is an increasing function of ρ, as ρ → 0, [log w(x, θ, ρ)]− increases. By the monotone convergence theorem, lim  ρ→0  [log w(x, θ, ρ)]− dP (x) =  [log f (x, θ)]− dP (x).  When ρ is sufficiently small, [log w(x, θ, ρ)]+ is dominated by [log w(x, θ, ρ0 )]+ for some ρ0 . And the latter is integrable by Lemma 5.3. By Lebesgue’s dominated convergence theorem, lim  ρ→0  [log w(x, θ, ρ)]+ dP (x) =  [log f (x, θ)]+ dP (x).  Notice that [log w(x, θ, ρ)]+ dP (x) and [log f (x, θ)]+ dP (x) are finite, the lemma is proved. Note that points in Θc are not identifiable for f (x; ·). The function f (x; ·) remains the same if we permutate the labels 1, . . ., K; pk1 and pk2 are not identifiable if (ak1 , bk1 , σk21 ) = (ak2 , bk2 , σk22 ). Thus the consistency result is in a quotient topology space. Let ∼ be an equivalent relation on Θc such that θ 1 ∼ θ 2 if and only if f (x; θ 1 ) = f (x; θ 2 ) almost surely in P . Denote by Θqc the quotient topological space consisting of all equivalent classes of ∼. For a point θ 0 that maximizes log f (x; θ)dP (x), its equivalent class in Θqc is denoted by θ q0 . 88  Chapter 5. Consistency and asymptotic normality Theorem 5.3 Let X1 , . . . , Xn be a sample from P which is absolutely continuous with finite second moments. Let B be a compact subset of Θc that ˆ (n) be a global maximizer of l(θ|X1 , . . . , Xn ) on B. Then contains C0 . Let θ ˆ (n) → θ q almost surely in the topological space B q . θ 0 Proof Let ω be a closed subset of B which does not intersect with C0 . For each point θ in ω, we associate a positive value ρθ such that E log w(X, θ, ρθ ) < E log f (X, θ 0 ), where θ 0 is a point in C0 . The existence of such a ρθ follows from Lemma 5.4 and the definition of C0 (C0 exists because P has finite second moments, ˆ (n) exists because P is absolutely continuous, by Theoby Theorem 5.1; θ rem 5.2). Since B is compact, its closed subset ω is also compact. There exists a finite number of points θ 1 , . . . , θ h in ω such that h  ω⊂  j=1  {θ : |θ − θ j | ≤ ρθj }.  Then h  0 ≤ sup f (x1 , θ) . . . f (xn , θ) ≤ θ∈ω  w(x1 , θ j , ρθj ) . . . w(xn , θ j , ρθj ). j=1  By the strong law of large numbers, n  Prob  lim  n→∞  i=1  [log w(Xi , θ j , ρθj ) − log f (Xi , θ 0 )] = −∞  = 1,  for j = 1, . . . , h. That is, Prob  w(X1 , θ j , ρθj ) . . . w(Xn , θ j , ρθj ) =0 n→∞ f (X1 , θ 0 ) . . . f (Xn , θ 0 ) lim  = 1,  for j = 1, . . . , h. Therefore, we have Prob  supθ∈ω f (X1 , θ) . . . f (Xn , θ) =0 n→∞ f (X1 , θ 0 ) . . . f (Xn , θ 0 ) lim  = 1.  (5.6)  Denote |θ − C0 | ≡ minθ0 ∈C0 |θ − θ 0 |. The minimum is attainable since C is a closed set in B and hence compact. We need only to prove that 89  Chapter 5. Consistency and asymptotic normality ˆ (n) are in C0 . If not, there exists a all limit points θ ∗ of the sequence θ limit point θ ∗ and an > 0 such that |θ ∗ − C0 | ≥ . This implies that ˆ (n) that lie in ω ≡ {θ : |θ − C0 | ≥ }. Thus there are infinitely many θ  ˆ f (x1 , θ  (n)  (n)  ˆ ) ≤ sup ) . . . f (xn , θ θ∈ω f (x1 , θ) . . . f (xn , θ) for infinitely many (n) (n) ˆ ) . . . f (xn , θ ˆ ) ≥ f (x1 , θ 0 ) . . . f (xn , θ 0 ). We have n. Since f (x1 , θ f (x1 , θ 0 ) . . . f (xn , θ 0 ) ≤ sup f (x1 , θ) . . . f (xn , θ), θ∈ω  for infinitely many n. Since ω is a closed set in B, this is an event with probability zero according to equation 5.6. This completes the proof. In next step, we shall show that limit points of the constrained global ˆ (n) over Θc are almost surely in a compact space and hence maximizer θ consistency follows from Theorem 5.3. Lemma 5.5 Let X1 , . . . , Xn be a sample of P which is absolutely continuˆ (n) is a constrained global maximizer ous with finite second moments. Let θ of l(θ|X1 , . . . , Xn ) over Θc . Then there exist positive numbers s1 and s2 such that (n) Prob{lim inf min(ˆ σk )2 ≥ s1 } = 1, (5.7) n→∞  k  and  (n)  Prob{lim sup max(ˆ σk )2 ≤ s2 } = 1. n→∞  (5.8)  k  (n)  Proof First we prove Equation (5.7). Let σ ˆk (n) Then τˆk ∈ [c, 1/c]. We write  (n) (n)  = τˆk σ ˆ1  for k = 1, . . . , K.  (n) (n) (n) (n) (n) (n) (n) (n) ˆ1 , . . . , a ˆ K , ˆb1 , . . . , ˆbK , τˆ1 σ, . . . , τˆK σ). h(σ) = l(ˆ p1 , . . . , pˆK , a  Then it satisfies  dh(σ) | (n) = 0. dσ σˆ1  This gives rise to n i=1  −  1 (n) σ ˆ1  ˆ f (Xi , θ  (n)  )+  1 (n) (ˆ σ 1 )3  (n) (n) ak ) Xi −ˆbk )2 (n) K ((ˆ ˆ (n) ) pˆk fk (Xi , θ (n) 2 k=1 (ˆ τk ) (n) ˆ ) f (Xi , θ  = 0.  90  Chapter 5. Consistency and asymptotic normality (n)  Solving for σ ˆ1  yields  (n)  (ˆ σ1 )2 =  (n) (n) ak ) Xi −ˆbk )2 (n) K ((ˆ ˆ (n) ) pˆk fk (Xi , θ (n) 2 k=1 (ˆ τk ) . ˆ (n) ) f (Xi , θ  n  1 n  i=1  Then, we have (n) (ˆ σ1 )2  ≥c  n  21  n  (n) (n) min ((ˆ ak ) Xi − ˆbk )2 .  (5.9)  1≤k≤K  i=1  We now show that for n0 ≥ Kd + 1,  n0  s(n0 ) = E{  inf  {(a1 ,...,aK ,b1 ,...,bK ):ak ak =1}  i=1  min (ak Xi − bk )2 } > 0.  1≤k≤K  (5.10)  In fact, n0  inf  {(a1 ,...,aK ,b1 ,...,bK ):ak ak =1}  i=1  min (ak xi − bk )2 dP n0 (x)  1≤k≤K  n0  =  inf  {(a1 ,...,aK ,b1 ,...,bK ):ak ak =1,|bk |≤r(x1 ,...,xn0 )}  =  i=1 n0  min  {(a1 ,...,aK ,b1 ,...,bK ):ak ak =1,|bk |≤r(x1 ,...,xn0 )}  i=1  min (ak xi − bk )2 dP n0 (x)  1≤k≤K  min (ak xi − bk )2 dP n0 (x)  1≤k≤K  s0 (x1 , . . . , xn0 )dP n0 (x) > 0,  ≡  since s0 (X1 , . . . , Xn0 ) > 0 almost surely from Corollary5.1 and the argument in Theorem 5.2. By the strong law of large numbers, we have 1 m  n0  m  inf  {(a1 ,...,aK ,b1 ,...,bK ):ak ak =1}  j=1  i=1  min (ak X(j−1)n0 +i − bk )2  1≤k≤K  → s(n0 ),  with probability 1. Let s1 = c2 s(n0 , P )/n0 . Then by equation (5.9), (n)  lim inf (ˆ σ1 )2 n→∞  1 ≥c lim inf { m→∞ mn0  mn0  2  ≥  c2  lim  n0 m→∞  1 m  inf  {(a1 ,...,aK ,b1 ,...,bK ):ak ak =1}  m j=1  {  i=1 n0  inf  {(a1 ,...,aK ,b1 ,...,bK ):ak ak =1}  min (ak Xi − bk )2 }  1≤k≤K  i=1  min (ak X(j−1)n0 +i − bk )2 }  1≤k≤K  =s1 , 91  Chapter 5. Consistency and asymptotic normality (n)  with probability 1. The same s1 serves as a lower bound of lim inf n (ˆ σk )2 for k = 2, . . . , K as well. Noticing that s1 does not depend on a set of observations x1 , . . . , xn , we have proved Equation (5.7). Now we prove Equation (5.8). Let Θsc = {θ : θ ∈ Θc , σk2 > s for some k}. Define q(x, s) = sup f (x, θ). θ∈Θsc 1 Then q(x, s) ≤ √2πsc . By Lemma 5.2, E| log f (X, θ 0 )| < ∞. There exists a positive number s2 , such that  E(log(q(X, s2 ))) < E(log f (X, θ 0 ). By the strong law of large numbers, n  Prob{ lim  n→∞  i=1  [log q(Xi , s2 ) − log f (Xi , θ 0 )] = −∞} = 1.  This implies that supθ∈Θsc2 f (X1 , θ) . . . f (Xn , θ) = 0} = 1. n→∞ f (X1 , θ 0 ) . . . f (Xn , θ 0 )  Prob{ lim  ˆ (n) always satisfies Equation 5.8 follows immediately, since θ ˆ (n) ) . . . f (Xn , θ ˆ (n) ) f (X1 , θ ≥ 1. f (X1 , θ 0 ) . . . f (Xn , θ 0 )  Now we prove the main consistency result. Theorem 5.4 Let X1 , . . . , Xn be a sample from an absolutely continuous ˆ (n) be a global maxprobability measure P with finite second moments. Let θ ˆ (n) → θ q almost surely in the imizer of l(θ |X1 , . . . , Xn ) over Θc . Then θ 0 topological space Θqc . Proof From Lemma 5.5, we need only to consider the subspace of Θc , Θsc = {θ ∈ Θc : s1 ≤ σk2 ≤ s2 , for all k = 1, . . . , K}, 92  Chapter 5. Consistency and asymptotic normality where s1 , s2 are positive numbers determined in Lemma 5.5. Since Θsc is not compact, Theorem 5.3 cannot be used directly. We shall use the compactification device in Hathaway (1985) and also in Kiefer and Wolfowitz (1956). In the space Θsc , define the metric δ(θ, θ ) = i  | arctan θ i − arctan θ i |,  where | · | is the Euclidean distance and θ i and θ i are components of θ and ¯ sc be the set of Θsc along with all its limit points. Then θ respectively. Let Θ θ : 0 ≤ pk ≤ 1, K k=1 pk = 1, mini,j σi /σj ≥ c > 0, ak ak = 1, −∞ ≤ bk ≤ ∞, s1 ≤ σk ≤ s2 , k = 1, . . . , K.  ¯ sc = Θ  ¯ sc as is compact. Since f (x, ·) is continuous on Θsc , it can be extended to Θ f (x, θ) = k  pk I(∞ < bk < ∞)fk (x, θ).  We have shown g(θ) = E(log f (x, θ)) is continuous on Θsc . It is continuous ¯ sc as well. To see this, let θ (n) be a sequence tending to θ. If all bk = ±∞, on Θ or all pk = 0 whenever bk = ±∞, then f (x, θ) = 0 and g(θ) = −∞. In this case, g(θ (n) )  = ≤ → =  log f (x, θ (n) )dP (x) max log fk (x, θ (n) )dP (x)  {k:pk >0}  − min −∞.  {k:pk >0}  (a x − bk )2 1 dP (x) log(2πσk2 ) + k 2 2 2σk  If there is some k such that pk = 0 and bk = ±∞. Then 1 log f (x, θ (n) ) ≤ log max fk (x, θ (n) ) ≤ − log(2πs1 ), 2 and log f (x, θ (n) ) ≥ log ≥ log  (n)  min  pk fk (x, θ (n) )  min  pk fk (x, θ) − ,  {k:pk >0,bk =±∞} {k:pk >0,bk =±∞}  93  Chapter 5. Consistency and asymptotic normality for some > 0 when n is sufficiently large. By Lemma 5.2, the latter term is integrable. By Lebesgue’s dominated convergence theorem, g is continuous at θ. Lemmas 5.3 and 5.4 are easily seen to hold and hence we can repeat literally the proof of Theorem 5.3. This completes the proof.  5.4  Asymptotical normality  ˆ (n) converges in the above quotient topological space. The global maximizer θ ˆ (n) , which converges to a There is hence a subsequence, still denoted by θ point θ 0 in C0 in the original space. If θ 0 is an interior point of Θc , then we ˆ (n) ) about θ 0 , can expand l (θ ˆ (n) ) l (θ ˆ (n) − θ 0 ) =l (θ 0 ) + l (θ 0 )(θ 1 ˆ (n) ˆ (n) − θ 0 )], . . . , (θ ˆ (n) − θ 0 )T l (θ ∗ )(θ ˆ (n) − θ 0 )]T , + [(θ − θ 0 )T l1 (θ ∗ )(θ s 2 (5.11) where s is the dimension of θ and θ ∗ is a point on the line segment connecting ˆ (n) and θ 0 ; l is an s-dimensional vector of first derivatives of l, l and l θ j  are s × s matrices of second and third derivatives of l. For ease of notation, we do not differentiate which components of p and ak are free parameters and which are not; it is assumed that the constraints pk = 1 and ak ak = 1 are taken care in the calculation of these derivatives. ˆ (n) satisfies the first order The left-hand side of (5.11) is 0, because θ conditions. Let v0 (θ 0 ) = E[((log f (X, θ 0 )) ]. Then v0 (θ 0 ) = g (θ 0 ) = 0 by Lebesgue’s dominated convergence theorem. Let v1 (θ 0 ) = E[((log f (X, θ 0 )) )2 ]. It is straightforward to verify that all the entries of v1 (θ 0 ) are finite if the underlying distribution P has finite fourth moments. Then by the central √ limit theorem, the first derivative l (θ 0 )/ n is asymptotical normal, 1 L √ l (θ 0 ) → N (0, v1 (θ 0 )). n 94  Chapter 5. Consistency and asymptotic normality Let v2 (θ 0 ) = E[(log f (X, θ 0 )) ], and v3 (θ 0 ) = E[(log f (X, θ 0 )) ]. Again, it is straightforward but tedious to verify that all the entries of v2 (θ 0 ) are finite if P has finite fourth moments and that all the entries of v3 (θ 0 ) are finite if P has finite sixth moments. By the strong law of large numbers, the second derivative l (θ 0 )/n tends to a constant matrix, 1 l (θ 0 ) → v2 (θ 0 ), n and the third derivative l is bounded entry-wise. So we have the following Theorem 5.5 Let X1 , . . . , Xn be a sample from an absolutely continuous ˆ (n) be a subsequence probability measure P with finite sixth moments. Let θ of global maximizer of l(θ|X1 , . . . , Xn ) over Θc , which tends to an interior point θ 0 . Then √ L ˆ (n) − θ 0 ) → n(θ N (0, v(θ 0 )),  where v(θ 0 ) = [v2 (θ 0 )]+ v1 (θ 0 )[v2 (θ 0 )]+ and A+ is the Moore-Penrose inverse of matrix A. ˆ Denote w ˆ In equation (5.3), w ˆik is a function of xi and θ. ˆik = hk (xi , θ), k = 1, . . . , K. By the Delta method, we have the following Corollary 5.2 Let X1 , . . . , Xn be a sample from an absolutely continuous ˆ (n) be a subsequence probability measure P with finite sixth moments. Let θ of global maximizer of L(θ|X1 , . . . , Xn ) over Θc , which tends to an interior point θ 0 . Let x be a data point. Then √  ˆ n(hk (x, θ  (n)  L  (0)  (0)  ) − hk (x, θ 0 )) → N (0, [hk (x, θ 0 )] v(θ 0 )[hk (x, θ 0 )]), (0)  for k = 1, . . . , K, where hk (x, θ 0 ) =  ∂hk (x,θ) |θ=θ0 . ∂θ  Using this corollary, we can build approximate confidence intervals for ˆ and hence evaluate the clustering of a data set. wik by replacing θ 0 with θ  95  Bibliography Eguchi, S. and J. Copas (2006). Interpreting Kullback–Leibler divergence with the Neyman–Pearson lemma. Journal of Multivariate Analysis 97 (9), 2034–2040. Garc´ıa-Escudero, L., A. Gordaliza, R. San Mart´ın, S. van Aelst, and R. Zamar (2007). Robust linear clustering. Preprint. Hathaway, R. (1983). Constrained maximum-likelihood estimation for a mixture of m univariate normal distributions. Ph. D. thesis, Rice University. Hathaway, R. (1985). A Constrained Formulation of Maximum-Likelihood Estimation for Normal Mixture Distributions. The Annals of Statistics 13 (2), 795–800. Kullback, S. and R. Leibler (1951). On Information and Sufficiency. The Annals of Mathematical Statistics 22 (1), 79–86. Redner, R. (1981). Note on the Consistency of the Maximum Likelihood Estimate for Nonidentifiable Distributions. The Annals of Statistics 9 (1), 225–228. Wald, A. (1949). Note on the Consistency of the Maximum Likelihood Estimate. The Annals of Mathematical Statistics 20 (4), 595–601. Yan, G., W. Welch, and R. Zamar (2008). A likelihood approach to linear clustering. Preprint.  96  Chapter 6  Future Work In this thesis work, we have developed a SNP genotype calling algorithm based on the linear grouping algorithm (Van Aelst et al., 2006), proposed a flexible model-based approach to linear clustering and introduced a Bayesian approach to linear clustering with specific relevance to the SNP genotyping problem. In this chapter, we briefly describe a few possible directions for future research.  6.1  Robustness consideration  Robustness to outliers is desirable in linear clustering, as the assumption of normal deviations around hyperplanes is sensitive to large deviations in the orthogonal direction. In addition to the inclusion of a uniform background cluster (Banfield and Raftery, 1993), one option would be to use a heavier tailed distribution, for example, Student’s t distribution with small degrees of freedom or with degrees of freedom depending on the data. In the partial likelihood approach, this would adapt Peel and McLachlan (2000)’s EM algorithm for t mixture models from the elliptical context to the linear clustering setting. The adaptation is straightforward but computationally more expensive. Further ideas include estimating the component covariance matrices in the M-step in a robust way, for example, trimming off some points as done by Garc´ıa-Escudero et al. (2007). In the Bayesian approach, we have already used a Student’s t distribution with small degrees of freedom.  6.2  Asymptotics  For the partial likelihood approach, we produced some asymptotic results in the case of normal orthogonal deviations. For more general cases, such as Student’s t distribution, these results have not been established. We shall study the asymptotic evaluation in more general cases.  97  Chapter 6. Future Work  6.3  Model Extension  With a x = b, we are specifying a hyperplane in d − 1 dimensions. With little effort, this could be generalized to a mixture of partial likelihoods, each of which specifies a hyperplane of dimension q < d, n  K  l(κ, p|x1:n ) = i=1 k=1  pk N (Ak xi − bk ; 0, Σk ),  (6.1)  where A is of dimension d × (d − q), b is a vector of dimension d − q, and Σk is a (d − q) × (d − q) covariance matrix for the deviation from the hyperplane. In the extreme case of a 0-dimension hyperplane, which is a point, we have the usual mixture of multivariate normal distributions. A mixture of components with various dimensions could be considered.  6.4  Variable/model selection  For a large, high dimensional dataset, efficient computation is essential in order to uncover linear patterns. In addition to develop more efficient algorithms, one idea is to use a subset of the data and/or a subset of variables. To this end, we are interested in methods to screen variables as well as to determine the number of linear structures.  6.5  Bayesian computation  In our examples in the Bayesian approach, label-switching is prevented either by a Gibbs sampler applied to a posterior distribution with isolated modes or by informative priors. In more general situations, we may need the ideas of tempering MCMC or Sequential Monte Carlo to explore the whole support of the posterior distribution and deal with the label-switching problem. In addition, the number of linear clusters are assumed known in our approach. In the situation of unknown number of clusters, our first thought is to investigate the feasibility of the Reversible Jump MCMC of Richardson and Green (1997). This may imply heavy computational burden. A related problem is again the scalability of the Bayesian approach to large datasets and high dimensions. We leave these problems for further research.  98  Bibliography Banfield, J. and A. Raftery (1993). Model-Based Gaussian and NonGaussian Clustering. Biometrics 49 (3), 803–821. Garc´ıa-Escudero, L., A. Gordaliza, R. San Mart´ın, S. van Aelst, and R. Zamar (2007). Robust linear clustering. Preprint. Peel, D. and G. McLachlan (2000). Robust mixture modelling using the t distribution. Statistics and Computing 10 (4), 339–348. Richardson, S. and P. Green (1997). On Bayesian Analysis of Mixtures with an Unknown Number of Components. Journal of the Royal Statistical Society. Series B (Methodological) 59 (4), 731–792. Van Aelst, S., X. Wang, R. Zamar, and R. Zhu (2006). Linear grouping using orthogonal regression. Computational Statistics and Data Analysis 50 (5), 1287–1312.  99  Appendix A  Glossary of some genetic terms (Extracted from the website http://www.genome.gov/glossary.cfm.) adenine (A) One of the four bases in DNA that make up the letters ATGC, adenine is the ”A”. The others are guanine, cytosine, and thymine. Adenine always pairs with thymine. allele One of the variant forms of a gene at a particular locus, or location, on a chromosome. Different alleles produce variation in inherited characteristics such as hair color or blood type. In an individual, one form of the allele (the dominant one) may be expressed more than another form (the recessive one). base pair Two bases which form a ”rung of the DNA ladder.” A DNA nucleotide is made of a molecule of sugar, a molecule of phosphoric acid, and a molecule called a base. The bases are the ”letters” that spell out the genetic code. In DNA, the code letters are A, T, G, and C, which stand for the chemicals adenine, thymine, guanine, and cytosine, respectively. In base pairing, adenine always pairs with thymine, and guanine always pairs with cytosine. chromosome One of the thread-like ”packages” of genes and other DNA in the nucleus of a cell. Different kinds of organisms have different numbers of chromosome. Humans have 23 pairs of chromosomes, 46 in all: 44 autosomes and two sex chromosomes. Each parent contributes one chromosome to each pair, so children get half of their chromosomes from their mothers and half from their fathers.  100  Appendix A. Glossary of some genetic terms cytosine One of the four bases in DNA that make up the letters ATGC, cytosine is the ”C”. The others are adenine, guanine, and thymine. Cytosine always pairs with guanine. deoxyribonucleic acid (DNA) The chemical inside the nucleus of a cell that carries the genetic instructions for making living organisms. diploid The number of chromosomes in most cells except the gametes. In humans, the diploid number is 46. dominant A gene that almost always results in a specific physical characteristic, for example, a disease, even though the patient’s genome possesses only one copy. With a dominant gene, the chance of passing on the gene (and therefore the disease) to children is 50-50 in each pregnancy. gene The functional and physical unit of heredity passed from parent to offspring. Genes are pieces of DNA, and most genes contain the information for making a specific protein. gene amplification An increase in the number of copies of any particular piece of DNA. A tumor cell amplifies, or copies, DNA segments naturally as a result of cell signals and sometimes environmental events. gene expression The process by which proteins are made from the instructions encoded in DNA. genetic code (ATCG) The instructions in a gene that tell the cell how to make a specific protein. A, T, G, and C are the ”letters” of the DNA code; they stand for the chemicals adenine, thymine, guanine, and cytosine, respectively, that make up the nucleotide bases of DNA. Each gene’s code combines the four chemicals in various ways to spell out 3-letter ”words” that specify which amino acid is needed at every step in making a protein. genetic marker A segment of DNA with an identifiable physical location on a chromosome and whose inheritance can be followed. A marker can be a gene, or it can be some section of DNA with no known function. Because DNA segments that lie near each other on a chromosome tend to be inherited together, markers are often used as indirect ways of tracking the inheritance 101  Appendix A. Glossary of some genetic terms pattern of a gene that has not yet been identified, but whose approximate location is known. genome All the DNA contained in an organism or a cell, which includes both the chromosomes within the nucleus and the DNA in mitochondria. genotype The genetic identity of an individual that does not show as outward characteristics. guanine One of the four bases in DNA that make up the letters ATGC, guanine is the ”G”. The others are adenine, cytosine, and thymine. Guanine always pairs with cytosine. haploid The number of chromosomes in a sperm or egg cell, half the diploid number. Haplotype It refers to a set of SNPs found to be statistically associated on a single chromatid. With this knowledge, the identification of a few alleles of a haplotype block unambiguously identifies all other polymorphic sites in this region. Such information is most valuable to investigate the genetics behind common diseases and is collected by the International HapMap Project (http://en.wikipedia.org/wiki/Haplotype). Hardy-Weinberg equilibrium (HWE) It states that, under certain conditions, after on generation of random mating, the genotype frequencies at a single gene locus will become fixed at a particular equilibrium value. It also specifies that those equilibrium frequencies can be represented as a simple function of the allele frequencies at that locus.(http://en.wikipedia.org /wiki /Hardy-Weinberg  equilibrium).  heterozygous Possessing two different forms of a particular gene, one inherited from each parent. homozygous Possessing two identical forms of a particular gene, one inherited from each parent. linkage The association of genes and/or markers that lie near each other on a chromosome. Linked genes and markers tend to be inherited together.  102  Appendix A. Glossary of some genetic terms linkage disequilibrium (LD) The non-random association of alleles at two or more loci on a chromosome. It describes a situation in which some combinations of alleles or genetic markers occur more or less frequently in a population than would be expected from a random formation of haplotypes from alleles based on their frequencies. In population genetics, linkage disequilibrium is said to characterize the haplotype distribution at two or more loci (http://en.wikipedia.org/wiki/Linkage disequilibrium). locus The place on a chromosome where a specific gene is located, a kind of address for the gene. The plural is ”loci,” not ”locuses.” non-coding DNA The strand of DNA that does not carry the information necessary to make a protein. The non-coding strand is the mirror image of the coding strand and is also known as the antisense strand. nucleotide One of the structural components, or building blocks, of DNA and RNA. A nucleotide consists of a base (one of four chemicals: adenine, thymine, guanine, and cytosine) plus a molecule of sugar and one of phosphoric acid. phenotype The observable traits or characteristics of an organism, for example hair color, weight, or the presence or absence of a disease. Phenotypic traits are not necessarily genetic. polymerase chain reaction (PCR) A fast, inexpensive technique for making an unlimited number of copies of any piece of DNA. Sometimes called ”molecular photocopying,” PCR has had an immense impact on biology and medicine, especially genetic research. polymorphism A common variation in the sequence of DNA among individuals. primer A short oligonucleotide sequence used in a polymerase chain reaction. probe A piece of labeled DNA or RNA or an antibody used to detect the function of a gene.  103  Appendix A. Glossary of some genetic terms single nucleotide polymorphism Common, but minute, variations that occur in human DNA at a frequency of one every 1,000 bases. These variations can be used to track inheritance in families. SNP is pronounced ”snip”. thymine One of the four bases in DNA that make up the letters ATGC, thymine is the ”T”. The others are adenine, guanine, and cytosine. Thymine always pairs with adenine. wild-type allele The allele designated as the standard (“normal”) for a strain of organism.  104  

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
United States 8 2
China 6 28
Japan 4 0
Russia 1 0
Canada 1 0
City Views Downloads
Beijing 6 0
Ashburn 5 0
Tokyo 4 0
Mountain View 2 2
Redmond 1 0
Unknown 1 0
Regina 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0066454/manifest

Comment

Related Items