UNIVERSIDADE FEDERAL DE SÃO CARLOS CENTRO DE CIÊNCIAS EXATAS E TECNOLOGIA PROGRAMA INTERINSTITUCIONAL DE PÓS-GRADUAÇÃO EM ESTATÍSTICA UFSCar-USP DANIELE CRISTINA TITA GRANZOTTO TRANSMUTATION MAPS: MODELING, STRUCTURAL PROPERTIES, ESTIMATION AND APPLICATIONS Thesis submitted to the Departamento de – Des/UFSCar and to the Instituto de Ciências Matemáticas e de Computação – ICMC-USP, in partial fulfillment for the PhD degree in Statistics - Interinstitucional Program of Pos-Graduation in Statistics UFSCar-USP. Advisor: Prof. Dr. Francisco Louzada Neto São Carlos January 2017 UNIVERSIDADE FEDERAL DE SÃO CARLOS CENTRO DE CIÊNCIAS EXATAS E TECNOLOGIA PROGRAMA INTERINSTITUCIONAL DE PÓS-GRADUAÇÃO EM ESTATÍSTICA UFSCar-USP DANIELE CRISTINA TITA GRANZOTTO MAPAS DA TRANSMUTAÇÃO: MODELAGEM, PROPRIEDADES ESTRUTURAIS, ESTIMAÇÃO E APLICAÇÕES Tese apresentada ao Departamento de Estatística – Des/UFSCar e ao Instituto de Ciências Matemáticas e de Computação – ICMC-USP, como parte dos requisitos para obtenção do título de Mestre ou Doutor em Estatística - Programa Interinstitucional de Pós-Graduação em Estatística UFSCar-USP. Orientador: Prof. Dr. Francisco Louzada Neto São Carlos Janeiro de 2017 Ficha catalográfica elaborada pelo DePT da Biblioteca Comunitária UFSCar Processamento Técnico com os dados fornecidos pelo(a) autor(a) G765t Granzotto, Daniele Cristina Tita Transmutation maps : modeling, structural properties, estimation and applications / Daniele Cristina Tita Granzotto. -- São Carlos : UFSCar, 2017. 101 p. Tese (Doutorado) -- Universidade Federal de São Carlos, 2016. 1. E-extension. 2. Log-logistic model. 3. Weibull model. 4. Cubic transmutation. 5. Quadratic transmutation. I. Título. Ao meu esposo e filhos Agradecimentos Acredito que todo caminho em que nos aventuramos, nos deparamos com pessoas que nos ajudam a aliviar as angústias,auxiliando e incentivando a seguir sempre em frente. Nestes quatro anos de doutoramento são muitos os que me auxiliaram e apoiaram, me ajudando a ter forças para concluir mais esta importante etapa. Primeiramente quero agradecer ao meu orientador e amigo Professor Francisco Louzada Neto que, além de compartilhar seus conhecimentos, acreditou no meu trabalhou e me incentivou desde o mestrado. Tenho muito a agradecer ao meu esposo, Carlos, que colaborou imensamente para a realização deste projeto. Aos meus filhos, Iasmim e João, que apesar da minha ausência demonstraram muita compreensão, sempre me motivando a continuar em busca desta conquista. Ao meu irmão e sua esposa, Marcos e Roberta, que muitas vezes me acudiram e auxiliaram. Muito obrigado ao professores do Deparamento de Estatística da Universidade Estadual de Maringá, departamento onde me graduei e hoje integro o corpo docente. Aos professores do Programa Interinstitucional de Pós-Graduação em Estatística da Universidades Federal de São Carlos e da Universidade de São Paulo por me acolherem e compartilharem comigo seus conhecimentos. E enfim, a todos que contribuíram para a realização deste trabalho, seja de forma direta ou indireta, fica registrado aqui, o meu muito obrigado! Resumo Inicialmente, usamos os mapas de transmutação quadráticos para compor um novo modelo de probabilidade: a distribuição log-logística transmutada. Mapas de transmutação são uma forma conveniente de construção de novas distribuições, em especial de sobrevivência/confiabilidade, e compreendem a composição funcional da função de distribuição acumulada e da função de distribuição acumulada inversa (quantil) de um outro modelo. Uma descrição detalhada de suas propriedades, tais como, momentos, quantis, estatística de ordem, dentre outras estatísticas, juntamente com o estudo de sobrevivência e métodos de estimação clássico e Bayesiano, também fazem parte deste tra- balho. Focando em análise sobrevivência, incluímos no estudo duas situações práticas comumente encontradas: a presença de variáveis regressoras, através do modelo de regressão transmutado log-logístico, e a presença de censura à direita. Em um segundo momento, buscando um modelo mais flexível que o transmutado, apresentamos uma generalização para esta classe de modelos, as distribuições transmutadas de rank cúbico. Usando a metodologia apresentada nesta primeira generalização, dois modelos foram considerados para compor as novas distribuições transmutadas cúbica: os modelos log-logístico e Weibull. Diante de problemas apresentados na classe transmutada de ordens quadrática e cúbica (tal como o espaço paramétrico restrito do parâmetro de transmu- tação λ), propomos neste trabalho, uma nova família de distribuição. Esta família, a qual chamamos e-transmutada ou e-extendida, é tão simples quanto o modelo transmutado, por incluir um único parâmetro ao modelo base, porém mais flexível do que a classe de modelos transmutados, sendo esta classe um caso particular da família proposta. Além disso, apresenta propriedades impor- tantes, como ortogonalidade entre os parâmetros do modelo base e o parâmetro de e-transmutação, e espaço paramétrico não restrito para o parâmetro de e- transmutação ω, que é definido em toda reta real. Estudos de simulação e aplicações a dados reais foram realizados para todos os modelos e generaliza- ções propostas. Palavras-Chave: E-extensão, modelo log-logístico, modelo Weibull, transmu- tação cúbica, transmutação quadrática. iii Abstract Initially, we use the quadratic transmutation maps to compose a new proba- bility model: the transmuted log-logistic distribution. Transmutation maps are a convenient way of constructing new distributions, in particular survival ones. It comprises the functional composition of the cumulative distribution function of one distribution with the inverse cumulative distribution (quantil) function of another. Its comprehensive description of properties, such as moments, quan- tiles, order statistics etc., along with its survival study and the classical and Bayesian estimation methods, are also part of this work. Focusing on analysis of survival, the study included two practical situations commonly found: the presence of regression variables, through the transmuted log-logistic regression model, and the presence of right censorship. In a second moment, searching for a more flexible model than the transmuted, we present its generalization, the transmuted distributions of cubic rank. Using the methodology presented in this first generalization, two models were considered to compose the new cu- bic transmuted distributions: the log-logistic and Weibull models. Faced with problems presented in the transmutated class of quadratic and cubic orders (such as the restricted parametric space of the transmutation parameter λ), we propose in this work, a new family of distribution. This family, which we call e-transmuted or e-extended, is as simple as the transmuted model, because it includes a single parameter to the base model, but more flexible than the class of transmuted models, once the transmuted is a particular case of the proposed family. In addition, the nem family presents important properties such as, or- thogonality between the baseline model parameters and the e-transmutation parameter, along with unrestricted parametric space for the ω e-transmutation parameter, which is defined on the real line. Simulation studies and real data applications were performed for all proposed models and generalizations. Keywords: E-extension, log-logistic model, Weibull model, cubic transmuta- tion, quadratic transmutation. iv Contents List of Figures vii List of Tables x 1 Introduction 1 1.1 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Organization of chapters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Transmuted model 5 2.1 The model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Properties of the model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.1 Moments and quantiles . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.2 Random number generation . . . . . . . . . . . . . . . . . . . . . . . 9 2.2.3 Survival analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.4 Order statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3 Parameter estimation of the regression model . . . . . . . . . . . . . . . . . 13 2.4 Simulation study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.5 Application: Tabapua cattle breed data . . . . . . . . . . . . . . . . . . . . . 15 2.5.1 Including covariates in the model . . . . . . . . . . . . . . . . . . . . 18 3 Bayesian and profile analysis 24 3.1 Hierarchical transmuted log-logistic model . . . . . . . . . . . . . . . . . . . 24 3.2 Tabapua cattle breed data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2.1 Influence analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.3 Transmuted log-logistic in the presence of censored lifetime . . . . . . . . . 38 3.3.1 Profile likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.3.2 Adjusted profile likelihood . . . . . . . . . . . . . . . . . . . . . . . . 40 3.3.3 Modified profile likelihood . . . . . . . . . . . . . . . . . . . . . . . . 41 3.3.4 Residual analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.3.5 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.3.6 Application to real dataset . . . . . . . . . . . . . . . . . . . . . . . 43 4 Cubic ranking transmuted model 51 4.1 The cubic ranking transmutation map . . . . . . . . . . . . . . . . . . . . . 51 4.2 Cubic rank transmuted weibull distribution . . . . . . . . . . . . . . . . . . 52 4.3 Cubic rank transmuted log-logistic distribution . . . . . . . . . . . . . . . . 55 4.4 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.5 Numerical studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 v 4.5.1 Simulation study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.5.2 Application: Carbon fibers data . . . . . . . . . . . . . . . . . . . . . 59 4.5.3 Application: Cattle sexual precocity data . . . . . . . . . . . . . . . 64 5 e-Transmuted family of distribution 67 5.1 Formulation of the model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.2 Hazard and related functions . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.2.1 Related distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.3 E-extended Weibull . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.4 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.4.1 Underlying exponential family structure . . . . . . . . . . . . . . . . 74 5.5 Weibull e-extended Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.5.1 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.5.2 Numerical experiments for the Weibull e-transmuted model . . . . . 78 5.5.3 Goodness-of-fit to the Weibull regression e-transmuted model . . . . 80 6 Conclusions and perspectives 86 A Hessian matrix 88 A.1 The hessian matrix of transmuted log-logistic model of order 1 . . . . . . . 88 A.2 The hessian matrix of transmuted Weibull of order 2 . . . . . . . . . . . . . 89 A.3 The hessian matrix of transmuted log-logistic of order 2 . . . . . . . . . . . 91 A.4 Hessian matrix of e-transmuted model . . . . . . . . . . . . . . . . . . . . . 93 B Tables 95 References 98 List of Figures 2.1 pdf of transmuted log-logistic distribution for µ = 1 and different values of λ and β. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Comparing the theoretical curve with generated values of transmuted log- logistic distribution for µ = 2, β = 5 and: (a) λ = −0.2; (b) λ = 0.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Some examples of hazard curves with the maximum points for λ = 0,−1 and 1, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4 Survival and hazard curves. . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.5 Coverage probability for γ1, γ2, β and λ, respectively, by considering the nominal value of 95%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.6 (a) TTT Plot, (b) histogram, (c) reliability curves and (d) hazard estimate curve, with the Tˆmax and the Tmax 95% confidence interval. . . . . . . . . . 19 2.7 TTT Plot by considering the covariates: prp (left panel) and period (right panel). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.8 Cross Validation results: Akaike weights for the complete sample (left panel) and without influential points (right panel). . . . . . . . . . . . . . . . . . . 21 2.9 Estimated residuals versus the estimated empirical survival for the residuals: log-logistic (left panel), transmuted log-logistic distributions (right panel). . 22 2.10 Hazard estimate curve, with the Tˆmax. . . . . . . . . . . . . . . . . . . . . . 23 3.1 Boxplot of times. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2 Traceplots and convergence plots, respectively, for: (a,f): µ; (b,g): β; (c,h): λ; (d,i): σ2 and; (e,j): θ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3 Marginal posterior densities for: (a) µ, (b) β, (c) λ, (d) σ2 and (e) θ . . . . 29 3.4 (a) hazard estimate curve, with the Tˆmax − 700 and the Tmax − 700 95% credible interval, (b) survival curves and (c) histogram. . . . . . . . . . . . . 30 3.5 Likelihood distance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.6 (a) Hazard estimate curve, with the Tˆmax − 700, (b) survival curves and (c) histogram. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.7 Mean for the parameters estimated: (a) µ, (b) β and (c) λ, by using × P, M AP and + MP methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.8 Mean Square Error (MSE) of the parameters estimated: (a) µ, (b) β and (c) λ, by using × P, M AP and + MP methods. . . . . . . . . . . . . . . . . . . 44 3.9 Coverage probability for the parameters estimated: (a) µ, (b) β and (c) λ, by using × P, M AP and + MP methods. . . . . . . . . . . . . . . . . . . . 45 3.10 (a) Boxplot of the observed times and (b) TTT plot. . . . . . . . . . . . . . 45 vii 3.11 (a) Survival curves: estimated versus empirical; (b) Hazard estimate curve, with the Tˆmax. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.12 (P) Profile, (A) Adjusted Profile and (M) Modified Profile relative log- likelihood for the Linezolid data. . . . . . . . . . . . . . . . . . . . . . . . . 47 3.13 (a) Cooks's distance and (b) Likelihood distance after response pertubation. 48 3.14 (a) Martingale residuals; (b) Deviance residuals. . . . . . . . . . . . . . . . . 49 3.15 (P) Profile, (A) Adjusted Profile and (M) Modified Profile relative log- likelihood for the Linezolid data. . . . . . . . . . . . . . . . . . . . . . . . . 50 4.1 Upper panels: Density and cumulative distribution functions; Lower panels: Survival, hazard and cumulative hazard functions, all of which are plotted for different values of model parameters. . . . . . . . . . . . . . . . . . . . . 54 4.2 Upper panels: Density and cumulative distribution functions; Lower panels: Survival, cumulative hazard and hazard functions, all of which are plotted for different values of model parameters. . . . . . . . . . . . . . . . . . . . . 56 4.3 Upper panels: TTT plot and hazard curves empirical and estimated by cu- bic rank transmuted Weibull distribution; Lower panels: Cuimulative and density functions estimated by transmuted Weibull and Weibull distributions. 61 4.4 P-P plot estimated for transmuted cubic and quadratic rank Weibull and Weibull distributions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.5 (a) Martingale residuals; (b) Deviance residuals. . . . . . . . . . . . . . . . . 62 4.6 (a) Cook's distance; (b) Likelihood distance. . . . . . . . . . . . . . . . . . . 63 4.7 Cumulative curves: empirical vs estimated for cubic rank transmutedWeibull distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.8 Upper Panels: histogram and survival curves; Lower Panel: hazard estimate curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.9 P-P plot estimated for transmuted cubic and quadratic rank log-logistic and log-logistic distributions, respectively. . . . . . . . . . . . . . . . . . . . . . . 65 5.1 The effect of ω in the cdf curves of the e-extended family with G(y) expo- nentially distributed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.2 The effect of ω in the hazard curves with G(y) exponential. . . . . . . . . . 70 5.3 The effect of ω in the density curves of the e-extended and transmuted families, with G(y|θ) exponential. . . . . . . . . . . . . . . . . . . . . . . . . 72 5.4 Probability density and cumulative distribution functionWeibull e-transmuted model, respectively, for different values of ω and φ = 0.7. . . . . . . . . . . . 73 5.5 The effect of ω in the hazard curves with G(y) Weibull(µ, β): (a) β > 1; (b) β < 1; and (c) β = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.6 Coverage probability of intervals for ω, µ and β parameters, respectively, generated by using the p-bootstrap method for φ = 2 fixed and different values of ω and β. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.7 Coverage probability of intervals for γ1, γ2, µ and β parameters, respectively, based on Monte Carlo experiment for φ = 2 fixed and different values of γ1, γ2 and β. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.8 Probability of error type I by using Anderson-Darling and Cramér-von-Mises statistics for different combinations of the parameters β and ω and fixed φ = 2 by considering α1 = 5% and 10% . . . . . . . . . . . . . . . . . . . . . 83 5.9 Power of the tests Anderson-Darling and Cramér-von Misses for n = 100 and different values of β. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.10 Cumulative distribution curves for the log-logistic andWeibull e-Transmuted models when β = 0.1 and 1.5, respectively; ω and φ were set in the same values then estimated in power of the test experiment. . . . . . . . . . . . . 84 5.11 Power of the tests Anderson-Darling and Cramér-von Mises with sample size n = 100 and for different combinations of the parameters β and ω and fixed φ = 2 when Y ∼Weibull(φ, β). . . . . . . . . . . . . . . . . . . . . . . . . . . 85 List of Tables 2.1 MLEs for the parameters of the transmuted log-logistic, log-logistic, Weibull and exponential distributions. . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 Computed −2 log likelihood and AIC values for the log-logistic, transmuted log-logistic, Weibull and Exponential distributions. . . . . . . . . . . . . . . 18 2.3 MLEs considering the log-logistic and transmuted log-logistic regression model. 20 2.4 MLEs, standard error and the 95% confidence interval considering the trans- muted log-logistic regression model in the sample without outliers. . . . . . 21 3.1 Posterior model summary of the hierarchical transmuted log-logistic model parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2 95% Credible Interval of parameters estimated. . . . . . . . . . . . . . . . . 27 3.3 Gelman and Rubin's criterion to verify the parameters convergence of the hierarchical transmuted log-logistic distribution. . . . . . . . . . . . . . . . . 27 3.4 RC (in %) and the corresponding TRC, MRC and LD(I). . . . . . . . . . . . 33 3.5 Posterior model summary of the hierarchical transmuted log-logistic model parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.6 95% Credible Interval of parameters estimated. . . . . . . . . . . . . . . . . 36 3.7 Estmates for the parameters of transmuted log-logistic model estimated by using profile methods, considering right censored lifetimes. . . . . . . . . . . 46 3.8 Wald confidence limits by considering 95% of confidence. . . . . . . . . . . . 46 3.9 RC (in %) and the corresponding TRC, MRC and LD(I). . . . . . . . . . . . 49 3.10 Estimates for the parameters of transmuted log-logistic model estimated by using profile methods after removing influential observations. . . . . . . . . 50 4.1 Means of estimates of all parameters and mean square errors (MSE), for different sample sizes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.2 Coverage probabilities of the confidence intervals for different sample sizes. 60 4.3 Selection criteria estimated for six different lifetime distributions. . . . . . . 61 4.4 MLEs of the parameters of the transmuted quadratic and cubic ranks log- logistic and log-logistic distributions. . . . . . . . . . . . . . . . . . . . . . . 64 4.5 Computed −2 log likelihood and AIC values for the log-logistic and trans- muted log-logistic of orders 1 and 2 distributions. . . . . . . . . . . . . . . . 65 x 5.1 Results of a Monte Carlo experiment for the Weibull e-transmuted model with φ = 2 and different combinations of ω and β values. The estimated val- ues presented are the median values of the 5000 generated samples. Within each horizontal block, the generated samples were subsampled from the larger one which was generated by using the inverse CDF method. . . . . . 79 5.2 Results of a Monte Carlo experiment for the Weibull regression e-transmuted model with φ = 2 and different combinations of γ1, γ2 and β values. The estimated values presented are the median values of the 5000 generated sam- ples. Within each horizontal block, the generated samples were subsampled from the larger one which was generated by using the inverse CDF method. 81 B.1 MLEs for the original sample generated with different parameter values and sample sizes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 B.2 Probability of coverage by considering 95% of confidence. . . . . . . . . . . 96 Chapter 1 Introduction Survival analysis researchers are usually concerned with the proposition of new survival probability models. The challenge is to derive statistical survival probability models or simply survival distributions of real world lifetime phenomena that can represent more consistently the random behavior of experimental observations. In this context, the literature on proposing new survival distributions is rich and grow- ing rapidly. There are many papers extending standard survival distributions designed to serve as statistical survival models for a wide range of real lifetime phenomena, which do not follow any of the standard survival distributions. Interested readers can refer to Ghi- tany (2001), Marshall and Olkin (2007a), and their references, in which they discuss some common approaches for constructing lifetime distributions. Generally transmutation maps are a convenient way of constructing new distributions, in particular, survival ones. According to Shaw and Buckley (2009), transmutation maps comprise the functional composition of the cumulative distribution function of one distribu- tion with the inverse cumulative distribution (quantile) function of another. Motivated by the need for parametric families of rich and yet tractable distributions in financial mathe- matics Shaw and Buckley (2009) used a transmutation map of a non-Gaussian distribution. After that, some studies involving quadratic rank transmutation map could be observed in other application areas such as survival analysis and reliability. Aryal and Tsokos (2009) proposed a generalization of the extreme value distribution by using quadratic rank trans- mutation maps and applied this new distribution to analyze snow fall data in Midway Airport in the state of Illinois, USA. Furthermore, Aryal and Tsokos (2011) proposed the transmuted Weibull distribution as a generalization of the Weibull probability distribution and its usefulness was illustrated using two published data sets. In this work, we propose for the first time the transmuted log-logistic model, obtained by considering a quadratic ranking transmutation map. A direct advantage of our modeling is the fact that the log-logistic distribution has a larger range of choices for the shape of the hazard function than Weibull or extreme value distributions do (see Ibrahim et al. (2005) 1 2 INTRODUCTION and Bennett (1983)). For instance, the log-logistic distribution can accommodate an uni- modal hazard function in detriment of its competitors. The log-logistic distribution is the probability distribution of a random variable whose logarithm has a logistic distribution, Lawless (2011). It is similar in shape to the log-normal distribution but has heavier tails and its cumulative distribution function can be written in closed form, unlike that of the log-normal. Moreover, we propose to apply the transmuted log-logistic model in a Bayesian context. In order to fit this new model, the subjective Bayesian analysis was used. To do so, we used the half-Cauchy prior distribution, cited by various authors (Polson and Scott (2012) and Gelman (2006)), as an alternative prior to an inverse Gamma distribution. Specially, Gelman (2006) made use of this particular prior for variance parameters in hierarchical models, which are our case. Although several studies involving quadratic rank transmutation maps can be seen in many areas, such as survival analysis and reliability, they do not provide a regression approach and do not consider the presence of censorship. In our work, we particularly consider the presence of right censorship and, in order to fit this new model, we used profile methods, such as pure, adjusted and modified profile; see, for example, Barndorff- Nielsen and McCullagh (1993), Severini (1998) and Ferrari et al. (2007). Note that the pure likelihood method cannot estimate the parameters in samples with high percentages of censored points, which is what our application case showed. A new order of transmuted models was introduced in Chapter 4, the cubic ranking transmutation map. This new order increases the flexibility of transmuted models and is able to analyze more complex data, for example, situations with bimodal hazard rates. In order to exemplify the applicability of a cubic ranking transmutation map, we used two of the best known distributions in reliability fields as a base: Weibull and log-logistic. Several mathematical properties of these new distributions, transmuted Weibull of order 2 and transmuted log-logistic of order 2, were derived. Although transmutation maps are a convenient way of constructing new distributions, the restricted parametric space of the extra parameter λ may be a problem in some sit- uations. As an alternative to this class of model, we present the regression e-transmuted family of model (or exponential transmuted), which has the property that the extra pa- rameter λ can take any real value. Without the restricted parameter space for λ, this new model continues to present characteristics of a good model: it is simple, more flexible than the transmuted model and continues being interpretable. Furthermore, an influential analysis was carried out in order to provide an indication of bad model fitting or influential observations; see Ortega et al. (2003) and Fachini et al. (2008). The bootstrap and Monte Carlo methods were used in order to validate the results presented. Moreover, two quadratic distance measures of goodness-of-fit: Anderson-Darling and Cramér-von Mises are used to verify the quality of the adjustment of the e-transmuted 3 INTRODUCTION model. 1.1 Objectives The main objective of this work was to study the transmuted models wich is a con- venient way of constructing new distributions. Particularly, we have worked with the log- logistic distribution, which is a continuous probability distribution for a non-negative ran- dom variable (it can be used in survival analysis) and presents all functions with closed forms. A second goal was to generalize the ranking of transmutation map. Thus, we introduced the cubic ranking transmutation map model in the last chapter; however, the kth ranking transmutation map can be seen as a continuation of this work. Finally, we proposed a new family of distributions to solve the problem of the re- stricted parametric space of λ. 1.2 Organization of chapters The dissertation was organized as a collection of papers and divided in five chapters besides the Introduction and Conclusion. It should be mentioned that these five chapters were accepted papers for publication, or submitted to specialized journals in the field. Chapter 2 comprises the paper entitled The Transmuted Log-Logistic Distribution: Modeling, Inference and Application to a Polled Tabapua Race Time up to First Calving Data, published by Communications in Statistics: Theory and Methods in 2014 and the paper entitled The transmuted log-logistic regression model: A new model for time up to first calving of cows, published by Statistical Paper in 2015. In these papers, we proposed a new lifetime model called Transmuted Log-Logistic, which is a generalization of the well- known survival log-logistic model. We also presented the main statistical properties as well as their main functions and characteristics in the area of survival analysis. The usefulness of the model was exemplified by using a real dataset. The third chapter presents the paper entitled Hierarchical Transmuted Log-Logistic Model: A Subjective Bayesian Analysis, written by us and co-authored by Vera L. D. Tomazella and the paper entitled Likelihood Based Inference for the Transmuted Log- Logistic Model in the Presence of Right Censored Lifetime, both submitted to specialized journals. There, we covered the transmuted log-logistic model by the Bayesian point of view and include censored observations in the transmuted model. Chapter 4 presents a new order of ranking transmutation map, order 2, which we called the cubic transmuted model. The study was carried out in partnership with Professor Narayanaswamy Balakrishnan, a senior lecturer at McMaster University and the paper entitled Cubic Rank Transmuted Distributions: Inferential Issues and Applications. 4 INTRODUCTION Finally, in Chapter 5, we present a new family of distribution called e-transmuted. The study was carried out in partnership with Professor Karim Anaya Izquierdo, a lecturer at University of Bath and the paper entitled Goodness-of-Fit Testing to the Regression e-Transmuted Family of Distribution was presented as a result. This study was done in a sandwich stage, in England. Finally, in the last section, we present some conclusions and propose further researcher topics. Thi is followed by the bibliography and appendix sections. Chapter 2 The transmuted log-logistic model In this Chapter, we introduce the generalization of the log-logistic model by using a convenient way of constructing new distributions, i.e, the rank transmutation map method- ology. Transmutation maps comprise the functional composition of the cumulative distribu- tion function with the inverse cumulative distribution (quantile) function. The new three parameters model is called transmuted log-logistic model and its main properties and statistics are presented in this Chapter as well as the model construction. 2.1 The model Let X be a nonnegative random variable denoting the lifetime of an individual in some population. The random variable X is said to be log-logistically distributed, with scale parameter µ and shape parameter β, if its probability density function (pdf) is given by g (x) = eµβxβ−1 (1 + eµxβ) 2 , (2.1) where β > 0 and −∞ < µ < +∞ (note that the shape parameter β determines the tail behavior of the distribution). The corresponding cumulative distribution function is given by G (x) = eµxβ (1 + eµxβ) . (2.2) According to Shaw and Buckley (2009), a ranking quadratic transmutation map has the following simple form: F (x) = (1 + λ)G(x)− λG2(x), (2.3) for | λ |≤ 1, where G(x) is the cumulative distribution function of the baseline distribu- tion. Observe that at λ = 0, we just obtain the baseline cumulative distribution function. Following this idea, several authors have considered extensions of some common survival 5 6 TRANSMUTED MODEL distributions; see Aryal and Tsokos (2009) and Ibrahim et al. (2005). The construction of the ranking quadratic transmutation map considered here is sim- ple and intuitive. Let X1 and X2 be independent and identically distributed nonnegative random variables with distribution G(x). Then, consider{ Y d = min(X1, X2), with probability pi, Y d = max(X1, X2), with probability 1− pi, where 0 ≤ pi ≤ 1. The distribution of Y is evidently FY (x) = piPr(min(X1, X2) ≤ x) + (1− pi)Pr(max(X1, X2) ≤ x). We know that F(min)(x) = 1 − [1−G(x)]n and F(max)(x) = [G(x)]n; see, for example, Arnold et al. (1992). Hence, FY (x) = pi [ 1− (1−G(x))2 ] + (1− pi)G2(x) = 2piG(x) + (1− 2pi)G2(x). (2.4) If we take 2pi = λ, the distribution in (2.4) is the well-known ranking quadratic transmu- tation (or simply the transmuted distribution) as presented in (2.3). Following this idea, several authors have considered extensions from usual survival distributions; see for example Aryal and Tsokos (2009) and Ibrahim et al. (2005). Hence, from (2.2) and (2.3), the cdf of a transmuted generalized log-logistic distribution is F (x) = eµxβ (1 + eµxβ) 2 ( 1 + eµxβ + λ ) . (2.5) The corresponding pdf is given by f (x) = eµβxβ−1 [( 1 + eµxβ )− λ (eµxβ − 1)] (1 + eµxβ) 3 . (2.6) Note that the transmuted log-logistic distribution is an extended model to analyze more complex data and it generalizes some important distributions in reliability analysis. The log-logistic distribution is clearly a special case for λ = 0. Some of the possible shapes of the transmuted log-logistic distribution were illustrated by Figure 2.1 for selected values of parameters λ and β and for µ = 1. The λ is responsible for introducing skewness into the log-logistic distribution. This is in full agreement with Shaw and Buckley (2009), who pointed out that the introduction of skewness into a distribution is a direct effect of the transmutation maps. The following results show that the random variable X has a probability density func- 7 TRANSMUTED MODEL 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 x D en si ty β = 0.5 β = 0.75 β = 1 β = 1.5 β = 2 β = 3 λ = − 1 λ = − 0.5 λ = 0 λ = 0.5 λ = 1 λ = 1 Figure 2.1: pdf of transmuted log-logistic distribution for µ = 1 and different values of λ and β. tion given by (2.6). Proposition 2.1.1 Let X be a non-negative random variable with log-logistic distribution and the quadratic transmutation is given, respectively, by Equations (2.1) and (2.3), then X is distributed according to a transmuted log-logistic with probability function given by (2.6). Proof According to Mood et al. (1974), any function f(. ) with domain the real line and conterdomain [0,∞) is defined as a probability density function if and only if (i) f(x) ≥ 0 for all x and; (ii) ∫∞ −∞ f(x)dx = 1. Therefore, the first property is satisfied for all x > 0. The second property is shown below: ∫ +∞ 0 f(x)dx = ∫ +∞ 0 eµβxβ−1 [( 1 + eµxβ )− λ (eµxβ − 1)] (1 + eµxβ) 3 dx =∫ +∞ 0 eµβxβ−1(1 + eµxβ)−3dx+ ∫ +∞ 0 e2µβx2β−1(1 + eµxβ)−3dx−∫ +∞ 0 λe2µβx2β−1(1 + eµxβ)−3dx+ ∫ +∞ 0 λeµβxβ−1(1 + eµxβ)−3dx Replacing y = xβeµ, x = e−µ/βy1/β and dy = eµβxβ−1dx in the equation above, it follows that ∫ +∞ 0 f(x)dx = (λ+ 1) ∫ +∞ 0 (1 + y)−3dy + (1− λ) ∫ +∞ 0 y(1 + y)−3dy. 8 TRANSMUTED MODEL We know that B(υ, ω) = ∫ +∞ 0 zω−1 (1 + z)υ+ω (2.7) in terms of Gamma function, we have B(υ, ω) = Γ(υ)Γ(ω) Γ(υ + ω) = (υ − 1)!(ω − 1)! (υ + ω − 1)! . Then, the integrate is given by (λ+ 1) ∫ +∞ 0 (1 + y)−3dy = λ+ 1 2 (2.8) and (1− λ) ∫ +∞ 0 y(1 + y)−3dy = 1− λ 2 . (2.9) Finally, adding the result of (2.8) with (2.9), we have λ+ 1 2 + 1− λ 2 = 1. Therefore, Equation (2.6) is a density function. 2.2 Properties of the model 2.2.1 Moments and quantiles In this section, we shall present the moments and quantiles for the transmuted log- logistic distribution. The rth order moments of a transmuted log-logistic random variable X is given by E (Xr) = ∫ +∞ 0 xr eµβxβ−1 [( 1 + eµxβ )− λ (eµxβ − 1)] (1 + eµxβ) 3 = (1 + λ) ∫ +∞ 0 eµβxr+β−1 (1 + eµxβ)3 dx+ (1− λ) ∫ +∞ 0 e2µβxr+2β−1 (1 + eµxβ)3 dx = (1 + λ)e −µr β ∫ +∞ 0 y r β (1 + y)3 dy + (1− λ)e−µrβ ∫ +∞ 0 y r β +1 (1 + y)3 dy. Then, the rth order moments of a transmuted log-logistic is E (Xr) = pir β2 e −µr β (β − rλ) csc [ pir β ] . 9 TRANSMUTED MODEL Therefore, the expected value E(X) and variance V ar(X) of a transmuted log-logistic random variable X are, respectively, given by E (X) = e −µ β pi β2 (β − λ) csc [ pi β ] (2.10) and V ar (X) = e −2µ β pi β2 [ 2(β − 2λ) csc ( 2pi β ) − pi β2 (β − λ)2 csc2 ( pi β )] . (2.11) The qth quantile xq of the transmuted log-logistic distribution can be obtained from (2.5) as xq = [√ (1 + λ)2 − 4qλ+ 2q − (1 + λ) 2eµ(1− q) ]1/β . (2.12) In particular, the median is given by x0.5 = e −µ/β [√ 1 + λ2 − λ ]1/β . (2.13) 2.2.2 Random number generation According to Ross (2009), we have the following Proposition 2.2.1. Proposition 2.2.1 Suppose that U is a random variable with a uniform distribution pat- tern, U ∼ Uniform(0, 1). Therefore, for a distribution function F , the random variable X defined by X = F−1(U) has distribution function F (x). Proposition 2.2.1 shows that we can generate a continuous random variable X from its cumulative distribution function. For that, we generate an uniform value (0, 1) and obtain X = F−1(U). Using this method, we generate random numbers of a transmuted log-logistic distribution when parameters µ, β and λ are known as follows, eµxβ (1 + eµxβ) 2 ( 1 + eµxβ + λ ) = u, (2.14) where u ∼ U(0, 1). Isolating x in the equation above, we obtain x = [√ (1 + λ)2 − 4uλ+ 2u− (1 + λ) 2eµ(1− u) ]1/β , (2.15) similar to 2.12. The histograms of Figures 2.2-a and 2.2-b were generated using Equation (2.15), from which we observe a full agreement to the theoretical curve with generated values of the transmuted log-logistic distribution. 10 TRANSMUTED MODEL t D en si ty 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0. 0 0. 5 1. 0 1. 5 2. 0 2. 5 (a) t D en si ty 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0. 0 0. 5 1. 0 1. 5 2. 0 2. 5 (b) Figure 2.2: Comparing the theoretical curve with generated values of transmuted log-logistic dis- tribution for µ = 2, β = 5 and: (a) λ = −0.2; (b) λ = 0.2. 2.2.3 Survival analysis The reliability function R(t), which is the probability of an item not failing prior to some time t, is defined by R(t) = 1− F (t). The reliability function of a transmuted log-logistic distribution is given by R(t) = 1 + eµtβ(1− λ) (1 + eµtβ)2 . (2.16) The other characteristic of interest of a random variable is the hazard rate function defined by h(t) = f(t) 1− F (t) = f(t) R(t) , (2.17) which is an important quantity characterizing a life phenomenon. It can be loosely inter- preted as the conditional probability of failure, given the fact that it has survived to time t; see Lawless (2011). The hazard rate function for a transmuted log-logistic random variable is given by h(t) = βeµtβ−1 [ 1 + eµtβ − λ(eµtβ − 1)] (1 + eµtβ)(1 + eµtβ − λeµtβ) . (2.18) The hazard function of the transmuted log-logistic distribution has the following prop- erties: 11 TRANSMUTED MODEL (i) If λ = 0, we have the log-logistic hazard as a particular case given by h(t) = βeµtβ−1 (1 + eµtβ) . It can be easily shown that h(t) is increasing for β ≤ 1 and, for β > 1, the hazard function initially increases to the maximum Tmax = e −µ/β(β − 1)1/β and tends to zero for t→∞, Figure 2.3-a. (ii) If λ = −1, we have h(t) = 2βe2µt2β−1 (1 + eµtβ)(1 + 2eµtβ) , which is increasing for β ≤ 12 and unimodal for β > 12 with the maximum in Tmax =[ 3(β−1)+ √ 9β2−2β+1 4eµ ]1/β , Figure 2.3-b. (iii) If λ = 1, we have h(t) = 2βeµtβ−1 1 + eµtβ , which is increasing for β ≤ 1 and unimodal for β > 1 with the maximum in Tmax = exp [(ln(β − 1)− µ) /β], Figure 2.3-c. 0.0 0.5 1.0 1.5 2.0 0. 4 0. 6 0. 8 1. 0 1. 2 t h(t ) (a) 0.0 0.5 1.0 1.5 2.0 0. 6 0. 8 1. 0 1. 2 1. 4 1. 6 1. 8 t h(t ) (b) 0.0 0.5 1.0 1.5 2.0 1. 0 1. 5 2. 0 2. 5 3. 0 3. 5 t h(t ) (c) Figure 2.3: Some examples of hazard curves with the maximum points for λ = 0,−1 and 1, respectively. Some properties and characterization of the risk of an event occurring at a given point, conditioned by the fact that the event has not already occurred can be seen in Chechile (2003). Figure 2.4-a and b illustrate, respectively, the reliability and hazard behavior of a transmuted log-logistic distribution as the value of the parameter λ ranges from −1 to 1. In the survival analysis literature, the main relationship between the hazard and the reliability function is given by the cumulative hazard rate function H(t), that is H(t) = 12 TRANSMUTED MODEL 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 time Su rv iva l β = 0.5 β = 0.75 β = 1 β = 1.5 β = 2 β = 3 λ = − 1 λ = − 0.5 λ = 0 λ = 0.5 λ = 1 λ = 1 (a) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0. 0 0. 5 1. 0 1. 5 2. 0 time H az ar d β = 0.5 β = 0.75 β = 1 β = 1.5 β = 2 β = 3 λ = − 1 λ = − 0.5 λ = 0 λ = 0.5 λ = 1 λ = 1 (b) Figure 2.4: Survival and hazard curves. − lnR(t). In our case, H(t) = 2 ln(1 + eµtβ)− ln(1 + eµtβ − λeµtβ), (2.19) where H(0) = 0, H(t) is nondecreasing for all t ≥ 0 and limt→∞H(t) =∞. 2.2.4 Order statistics According to Aryal and Tsokos (2011), suppose that we have a system containing two components with each of them having independent and identical "base" distribution, for example, a log-logistic. If the components are connected in series, then the overall system will have the transmuted baseline distribution with λ = 1 whereas if the components are parallel, then the overall system will have a transmuted baseline. It has been observed that a transmuted log-logistic distribution with λ = 1 is the distribution of min(X1, X2) and a transmuted log-logistic distribution with λ = −1 is the distribution of the max(X1, X2) where X1 and X2 are independent and identically distributed 2-parameter log-logistic random variables. Now, suppose that X(1), X(2), . . . , X(n) denotes the order statistics of a random sample X1, X2, . . . , Xn from a continuous population with cumulative density function FX(x) and probability density function fX(x). Therefore, the pdf of X(j) is given by fX(j)(x) = n! (j − 1)!(n− j)!fX(x)[FX(x)] j−1[1− FX(x)]n−j , (2.20) for j = 1, 2, . . . , n. Then, the pdf of the jth order transmuted log-logistic random variable X(j) is given by 13 TRANSMUTED MODEL fX(j)(x) = n! (j − 1)!(n− j)! eµβxβ−1 (1 + eµxβ)2n+1 (2.21){[ (1 + eµxβ)− λ(eµxβ − 1) ] [eµxβ(1 + eµxβ + λ)]j−1 × [(1 + eµxβ)2 + eµxβ(1 + eµxβ + λ)]n−j } . Using Equation (2.21), the pdf of the nth order transmuted log-logistic statistics X(n) is given by fX(n)(x) = neµβxβ−1 (1 + eµxβ)2n+1 {[ (1 + eµxβ)− λ(eµxβ − 1) ] [ eµxβ(1 + eµxβ + λ) ]n−1} , (2.22) while the 1st order is given by fX(1)(x) = neµβxβ−1 (1 + eµxβ)2n+1 {[ (1 + eµxβ)− λ(eµxβ − 1) ] × (2.23)[ (1 + eµxβ)2 + eµxβ(1 + eµxβ + λ) ]n−1} . Note that λ = 0 yields the order statistics of the two parameter log-logistic distribution. 2.3 Parameter estimation of the regression model In order to present the parameter estimation, let us first define the transmuted log- logistic regression model. For that, let Y be a random variable denoting the lifetimes with pdf (2.6) and µ = γ(x) as a parameter depending on a covariate vector x = (1, x1, . . . , xp) ′ such as γ(x) = γ0 + γ1x1 + . . .+ γpxp. Then, the pdf (2.6) may be written as f (y | γ(x)) = e γ(x)βyβ−1 [( 1 + eγ(x)yβ )− λ (eγ(x)yβ − 1)]( 1 + eγ(x)yβ )3 , (2.24) where λ > 0, β > 0 and γ(x) is a regression defined above. The maximum likelihood estimation of the parameters that are inherent within the transmuted log-logistic regression probability distribution function is given by the follow- ing: Let y1, y2, . . . , yn be a sample of size n from a transmuted log-logistic distribution and x = (1, x1, x2, . . . , xp) ′ be a vector of covariates. Moreover, consider the following relation- ship between the vector of covariates and the parameters γ(x) = γ0 + γ1x1 + . . .+ γpxp. 14 TRANSMUTED MODEL Hence, the log-likelihood function is given by l = n lnβ + n∑ i=1 p∑ j=0 γjxij + n∑ i=1 ln [ 1 + e ∑p j=0 γjxijyβi − λ ( e ∑p j=0 γjxijyβi − 1 )] +(β − 1) n∑ i=1 ln yi − 3 n∑ i=1 ln ( 1 + e ∑p j=0 γjxijyβi ) . (2.25) Therefore, the normal equations are given by ∂l ∂γj = n∑ i=1  (1− λ)yβi xij exp [∑p j=0 γjxij ] 1 + exp [∑p j=0 γjxij ] yβi − λ [ exp [∑p j=0 γjxij ] yβi − 1 ]  + n∑ i=1 xij − 3 n∑ i=1  xijyβi exp [∑p j=0 γjxij ] 1 + yβi exp [∑p j=0 γjxij ]  , j = 0, 1, . . . , p ∂l ∂β = n β + n∑ i=1  (1− λ)yβi ln(yi) exp [∑p j=0 γjxij ] 1 + yβi exp [∑p j=0 γjxij ] − λ [ yβi exp [∑p j=0 γjxij ] − 1 ]  + n∑ i=1 ln(yi)− 3 + n∑ i=1 yβi ln(yi) exp [∑p j=0 γjxij ] 1 + yβi exp [∑p j=0 γjxij ]  , ∂l ∂λ = n∑ i=1 1− yβi exp [∑p j=0 γjxij ] 1 + yβi exp [∑p j=0 γjxij ] − λ [ yβi exp [∑p j=0 γjxij ] − 1 ] . The maximum likelihood estimator (MLE) θˆ = (γˆ0, . . . , γˆp, βˆ, λˆ) ′ can be obtained by solving the above nonlinear system of equations. It is usually more convenient to use nonlinear optimization algorithms such as quasi-Newton or Newton-Raphson algorithms to numerically maximize the log-likelihood function given in (2.25). Following Aryal and Tsokos (2011), in order to compute the standard error and asymp- totic confidence interval, we use the usual large sample approximation in which the MLEs of θ can be treated as being approximately a (p+ 3)-variate normal distribution. Hence as n→∞, the asymptotic distribution of the MLE (γˆ0, . . . , γˆp, βˆ, λˆ) is given by γˆ0 γˆ1 ... γˆp βˆ λˆ  ∼ N   γˆ0 γˆ1 ... γˆp βˆ λˆ  ,  Vˆ11 Vˆ12 . . . Vˆ1(p+3) Vˆ21 Vˆ22 . . . Vˆ2(p+3) ... ... . . . ... Vˆ(p+3)1 Vˆ(p+3)2 . . . Vˆ(p+3)(p+3)   , (2.26) 15 TRANSMUTED MODEL where, Vˆij = Vij |θ=θˆ is determined by the inverse of the Hessian matrix. Therefore, the approximate 100(1−α)% two sided confidence intervals for γj , β and λ are, respectively, given by γˆj ± zα/2 √ Vˆ(j+1),(j+1), βˆ ± zα/2 √ Vˆ(p+2),(p+2) and λˆ± zα/2 √ Vˆ(p+3),(p+3), where j = 0, . . . , p, zα is the upper α th percentile of the standard normal distribution. The needed Hessian matrix is shown in Appendix B. Different models can be compared by penalizing over-fitting by using the Akaike infor- mation criterion (Akaike, 1973), which intends to minimize the Kullback-Leibler divergence between the true distribution and the estimate from a candidate model. It is given by AIC = −2l(θˆ) + 2size(θ), where l(θˆ) denotes the log likelihood function evaluated at the max- imum and size(θ) is the number of model parameters. The model with the lowest value of this criterion (among all considered models) is regarded as the preferred model for describing the given dataset. 2.4 Simulation study This section presents the results of a Monte Carlo experiment on the finite sample behavior of the MLEs. For that, we generated according to a transmuted log-logistic re- gression distribution in the presence of two covariates. All results were obtained from 1000 Monte Carlo replications and the results are summarized in two tables, presented in Ap- pendix B. Table B.1 shows the relative difference of generated and estimated parameters and their respective standard errors over the 1000 MLEs, which decreased as the sample size increased. Table B.2 shows the coverage probability of 95% two sided confidence in- tervals for the model parameters, which are close to the nominal coverage for large sample sizes, though they usually differ less than 5% from the nominal coverage probability for the smallest sample size considered. In order to summarize the results of the simulation study, Figure 2.5 presents the coverage probability for γ1, γ2, β and λ, respectively. 2.5 Application: Tabapua cattle breed data Economic results related to beef cattle are directly related to their genetic prepotency, aiming at increasing the production of kilograms of meat per hectare, at a certain time and at less cost. Therefore, the sexual precocity is an important factor in the choice of cattle breed. According to the literature, heifers from temperate regions have sexual precocity, having their first calving when they are little more than two years old, half the time 16 TRANSMUTED MODEL l l l l l n Pr ob ab ilid ad e de C ob er tu ra 100 150 300 500 1000 0. 75 0. 80 0. 85 0. 90 0. 95 1. 00 l l l l l (a) l l l l l n Pr ob ab ilid ad e de C ob er tu ra 100 150 300 500 1000 0. 75 0. 80 0. 85 0. 90 0. 95 1. 00 l l l l l (b) l l l l l n Pr ob ab ilid ad e de C ob er tu ra 100 150 300 500 1000 0. 75 0. 80 0. 85 0. 90 0. 95 1. 00 l l l l l (c) l l l l l n Pr ob ab ilid ad e de C ob er tu ra 100 150 300 500 1000 0. 5 0. 6 0. 7 0. 8 0. 9 1. 0 l l l l (d) Figure 2.5: Coverage probability for γ1, γ2, β and λ, respectively, by considering the nominal value of 95%. observed in heifers from tropical regions, which have their first calving when they are approximately 4 years old on average (Pereira, 2000). In this context, a long-term study conducted on the Tabapua cattle breed was held at EMBRAPA, a Brazilian agricultural research institute, in order to infer the time up to the first calving in polled Tabapua cows. The Tabapua breed is strong, rustic and chunky, with excellent mothering ability and potential to gain weight. It originated from a crossover made in the municipality of Tabapua, in the state of São Paulo, Brazil, with crossbred zebu, which showed characteristics of spontaneous polled cattle coming from European temperate regions with high genetic potential for meat production and milk (Paro et al., 2013). The data consist of the time up to first calving, in days, of 17026 animals observed from 1983 to 2007. 17 TRANSMUTED MODEL This breed is distinguished by meekness, sexual precocity, high milk production, fer- tility, good meat quality, and adaptability to various regions and climate. However, it is mainly bred for meat purposes because it also has very favorable characteristics for this. It is also widely used in crosses with other breeds, generating hardy animals and good productivity (Pereira, 2000). According to the characteristics mentioned and especially its sexual precocity, the polled Tabapua breed provides an important economic result related to beef cattle since it aims at increasing production of kilograms of meat per hectare, at a certain time and at less cost. The resulting heifers from tropical regions have their first calving when they are approximately 4 years old, twice the time of heifers from temperate regions (Paro et al., 2013; Pereira, 2000). In this section, we illustrated the usefulness of the transmuted log-logistic distribution modeling the polled Tabapua breed data. Moreover, we provided a comparison model with the log-logistic, Weibull, and exponential distributions. Firstly, a brief descriptive analysis was made. The minimum observed time was 721 days, or approximately 24 months; the maximum observed was 6184 days or 206 months. The median time for the first calving is 1140 days (38 months) and even the first and third quantiles are 1068 and 1365 days, respectively. Figure 2.6-a shows the TTT plot (for example see Barlow and Campo (1975)), in order to verify the possible shape for the hazard function. If the TTT plot is concave, convex and then concave again, it indicates unimodal hazard, which is our case. The transmuted log-logistic, log-logistic, Weibull and exponential distributions were fitted to the data. Table 2.1 provides the MLEs, their corresponding standard errors and 95% confidence intervals. For these models, the computed −2 log likelihood and AIC are shown in Table 2.2. Both criteria provide evidence in favor of the transmuted log-logistic distribution. This result is corroborated by the fitted distribution under the data histogram in Figure 2.6-b and by the fitted survival functions under the Kaplan-Meier estimator, Figure 2.6-c. Table 2.1: MLEs for the parameters of the transmuted log-logistic, log-logistic, Weibull and expo- nential distributions. Model Parameter Estimate Standard 95% Confidence Interval Error Lower Upper Transmuted α −17.9079 0.1369 −18.1763 −17.6395 Log-Logistic β 3.0503 0.0220 3.0073 3.0933 λ −0.8139 0.0120 −0.8374 −0.7904 Log-Logistic α −16.7482 0.1003 −16.9448 −16.5516 β 2.5832 0.0160 2.5519 2.6145 Weibull µ 615.7300 3.1161 609.6300 621.8400 β 1.6078 0.0082 1.5916 1.6239 Exponential µ 547.2500 4.1940 539.0300 555.4700 18 TRANSMUTED MODEL After selecting the most appropriate model for describing the time up to first calv- ing (the transmuted log-logistic distribution), our interest turned to estimating the most probable time the first calving would occur, defined as the Tmax according to (2.18). The Tˆmax is equal to 1267 days (42.2 months). Its 95% confidence interval is given by IC[Tmax, 95%] = (1182.05; 1370.83) days. Figure 2.6-d shows the hazard estimate curve, with the Tˆmax and the Tmax 95% confidence interval. Furthermore, the median time up to first calving is equal to 1173.42 days. Therefore, half of the cows experiment the event of in- terest up to approximately 39 months and the mean time is 1259.63 days (or approximately 42 weeks), with a standard deviation of 398.33 days. Table 2.2: Computed −2 log likelihood and AIC values for the log-logistic, transmuted log-logistic, Weibull and Exponential distributions. Model Criteria −2 log likelihood AIC Transmuted 236201 236207 Log-Logistic 237996 238000 Weibull 242588 242592 Exponential 248747 248749 From the practical point of view, the use of a misspecificated distribution to represent the random behavior of the time up to first calving would, of course, lead to a misleading conclusion concerning the parameters of interest. For the sake of argument, by considering the log-logistic distribution instead of the transmuted log-logistic one, the Tˆmax is equal to 1501.62 days (approximately 50.0 months); almost 8 months longer than the Tˆmax obtained by considering the transmuted log-logistic one. This is another drawback that corroborates in favor of the proposed distribution. 2.5.1 Including covariates in the model In this section, we used a random sample of 500 animals observed and two covariates: the time the calf was born, before or after 2000 (period) and the age the first oestrus occurred (prp), up to or after one year of age. By using the nonparametric Log-Rank test to compare the survival distributions of two samples, we observed p-values < 0.0001 and 0.0248, respectively, for covariates prp and period (see the survival curves in Figure 2.10). These p-values showed us the significance of these covariates to describe the time of the first calf. Another nonparametric test was carried out (Peto-Peto test; see for example Lawless (2011)), showing almost the same 19 TRANSMUTED MODEL 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 TTT plot r/n f (a) time (days) D en si ty 0 500 1000 1500 2000 2500 3000 0. 00 00 0. 00 05 0. 00 10 0. 00 15 0. 00 20 Transmuted Log−Logistic Weibull Exponential (b) 0 500 1000 1500 2000 2500 3000 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 time (days) Su rv iva l Kaplan−Meier Transmuted Log−Logistic Weibull Exponencial (c) 1000 1500 2000 2500 3000 0. 00 0 0. 00 1 0. 00 2 0. 00 3 0. 00 4 time (days) h(t ) (d) Figure 2.6: (a) TTT Plot, (b) histogram, (c) reliability curves and (d) hazard estimate curve, with the Tˆmax and the Tmax 95% confidence interval. results (p-values < 0.0001 and 0.0006, respectively). In order to verify the behavior of the hazard function, Figures 2.7 show the TTT plots according to two covariates: the time the calf was born (after or before 2000) (period) and the age the first oestrus occurred (prp), up to or after one year. Interested readers can refer to Barlow and Campo (1975) for more information on TTT plotting. Overall, if the TTT plot is concave, convex and then concave again, it indicates an unimodal hazard, which is our case. Then, the transmuted log-logistic and the simple log-logistic distributions were fitted to the data. Table 2.3 provides the MLEs, their corresponding standard errors and 95% confidence intervals, as well as the computed −2 log likelihood (−2l) and the AIC. Both criteria provide evidence in favor of the transmuted log-logistic regression distribution. Ob- serve that, as the transmutation map presented skewness and the log-logistic distribution 20 TRANSMUTED MODEL 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 r/n f prp until one year prp after one year 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 r/n f calves after 2000 calves before 2000 Figure 2.7: TTT Plot by considering the covariates: prp (left panel) and period (right panel). is an asymmetric model, the mode point is higher, which is totally interpreted in our ap- plication. In general, the crossbreeds are at similar periods, increasing the chance of a first calf around a certain period of time and, thus, increasing the mode point. Table 2.3: MLEs considering the log-logistic and transmuted log-logistic regression model. Model Parameter Estimate Standard IC 95% Error Lower Upper −2l AIC γ0 −17.092 0.801 −18.667 −17.860 6925.3 6935.3 Transmuted γ1 −0.364 0.137 −0.632 −0.095 γ2 3.281 0.794 1.722 4.840 β 2.954 0.128 2.703 3.205 λ −0.764 0.084 −0.930 −0.598 γ0 −19.319 0.742 −20.777 −17.860 6936.8 6944.8 Log-Logistic γ1 −0.419 0.154 −0.722 −0.117 γ2 3.213 0.809 1.623 4.804 β 3.208 0.122 2.968 3.448 Furthermore, according to Burnham and Anderson (2002), we verified that if the MLE values are influenced by a specific experimental observation, they can be regarded as an out- lier. We considered the one-leaved-out approach on the basis of a cross validation scheme, in the sense that the parameters of the transmuted log-logistic regression were re-estimated 500 times. Each time, one specific observed time was withdrawn. We also estimated the standard errors, confidence intervals and the −2 log likelihood and AIC criteria. Moreover, the differences, ∆i, given by ∆i = AICi − AICmin, were calculated. After, the Akaike weights given by ωi = exp(−∆i/2) were obtained and plotted as shown on the left panel of Figure 2.8. The presence of outliers is clear. We decide to remove all the experimental points with ωi > 0.20, corresponding to observations 15, 143, 231, 242, 289 and 360. The Akaike weights in the sample, without the outliers, are shown on the right 21 TRANSMUTED MODEL panel of Figure 2.8. 0 100 200 300 400 500 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Index Ak ai ke W e ig ht s 0 100 200 300 400 500 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Index w i Figure 2.8: Cross Validation results: Akaike weights for the complete sample (left panel) and without influential points (right panel). Table 2.4 shows the MLEs for the parameters of the transmuted log-logistic regression model in the sample without the outliers. There seem to be no important changes in the MLEs, althought the −2l and AIC values are much smaller than those obtained when the complete sample is considered. Table 2.4: MLEs, standard error and the 95% confidence interval considering the transmuted log-logistic regression model in the sample without outliers. Model Parameter Estimate Standard IC 95% Error Lower Upper −2l AIC γ0 −20.611 1.065 −20.704 −18.518 6366.9 6376.9 Transmuted γ1 −0.430 0.145 −0.552 −0.145 without γ2 3.646 0.913 1.851 5.440 outliers β 3.549 0.168 3.218 3.880 λ −0.710 0.118 −0.943 −0.478 Moreover, as a goodness-of-fit procedure, we performed a residual analysis for both models (transmuted log-logistic and simple log-logistic regression ones) using the Cox- Snell residuals (Cox and Snell, 1968). The Cox-Snell residuals are defined as eˆi = Λˆ(ti | xi), where Λˆ(. ) is the cumulative hazard function of the adjusted model. If we consider the log-logistic and transmuted log-logistic models, the Cox-Snell residuals are, respectively, given by eˆiLL = ln [ 1 + exp(γ(x))tβi ] (2.27) and eˆiTLL = 2 ln [ 1 + exp(γ(x))tβi ] − ln [ 1 + exp(γ(x))tβi − λ exp(γ(x))tβi ] , (2.28) 22 TRANSMUTED MODEL where γ(x) = γ0 + γ1x1 + . . .+ γpxp. Figure 2.9 shows the estimated residuals versus the estimated empirical survival for the residuals. The criteria provide clear evidence in favor of the transmuted log-logistic regression model. Figure 2.10 shows the estimated hazard curve and the estimated most probable time for the first calf (Tmax) according to both covariates. The Tˆmax is equal to 29.52 months with a 95% confidence interval equal to (26.56, 32.49) if the prp occurs before the first year. For calves born until the year 2000, the Tˆmax is equal to 42.97 months with a 95% confidence interval equal to (41.43, 44.50). For calves born after the year 2000 or prp occurring after the first year, the Tˆmax is equal to 40.77 months with a 95% confidence interval equal to (39.35, 42.18). Considering the period prior to the year 2000, we observed an older age at the first calving, occurring close to 4 years. This time decreased to less than 3.5 years showing that the current reproductive techniques of livestock are more efficient than in the previous period. As the reproductive life begins earlier, cows have an increase number of conceptions and consequently increased their reproductive life by reducing the number of cull cows. lllllll lllll lllllllll lllllllllll lllllll llllll lllll lllll lll llllll llll lllll llll llll ll l ll l llll ll lll ll ll l l ll l l l l l 0 2 4 6 8 10 0 2 4 6 8 10 ei − lo g(S (ei )) llll llll llll lll lll lll lll lll llll llll lllll llll lll llll lll llll lll lll llll lll lll ll ll ll ll lll lll ll ll ll ll ll ll l l l l l l l l l 0 1 2 3 4 5 6 0 1 2 3 4 5 6 ei − lo g(S (ei )) Figure 2.9: Estimated residuals versus the estimated empirical survival for the residuals: log- logistic (left panel), transmuted log-logistic distributions (right panel). 23 TRANSMUTED MODEL 1000 1500 2000 2500 3000 0. 00 0 0. 00 2 0. 00 4 0. 00 6 0. 00 8 0. 01 0 0. 01 2 time (days) h(t ) prp until one year prp after one year 1000 1500 2000 2500 3000 0. 00 0 0. 00 1 0. 00 2 0. 00 3 0. 00 4 time (days) h(t ) calves after 2000 calves before 2000 Figure 2.10: Hazard estimate curve, with the Tˆmax. Chapter 3 Bayesian and profile analysis for the transmuted log-logistic model In this Chapter, we present the transmuted log-logistic model in a Bayesian context. For that, we used the hierarchical structure using the half-Cauchy prior proposed by Gelman (2006). Furthermore, we presented the transmuted log-logistic model in the presence of the censored lifetime and the profile method was used in the estimation process. 3.1 Hierarchical transmuted log-logistic model According to Chen and Ibrahim (2006), one of most common ways of combining several sources of information is though hierarchical modeling. Thus, the authors show us the relationship between the power prior and hierarchical models using regression models as an example. Moreover, Gelman (2006) points out that several studies using multilevel models are central to modern Bayesian statistics for both conceptual and practical reasons. The au- thors suggested the use of half-t family as a prior distribution for variance parameters such as the half-Cauchy distribution, which is a special conditionally-conjugate folded-non central-t family case of prior distributions for parameters representing the discrepancy. Even though several studies use the half-Cauchy prior for scale parameters (see for exam- ple Polson and Scott (2012)), Gelman (2006) used this prior not for scale, but for variance parameters and showed serious problems with the inverse-Gamma prior, the most com- monly used prior for variance component, Daniels (1999). In this study, we proposed to use the hierarchical models in two levels; to do so, suppose the hierarchical model is given as [X|µ, β, λ] ∼ f(x|µ, β, λ), µ|σ2 ∼ piµ(µ|σ2), β|θ ∼ piβ(β|θ), λ ∼ piλ(λ), σ2 ∼ ψσ(σ2) and θ ∼ ψθ(θ). The posterior distribution can be constructed as follows. Proposition 3.1.1 Let us suppose that, in the first stage, we considered a class Γ of priors 24 25 BAYESIAN AND PROFILE ANALYSIS that led to the following Γ = { pi(µ, β, λ|σ2, θ) : pi(µ, β, λ|σ2, θ) = piµ(µ|σ2)piβ(β|θ)piλ(λ) ∣∣ piµ being N(τ, σ 2), (τ, σ2) ∈ R× R+; piβ being HC(θ), θ ∈ R+; piλ being U(a, b), (a, b) ∈ R× R, a < b} . Moreover, in the second stage (sometimes called a hyperprior), that would consist of putting a prior distribution ψk(·) on the hyperparameters σ2 and θ where Ψ = { ψ(σ2, θ) : ψ(σ2, θ) = ψσ2(σ 2)ψθ(θ); ψσ2 is Gamma(α, ζ), (α, ζ) ∈ R+ × R+; ψβ is Gamma(η, ϑ), (η, ϑ) ∈ R+ × R+; α, ζ, η, ϑ are known and do not depend on any other hyperparameter} . Thus, the hierarchical log-logistic posterior distribution is written as pi(µ, β, λ|x) ∝ e µβθ σ [ xβ+α+η−3 [ (1 + eµxβ)− λ(eµxβ − 1)] (1 + eµxβ)3 ] × exp { − ( ζ + ϑ+ x 2σ2 )} . (3.1) Proof: The demonstration is direct. Note that the β parameter is supposed to be a half-Cauchy distribution whose proba- bility density function is given by f(x) = 2θ pi (x2 + θ2) , x > 0, θ > 0, (3.2) where θ is a scale parameter which has a broad peak at zero and, in limit, θ → ∞ this becomes a uniform prior density. However, large finite values for θ represent prior distribu- tions we call weakly informative. For example, Gelman (2006) shows us that, for θ = 25, the half-Cauchy is nearly flat, although not completely. 3.2 Tabapua cattle breed data In this section, we showed the usefulness of the hierarchical transmuted log-logistic model on modeling of the polled Tabapua breed data using Bayesian methods. The data consist of the time up to first calving, in days, of 17026 animals observed from 1983 to 2007, which was presented in Chapter 2. Figure 3.1-b shows the distribution of the time when the first calving occurred. The median time for the first calving was 1140 26 BAYESIAN AND PROFILE ANALYSIS days (38 months) and the first and third quantiles are 1068 and 1365 days, respectively. l ll ll l l l l 0 50 10 0 15 0 Ti m e (a) Figure 3.1: Boxplot of times. After an initial analysis, the hierarchical transmuted log-logistic was fitted to the data, as we specified in Section 3.1, The posterior samples were generated by the Metropolis- Hastings technique. Three chains of the dimension 100, 000 were considered for each pa- rameter, discarding the first 10, 000 iterations in order to eliminate the effect of the initial values and to avoid correlation problems. A lag size 10 was used, resulting in a final sample size of 10, 000. Tables 3.1 and 3.2 show the posterior summaries for the parameters and the 95% credible intervals considering the mentioned priors. Table 3.1: Posterior model summary of the hierarchical transmuted log-logistic model parameters. Parameter Mean Standard Percentiles Deviation 25% 50% 75% α −17.865 0.138 −17.958 −17.871 −17.775 β 3.043 0.022 3.029 3.044 3.058 λ −0.815 0.012 −0.823 −0.815 −0.807 σ2 900.100 822.500 333.800 640.400 1208.100 θ 198.800 195.100 60.171 139.800 273.400 The chain convergence was verified by the Gelman and Rubin's convergence diagnostic criterion, see for example Gelman and Rubin (1992), which demonstrated that these criteria are satisfied (Table 3.3). The convergence can also be seen in Figures 3.2-a to j. Furthermore, the marginal posterior densities for µ, β and λ, respectively, can be ana- lyzed in Figures 3.3-a to e. After estimating and analyzing the model convergence, Figure 3.4-a and b show, re- spectively, the hazard estimate curve, with the Tˆmax and the Tmax 95% credible interval; 27 BAYESIAN AND PROFILE ANALYSIS Table 3.2: 95% Credible Interval of parameters estimated. Parameter Equal-Tail Interval HPD Interval α −18.123 −17.576 −18.118 −17.571 β 2.997 3.085 2.997 3.084 λ −0.838 −0.791 −0.838 −0.790 σ2 96.474 3044.600 37.837 2516.300 θ 5.395 733.800 0.008 594.800 Table 3.3: Gelman and Rubin's criterion to verify the parameters convergence of the hierarchical transmuted log-logistic distribution. Parameter Estimate Upper Bound µ 1.0085 1.0060 β 1.0082 1.0057 λ 1.0020 1.0017 σ2 1.0016 1.0019 θ 1.0004 1.0009 the estimated vs empirical survival curves and the histograms allow us to see how well it fits a set of observations. Considering the hierarchical transmuted log-logistic fitting, the Tˆmax is equal to 1246.77 days (41.56 months) and its 95% credible interval is given by IC[Tmax, 95%] = (1160.04; 1352.86) days (see Figure 3.4-a). Furthermore, the median time up to first calving is equal to 1152.48 days (or approximately 38.42 months), and the mean time is 1240.13 days (or approximately 41.34 months), with a standard deviation equal to 13.34 months. 28 BAYESIAN AND PROFILE ANALYSIS (a) (b) (c) (d) (e) 0 5000 10000 15000 20000 1. 0 1. 1 1. 2 1. 3 1. 4 1. 5 Iteration µ Estimate Upper Bound (f) 0 5000 10000 15000 20000 1. 0 1. 1 1. 2 1. 3 1. 4 1. 5 Iteration β Estimate Upper Bound (g) 0 5000 10000 15000 20000 1. 0 1. 1 1. 2 1. 3 1. 4 1. 5 Iteration λ Estimate Upper Bound (h) 0 5000 10000 15000 20000 1. 0 1. 1 1. 2 1. 3 1. 4 1. 5 Iteration σ 2 Estimate Upper Bound (i) 0 5000 10000 15000 20000 1. 0 1. 1 1. 2 1. 3 1. 4 1. 5 Iteration θ Estimate Upper Bound (j) Figure 3.2: Traceplots and convergence plots, respectively, for: (a,f): µ; (b,g): β; (c,h): λ; (d,i): σ2 and; (e,j): θ 29 BAYESIAN AND PROFILE ANALYSIS −18.4 −18.2 −18.0 −17.8 −17.6 −17.4 0. 0 0. 5 1. 0 1. 5 2. 0 2. 5 D en si ty (a) 2.95 3.00 3.05 3.10 0 5 10 15 D en si ty (b) −0.86 −0.84 −0.82 −0.80 −0.78 −0.76 0 5 10 15 20 25 30 D en si ty (c) 0 2000 4000 6000 8000 10000 0e +0 0 2e −0 4 4e −0 4 6e −0 4 8e −0 4 1e −0 3 D en si ty (d) 0 200 400 600 800 1000 0. 00 0 0. 00 2 0. 00 4 0. 00 6 0. 00 8 D en si ty (e) Figure 3.3: Marginal posterior densities for: (a) µ, (b) β, (c) λ, (d) σ2 and (e) θ 30 BAYESIAN AND PROFILE ANALYSIS 0 500 1000 1500 2000 2500 3000 0. 00 0 0. 00 1 0. 00 2 0. 00 3 0. 00 4 time (days) h(t ) (a) 0 1000 2000 3000 4000 5000 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 time (days) Su rv iva l Kaplan−Meier Transmuted (b) time (days) D en si ty 0 1000 2000 3000 4000 5000 0. 00 00 0. 00 05 0. 00 10 0. 00 15 0. 00 20 Transmuted (c) Figure 3.4: (a) hazard estimate curve, with the Tˆmax − 700 and the Tmax − 700 95% credible interval, (b) survival curves and (c) histogram. 31 BAYESIAN AND PROFILE ANALYSIS 3.2.1 Influence analysis In this section, we will make an analysis of global influence for the data set given, using the transmuted log-logistical model in a Bayesian context. The first tool to assess the sensitivity analysis is measuring the global influence. We started with the case-deletion, in which we study the effect of withdrawing the ith ele- ment sampled. The first measure of global influence analysis is known as the generalized Cook's distance, which is defined as the standard norm of ζi = (αi, βi, λi, σ 2 i , θi) and ζˆ = (αˆ, βˆ, λˆ, σˆ2, θˆ) and is given by CDi(ζ) = [ ζi − ζˆ ]T [−L¨(ζ)] [ ζi − ζˆ ] (3.3) where L¨(ζ) can be approximated by the estimated covariance and variance matrix. Another way to measure the global influence is by the difference in likelihoods given by LDi(ζ) = 2 { l(ζˆ)− l(ζi) } . (3.4) Figures 3.5 show the likelihood distances where some possible influence points can be observed. In order to reveal the impact of these, the relative changes were measured as RCζj = ∣∣∣∣∣ ζˆj − ζˆj(I)ζˆj ∣∣∣∣∣× 100%, j = 1, . . . , p+ 1 where ζˆj(I) denotes the MLE of ζj after set I of observations was removed. As suggested by Lee et al. (2006), we use the total and maximum relative changes and the likelihood displacement. To reveal the impact of the detected influential observations, we used three measures as defined by Lee et al. (2006), TRC = np∑ i=1 ∣∣∣∣∣ ζˆi − ζˆ0iζˆi ∣∣∣∣∣ , MRC = maxi ∣∣∣∣∣ ζˆi − ζˆ0iζˆi ∣∣∣∣∣ and LD(l)(ζ) = 2{l(ζˆ − l(ζˆ0)} where TRC is the total relative changes, MRC the maximum relative changes and LD the likelihood displacement, with np (the number of parameters) and ζˆ 0 denoting MLE of ζ after set I of observations was removed. In order to analyze these influential points, we withdrew the identified points in Figure 3.5. Moreover, a percentage of points (0.01% to 5%) that stood out was also withdrawn. The results can be seen in Table 3.4. 32 BAYESIAN AND PROFILE ANALYSIS Figure 3.5: Likelihood distance. 33 BAYESIAN AND PROFILE ANALYSIS Table 3.4: RC (in %) and the corresponding TRC, MRC and LD(I). Removed Case Parameter RC TRC MRC LD(I) µ 8.271 86.438 68.336 142 Identify β 5.142 Points λ 68.336 σ2 2.877 θ 1.811 µ 0.844 3.598 1.358 217 β 0.818 0.1% λ 0.123 σ2 0.456 θ 1.358 µ 2.098 8.634 3.219 949 β 2.073 0.5% λ 0.098 σ2 1.144 θ 3.219 µ 3.257 9.782 3.257 1821 β 3.250 1% λ 0.368 σ2 1.700 θ 1.207 µ 5.843 17.064 5.855 3519 β 5.855 2% λ 0.393 σ2 2.911 θ 2.062 µ 8.096 22.086 8.162 5178 β 8.162 3% λ 1.460 σ2 3.111 θ 1.258 µ 10.458 28.558 10.547 6775 β 10.547 4% λ 0.785 σ2 6.566 θ 0.201 µ 12.463 34.232 12.627 8383 34 BAYESIAN AND PROFILE ANALYSIS β 12.627 5% λ 2.160 σ2 6.277 θ 0.704 Notice that when we withdrew the 10 most influential points, λ parameter was the most affected, RC ' 68%. Then, the hierarchical transmuted log-logistic was re-fitted to the data. The posterior samples were generated by the Metropolis-Hastings technique with one chain of the dimension 100000, discarding the first 10000 iterations in order to eliminate the effect of the initial values. To avoid the correlation problems, a lag size 10 was used resulting in a final sample size of 10000. Tables 3.5 and 3.6 show the posterior summaries for the parameters and the 95% credible intervals considering the mentioned priors for all removed cases, respectively. Table 3.5: Posterior model summary of the hierarchical transmuted log-logistic model parameters. Removed Parameter Mean Standard Percentiles Case Deviation 25% 50% 75% α −19.343 1.269 −20.522 −20.327 −17.975 Identify β 3.200 0.130 3.062 3.285 3.319 Points λ −0.258 0.519 −0.813 0.153 0.232 σ2 926.00 840.30 348.70 664.10 1224.00 θ 202.40 204.70 59.08 139.20 282.00 α −18.016 0.129 −18.098 −18.015 −17.933 β 3.068 0.021 3.055 3.068 3.082 0.1% λ −0.814 0.012 −0.822 −0.814 −0.806 σ2 904.20 826.40 334.50 651.90 1207.20 θ 196.10 194.80 56.84 137.20 270.50 α −18.240 0.135 −18.336 −18.238 −18.146 β 3.107 0.022 3.091 3.106 3.122 0.5% λ −0.816 0.012 −0.824 −0.816 −0.808 σ2 910.40 842.80 339.00 643.80 1203.30 θ 205.20 203.70 61.83 142.40 284.10 α −18.447 0.131 −18.536 −18.441 −18.352 β 3.142 0.021 3.127 3.142 3.157 1% λ −0.818 0.012 −0.826 −0.818 −0.810 σ2 915.40 853.70 337.50 647.80 1221.70 θ 196.40 195.20 56.19 136.90 273.10 α −18.909 0.144 −19.007 −18.907 −18.811 β 3.222 0.023 3.206 3.221 3.238 35 BAYESIAN AND PROFILE ANALYSIS 2% λ −0.818 0.013 −0.827 −0.818 −0.810 σ2 926.30 826.60 353.60 667.40 1235.00 θ 202.90 199.40 60.53 143.40 279.70 α −19.312 0.142 −19.410 −19.304 −19.222 β 3.292 0.023 3.277 3.291 3.308 3% λ −0.827 0.012 −0.835 −0.827 −0.819 σ2 928.10 849.10 352.20 668.20 1215.20 θ 196.30 193.30 56.37 139.00 275.50 α −19.734 0.154 −19.834 −19.737 −19.641 β 3.364 0.025 3.350 3.365 3.381 4% λ −0.821 0.013 −0.830 −0.821 −0.813 σ2 959.20 858.50 365.20 698.70 1267.90 θ 198.40 195.60 58.00 137.50 277.00 α −20.092 0.175 −20.221 −20.093 −19.967 β 3.428 0.028 3.407 3.428 3.448 5% λ −0.832 0.013 −0.842 −0.833 −0.824 σ2 956.60 854.00 372.90 704.10 1259.10 θ 200.20 198.30 59.91 141.00 273.20 It could be clearly seen that parameter λ was estimated in the first moment −0.838. It changed abruptly only when the identified points were removed, remaining close to −0.838 in all other cases, see in Tables 3.1 and 3.5. Furthermore, when we withdrew 0.1% from the sample, i.e., just 17 observations, we did not lose much information and we improved the fitted model. Figures 3.6-a, b and c show the fitted model, when using this number of removed cases. Finally, considering the hierarchical transmuted log-logistic fitting, the Tˆmax changes to 1267.71 compared to 1266.77 days (42.22 months) and its 95% credible interval is given by CI[Tmax, 95%] = (39.56; 45.47) months. Furthermore, the median time up to the first calv- ing is equal to 1172.41 days (or approximately 39.08 months), and the mean time is 1258.54 days (or approximately 41.95 months), with standard deviation equal to 13.11 months. 36 BAYESIAN AND PROFILE ANALYSIS Table 3.6: 95% Credible Interval of parameters estimated. Removed Case Parameter Equal-Tail Interval HPD Interval α −20.709 −17.754 −20.713 −17.760 Identify β 3.026 3.351 3.027 3.352 Points λ −0.833 0.293 −0.837 0.284 σ2 104.300 3096.600 43.226 2603.000 θ 5.211 746.000 0.143 601.600 α −18.2762 −17.7678 −18.2845 −17.7776 β 3.0289 3.11 3.0283 3.1093 0.1% λ −0.8365 −0.7896 −0.8372 −0.7906 σ2 94.0797 3151.6 40.6435 2574.8 θ 5.2045 728.6 0.1098 586.3 α −18.4967 −17.9843 −18.4981 −17.9863 β 3.0654 3.1476 3.0664 3.1483 0.5% λ −0.8391 −0.7905 −0.84 −0.7918 σ2 98.81 3143.7 28.4125 2605.9 θ 5.1449 745.1 0.0156 613.7 α −18.7121 −18.205 −18.7066 −18.2029 β 3.103 3.1855 3.1015 3.1829 1% λ −0.8407 −0.793 −0.8421 −0.7948 σ2 102.2 3172.8 34.7775 2578.2 θ 4.6184 727.6 0.0235 591.8 α −19.1832 −18.6319 −19.184 −18.6338 β 3.1769 3.266 3.1779 3.2663 2% λ −0.8415 −0.7927 −0.8428 −0.7944 σ2 106.2 3111.3 39.5096 2606.1 θ 5.4845 729.9 0.0306 600.2 α −19.5917 −19.0343 −19.5817 −19.0259 β 3.247 3.3364 3.2482 3.3373 3% λ −0.8504 −0.8013 −0.851 −0.8023 σ2 107.5 3224.1 30.9968 2597.3 θ 4.8582 710.8 0.0153 578.9 α −20.0252 −19.4068 −20.0321 −19.4161 β 3.312 3.4113 3.3143 3.4124 4% λ −0.8458 −0.7955 −0.8459 −0.7958 σ2 111 3307.1 50.2358 2671.5 θ 5.2325 734.3 0.0631 604.8 α −20.4078 −19.7676 −20.418 −19.7839 β 3.3762 3.4789 3.3778 3.4802 5% λ −0.8572 −0.8055 −0.8581 −0.8066 σ2 119 3273.2 56.2007 2652.8 θ 5.0109 743.8 0.0992 605.8 37 BAYESIAN AND PROFILE ANALYSIS 0 1000 2000 3000 4000 5000 0. 00 0 0. 00 1 0. 00 2 0. 00 3 0. 00 4 time (days) h(t ) (a) 0 1000 2000 3000 4000 5000 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 time (days) Su rv iva l Kaplan−Meier Transmuted (b) time (days) D en si ty 0 1000 2000 3000 4000 5000 0. 00 00 0. 00 05 0. 00 10 0. 00 15 0. 00 20 Transmuted (c) Figure 3.6: (a) Hazard estimate curve, with the Tˆmax−700, (b) survival curves and (c) histogram. 38 BAYESIAN AND PROFILE ANALYSIS 3.3 Transmuted log-logistic in the presence of censored life- time Let us consider the transmuted log-logistic model presented in Chapter 2. Further- more, consider the n lifetimes observed t1, t2, . . . , tn and the observed censor indicators c1, c2, . . . , cn, where ci = 1 if ti is exactly observed, or ci = 0 otherwise. For the right censoring, we are currently considering in this study, we have ci = 1 if ti ≤ t(r), and ci = 0 otherwise. Then, the likelihood function can be written as Ln(θ|t1, t2, . . . , tn) = k n∏ i=1 f(ti|θ)ci (1− F (ti|θ))(1−ci) , where K is a constant factor not depending on ti, ci or any unknown parameters, see for example Lawless (2011). By considering the transmuted log-logistic model, the likelihood function is written as Ln(µ, β, λ|t1, . . . , tn, c1, . . . , cn) = (βeµ) ∑n i=1 ci (3.5) × n∏ i=1 tci(β−1)i [ (1 + eµtβi )− λ(eµtβi − 1) (1 + eµtβi ) 3 ]ci . [ 1 + eµtβi (1− λ) (1 + eµtβi ) 2 ]1−ci . Hence, the log-likelihood function ln(·) = lnLn(·) becomes ln(µ, β, λ | t1, . . . , tn, c1, . . . , cn) = n∑ i=1 ci lnβ + µ n∑ i=1 ci + (β − 1) n∑ i=1 ci ln ti + n∑ i=1 ci ln [ (1 + eµtβi )− λ(eµtβi − 1) ] − 3 n∑ i=1 ci ln(1 + e µtβi ) + n∑ i=1 (1− ci) ln [ 1 + (1− λ)eµtβi ] − 2 n∑ i=1 (1− ci) ln(1 + eµtβi ). 3.3.1 Profile likelihood Estimating the parameters of a model using the likelihood method can be quite difficult in some special cases. In these cases, computational alternatives, such as the profile likeli- hood function, are required. They are used to calculate of classical estimates that allow us to estimate parameters by eliminating the numerical problems of the classical likelihood function; see for example Barndorfe-Nielsen and Cox (1994). Following Venzon and Moolgavkar (1988), suppose there are two sets of parameters, ψ and φ, where ψ represents the vector of factor loadings and φ represents all the other model parameters (in some cases, φ is a nuisance parameter or a vector of nuisance parameters). The profile method estimating parameters (ψ, φ) in two stages is given below. 39 BAYESIAN AND PROFILE ANALYSIS Let us suppose that X1, X2, . . . , Xn are iid random variables with density f(x;ψ, φ). As our objective is to estimate ψ and φ, the log-likelihood function is given by lP (ψ, φ|x1, . . . , xn) = n∑ i=1 log f(x;ψ, φ). (3.6) Then, to estimate ψ and φ the following can be used (ψˆP , φˆP ) = arg maxψ,φ lP (ψ, φ|x1, . . . , xn) what can be quite difficult and can lead to expression which are hard to maximize. Instead let us consider a different method where suppose for now, ψ is known, then we rewrite the likelihood as lP (ψ, φ|x1, . . . , xn) = lψ(φ|x1, . . . , xn) (to show that ψ is fixed and φ varies). To estimate φ, we maximize lψ(φ|x1, . . . , xn) with respect to φ, ie, φˆψ = arg max φ lψ(φ|x1, . . . , xn). (3.7) In reality ψ is unknown, hence for each ψ we have a new curve lψ(φ|x1, . . . , xn) over φ. Now, to estimate ψ, we evaluate the maximum lψ(φ|x1, . . . , xn), over φ, and choose the ψ, which is the maximum over all these curves, i.e., ψˆP = arg max ψ lψ(φˆψ|x1, . . . , xn) = arg max ψ lP (ψ, φˆψ|x1, . . . , xn), (3.8) and, logically, ψˆ and φψˆP are the maximum likelihood estimators (ψˆ, φˆ) = arg maxψ,φ lP (ψ, φ|x1, . . . , xn). For the transmuted log-logistic model, consider the n lifetimes observed (t1, t2, . . . , tn), with observed right censorship, i.e., δi = 1 if the lifetime is observed or 0 otherwise, i = 1 . . . , n. Notice that in the presence of a high number of censored lifetimes, it is quite difficult to estimate the model parameter using the pure likelihood. In addition, there is the presence of the nuisance parameter in the model, given by parameter λ, and a vector of parameters (µ, β). Thus, in our case, (µˆP , βˆP ) = arg max µ,β lµ,β(λˆµ,β|x1, . . . , xn) = arg max µ,β lP (µ, β, λˆµ,β|x1, . . . , xn), (3.9) where lP (µ, β, λˆµ,β|t) =  (β − 1)∑ni=1 ln ti − 3∑ni=1 ln(1 + eµtβi ) + k(lnβ + µ) + ∑n i=1 ln [ (1 + eµtβi )− λˆµ,β(eµtβi − 1) ] , if i : δi = 1, −2∑ni=1 ln(1 + eµtβi ) +∑ni=1 ln [1 + (1− λˆµ,β)eµtβi ] , if i : δi = 0, 40 BAYESIAN AND PROFILE ANALYSIS where k is the number of censored lifetimes. 3.3.2 Adjusted profile likelihood We shall now consider a different adjustment to the profile likelihood function. Suppose that the parameters that index the model are orthogonal, i.e., the elements of the score vector, ∂l/∂µ, ∂l/∂β and ∂l/∂λ , are uncorrelated. As proposed by Cox and Reid (1987), the adjustment can be applied to the profile likelihood function in this setting and it is an approximation to a conditional probability density function of the observations given, the maximum likelihood estimator of µ and β and can be written as LAP (µ, β, λ|x) = LP (µ, β, λˆµ,β|x) ∣∣∣Jλλ(µ, β, λˆµ,β)∣∣∣−1/2 (3.10) where Jλλ(µ, β, λˆµ,β) is the (λ, λ) element of the observed Fischer information J(µ, β, λ), see for example Barndorff-Nielsen and McCullagh (1993). In order to obtain the Jλλ(µ, β, λ) element, the Hessian matrix is presented in Appendix A. Then, the expected value is given by Jλλ(µ, β, λ) = ∫ +∞ 0 ( 1− eµxβ 1 + eµxβ − λeµxβ + λ )2 f(x|µ, β, λ)dx ' pie −µ/β 3β2 (β + 1) [ 2(β + 1)2 + 1 ] csc ( pi(β + 1) β ) . (3.11) By considering n lifetimes observed (t1, t2, . . . , tn), with observed right censored, i.e., δi = 1 if the lifetime is observed or 0 otherwise, i = 1 . . . , n, the adjusted profile log- likelihood function can be written as lAP (µ, β, λ|t) =  (β − 1)∑ni=1 ln ti − 3∑ni=1 ln(1 + eµtβi )+ +k(lnβ + µ) + ∑n i=1 ln [ (1 + eµtβi )− λˆµ,β(eµtβi − 1) ] + ζ, if i : δi = 1, −2∑ni=1 ln(1 + eµtβi ) +∑ni=1 ln [1 + (1− λˆµ,β)eµtβi ]+ ζ, if i : δi = 0, (3.12) where ζ = β e−µ/2β [ − 3 pi [2(β + 1)2 + 1] sin ( pi(β + 1) β )]1/2 . 41 BAYESIAN AND PROFILE ANALYSIS 3.3.3 Modified profile likelihood The modified profile likelihood was proposed by Barndorff-Nielsen (1993) as an im- provement to the profile likelihood and it is defined as LMP = LP (θ,x) ∣∣∣Jφφ(θ, φˆ(θ))∣∣∣−1/2 ∣∣∣∣∣ ∂φˆ∂φˆ(θ) ∣∣∣∣∣ . (3.13) The main difficulty in computing the modified profile likelihood function lies in obtain- ing |∂φˆ/∂φˆ(θ)| and several approximations were proposed in order to simplify its evaluation. Severini (1998) proposed an approximation based on empirical covariates and LMP can be rewritten as LMP = LP (θ,x) ∣∣∣Jφφ(θ, φˆ(θ))∣∣∣1/2∣∣∣I(θ, φˆ(θ); θˆ, φˆ)∣∣∣ , (3.14) and replacing I(θ, φˆ(θ); θˆ, φˆ) by I˘(θ, φˆ(θ); θˆ, φˆ) = n∑ j=1 ` (j) θ (θ, φˆ(θ))` (j) θ (θˆ, φˆ) T (3.15) where ` (j) θ (·) is the score function for the jth observation. This approximation is particularly useful when the computation of expected values of products of log-likelihood derivatives is cumbersome. Thus, for the transmuted log-logistic model, I˘(µ, β, λˆ(µ, β); µˆ, βˆ, λˆ) = n∑ j=1 ` (j) λ (µ, β, λˆ(µ, β))` (j) λ (µˆ, βˆ, λˆ) T and ` (j) λ (µ, β, λ) = (1− eµxβi )( 1 + eµxβi − λ(eµxβi − 1) ) . 3.3.4 Residual analysis In order to study the presence of atypical observations, or departures from the error assumptions, we consider two kinds of residuals: martingale-type and deviance. (see for example McCullagh and Nelder (1989), Barlow and Prentice (1988) and Therneau et al. (1990)). The first one, martingale-type residual, was introduced by Therneau et al. (1990) and was firstly used in a counting process, which is skewed and has a maximum value at +1 and a minimum value at −∞. By considering the transmuted log-logistic model, the martingale-type residual can be written as 42 BAYESIAN AND PROFILE ANALYSIS rMi =  1 + log [ 1 + eµtβi (1− λ) ] − 2 log(1 + eµtβi ), if i : δi = 1, log [ 1 + eµtβi (1− λ) ] − 2 log(1 + eµtβi ), if i : δi = 0, (3.16) where i = 1, . . . , n. In addition, the deviance residual can be used, widely used in GLMs (generalized linear models). This was proposed by the same authors (Therneau et al. (1990)) and it is a trans- formation of the martingale residual to attenuate the skewness. In our case, the deviance residuals are given by rDi =  sign(rMi) [ −2 [ 1 + log(1 + eµtβi (1− λ))− 2 log(1 + eµtβi ) + log ( − log(1 + eµtβi (1− λ)) ) + 2 log(1 + eµtβi ) ]1/2] , if i : δi = 1, sign(rMi) [ −2 [ 1 + log(1 + eµtβi (1− λ))− 2 log(1 + eµtβi ) ]] , if i : δi = 0, (3.17) i = 1, . . . , n. 3.3.5 Simulation In this section, we describe a simulation study designed to assess the frequentist prop- erties of the model in the presence of right censored observations in order to validate the results that will be shown in Section 3.3.6. We must observe that transmuted models have never been built considering the presence of any censor type. Consider the sample size 148 of patients treated with the drug Linezolid. A bootstrap study was made considering 6 different sizes of re-samples: 50, 80, 100, 200, 300 and 500 (index 1 to 6 respectively). For each sample size, we obtained 1000 re-sampling. The transmuted log-logistic model was fitted considering three estimation methods, presented before: Profile (P), Adjusted Profile (AP) and Modified Profile (MP). The results of the bootstrap simulations were summarized and are presented in Figures 3.7, 3.8 and 3.9. In the first figure, we present the mean value of estimates of parameters µ, β and λ. Notice that for µ and β parameters, the three profile methods behave similarly. However, in Figure 3.7 the adjusted profile method is better than the other methods for samples less than 200. As for the mean square error (MSE, in Figures 3.8) it is clear that the adjusted method presents smaller values even though the MSEs are close. Furthermore, with respect to the coverage probability, we expected it to be around 95% in all cases (we can observe that 43 BAYESIAN AND PROFILE ANALYSIS 1 2 3 4 5 6 − 5. 5 − 5. 0 − 4. 5 − 4. 0 µ (a) 1 2 3 4 5 6 1. 3 1. 4 1. 5 1. 6 1. 7 1. 8 1. 9 2. 0 β (b) 1 2 3 4 5 6 − 1. 0 − 0. 9 − 0. 8 − 0. 7 − 0. 6 λ (c) Figure 3.7: Mean for the parameters estimated: (a) µ, (b) β and (c) λ, by using × P, M AP and + MP methods. the points are close to the dashed line in Figures 3.9 a, b and c). 3.3.6 Application to real dataset Linezolid is a synthetic antibiotic developed by a team at the Pharmacia and Upjohn Company, used to treat of serious infections, Brickner (1996). Discovered in the 1990s and first approved for use in 2000, linezolid was the first commercially available 1, 3-oxazolidinone antibiotic. The main indication of this substance is the treatment of severe infections caused by Gram-positive bacteria that are resistant to other antibiotics. In both the popular press and the scientific literature, linezolid has been called a reserve antibiotic - one that should be used sparingly so that it will remain effective as a drug of last resort against potentially intractable infections, see Wilson et al. (2006). In Brazil, linezolid has been used since 2007 and researchers are concerned about their indiscriminate use. It is known that its use, for short periods of time, is safe but the effects 44 BAYESIAN AND PROFILE ANALYSIS 1 2 3 4 5 6 0. 0 0. 5 1. 0 1. 5 µ (a) 1 2 3 4 5 6 0. 0 0. 1 0. 2 0. 3 0. 4 0. 5 β (b) 1 2 3 4 5 6 0. 0 0. 1 0. 2 0. 3 0. 4 0. 5 λ (c) Figure 3.8: Mean Square Error (MSE) of the parameters estimated: (a) µ, (b) β and (c) λ, by using × P, M AP and + MP methods. of its prolonged use are not known. Given this situation, we conducted a study related to the time of using this drug in patients at the ICU (Intensive Care Unit) at the University Hospital in Maringá city, from 2008 and 2012. Total of 148 patients were treated with the drug, which is only used in extreme cases of bacterial resistance against other drugs, such as Vancomycin (drug used before linezolid in those cases). An initial analysis was made, in which we noticed that the mean time was 28.4 days and the median time was 22.5 days. Moreover, Figure 3.10-a show the distribution of the time and, it can be seen that the first and third observed quantiles of time of treatment was 13.0 and 36.3 days, respectively. In order to know the behavior of the failure rate, the TTT plot was plotted in Figure 3.10-b. It is possible to see that the hazard is unimodal, since the TTT plot is initially concave and then convex. After that, the transmuted log-logistic model was fitted, using the pure likelihood and the profile method. Both results are shown in Table 3.7 below. It can be clearly seen 45 BAYESIAN AND PROFILE ANALYSIS 1 2 3 4 5 6 0. 5 0. 6 0. 7 0. 8 0. 9 1. 0 µ (a) 1 2 3 4 5 6 0. 5 0. 6 0. 7 0. 8 0. 9 1. 0 β (b) 1 2 3 4 5 6 0. 5 0. 6 0. 7 0. 8 0. 9 1. 0 λ (c) Figure 3.9: Coverage probability for the parameters estimated: (a) µ, (b) β and (c) λ, by using × P, M AP and + MP methods. l l l l l l l l l 0 50 10 0 15 0 Te m po (a) 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 r/n f (b) Figure 3.10: (a) Boxplot of the observed times and (b) TTT plot. 46 BAYESIAN AND PROFILE ANALYSIS that the estimated parameters are very close in both likelihood methods. However, the profile confidence intervals are shorter when using the pure likelihood. Furthermore, if we consider the Wald confidence intervals for the profile method, the intervals are even smaller, as shown in Table 3.8. Table 3.7: Estmates for the parameters of transmuted log-logistic model estimated by using profile methods, considering right censored lifetimes. Method Parameter Estimate Error Confidence Intervals 95% Lower Upper α −4.4446 0.8072 −5.0528 −3.9723 Profile β 1.4804 0.2089 1.3526 1.6301 λ −0.9483 0.1685 −1.0000 −0.7928 α −4.5909 0.6548 −5.0649 −4.1828 Adjusted β 1.5415 0.1702 1.4319 1.6614 λ −0.9473 0.1382 −1.0000 −0.8296 α −4.3460 0.8980 −5.0183 −3.8702 Modified β 1.4480 0.2344 1.3142 1.6135 λ −0.9577 0.1796 −1.0000 −0.7888 Table 3.8: Wald confidence limits by considering 95% of confidence. Method Lower Upper −4.9890 −3.9001 Profile 1.3394 1.6213 −1.0000 −0.8346 −5.0649 −4.1828 Adjusted 1.4319 1.6614 −1.0000 −0.8296 −4.9517 −3.7403 Modified 1.2899 1.6060 −1.0000 −0.8366 Figure 3.11-a and b, shows the estimated survival curves versus empirical (by Kaplan- Meier). It can be clearly observed how close the curves are. Moreover the hazard curve can also be seen. Notice that the most probable time for patient discharge is 25.44 days, with 95% confidence interval given by (15.47; 48.69) days. Furthermore, Figure 3.12 shows the graphics of the Profile, Adjusted and Modified relative log-likelihoods (divided by the maximum absolute value of the log-likelihood). It can be observed that the Profile and Adjusted Profile showed a lower decrease (after the maximum value) than the Modified Profile log-likelihoods. 47 BAYESIAN AND PROFILE ANALYSIS 0 50 100 150 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Time Su rv iva l (a) 0 50 100 150 0. 00 5 0. 01 0 0. 01 5 0. 02 0 0. 02 5 Time H az ar d (b) Figure 3.11: (a) Survival curves: estimated versus empirical; (b) Hazard estimate curve, with the Tˆmax. 0 50 100 150 − 1. 0 − 0. 9 − 0. 8 − 0. 7 − 0. 6 − 0. 5 time P A M Figure 3.12: (P) Profile, (A) Adjusted Profile and (M) Modified Profile relative log-likelihood for the Linezolid data. Global and local influence In this section, we will make an analysis of the global and local influence for the data set given, using the transmuted log-logistical model, as presented in Section 3.2.1. 48 BAYESIAN AND PROFILE ANALYSIS In order to analyze the local influence, here we consider the response variable pertur- bation, i.e., we will consider that each ti is perturbed as tim = ti + miSt , where St is a scale factor that may be the estimated standard deviation of T and mi ∈ R. Then, the perturbed log-likelihood function becomes ln(µ, β, λ | t, c,m) = n∑ i=1 ci lnβ + µ n∑ i=1 ci + (β − 1) n∑ i=1 ci ln tim + n∑ i=1 ci ln [ (1 + eµtβim)− λ(eµtβim − 1) ] − 3 n∑ i=1 ci ln(1 + e µtβim) + n∑ i=1 (1− ci) ln [ 1 + (1− λ)eµtβim ] − 2 n∑ i=1 (1− ci) ln(1 + eµtβim). Figure 3.13-a and b, show, respectively, the Cook's generalized and likelihood distances and it can be observed that the pertubation causes some disproportionate effects. (a) (b) Figure 3.13: (a) Cooks's distance and (b) Likelihood distance after response pertubation. After analyzing Figure 3.13-a and b, we can see the distinction of some observations in relation to others. Furthermore, we made a residual analysis using the Martingale-type and deviance presented in Section 3.3.4. Figures 3.14 show the results of these analyses. To reveal the impact of the detected influential observations, we use three measures defined by Lee et al. (2006): TRC, MRC and LD. In our case, we find TRC = 22.78, MRC= 0.0873 and LD(I) = −32.3, for I = {50, 53, 130, 144}. Hence, the results are more sensitive for the influential observation group. After the influence and residual analysis, the possible influential observations were identified: case 50, 53, 130 and 144. Then, we removed each of these points and fitted the transmuted log-logistic model again for each case. In order to reveal the impact of them, the relative changes were measured RCθj and all these results are shown in Table 3.9. As a complementary analysis, the possible influential observations were removed. Aided by the analysis of local influence and residual analysis, observations 50, 53, 130 and 144 were 49 BAYESIAN AND PROFILE ANALYSIS (a) (b) Figure 3.14: (a) Martingale residuals; (b) Deviance residuals. Table 3.9: RC (in %) and the corresponding TRC, MRC and LD(I). Case removed Parameter RC TRC MRC LD(I) {50} µ 0.682 7.336 5.452 4.300 β 1.202 λ 5.452 {53} µ 1.442 3.831 4.144 14.300 β 1.513 λ 0.875 {130} µ 4.144 10.704 4.144 13.400 β 3.702 λ 2.858 {144} µ 2.187 5.707 2.187 12.600 β 2.128 λ 1.392 {50; 53; 130; 144} µ 8.734 22.777 8.734 33.300 β 8.127 λ 5.916 then removed. The results can be found in Table 3.10. The graphics of the Profile, Adjusted and Modified relative log-likelihoods (divided by the maximum absolute value of the log-likelihood) were plotted again and they are presented in Figure 3.15. It can be observed that the Profile and Adjusted Profile had a lower decrease (after the maximum value) than the Modified Profile log-likelihoods, as seen in Figure 3.12 showing some advantage, which are not significative, for this method. 50 BAYESIAN AND PROFILE ANALYSIS Table 3.10: Estimates for the parameters of transmuted log-logistic model estimated by using profile methods after removing influential observations. Method Parameter Estimate Error Confidence Intervals 95% Lower Upper α −4.8327 0.9350 −5.4634 −4.2021 Profile β 1.6006 0.2374 1.4405 1.7607 λ −0.8922 0.2139 −1.0000 −0.7479 α −4.8334 0.7183 −5.3179 −4.3489 Adjusted β 1.6248 0.1872 1.4985 1.7511 λ −0.9192 0.1542 −1.0000 −0.8152 α −4.8259 1.0426 −5.5292 −4.1227 Modified β 1.5922 0.2624 1.4153 1.7692 λ −0.8848 0.2419 −1.0000 −0.7216 0 50 100 150 − 1.0 − 0.9 − 0.8 − 0.7 − 0.6 − 0.5 time P A M Figure 3.15: (P) Profile, (A) Adjusted Profile and (M) Modified Profile relative log-likelihood for the Linezolid data. Chapter 4 Cubic ranking transmuted model This Chapter introduces a new order of transmuted distribution, the cubic rank trans- mutation map distribution. This new order increases the flexibility of transmuted distri- butions providing the accommodation of more complex data; for instance, data in the presence of bimodal hazard rates. In order to exemplify the applicability of the cubic rank transmutation map, we used two of the best known distributions maps in the reliability field: the Weibull and the log-logistic. Several mathematical properties of these new distri- butions, namely the transmuted Weibull distribution of order 2 and transmuted log-logistic distribution of order 2, are derived (either cubic rank transmuted Weibull or cubic rank transmuted log-logistic). Inference is based on maximum likelihood and applications were made using real datasets. Analysis of diagnostic and a bootstrap simulation study are part of this work. 4.1 The cubic ranking transmutation map In this section we propose a cubic ranking transmutation map, following the transmu- tation described in Section 2.1, but now considering order 2. Theorem 4.1.1 Let X1, X2 and X3 be independent and identically random variables with distribution G(x). Then, the ranking cubic transmutation map is given by F (x) = λ1G(x) + (λ2 − λ1)G(x) + (1− λ2)G3(x), (4.1) with λ1 ∈ [0, 1] and λ2 ∈ [−1, 1]. Proof: Let X1, X2 and X3 be independent and identically random variables distributed with distribution G(x). Now, consider the following order X1:3 = min(X1, X2, X3), X2:3 and X3:3 = max(X1, X2, X3), 51 52 CUBIC RANKING TRANSMUTED MODEL and let  Y d = X1:3, with prob pi1, Y d = X2:3, with prob pi2, Y d = X3:3, with prob pi3, where ∑3 i=1 pii = 1⇒ pi3 = 1− pi1 − pi2. Evidently, FY (x) is given by FY (x) = pi1Pr(min(X1, X2, X3) ≤ x) + pi2Pr(X2:3 ≤ x) + pi3Pr(max(X1, X2, X3) ≤ x) = 3pi1G(x) + 3(pi2 − pi1)G2(x) + (1− pi2)G3(x), and if 3pi1 = λ1 and 3pi2 = λ2, the cubic ranking transmuted distribution (or transmuted of order 2) becomes the one given in (4.1).  Definition 4.1.1 The density function f(x) of the cubic ranking transmutation distribu- tion is given by f(x) = g(x) [ λ1 + 2(λ2 − λ1)G(x) + 3(1− λ2)G2(x) ] . (4.2)  In what follows, we focus on two particular cases based on the well-known Weibull and log-logistic lifetime distributions. 4.2 Cubic rank transmuted weibull distribution The Weibull distribution, originally proposed by Weibull (1951), is being used effec- tively in many applied problems; see for example, Johnson et al. (1996), Nichols and Pad- gett (2006), Pham and Lai (2007), and Aryal and Tsokos (2011). Let X be a positive random variable with Weibull distribution, then, the cumulative function is given by G(x) = 1− exp [ − ( x β )µ] , x > 0, (4.3) where µ > 0 and β > 0 are shape and scale parameters, respectively. It is important to note that for µ < 1 the failure rate decreases over time (which happens, for instance, in data where there is significant infant mortality), µ = 1 indicates that the failure rate is constant over time (i.e., it is an exponential distribution as a special case), and µ > 1 possesses the failure rate increasing over time, (Lawless, 2011). Proposed by Aryal and Tsokos (2011), the transmuted Weibull distribution was for- mulated as a generalization of the Weibull distribution. The transmuted quadratic rank Weibull cumulative distribution function has the form F (x) = [ 1− exp ( − ( x β )µ)][ 1 + λ exp ( − ( x β )µ)] , 53 CUBIC RANKING TRANSMUTED MODEL where λ ∈ [−1, 1]. The cubic rank transmuted Weibull distribution introduced here has the following form. Proposition 4.2.1 Let X have a Weibull distribution with parameters µ and β. Then, the density function of the cubic rank transmuted Weibull distribution is given by f(x) = µ βµ xµ−1e−(x/β) µ [ (3− λ1 − λ2) + 2(λ1 + 2λ2 − 3)e−(x/β)µ + 3(1− λ2)e−2(x/β)µ ] , (4.4) for x > 0. Proof: With X having a Weibull distribution with parameters µ and β, we obtain from equation (4.1) that F (x) = λ1G(x) + (λ2 − λ1)G2(x) + (1− λ2)G3(x) = 1 + (λ1 + λ2 − 3)e−(x/β)µ + (3− λ1 − 2λ2)e−2(x/β)µ + (λ2 − 1)e−3(x/β)µ , for x > 0. Upon differentiating the above expression with respect to x, we obtain the density function presented in (4.1).  Note that the cubic rank transmuted Weibull is an extended distribution to analyze more complex data, generalizing theWeibull distribution (λ1 = λ2 = 1) and the transmuted Weibull distribution (λ2 = 1). Some of the possible shapes that the cubic rank transmuted Weibull distribution can take on are illustrated in the upper panels of Figure 4.1 for some selected values of the parameters. The λ's are responsible for introducing extra skewness into the Weibull distribution. This is in full agreement with Shaw and Buckley (2009) who pointed out that the introduction of extra skewness into a distribution is a direct effect of transmutation maps of order 1, or transmuted quadratic rank ones. Now, let T be a random variable representing the lifetime of an unit. Then, the survival and hazard functions of T are given, respectively, by R(t) = (3− λ1 − λ2)ζ(t) + (λ1 + 2λ2 − 3)(ζ(t))2 + (1− λ2)(ζ(t))3, t > 0, (4.5) and h(t) = µ βµ tµ−1 [ (3− λ1 − λ2) + 2(λ1 + 2λ2 − 3)ζ(t) + 3(1− λ2)(ζ(t))2 (3− λ1 − λ2) + (λ1 + 2λ2 − 3)ζ(t) + (1− λ2)(ζ(t))2 ] , t > 0, (4.6) where ζ(t) = e−(t/β)µ . Several plots of the above survival and hazard functions are presented in the lower panels of Figure 4.1. Moments are essential for inferential proposes. Moreover, they are useful to study some important features and characteristics of a distribution such as tendency, dispersion, skew- ness and kurtosis. We derive the rth moment for the cubic rank transmuted Weibull dis- 54 CUBIC RANKING TRANSMUTED MODEL 0 100 200 300 400 0 1 2 3 4 Time D en si ty F un ct io n c1 c2 c3 c4 c5 c6 0 100 200 300 400 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Time Cu m u la tiv e F un ct io n c1 c2 c3 c4 c5 c6 0 100 200 300 400 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Time Cu m u la tiv e F un ct io n c1 c2 c3 c4 c5 c6 0 50 100 150 0. 0 0. 5 1. 0 1. 5 2. 0 Time H az ar d Fu nc tio n c1 c2 c3 c4 c5 c6 0 100 200 300 400 0 10 20 30 Time Cu m u la tiv e H az ar d Fu nc tio n c1 c2 c3 c4 c5 c6 Figure 4.1: Upper panels: Density and cumulative distribution functions; Lower panels: Survival, hazard and cumulative hazard functions, all of which are plotted for different values of model parameters. tribution as follows. The rth moment is given by E (Xr) = ∫ +∞ 0 xr µ βµ xµ−1e−(x/β) µ [ (3− λ1 − λ2) + 2(λ1 + 2λ2 − 3)e−(x/β)µ +3(1− λ2)e−2(x/β)µ ] dx. With the use of gamma function, we can obtain an expression for the above integral as E (Xr) = µ β1−r Γ ( r − 1 µ + 2 )[ (3− λ1 − λ2) + (λ1 + 2λ2 − 3) 2 r−1 µ +1 + (1− λ2) 3 r−1 µ +1 ] , where Γ(·) denotes a complete gamma function. In particular, the mean and variance are given, respectively, by E (X) = µ 6 (11− 3λ1 − 2λ2) (4.7) and V (X) = Γ ( 1 µ + 2 ) βµ [ (3− λ1 − λ2) + (λ1 + 2λ2 − 3) 2 1 µ +1 + (1− λ2) 3 1 µ +1 ] −µ 2 36 (11− 3λ1 − 2λ2)2. 55 CUBIC RANKING TRANSMUTED MODEL 4.3 Cubic rank transmuted log-logistic distribution Let X be a nonnegative random variable representing the lifetime of an individual in a population log-logistically distributed as presented in Section 2.1. Then, the cubic rank transmuted log-logistic distribution can be defined as follows. Proposition 4.3.1 Let X be random variables having a log-logistic distribution with pa- rameters β and µ. Then, the density function of the cubic rank transmuted log-logistic distribution is given by f(x) = µxµ−1eβ (1 + eβxµ)4 [ λ1(1− e2βx2µ) + λ2eβxµ(2− eβxµ) + 3e2βx2µ ] , x > 0.(4.8) Proof: Let X be a log-logistic random variable with parameters β and µ. Then, from equation (4.1), we obtain F (x) = λ1G(x) + (λ2 − λ1)G2(x) + (1− λ2)G3(x) = xβeµ (1 + eµxβ) [ λ1 + (λ2 − λ1) x βeµ (1 + eµxβ) + (1− λ2) x 2βe2µ (1 + eµxβ)2 ] , x > 0. Upon differentiating the above expression with respect to x, we obtain the density function as presented in (4.8).  Note that the cubic rank transmuted log-logistic model is also an extended distribution that is quite useful for analysing more complex data, as presented for the Weibull model in Section 4.2. The log-logistic distribution is clearly a special case when λ1 = λ2 = 1. Some of the possible shapes that the cubic rank transmuted log-logistic distribution can take on are illustrated in the upper panels of Figure 4.2 for some selected values of the parameters. The λ's are responsible for introducing extra skewness into the log-logistic distribution. Furthermore, the reliability and hazard functions of this distribution are given by R(t) = 1− κ(t) [λ1 + (λ2 − λ1)κ(t) + (1− λ2)κ2(t)] , t > 0, (4.9) and h(t) = µtµ−1eβ (1 + eβtµ)3 [ λ1(1− e2βt2µ) + λ2eβtµ(2− eβtµ) + 3e2βt2µ 1 + tµeβ − tµeβ [λ1 + (λ2 − λ1)κ(t) + (1− λ2)κ2(t)] ] , (4.10) t > 0, respectively, where κ(t) = e βtµ (1+eβtµ) . Several plots of survival and hazard functions are presented in the lower panels of Figure 4.2. The rth moment of the cubic rank transmuted log-logistic distribution is given by E (Xr) = ∫ +∞ 0 µxr+µ−1eβ (1 + eβxµ)4 [ λ1(1− e2βx2µ) + λ2eβxµ(2− eβxµ) + 3e2βx2µ ] dx. 56 CUBIC RANKING TRANSMUTED MODEL 0.0 0.5 1.0 1.5 2.0 0 2 4 6 8 Time D en si ty F un ct io n c1 c2 c3 c4 c5 c6 0 100 200 300 400 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Time Cu m u la tiv e F un ct io n c1 c2 c3 c4 c5 c6 0 100 200 300 400 0. 2 0. 4 0. 6 0. 8 1. 0 Time Su rv iva l F un ct io n c1 c2 c3 c4 c5 c6 0 50 100 150 200 250 0 1 2 3 4 5 Time H az ar d Fu nc tio n c1 c2 c3 c4 c5 c6 0 50 100 150 200 250 0 1 2 3 4 5 Time Cu m u la tiv e H az ar d Fu nc tio n c1 c2 c3 c4 c5 c6 Figure 4.2: Upper panels: Density and cumulative distribution functions; Lower panels: Survival, cumulative hazard and hazard functions, all of which are plotted for different values of model parameters. Upon using beta function, we obtain an expression for the rth moment as E (Xr) = e−βr/µ 6 B ( 1− r µ , r µ + 1 ){ r2 µ2 (3− λ2 + 2λ1) + 3 r µ (3 + λ2) + 2(2λ1 + λ2 + 3) } , where B(·, ·) denotes the complete beta function. In particular, the mean and variance are given by E (X) = e−β/µ 6µ2 B ( µ− 1 µ , µ+ 1 µ ){ 2λ1(1 + 2µ 2) + λ2(µ 2 + 3µ− 1) + 7µ2 + 9µ} (4.11) and V (X) = e−2β/µ 6µ2 [ B ( µ− 2 µ , µ+ 2 µ ){ 4λ1(2 + µ 2) + 2λ2(µ 2 + 3µ− 2) +6(µ2 + 3µ+ 2) }−B2(µ− 1 µ , µ+ 1 µ ){ 2λ1(1 + 2µ 2) + λ2(µ 2 + 3µ− 1) +7µ2 + 9µ }2] . (4.12) 57 CUBIC RANKING TRANSMUTED MODEL 4.4 Parameter estimation Maximum likelihood approach can be used for the estimation of model parameters. Let X1, . . . , Xn be a random sample of size n so that the likelihood function is given by L = L(θ|X1, X2, . . . , Xn) = n∏ i=1 f(Xi|θ), where θ is the parameter vector. For the cubic rank transmuted Weibull distribution, the log-likelihood function becomes lW = n log(µ)− nµ log(β) + (µ− 1) n∑ i=1 log xi − n∑ i=1 ( xi β )µ (4.13) + n∑ i=1 log [ (3− λ1 − λ2) + 2(λ1 + 2λ2 − 3)e−(xi/β)µ + 3(1− λ2)e−2(xi/β)µ ] . Therefore, the maximum likelihood estimates (MLEs) of µ, β, λ1 and λ2, which maxi- mize (4.13), must satisfy the following normal equations: ∂lW ∂µ = n µ − n∑ i=1 log xi − n∑ i=1 ( xi β )µ log ( xi β ) +2 n∑ i=1 ( xi β )µ log ( xi β ) e−(xi/β) µ [ λ1 + 2λ2 − 3 ϑi − 3(1− λ2)e−(xi/β)µ ] = 0, ∂lW ∂β = 2µ β n∑ i=1 ( xi β )µ e−(xi/β) µ [ λ1 + 2λ2 − 3 ϑi + 3(1− λ2)e−(xi/β)µ ] −nµ β + µ β n∑ i=1 ( xi β )µ = 0, ∂lW ∂λ1 = n∑ i=1 2e−(xi/β)µ − 1 ϑi = 0, ∂lW ∂λ2 = n∑ i=1 4e−(xi/β)µ − 1 ϑi − 3e−2(xi/β)µ = 0, where ϑi = (3 − λ1 − λ2) + 2(λ1 − 2λ2 − 3)e−(xi/β)µ . The MLE θˆ = (µˆ, βˆ, λˆ1, λˆ2)′ is obtained by solving the above nonlinear system of equations. It is usually more convenient to use nonlinear optimization algorithms such as quasi-Newton or Newton-Raphson to numerically maximize the log-likelihood function in (4.13). Similarly, for the cubic rank transmuted log-logistic distribution, the log-likelihood 58 CUBIC RANKING TRANSMUTED MODEL function is given by lLL = nβ + n logµ− 4 n∑ i=1 log(1 + eβxµi ) + (µ− 1) n∑ i=1 log xi + + n∑ i=1 log(λ1 − λ1e2βx2µi + 2λ2eβxµi − λ2e2βx2µi + 3e2βx2µi ). (4.14) Therefore, the MLE of µ, β, λ1 and λ2 which maximize (4.14) must satisfy the following normal equations: ∂lLL ∂µ = n µ − 4 n∑ i=1 ηi log xi 1 + ηi + n∑ i=1 log xi + n∑ i=1 2ηi log xi [ 3ηi + λ2(1− ηi)− λ1ηi η2i (3− λ1 − λ2) + 2ηi + λ1 ] = 0, ∂lLL ∂β = n− 4 n∑ i=1 ηi 1 + ηi + 2 n∑ i=1 ηi [ ηi(3− λ1 − λ2) + λ2 η2i (3− λ1 − λ2) + 2ηi + λ1 ] = 0, ∂lLL ∂λ1 = n∑ i=1 1− η2i η2i (3− λ1 − λ2) + 2ηi + λ1 = 0, ∂lLL ∂λ2 = n∑ i=1 ηi(2− ηi) η2i (3− λ1 − λ2) + 2ηi + λ1 = 0, where ηi = e βxµi . The MLE θˆ = (µˆ, βˆ, λˆ1, λˆ2) ′ is obtained by solving the above nonlinear system of equations. In order to compute the asymptotic confidence intervals, we use the usual large-sample approximation for the distribution of the maximum likelihood estimator of θ which is four- variate normal (Migon et al., 2014). Thus, as n → ∞, the asymptotic distribution of the MLE, for both cubic rank transmuted Weibull and log-logistic distributions is given by, µˆ βˆ λˆ1 λˆ2  ∼ N   µˆ βˆ λˆ1 λˆ2  ,  Vˆ11, Vˆ12, Vˆ13, Vˆ14 Vˆ21, Vˆ22, Vˆ23, Vˆ24 Vˆ31, Vˆ32, Vˆ33, Vˆ34 Vˆ41, Vˆ42, Vˆ43, Vˆ44   , (4.15) where Vˆij = Vij |θ=θˆ. The asymptotic variance-covariance matrix V is determined by the inverse of Hessian matrix; see A for the corresponding details for the cubic rank transmuted Weibull and log-logistic distributions. Then, approximate 100(1−α)% two-sided confidence intervals for µ, β, λ1 and λ2 are, respectively, given by µˆ± zα/2 √ Vˆ11, βˆ ± zα/2 √ Vˆ22, λˆ1± zα/2 √ Vˆ33 and λˆ2± zα/2 √ Vˆ44, where zα is the upper α− th percentile of the standard normal distribution. 59 CUBIC RANKING TRANSMUTED MODEL Fits of different distributions can be compared by penalizing over-fitting by using the Akaike information criterion (Akaike, 1973), which minimizes the Kullback-Leibler diver- gence between the true distribution and the estimate from a candidate distribution. It is given by AIC = −2l(θˆ)+2size(θ), where l(θˆ) denotes the log-likelihood function evaluated at the maximum value and size(θ) is the number of model parameters. The distribution with the lowest value of this criterion (among all considered distributions) is regarded as the preferred distribution for describing a given dataset. 4.5 Numerical studies In this section, we present the results of a simulation study designed to assess the properties of the proposed estimation procedure. Moreover, we illustrate the usefulness of the proposed models with two real data sets. An analysis of global and local influence is also carried out. 4.5.1 Simulation study A simulation study was performed by considering samples of size 50, 80, 100, 150, 300 and 500 from the cubic rank transmuted Weibull and cubic rank transmuted log-logistic distributions. A total of 1000 random samples were generated for each set up with the parameters fixed as µ = 2.0, β = 2.5, λ1 = 0.2, and λ2 = 1.0 for the cubic rank transmuted Weibull distribution, and µ = 3.0, β = −16.0, λ1 = 0.2, and λ2 = −0.9 for the cubic rank transmuted log-logistic distribution. Tables 4.1 and 4.2 present the means of the estimates and the mean square errors (MSE) as well as the coverage probabilities of 95% two-sided confidence intervals for the model parameters. The MSE decreases with increasing sample size, showing the consistency of the estimators. Moreover, the coverage probabilities become closer to the nominal level as the sample size increases. 4.5.2 Application: Carbon fibers data In this section, we provide an application of the cubic rank transmuted Weibull distri- bution. For this purpose, we consider a data set from a study on breaking stress of carbon fibers (in Gba), taken from Nichols and Padgett (2006). The lifetimes are times until the break of the fibers. TTT plot of the lifetimes is presented in the upper left panel of Figure 4.3, indicating a possible increasing hazard function. Also, for the sake of comparison, we have considered six alternative distributions: the transmuted Weibull distributions of orders 2 and 1, and the prominent Weibull, exponen- tial, gamma and log-normal distributions. 60 CUBIC RANKING TRANSMUTED MODEL Table 4.1: Means of estimates of all parameters and mean square errors (MSE), for different sample sizes. Distribution Sample Estimate MSE Size µ β λ1 λ2 µ β λ1 λ2 50 2.226 2.338 0.317 0.492 0.895 0.413 0.821 1.787 Trans 80 2.2 2.35 0.308 0.569 0.714 0.416 0.667 1.376 Weibull 100 2.179 2.346 0.289 0.597 0.63 0.397 0.563 1.258 Order2 150 2.176 2.363 0.29 0.652 0.53 0.374 0.457 1.028 300 2.15 2.364 0.245 0.746 0.402 0.336 0.296 0.771 500 2.142 2.376 0.228 0.807 0.314 0.281 0.23 0.568 50 2.893 −16.138 0.293 −0.822 0.304 4.559 0.120 0.257 Trans 80 2.889 −16.115 0.299 −0.823 0.253 3.507 0.107 0.244 Log-Logistic 100 2.890 −16.122 0.313 −0.838 0.166 1.929 0.080 0.187 Order2 150 2.880 −16.061 0.317 −0.845 0.110 1.083 0.059 0.138 300 2.876 −16.036 0.318 −0.848 0.096 0.915 0.053 0.124 500 2.876 −16.030 0.321 −0.853 0.040 0.309 0.025 0.059 Table 4.2: Coverage probabilities of the confidence intervals for different sample sizes. Distribution Sample Coverage Probability Size µ β λ1 λ2 50 0.577 0.671 0.551 0.732 Transmuted 80 0.662 0.72 0.642 0.781 Cubic Rank 100 0.702 0.734 0.690 0.787 Weibull 150 0.783 0.808 0.777 0.815 300 0.894 0.897 0.903 0.898 500 0.959 0.959 0.956 0.947 50 0.830 0.824 0.864 0.521 Transmuted 80 0.852 0.855 0.887 0.624 Cubic Rank 100 0.879 0.882 0.906 0.749 Log-Logistic 150 0.893 0.898 0.902 0.835 300 0.893 0.894 0.897 0.857 500 0.901 0.904 0.940 0.950 In Table 4.3, we have presented, for all these distributions, seven different compari- son measures used as selection criteria: −2× log-likelihood (-2log), Akaike's information criterion (AIC), corrected Akaike's information criterion (AICC), Schwarz's Bayesian in- formation criterion (BIC), Kolmogorov-Smirnov statistic (KS), Anderson-Darling statistic (A) and Cramér-von-Mises statistic (W). The calculated values of theses statistics (the smaller the better) all reveal that the cubic rank transmuted Weibull distribution is the most appropriate model according to four different criteria and the Weibull distribution is the most appropriate one according to AIC, AICC and BIC (it must be noted that, in these criteria, the number of estimated parameters has a huge impact). Figure 4.4 shows P-P plots that represent the empirical cumulative versus estimated cumulative functions of the rank transmuted Weibull distributions of orders 2 and 1 and the Weibull distribution, respectively. In order to examine the global fit of the distribution, we carried out a residual anal- 61 CUBIC RANKING TRANSMUTED MODEL 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 r/n f 0 1 2 3 0. 0 0. 5 1. 0 1. 5 2. 0 2. 5 x H az ar d Ra te 1 2 3 4 5 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 x Cu m u la tiv e F un ct io n Empirical Trans Order 2 Trans Order 1 Weibull tempo D en si ty 0 1 2 3 4 5 6 0. 0 0. 1 0. 2 0. 3 0. 4 Empirical Trans Order 2 Trans Order 1 Weibull Figure 4.3: Upper panels: TTT plot and hazard curves empirical and estimated by cubic rank transmuted Weibull distribution; Lower panels: Cuimulative and density functions estimated by transmuted Weibull and Weibull distributions. Table 4.3: Selection criteria estimated for six different lifetime distributions. Distribution −2 log AIC AICC BIC KS A∗ W∗ Cubic rank transmuted 282.694 290.694 291.115 301.115 0.647 0.371 0.064 Transmuted quadratic rank 282.702 288.702 288.952 296.518 0.688 0.372 0.074 Weibull 283.059 287.059 287.182 292.269 0.651 0.390 0.068 Exponential 392.742 394.742 394.783 397.347 3.225 17.076 3.392 Gamma 286.467 290.467 290.591 295.678 1.027 0.695 0.147 Log-Normal 296.840 300.840 300.963 306.050 1.313 1.380 0.267 l ll ll lll ll ll l l l l lll l l ll l l l l l l ll ll ll l ll l l llll l lll lll lll l l l l l ll l ll l ll ll lll l ll llll l 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Estimate of EDF Es tim at e CD F fo r Tr a n s W e ib u ll 2 M od el l ll ll lll ll ll l l l l lll l l ll l l l l l l ll ll ll l ll l l llll l lll lll lll l l l l l ll l ll l ll ll lll l ll llll l 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Estimate of EDF Es tim at e CD F fo r Tr a n s W e ib u ll 1 M od el l ll ll lll ll ll l l l l lll l l ll l l l l l l ll ll ll l ll l l llll l lll lll lll l l l l l ll l ll l ll ll lll l ll llll l 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Estimate of EDF Es tim at e CD F fo r W e ib u ll M od el Figure 4.4: P-P plot estimated for transmuted cubic and quadratic rank Weibull and Weibull distributions. 62 CUBIC RANKING TRANSMUTED MODEL yses by using the Martingale-type residual and the deviance residual; see, for example, McCullagh and Nelder (1989), Barlow and Prentice (1988) and Therneau et al. (1990) for pertinent details. The first one, the martingale-type residual, was introduced by Therneau et al. (1990) and is based on a counting process. For the cubic rank transmuted Weibull distribution, the martingale-type residual can be expressed as rMi = (3− λ1 − λ2)ζ(ti) + (λ1 + 2λ2 − 3)(ζ(ti))2 + (1− λ2)(ζ(ti))3, ti > 0, (4.16) where ζ(ti) = e −(ti/β)µ , i = 1, . . . , n. In addition, it is possible to use the deviance residual that has been widely used in generalized linear models. This was also proposed by the same authors (Therneau et al. (1990)) and it is a transformation of the martingale residual to attenuate the skewness. In our case, the deviance residuals are given by rDi = sign(rˆMi) [−2 (rˆMi + log (1− rˆMi))] (4.17) for i = 1, . . . , n. Figure 4.5 shows the Martingale and deviance residuals. l ll l l l l l l l l l l l l ll l l l l l l l l ll l l l l l l l l l ll l l l llll l l l l l l l l l l l l ll l l l l l l l l l l l l l l l l l l ll l l l l l lll l l l llll lll l l l 0 20 40 60 80 100 − 4 − 2 0 2 4 Index M ar tin ga le R es id ua ls 25 19 (a) l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l llll l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l ll l l l l l l 0 20 40 60 80 100 − 3 − 2 − 1 0 1 2 3 Index D ev ia nc e R es id ua ls (b) Figure 4.5: (a) Martingale residuals; (b) Deviance residuals. Global and local influence In this section, we perform global and local influence analyses by considering the cubic rank transmuted Weibull distribution as presented in Section 3.2.1. For analysing the local influence, we consider the response variable perturbation, i.e., we consider that each ti is peturbed as tim = ti+miVt , where Vt is a scale factor that may be the estimated standard deviation of T and mi ∈ R. Then, the perturbed log-likelihood 63 CUBIC RANKING TRANSMUTED MODEL function can be expressed as ` (ψ|t,m) = n log(µ)− nµ log(β) + (µ− 1) n∑ i=1 log tim − n∑ i=1 ( tim β )µ (4.18) + n∑ i=1 log [ (3− λ1 − λ2) + 2(λ1 + 2λ2 − 3)e−(tim/β)µ + 3(1− λ2)e−2(tim/β)µ ] . We can see possible presence of influential points in panels a and b of Figure 4.6, which show the generalized Cook's and likelihood distances. 0 20 40 60 80 100 0. 0 0. 2 0. 4 0. 6 Index Co ok 's D is ta nc e 25 (a) 0 20 40 60 80 100 0 2 4 6 8 10 Index Li ke lih oo d Di st an ce 25 13 19 56 98 (b) Figure 4.6: (a) Cook's distance; (b) Likelihood distance. Note that observation 25 appears to be a discrepant point under both the residual analysis of adjusted distribution and Cook's and likelihood distances for perturbed and non-perturbed data. To assess the impact of the detected influential observations, the RCθj was calculated and the impact was measured by using the total and maximum relative changes and the likelihood displacement given by TRC = ∑np i=1 ∣∣RCθj ∣∣ = 4.18, MRC = maxj ∣∣RCθj ∣∣ = 3.65 and LD(I)(θ) = 2{l(θˆ− l(θˆI)} = 9.79, with np = 4 (the num- ber of parameters). In fact, we consider this point as an outlier and the cubic rank trans- muted Weibull distribution was re-fitted. The MLEs of the parameters and their respective standard deviations (in brackets) are found to be µˆ = 2.702 (1.119), βˆ = 2.731 (1.161), λˆ1 = 0.721 (1.289), and λˆ2 = 0.968 (1.617). The re-fitted distribution can be seen in Figure 4.7. Also, the estimated mean of breaking stress to the carbon fibres turns out to be 3.110 and the corresponding standard error to be 1.025. 64 CUBIC RANKING TRANSMUTED MODEL 1 2 3 4 5 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 x Cu m u la tiv e F un ct io n Empirical Trans Order 2 Figure 4.7: Cumulative curves: empirical vs estimated for cubic rank transmuted Weibull distri- bution. 4.5.3 Application: Cattle sexual precocity data In order to illustrate the usefulness of the cubic rank transmuted log-logistic distribu- tion for modelling, in this section we consider a data set concerned with a study on the economic results concerning beef cattle which are directly related to their genetic prepo- tency, presented in Section 2.5. The transmuted quadratic and cubic rank log-logistic distributions and the log-logistic distribution were fitted to these data. Table 4.4 provides the MLEs, their corresponding standard errors and 95% confidence intervals for the model parameters. For these distri- butions, the computed −2 log -likelihood and AIC values are presented in Table 4.5. Table 4.4: MLEs of the parameters of the transmuted quadratic and cubic ranks log-logistic and log-logistic distributions. Distribution Parameter Estimate Standard 95% Confidence Interval Error Lower Upper µ 2.8758 0.0208 2.8351 2.9165 Transmuted β −16.0305 0.1263 −16.2781 −15.7829 Cubic Rank λ1 0.3201 0.0131 0.2942 0.3457 λ2 −0.8505 0.0304 −0.9101 −0.7908 Transmuted µ 3.0503 0.0220 3.0073 3.0933 Quadratic Rank β −17.9079 0.1369 −18.1763 −17.6395 λ −0.8139 0.0120 −0.8374 −0.7904 Log-Logistic µ 2.5832 0.0160 2.5519 2.6145 β −16.7482 0.1003 −16.9448 −16.5516 Both criteria provide evidence in favor of the cubic rank transmuted log-logistic dis- tribution. This result is corroborated further by the fitted distribution over the histogram of the data presented in the upper right panel of Figure 4.8, and by the fitted survival functions under the Kaplan-Meier estimator on the lower left panel of Figure 4.8. In order to examine the global fit of the distribution, the P-P plot is presented in Figure 4.9. 65 CUBIC RANKING TRANSMUTED MODEL Table 4.5: Computed −2 log likelihood and AIC values for the log-logistic and transmuted log- logistic of orders 1 and 2 distributions. Distribution −2 log AIC AICC BIC A∗ W ∗ Transmuted cubic 235160 235168 235168 235199 94.15 19.30 Transmuted quadratic 236201 236207 236207 236230 142.53 26.63 Log-Logistic 236745 236749 236749 236764 3113.52 662.43 time (days) D en si ty 0 500 1000 1500 2000 2500 3000 0. 00 00 0. 00 05 0. 00 10 0. 00 15 0. 00 20 0. 00 25 0. 00 30 Transmuted 2 Transmuted 1 Log−Logistic 0 500 1000 1500 2000 2500 3000 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 time (days) Su rv iva l Kaplan−Meier Transmuted 2 Transmuted 1 Log−Logistic 0 1000 2000 3000 4000 5000 0. 00 0 0. 00 1 0. 00 2 0. 00 3 0. 00 4 time (days) h(t ) Figure 4.8: Upper Panels: histogram and survival curves; Lower Panel: hazard estimate curve. llll llll lll llll lllll llll llll llll lll lll lll lll lll lll lll lll lll lll lll lll lll lll lll lll lll ll lll lll lll llll llll lllll llll lllll llllll lllll lllll lllll llll llll lll lll lll lll lll llll lllll lll lll lll lll lll ll 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Estimate of EDF Es tim at e CD F fo r Tr a n s Lo g− Lo g 2 M od el lll llll llllll lllll lllll llll llll llll lll lll lll lll lll lll lll lll lll lll ll lll lll lll lll lll ll lll ll lll lll lll llll llll lllll llll lllll llllll llllll lllll lllll llll llll lll lll lll lll lll llll lllll lll lll lll lll lll ll 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Estimate of EDF Es tim at e CD F fo r Tr a n s Lo g− Lo g 1 M od el lll llll llll llll lll lll lll lll lll lll lll lll lll lll lll lll lll lll lll lll lll lll lll lll lll ll lll lll lll lll lll lll lll llll lllll lllll lllll lllll llll llll llll llll llll llll lllll llllll llllllll lllllllll llllll lllll lllllllll lllllll llllllllll lllllll 0.0 0.2 0.4 0.6 0.8 1.0 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Estimate of EDF Es tim at e CD F fo r Lo g− Lo gi st ic M od el Figure 4.9: P-P plot estimated for transmuted cubic and quadratic rank log-logistic and log-logistic distributions, respectively. 66 CUBIC RANKING TRANSMUTED MODEL Upon using the most appropriate distribution for describing the time up to first calving, viz., the cubic rank transmuted log-logistic distribution, we turn our attention to estimating the mean time for the first calf to occur and it is determined to be 1252.54 days (or approximately 41.75 months), with a standard error of 308.32 days. Also, the hazard curve possesses two modes: the first one occuring at 814.5 days and the second one occuring at 1205 days. Chapter 5 e-Transmuted family of distribution Although transmutation maps are a convenient way to construct new distributions, the restricted parametric space of the extra parameter λ may be a problem in some situations. As an alternative to this class of model, in this Chapter we present the regression e- Transmuted family of model (or Exponential Transmuted), which has the property that the extra parameter λ can take any real value. 5.1 Formulation of the model Let Y be a continuous and positive random variable that represents the lifetime of an individual in some population of interest. Consider the parametric family of cumula- tive distribution functions G = {G(y|θ) : θ ∈ Θ} where θ contains a set of k unknown parameters with Θ ⊆ Rk. The first objective of this work is to construct a larger family of distributions that contains the baseline family G. More specifically, we would like to find a family of distributions F (y|θ, ω) where ω ∈ R that satisfies F (y|θ, ω = 0) = G(y|θ) , for all y , θ ∈ Θ . (5.1) In this sense, the new family FG = {F (·|θ, ω) : θ ∈ Θ , ω ∈ R} not only generalizes G but is also centred around G in the scale defined by the parameter ω. The proposed generalized family is defined next. Definition 5.1.1 Given a parametric family of distributions G, the cumulative distribution functions in the extended family FG are defined as follows: F (y|θ, ω) = 1− e −ωG(y|θ) 1− e−ω , ω ∈ R ,θ ∈ Θ. (5.2) We will call the generated family FG, the exponential extension (or e-extension) of the baseline family G and we will use the notation Y ∼ e-extension(ω,G(·|θ)). 67 68 E-TRANSMUTED FAMILY OF DISTRIBUTION It is trivial to see that F (y|θ, ω) is in fact, a cumulative distribution function. In order to graphically see that the e-extended family is centred at G, we can take, for example, G to be an exponential distribution with unit mean, that is G(y|θ) = 1− exp(−y). Cumulative distribution functions in FG have been plotted in Figure 5.1 where we can clearly see that the functions are centred around the baseline distribution G either above or below, depending on the sign of ω. 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 y F ω > 0 ω = 0 ω < 0 Figure 5.1: The effect of ω in the cdf curves of the e-extended family with G(y) exponentially distributed. In more generality, that is, for any distributionG(y|θ), a simple application of L'Hospital's rule shows that F (y|θ, ω = 0) = G(y|θ) since lim ω→0 F (y|θ, ω) = lim ω→0 ∂/∂ω [ ωe−ωG(y|θ) ] ∂/∂ω [ 1− e−ω] = G(y|θ) for all y,θ In this Chapter, we focus on distributions defined on the positive real line but clearly our construction is valid for any support of the random variable Y . Also, we will assume from now on that the functional form of the distributions in G is known. Now, the probability density functions in an e-extended family can be written as: f(y|θ, ω) = ωg(y|θ)e −ωG(y|θ) 1− e−ω , θ ∈ Θ , ω ∈ R, (5.3) where g(y|θ) is the density corresponding to G(y|θ). Applying L'Hospital rule again, shows that f(y|θ, ω = 0) = g(y|θ). The following result, provides a simple probabilistic interpretation of our construction (5.2) and (5.3). More explicitly, a random variable Y that follows an e-extension of a given G can be constructed via a truncated (to lie in the unit interval (0, 1)) and then transformed 69 E-TRANSMUTED FAMILY OF DISTRIBUTION (by the quantile function G−1) exponential random variable. Result 1 Let Z be an exponential random variable with rate ω which has been truncated to lie in the unit interval (0, 1). Then Y = G−1(Z|θ) ∼e-Transmuted(ω,G). We clearly have that P (Z ≤ z) = 1− exp(−ω z) 1− exp(−ω) for any z ∈ (0, 1) so that for any y > 0 we have P (Y ≤ y) = P (G−1(Z|θ) ≤ y) = P (Z ≤ G(y|θ)) = 1− exp(−ωG(y|θ)) 1− exp(−ω) . This result can be seen as a particular case of the construction of extended parametric models presented by Verdinelli et al. (1998). Clearly, the result also provides a simple method to generate random samples from an e-extended family. In particular, the lower qth quantile of e-extended model is given by yq = G −1(zq(ω)|θ) where zq(ω) = − log [1− q(1− e −ω)] ω that is, the composition of two quantiles, first for the truncated exponential (zq(ω)) and then the quantile G−1(·|θ) of the baseline family G . Generation a random samples from an e-extended family is then trivial by using the inverse method. 5.2 Hazard and related functions Given that we are interested in distributions with positive support we develop here ex- pressions for survival and hazard functions. The survival function, which is the probability of the survival time Y being larger than or equal to some time y, has a simple expression in terms of the baseline survival function S(y|θ) = 1−G(y|θ), namely S(y|θ, ω) = e ω S(y|θ) − 1 eω − 1 , (5.4) so that clearly S(y|θ, ω = 0) = S(y|θ). The other function of interest is the hazard rate function, which for an e-extended model can be written as h(y|θ, ω) = ω h(y|θ)S(y|θ) 1− e−ωS(y|θ) , (5.5) where h(y|θ) = g(y|θ)/S(y|θ) is the baseline hazard function. Figure 5.2 shows the be- haviour of the hazard rate functions for different values of ω when the baseline distribution 70 E-TRANSMUTED FAMILY OF DISTRIBUTION is exponential with mean one. y h ω > 0 ω = 0 ω < 0 Figure 5.2: The effect of ω in the hazard curves with G(y) exponential. We can see from the plot that the hazard rate curves are all asymptote to the constant baseline when y is large. This result can be generalized as follows. Proposition 5.2.1 For all θ and ω, we have that h(y|ω,θ) ∼ h(y|θ) when y →∞. Proof: The hazard rate function (5.5) can be written as h(y|θ, ω) = h(y|θ)×R(y,θ) where R(y,θ) = ωS(y|θ) 1− e−ωS(y|θ) . For any θ ∈ Θ and ω ∈ R, we have that S(y|θ) ∼ 0 and 1− exp(−ω S(y|θ)) ∼ 0 when y is large, so an application of L'Hospital rule shows that R(y,θ) ∼ 1, which proves the result.  5.2.1 Related distributions The modified negative Gompertz distribution (see Marshall and Olkin (2007b) page 390 and Dahiya and Hossain (1996)) can be obtained from the construction (5.2) in the particular case where G is the exponential distribution family with rate θ, that is, when G(y|θ) = 1−exp(−θ y). It is interesting to note that we naturally allow the shape parameter ω to be negative but this does not seem to be acknowledged in the work of Marshall and Olkin (2007b) and Dahiya and Hossain (1996). The corresponding density can be written as f(y|ω, θ) = sign(ω)e −ω σ(θ)(1− e−ω) exp { −y − κ(ω, θ) σ(θ) + sign(ω) exp ( −y − κ(ω, θ) σ(θ) )} , (5.6) 71 E-TRANSMUTED FAMILY OF DISTRIBUTION with ω ∈ R, κ(ω, θ) = log |ω|/θ and σ(ω) = 1/θ. Therefore, if ω < 0, expression (5.6) shows that the e-extended model corresponds to a subset of truncated (to the positive real line) Gumbel distribution. For a given value of θ ∈ Θ, the densities (5.3) form an exponential family with natural parameter ω and sufficient statistic G(y|θ). This will clearly have important implications for making inference on the family F , namely, for maximum likelihood. In fact, our con- struction can be seen as a particular case of a more general construction of the form F (y|ω,θ) = g(y|θ) exp (ωR(y|θ)− ψ(ω,θ)) (see, for example, Rayner et al. (2009)), where our construction corresponds to the choice R(y|θ) = −G(y|θ) and as a consequence of that choice, the function ψ(ω,θ) = log((1 − e−ω)/ω) does not depend on θ. Another related construction is the so-called rank transmutation mapping introduced by Shaw and Buckley (2009). The authors define a rank transmutation map, in the con- tinuous case, as T (u) = F (G−1(u)) where F and G are continuous cumulative distribution functions with the same support and varying in different parametric families. The authors use this, and other related mappings, with copula-base simulation in mind, which we do not pursue here, but we can think of their construction as similar to ours in the following sense. Assume, for illustration purposes, that the functions F and G vary in parametric families so that the mapping T can be indexed by a parameter ω so that we have T (u|ω). If the identity mapping is included into consideration, that is, T (u|ω = 0) = u then clearly T (G(y|θ)|ω) will generate a family that includes G as a particular case when ω = 0. The simplest example is the quadratic transmutation T (u|ω) = u+ ω 2 u(1− u). Then we can construct a larger family via F (y|θ, ω) := T (G(y|θ)|ω) = G(y|θ) + ω 2 G(y|θ)S(y|θ), (5.7) with the only difference that this instance requires that ω ∈ [−2, 2] since any cumulative distribution function is bounded between 0 and 1. There is a simple heuristic relationship between the two approaches which helps to understand the parametrizations used. We can write the densities corresponding to (5.7) as: f(y|θ, ω) = f(y|θ) [ 1 + ω ( 1 2 −G(y|θ) )] 72 E-TRANSMUTED FAMILY OF DISTRIBUTION For small values of ω(1−G(y|θ)) we have that 1 + ω ( 1 2 −G(y|θ) ) ≈ exp(1/2−G(y|θ)) and a renormalization (in order to integrate to one) will arrive at the same expression as (5.3), hence the use of the same parameter ω in both (5.2) and (5.7). There also is a relationship with the so-called local mixture models of Anaya-Izquierdo and Marriott (2007). In this case the extension of the densities is as follows f(y|θ, ω) = f(y|θ) [1 + ωR(y,θ)] , where the function R(y,θ) requires that the expectation with respect to f(y|θ) is zero. In this case the choice R(y,θ) = 1/2−G(y|θ) gives clearly the same distribution as (5.7). Our construction has some advantages over the alternative construction in the sense of having a wider range of behaviours. Figure 5.3 shows the comparison between the e- extended and the quadratic rank transmuted model, respectively, by assuming that G(y|θ) is exponentially distributed. We can observe that the e-extended model has a wider range of shapes when the parameter |ω| ≥ 2 and a similar behaviour to the rank transmuted one when the parameter |ω| < 2. e−extended y f − 2 < ω < 0 ω = 0 0 < ω < 2 transmuted y f − 2 < ω < 0 ω = 0 0 < ω < 2 e−extended y f ω ≤ − 2 ω = 0 ω ≥ 2 Figure 5.3: The effect of ω in the density curves of the e-extended and transmuted families, with G(y|θ) exponential. 5.3 E-extended Weibull Given its central place in survival analysis, we focus here on e-extensions of the Weibull distribution. If Y is a positive random variable following a Weibull distribution then the 73 E-TRANSMUTED FAMILY OF DISTRIBUTION probability density and cumulative distribution functions are given by g(y|φ, β) = β φβ yβ−1e−(y/φ) β , G(y|φ, β) = 1− e−(y/φ)β , (5.8) where φ > 0 and β > 0 are scale and shape parameters, respectively, so in this case we have that θ = (φ, β). The shape parameter β is linked to the hazard behavior: if β ∈ [0, 1) the hazard is decreasing; if β > 1 the hazard is increasing and; the exponential distribution is a particular case if β = 1 and the hazard is constant. By using the expressions above the Weibull e-extended model density can be written as: f(y|φ, β, ω) = ωβe −ω (1− e−ω)φβ y β−1 exp { − ( y φ )β + ω exp [ − ( y φ )β]} , (5.9) with ω ∈ R and (φ, β) ∈ R2+. Figure 5.4 shows some examples of probability density Weibull e-extended curves for different values of parameters β and ω, with fixed φ = 0.7. 0 1 2 3 4 5 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 x f*( x) ω = 0.5 ω = −5.0 ω = 2.5 ω = −2.5 ω = 5.0 ω = 5.0 β = 0.5 β = 2.0 β = 1.5 β = 1.5 β = 3.5 β = 1.0 Figure 5.4: Probability density and cumulative distribution function Weibull e-transmuted model, respectively, for different values of ω and φ = 0.7. The corresponding hazard rate function is given by h(y|ω, φ, β) = ωβy β−1e−(y/φ)β eφβ [ e−γ(y)e−(y/φ) β − 1 ] . (5.10) There are three different behaviours of the hazard rate function, which are described next: Case I: if the baseline model is Weibull with parameter β > 1, i.e., the base distribution 74 E-TRANSMUTED FAMILY OF DISTRIBUTION hazard is increasing, see for example Figure 5.5-(a). • For ω ≤ 0 the hazard curves continue being increasing but they are dominated by the Weibull baseline model; • For ω > 0 the hazard curves present two inflexion points, i.e., initially increase until the point t1 when start to decrease up to t2 and increase again for times greater than t2 (the changes of curves behavior are evident for greater values of ω). Case II: if the baseline model is Weibull with parameter β < 1, i.e. the baseline distribu- tion hazard is decreasing, then the behavior of the Weibull e-extended hazard is the same for all ω, see for example Figure 5.5-(b). Case III: if the baseline model is Weibull with parameter β = 1, i.e the baseline distri- bution is exponential and the hazard is constant, see for example Figure 5.5-(c). • For ω > 0 the hazard of Weibull e-extended model is decreasing; • For ω < 0 the hazard increasing. Note that in all cases, as explained by Proposition 5.2.1, the e-extended hazard converges to the baseline hazard for large values of y. 0 1 2 3 4 5 x h* (x) 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 1. 2 x h* (x) 0 1 2 3 4 5 x h* (x) (a) (b) (c) Figure 5.5: The effect of ω in the hazard curves with G(y) Weibull(µ, β): (a) β > 1; (b) β < 1; and (c) β = 1. 5.4 Parameter estimation 5.4.1 Underlying exponential family structure One of the advantages of this choice is that we can found easily expressions for entries of the expected Fisher information. 75 E-TRANSMUTED FAMILY OF DISTRIBUTION This particular choice of R (R(y|θ) = −G(y|θ)), presented in Section 5.2.1 will give orthogonality between the parameters θ and ω at the underlying base model G. Note that, −∂ 2 log f(y|θ, ω) ∂θ ∂ω = ∂ ∂θ G(y|θ). Assuming the standard regularity conditions on the parametric family G we can differen- tiate under the integral sign so that E [ −∂ 2 log f(Y |θ, ω) ∂θ ∂ω ] = E [ ∂ ∂θ G(Y |θ) ] = ∂ ∂θ E [G(Y |θ)] , where the expectation is taken with respect to f(y|θ, ω). At the baseline model (ω = 0) we have that f(y|θ, ω) = f(y|θ) so that clearly E[G(y|θ)] = 1/2 since, in that case, G(Y |θ) follows a uniform distribution on the unit interval. Therefore it follows that ∂ ∂θ E [G(Y |θ)] = 0. The above result implies that maximum likelihood can proceed iteratively and also using the so-called profile likelihood. Fixing a value of θ the MLE of ω is easily found since the resultant family is of exponential form. Let y1, . . . , yn be a random sample of size n from an e-extended distribution. We will allow for unknown parameters in the baseline distribution G, namely G(yi|θ). Given the exponential family structure defined above, if we fix the values of θ then we can find the corresponding maximum value of ω trivially by solving the non-linear equation k(ω) = ∑n i=1G(yi|θ) n := G¯(θ). (5.11) Now, suppose X = (X1, · · · , Xp)′ a vector of covariates and consider the following relationship between the vector of covariates and the parameters γ(x) = γ1x1 + . . .+γpxp. Then the likelihood function is given by L(ω,θ|y) = n∏ i=1 p∏ j=1 γjxij (1− e−γjxij )g(yi|θ)e −γjxijG(yi|θ). (5.12) Hence, the log-likelihood function becomes l = logL(ω,θ|y) = n∑ i=1 p∑ j=1 log ( γjxij 1− e−γjxij ) + p n∑ i=1 log g(yi|θ) − n∑ i=1 p∑ j=1 γjxijG(yi|θ). (5.13) 76 E-TRANSMUTED FAMILY OF DISTRIBUTION Therefore, the maximum likelihood estimates of γ and the vector of parameters θ, inher- ent to the baseline distribution, which maximize (5.13) must satisfy the following normal equations ∂l ∂γj = n∑ i=1 x∗i p∑ j=1 (1− e−γjxij − γjxije−γjxij ) γjxij (1− e−γjxij ) − n∑ i=1 x∗iG(yi|θ) = 0, ∂l ∂θj = p n∑ i=1 g′θj (yi|θ) g(yi|θ) − n∑ i=1 p∑ j=1 γjxij n∑ i=1 G′θj (yi|θ) = 0, with x∗i the values of xi fixed at the point j. Note that, the maximum likelihood estimator (γˆ, θˆ) is obtained by solving the above nonlinear system of equations. Furthermore, the observed Fisher Information matrix is given in Appendix A.4. In order to compute the standard error and asymptotic confidence interval we use the usual large sample approximation in which the maximum likelihood estimators of γ,θ can be treated as being approximately (p + k)−variate normal, p is the number of covariates and k the number of parameters of the baseline distribution. For example, as n→∞ the asymptotic distribution of the MLE (γˆp, θˆk), for p = k = 1, is given by,( γˆ θˆ ) ∼ N [( γˆ θˆ ) , ( Vˆ11, Vˆ12 Vˆ21, Vˆ22 )] , (5.14) with, Vˆij = Vij |θ=θˆ and it is determined by the inverse of Fisher information matrix. Thereby, an approximate 100(1 − α)% two sided confidence intervals for γ and θ are, respectively, given by γˆ ± zα/2 √ Vˆ11 and θˆ ± zα/2 √ Vˆ22 , where zα is the upper α− th percentile of the standard normal distribution. 5.5 Weibull e-extended Model 5.5.1 Parameter estimation The maximum likelihood estimates, MLEs, of the parameters of the Weibull regression e-transmuted probability distribution function are given by the following: Let Y1, . . . , Yn be a sample of size n from a Weibull e-transmuted distribution and x = (x1, . . . , xp) ′ the vector of covariates with γ(x) = γ1x1 + . . . + γpxp. Then the likelihood function of the Weibull regression e-transmuted model is given by 77 E-TRANSMUTED FAMILY OF DISTRIBUTION L(ω, φ, β|y) = β np φnpβ n∏ i=1 p∏ j=1 {( γjxij 1− e−γjxij ) e−γjxijyβ−1i × exp { − ( yi φ )β + γjxij exp [ − ( yi φ )β]}} . (5.15) Hence, the log-likelihood function becomes l = logL(ω, φ, β|y) = np log β − npβ log φ+ n∑ i=1 p∑ j=1 log ( γjxij 1− e−γjxij ) − n∑ i=1 p∑ j=1 γjxij +(β − 1)p n∑ i=1 log yi − n∑ i=1 p∑ j=1 [( yi φ )β − γjxije−(yi/φ)β ] . (5.16) Therefore, the maximum likelihood estimates of γj , φ and β which maximize (5.16) must satisfy the following normal equations: ∂l ∂γj = n∑ i=1 x∗i p∑ j=1 (1− e−γjxij − γjxije−γjxij ) γjxij (1− e−γjxij ) − n∑ i=1 x∗i + p n∑ i=1 x∗i e −(yi/φ)β = 0, ∂l ∂φ = npβ φ + p φ n∑ i=1 ϕi + n∑ i=1 p∑ j=1 γjxij φ ϕie −ϕi = 0, ∂l ∂β = np β − npφ+ p n∑ i=1 log yi − n∑ i=1 p∑ j=1 ϕi log ( yi φ ) − n∑ i=1 p∑ j=1 γjxijϕi log ( yi φ ) e−ϕi = 0, with ϕi = ( yi φ )β and x∗i the values of the covariate xi fixed at the point j. We denote, the maximum likelihood estimator (γˆj , φˆ, βˆ) and the observed Fisher Infor- mation is indicated in Appendix A.4. In order to compute the standard error and asymptotic confidence intervals we use the usual large sample approximation in which the maximum likelihood estimators of γj , φ, β can be treated as being approximately (p + 2)−variate normal. For example, hence as n→∞, fixing j = 1, the asymptotic distribution of the MLE (γˆ1, φˆ, βˆ), are given by, γˆ1φˆ βˆ  ∼ N   γ1φ β   V11 V12 V13V21 V22 V23 V31 V32 V33   , (5.17) with Vˆij = Vij |θ=θˆ and it is determined by the inverse of observed Fisher information 78 E-TRANSMUTED FAMILY OF DISTRIBUTION matrix. Therefore, an approximate 100(1− α)% two sided confidence intervals for γ1, φ and θ are, respectively, given by γˆ1 ± zα/2 √ Vˆ11, φˆ± zα/2 √ Vˆ22 and βˆ ± zα/2 √ Vˆ33, where zα is the upper α− th percentile of the standard normal distribution. 5.5.2 Numerical experiments for the Weibull e-transmuted model This section presents the results of a Monte Carlo experiment to investigate the finite sample behavior of the MLEs. First of all, we consider the Weibull e-transmuted model without covariates and all results were obtained from 5000 Monte Carlo replications. The sample sizes n range from 50 to 1000, generated according to a Weibull e-transmuted distribution as presented in Section 5.3. For that, since φ is a location parameter the results of the Monte Carlo experiment are unaffected by the choice of φ. We fixed the value of parameter φ = 2 and we consider different combinations of ω and β parameters values. We consider a wide range of hazard behaviors, namely β = 0.5, 1.0 and 1.5 (decreasing, constant and increasing ones, respectively). Also, we considered positive and negative ω's. In order to check the coverage probability of the confidence intervals a boostrap study was made. For that, we generated 200 samples with sample sizes ranging from 50 to 500, generated according to a Weibull e-transmuted distribution, without covariates, and con- sidering several values of parameters. We fixed the value of parameter φ = 2 and considered different combinations of ω and β. For each of the 200 samples, we resampled B = 1000 and constructed the intervals by using the p-boostrap method; Gentle (2009). The coverage probability of these intervals is presented in Figure 5.6, for the parameters ω, µ and β, left, middle and right panels respectively, where we can see the closeness to the nominal value, 95%, in all simulated situations. Sample Size co ve ra ge P ro ba bi lity l l l l l 50 80 100 200 500 0. 5 0. 6 0. 7 0. 8 0. 9 1. 0 l ω = −2.0 ω = 2.0 ω = −2.0 ω = 2.0 β = 0.5 β = 0.5 β = 1.5 β = 1.5 Sample Size co ve ra ge P ro ba bi lity l l l l l 50 80 100 200 500 0. 5 0. 6 0. 7 0. 8 0. 9 1. 0 l ω = −2.0 ω = 2.0 ω = −2.0 ω = 2.0 β = 0.5 β = 0.5 β = 1.5 β = 1.5 Sample Size co ve ra ge P ro ba bi lity l l l l l 50 80 100 200 500 0. 5 0. 6 0. 7 0. 8 0. 9 1. 0 l ω = −2.0 ω = 2.0 ω = −2.0 ω = 2.0 β = 0.5 β = 0.5 β = 1.5 β = 1.5 Figure 5.6: Coverage probability of intervals for ω, µ and β parameters, respectively, generated by using the p-bootstrap method for φ = 2 fixed and different values of ω and β. Now, we consider the Weibull regression e-transmuted model in the presence of two continuous covariates and all results were obtained from 5000 Monte Carlo replications. 79 E-TRANSMUTED FAMILY OF DISTRIBUTION Table 5.1: Results of a Monte Carlo experiment for the Weibull e-transmuted model with φ = 2 and different combinations of ω and β values. The estimated values presented are the median values of the 5000 generated samples. Within each horizontal block, the generated samples were subsampled from the larger one which was generated by using the inverse CDF method. Generated Sample Estimated Standard Error MSE (ω;β) Size ω φ β −2 log ω φ β ω φ β (−0.5; 0.5) 50 −0.409 2.047 0.513 338.2 1.900 1.057 0.104 3.556 1.696 0.116 80 −0.451 2.027 0.506 542.3 1.624 0.923 0.088 3.029 1.461 0.098 100 −0.444 2.033 0.507 679.5 1.516 0.858 0.082 2.785 1.338 0.091 200 −0.467 2.016 0.502 1359.2 1.180 0.680 0.066 2.057 1.010 0.072 500 −0.468 2.014 0.502 3404.0 0.832 0.483 0.049 1.311 0.661 0.052 1000 −0.490 2.002 0.501 6803.6 0.629 0.368 0.038 0.904 0.472 0.040 (0.5; 0.5) 50 0.390 1.917 0.501 280.2 1.911 1.131 0.090 3.599 1.770 0.100 80 0.420 1.948 0.498 447.2 1.659 0.988 0.075 2.985 1.500 0.083 100 0.468 1.961 0.499 559.4 1.540 0.925 0.070 2.710 1.380 0.076 200 0.465 1.979 0.498 1121.7 1.185 0.711 0.054 1.970 1.003 0.057 500 0.491 1.992 0.500 2805.6 0.800 0.479 0.037 1.125 0.609 0.038 1000 0.477 1.987 0.499 5612.0 0.575 0.344 0.027 0.756 0.414 0.028 (−2; 1) 50 −1.577 2.103 1.074 334.6 1.994 0.481 0.253 4.574 0.683 0.327 80 −1.699 2.076 1.048 536.6 1.719 0.421 0.218 3.678 0.586 0.275 100 −1.778 2.062 1.040 671.6 1.603 0.397 0.202 3.356 0.552 0.254 200 −1.845 2.038 1.029 1346.2 1.258 0.317 0.162 2.493 0.434 0.196 500 −1.904 2.024 1.015 3369.9 0.891 0.230 0.119 1.576 0.298 0.138 1000 −1.953 2.012 1.009 6744.2 0.682 0.179 0.093 1.061 0.220 0.103 (2; 1) 50 1.554 1.843 0.975 240.3 2.011 0.612 0.143 3.912 0.810 0.163 80 1.671 1.887 0.978 385.2 1.748 0.533 0.115 3.328 0.693 0.126 100 1.708 1.905 0.977 482.4 1.628 0.495 0.102 2.977 0.637 0.112 200 1.854 1.946 0.985 966.8 1.274 0.386 0.073 2.128 0.478 0.077 500 1.904 1.973 0.992 2422.8 0.865 0.262 0.046 1.266 0.305 0.048 1000 1.943 1.984 0.996 4849.0 0.627 0.188 0.032 0.827 0.211 0.033 (−0.5; 1.5) 50 −0.367 2.022 1.546 282.1 1.881 0.349 0.309 3.477 0.430 0.403 80 −0.393 2.025 1.535 452.9 1.614 0.303 0.266 2.984 0.373 0.343 100 −0.391 2.020 1.528 566.9 1.496 0.284 0.246 2.785 0.347 0.316 200 −0.429 2.011 1.515 1136.4 1.185 0.227 0.199 2.095 0.273 0.248 500 −0.461 2.005 1.504 2845.3 0.840 0.162 0.147 1.341 0.187 0.173 1000 −0.486 2.001 1.499 5692.9 0.633 0.124 0.115 0.909 0.138 0.129 (0.5; 1.5) 50 0.355 1.970 1.489 269.1 1.908 0.379 0.269 3.538 0.458 0.335 80 0.406 1.978 1.489 431.5 1.645 0.329 0.227 3.033 0.395 0.271 100 0.429 1.983 1.491 540.2 1.523 0.305 0.208 2.678 0.360 0.246 200 0.452 1.988 1.493 1082.4 1.173 0.236 0.162 1.886 0.272 0.182 500 0.480 1.995 1.498 2709.7 0.794 0.159 0.112 1.129 0.176 0.120 1000 0.493 1.999 1.501 5424.6 0.570 0.114 0.081 0.739 0.123 0.085 Similar to the first numerical example, the sample sizes n range from 50 to 1000, generated according to a Weibull regression e-transmuted distribution as presented in Section 5.3. Again, we fixed the value of parameter φ = 2 and we consider different combinations of γ1, γ2 (the regressors) and β parameters values. We consider a wide range of hazard behaviors, namely β = 0.5, 1.0 and 1.5 (decreasing, constant and increasing ones, respectively). We checked the coverage probability of the confidence intervals for sample sizes ranging from 50 to 2000, generated according to a Weibull e-transmuted distribution, in the pres- ence of two covariates, and considering several values of parameters as presented in Table 5.2. Similar to which was made in the first simulation study, we fixed the value of parame- 80 E-TRANSMUTED FAMILY OF DISTRIBUTION ter φ = 2 and considered different combinations of γ1, γ2 and β. The coverage probability of these intervals is presented in Figure 5.7, for the parameters γ1 and γ2 (upper panels), µ and β (lower panels), respectively, where we can see the closeness to the nominal value, 95%, in all simulated situations. l l l l l l l Sample Size Co ve ra ge P ro ba bi lity l l l l l l l 50 80 100 200 500 1000 2000 0. 5 0. 6 0. 7 0. 8 0. 9 1. 0 l l γ1 = 0.5 γ1 = 2.0 γ1 = 2.0 γ1 = 0.5 γ2 = 2.0 γ2 = 0.5 γ2 = 0.5 γ2 = 2.0 β = 0.5 β = 0.5 β = 1.5 β = 1.5 l l l l l l l Sample Size Co ve ra ge P ro ba bi lity l l l l 50 80 100 200 500 1000 2000 0. 5 0. 6 0. 7 0. 8 0. 9 1. 0 l l γ1 = 0.5 γ1 = 2.0 γ1 = 2.0 γ1 = 0.5 γ2 = 2.0 γ2 = 0.5 γ2 = 0.5 γ2 = 2.0 β = 0.5 β = 0.5 β = 1.5 β = 1.5 l l l l l l l Sample Size Co ve ra ge P ro ba bi lity l l l l l l l 50 80 100 200 500 1000 2000 0. 5 0. 6 0. 7 0. 8 0. 9 1. 0 l l γ1 = 0.5 γ1 = 2.0 γ1 = 2.0 γ1 = 0.5 γ2 = 2.0 γ2 = 0.5 γ2 = 0.5 γ2 = 2.0 β = 0.5 β = 0.5 β = 1.5 β = 1.5 l l l l l l l Sample Size Co ve ra ge P ro ba bi lity l l l l 50 80 100 200 500 1000 2000 0. 5 0. 6 0. 7 0. 8 0. 9 1. 0 l l γ1 = 0.5 γ1 = 2.0 γ1 = 2.0 γ1 = 0.5 γ2 = 2.0 γ2 = 0.5 γ2 = 0.5 γ2 = 2.0 β = 0.5 β = 0.5 β = 1.5 β = 1.5 Figure 5.7: Coverage probability of intervals for γ1, γ2, µ and β parameters, respectively, based on Monte Carlo experiment for φ = 2 fixed and different values of γ1, γ2 and β. 5.5.3 Goodness-of-fit to the Weibull regression e-transmuted model One of the most complex problems in parametric statistical methods is the identification of the distribution that can best fit the data. In general words, the goodness-of-fit of a statistical model describes how well it fits into a set of observations. The goodness-of-fit indices summarize the discrepancy between the observed values and the values expected under a statistical model. Also, the goodness-of-fit statistics are usually obtained using asymptotic methods, that are used in statistical hypothesis testing. As large sample approximations may behave poorly in small samples, a great deal of research using simulation studies has been devoted to investigate under which conditions the asymptotic p-values of goodness-of-fit statistics are accurate (i.e., how large the sample size must be for models of different sizes), Tollenaar and Mooijaart (2003). Some general indices and statistics can be used in goodness-of-fit tests like Akaike in- formation criteria, bayesian factor, Anderson-Darling, Cramér-von Mises, among others. 81 E-TRANSMUTED FAMILY OF DISTRIBUTION Table 5.2: Results of a Monte Carlo experiment for the Weibull regression e-transmuted model with φ = 2 and different combinations of γ1, γ2 and β values. The estimated values presented are the median values of the 5000 generated samples. Within each horizontal block, the generated samples were subsampled from the larger one which was generated by using the inverse CDF method. Generated Sample Estimated Estimated Error MSE (γ1; γ2;β) Size γ1 γ2 φ β γ1 γ2 φ β γ1 γ2 φ β (0.5, 2, 0.5) 50 0.524 2.018 1.924 0.509 1.169 0.972 1.579 0.055 1.969 3.346 2.450 0.056 80 0.513 2.026 1.987 0.505 0.882 0.799 1.347 0.043 1.301 3.164 1.968 0.043 100 0.505 2.018 1.986 0.505 0.772 0.721 1.224 0.038 1.085 3.061 1.747 0.039 500 0.505 1.994 1.990 0.501 0.330 0.329 0.574 0.017 0.383 2.558 0.695 0.017 1000 0.500 1.995 1.993 0.500 0.231 0.233 0.408 0.012 0.257 2.470 0.474 0.012 (2, 0.5, 0.5) 50 2.207 0.508 2.047 0.514 1.354 0.450 1.284 0.060 2.089 2.726 1.833 0.062 80 2.095 0.512 2.059 0.509 1.035 0.349 1.048 0.047 1.490 2.573 1.386 0.048 100 2.078 0.505 2.035 0.507 0.914 0.309 0.923 0.042 1.270 2.548 1.215 0.043 500 2.009 0.499 2.000 0.501 0.398 0.136 0.409 0.018 0.470 2.389 0.478 0.019 1000 2.001 0.502 1.997 0.501 0.280 0.096 0.291 0.013 0.315 2.340 0.324 0.013 (0.5, 2, 1) 50 0.525 2.020 1.963 1.018 1.106 0.908 0.764 0.109 1.971 3.350 1.170 0.115 80 0.514 2.027 1.994 1.011 0.866 0.773 0.667 0.086 1.301 3.165 0.906 0.089 100 0.505 2.018 1.993 1.010 0.762 0.707 0.608 0.076 1.084 3.061 0.801 0.079 500 0.505 1.994 1.995 1.002 0.330 0.329 0.288 0.034 0.383 2.558 0.324 0.034 1000 0.500 1.995 1.997 1.001 0.231 0.233 0.204 0.024 0.257 2.470 0.223 0.024 (2, 0.5, 1) 50 2.207 0.508 2.023 1.029 1.344 0.447 0.638 0.120 2.089 2.726 0.836 0.127 80 2.095 0.512 2.029 1.018 1.033 0.349 0.513 0.094 1.490 2.573 0.630 0.098 100 2.078 0.505 2.017 1.014 0.914 0.309 0.457 0.084 1.270 2.548 0.553 0.087 500 2.009 0.499 2.000 1.001 0.398 0.136 0.205 0.037 0.470 2.389 0.225 0.038 1000 2.001 0.502 1.999 1.002 0.280 0.096 0.145 0.026 0.315 2.340 0.155 0.026 (0.5, 2, 1.5) 50 0.527 2.019 1.975 1.527 1.109 0.911 0.528 0.164 1.970 3.348 0.738 0.176 80 0.514 2.027 1.996 1.516 0.864 0.775 0.449 0.128 1.299 3.165 0.572 0.136 100 0.505 2.018 1.995 1.515 0.764 0.708 0.409 0.114 1.084 3.061 0.509 0.120 500 0.505 1.994 1.997 1.502 0.330 0.329 0.192 0.050 0.383 2.558 0.210 0.052 1000 0.500 1.995 1.998 1.501 0.231 0.233 0.136 0.035 0.257 2.470 0.145 0.036 (2, 0.5, 1.5) 50 2.207 0.508 2.015 1.543 1.345 0.448 0.426 0.179 2.089 2.726 0.528 0.196 80 2.095 0.512 2.020 1.527 1.034 0.349 0.341 0.140 1.490 2.573 0.402 0.150 100 2.078 0.505 2.011 1.521 0.914 0.309 0.306 0.125 1.270 2.548 0.352 0.133 500 2.009 0.499 2.000 1.502 0.398 0.136 0.137 0.055 0.470 2.389 0.147 0.057 1000 2.001 0.502 1.999 1.502 0.280 0.096 0.097 0.039 0.315 2.340 0.102 0.040 Another ones are graphical based, for example the probability-probability and QQ plots; or quantify a distance between the empirical distribution function of the sample and the cumulative distribution function of the reference distribution, for example, the widely used Kolmogorov-Smirnov test. In this section, we will focused in two quadratic distance meth- ods: Anderson-Darling and Cramér-von Mises and consider the following description of the test. First, let Y1, . . . , Yn be random variables from a family of distribution and F (·) the cumulative distribution function. Then, consider the following hypothesis: 82 E-TRANSMUTED FAMILY OF DISTRIBUTION { H0 : F (y|θ) ∈ F = { F (·|ω,θ) : ω ∈ R,θ ∈ Θk} H1 : F (y|θ) /∈ F = { F (·|ω,θ) : ω ∈ R,θ ∈ Θk} . (5.18) It is equivalent to test the alternative hypothesis as being H1 : F (y) 6= F (x|θ),∀θ ∈ Θ. Let ψ(Y ) be a test statistic, in this particular case two that belong to the quadratic class of empirical distribution function statistics, Anderson-Darling and Cramér-von Mises statistics (see, for example, Anderson (1962) and Duchesne et al. (1997)). This class of statistics is based on the squared difference (Fn(y)− F0(y))2 and have the following general form: ψ(Y ) = n ∫ +∞ −∞ ( Fn(y)− F(y|θˆ) )2 κ(y)dF (y) where κ(·) weights the squared difference and Fn(·) is the empirical distribution function and F0(·) is the cumulative distribution function under H0. The Anderson-Darling and the Cramér-von Mises statistics are defined, respectively, by ψAD(Y ) = n ∫ +∞ −∞ ( Fn(y)− F(y|θˆ) )2 [ F(y|θˆ) ( 1− F(y|θˆ) )] dF (y) (5.19) and ψCvM (Y ) = n ∫ +∞ −∞ ( Fn(y)− F(y|θˆ) )2 dF (y). (5.20) The first interest is to investigate the type I error, which is the incorrect rejection of a true null hypothesis, i.e. P (Reject H0|H0 is true) = α1. For that, we need the sampling distribution of ψAD(Y ) and ψCvM (Y ) under the null hypothesis F (x|θ) for some θ ∈ Θ. Note that, the sampling distribution of ψAD(Y ) and ψCvM (Y ) depending on θ. In this case, we use the θˆ0 and use the sampling distribution of ψAD(Y ) and ψCvM (Y ) when θ = θˆ0. Then, we simulated from F (y|θˆ), i.e., from the Weibull e-transmuted model and the quantile 100(1 − α1)% of the Anderson-Darling and the Cramér-von misses statistics were calculated and given by ψ1−α1AD (Y ) and ψ 1−α1 CvM (Y ), respectively. Then, the initial Monte Carlo experiment with 1000 samples was used to set the corre- sponding nominal values when α1 = 5% and 10% for different combinations of the param- eters. In order to check the rejection rate behaviors, we use the non-parametric bootstrap method to resample B = 1000 times from Monte Carlo samples and the rejection rates are showed in the Figure 5.8. Note that, we analysed for different situations: decreasing and increasing hazards (β = 0.5 and 1.5, respectively) and ω parameter being negative and positive. From Figure 5.8, it is clear that in all situations the rejection rates are close to the nominal values. Also, extensions to the theory of hypothesis testing include the study of the power of tests. Then, let's define the type II error as the probability of failure to reject a false null hypothesis, i.e., when the null hypothesis is false but erroneously fails to be rejected, 83 E-TRANSMUTED FAMILY OF DISTRIBUTION Sample size Reje ction Rat e 50 80 100 150 300 500 1000 0.0 0 0.0 5 0.1 0 0.1 5 Anderson Darling Cramér−von Misses β = 0.5 ω = −2.0 Sample size Reje ction Rat e 50 80 100 150 300 500 1000 0.0 0 0.0 5 0.1 0 0.1 5 Anderson Darling Cramér−von Misses β = 0.5 ω = 2.0 Sample size Reje ction Rat e 50 80 100 150 300 500 1000 0.0 0 0.0 5 0.1 0 0.1 5 Anderson Darling Cramér−von Misses β = 1.5 ω = −2.0 Sample size Reje ction Rat e 50 80 100 150 300 500 1000 0.0 0 0.0 5 0.1 0 0.1 5 Anderson Darling Cramér−von Misses β = 1.5 ω = 2.0 Figure 5.8: Probability of error type I by using Anderson-Darling and Cramér-von-Mises statistics for different combinations of the parameters β and ω and fixed φ = 2 by considering α1 = 5% and 10% . P(Not reject H0|H0 is false) = 1− α2. Then the probability of correctly rejecting the null hypothesis given that it is false is the complement of the false negative rate, α2. In order to check the power of the test, initially, let's assume a more general test, pre- 84 E-TRANSMUTED FAMILY OF DISTRIBUTION sented in the composite hypothesis (5.18). As described before, a Monte Carlo experiment with 1000 samples was used to verify the sampling distribution of the tests statistics, for dif- ferent combinations of the parameters. After, we used the parametric bootstrap method to resample B = 1000 times from Monte Carlo samples and the critical regions were defined. Now, we simulated by using the Monte Carlo method from the log-logistic distribution with location parameter φ = φˆ, the MLE under the null hypothesis and β = β∗ varying as shown in Figure 5.9. Finally, 1000 ψAD(Y ) and ψCvM (Y ) were calculated for each β ∗ and the rejection rate was used to compute the power of the test as presented in Figure 5.9. 0.5 1.0 1.5 0.2 0.4 0.6 0.8 1.0 β Po we r o f Te st Anderson Darling Cramér−von Misses Figure 5.9: Power of the tests Anderson-Darling and Cramér-von Misses for n = 100 and different values of β. 0 200 400 600 800 1000 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 x Cu m u la tiv e F un ct io n Log−Logistic Weibull e−Transmuted 0 10 20 30 40 50 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 x Cu m u la tiv e F un ct io n Log−Logistic Weibull e−Transmuted Figure 5.10: Cumulative distribution curves for the log-logistic and Weibull e-Transmuted models when β = 0.1 and 1.5, respectively; ω and φ were set in the same values then estimated in power of the test experiment. In Figure 5.10 we show the cumulative distribution curves for the models log-logistic and Weibull e-transmuted. It is important to note that the left panel, when β is close to zero, the curves are close which explain the low Anderson-Darling and Cramér-von Mises 85 E-TRANSMUTED FAMILY OF DISTRIBUTION power of tests. However, when β is more then 1, in left panel we set β = 1.3, the curves start having different shapes, which is consistent with the observed test power, close to 1. Furthermore, a second study of power of the test was carried out. For this study, we consider the following directional simple hypothesis:{ H0 : Y ∼Weibull(φ, β) H1 : Y ∼ e-Transmuted Weibull(ω, φ, β) ⇒ { H0 : ω = 0 H1 : ω 6= 0 In order to test these hypothesis, a Monte Carlo experiment with 1000 samples with size n = 100 was used to check the sampling distribution of the statistics Anderson-Darling and Cramér-von Mises, i.e. to obtain the points ψ1−α1AD (Y ) and ψ 1−α1 CvM (Y ). Then, for each sample, we used the bootstrap non-parametric method resampling B = 1000 times and, by using this estimates, we calculated the quadratic distances presented in equations (5.20) and (5.19). The results of this experiment are showed in Figure 5.11. 0 5 10 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ω Po w e r o f T e st Anderson Darling Cramér−von Misses β = 0.5 −4 −2 0 2 4 6 8 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 ω Po w e r o f T e st Anderson Darling Cramér−von Misses β = 1.5 Figure 5.11: Power of the tests Anderson-Darling and Cramér-von Mises with sample size n = 100 and for different combinations of the parameters β and ω and fixed φ = 2 when Y ∼Weibull(φ, β). It is important to note that, when the parameter ω is close to zero, the power of the test is about 5%. The power of the test increases when ω moves away from zero positively or negatively. Also, note that the test is more sensible for negative ω's than for positive ones since the curves are asymmetric. Chapter 6 Conclusions and perspectives In this research, we propose a new generalization of the log-logistic, the transmuted log-logistic distribution. The proposed distribution is constructed by using a quadratic rank transmutation map and taking the log-logistic distribution with two parameters as the baseline distribution. We also propose the use of a new generalization of the log-logistic distribution, the transmuted log-logistic distribution, in a Bayesian context and using a regression approach. We have provided closed expressions for several probabilistic measures including the probability density function, function hazard, moments, quantile function, mean, variance, and median. For each presented model, a simulation study was performed, from which we learned that the maximum likelihood estimate bias and the standard error decrease when the sample size is increased. On the other hand, the coverage probability of 95% two sided confidence intervals for the model parameters becomes closer to the nominal ones as the sample size increases. Furthermore, we proposed to analyze the behavior of the transmuted log-logistic model in the presence of censored observations. Note that this is the first time that this class of model (transmuted) was applied in the presence of censored data. In order to fit the parameters of the model, profile methods were used and compared to each other: Pure, Adjusted and Modified profile. All profile methods can be used to fit this class of model in the presence of censors, but the Adjusted Profile was better when compared to others. The applicability of the model was shown using three real datasets. The first one refers to a long term study involving 17026 cows of Tabapua breed studied by EMBRAPA. This dataset was used in order to show the usefulness of the proposed log-logistic transmuted model in a classical and Bayesian context, in the presence of covariates (regression ap- proach) and by considering the cubic log-logistic model. The second one refers to 148 patients treated with the drug Linezolid. This drug is 86 87 CONCLUSIONS AND PERSPECTIVES used to treat serious infections, therefore, all patients were treated with Linezolid at an ICU (Intensive Care Unit), in the city of Maringá. Furthermore, an analysis of influential points, residuals and a bootstrap study were made in order to improve the fitting quality of the model and validate our results. This dataset was used in order to show the usefulness of the model in the presence of right censor times. The last one corresponds to an uncensored study on the breaking stress of carbon fibers (in Gba), from Nichols and Padgett (2006) and was used to provide an example of applicability of the cubic Weibull transmuted model. In all of these applications, we observed that the most probable time of the event of interest to occur could not be fitted by an usual model in survival analysis. This happens because the risk function has a high peak at its point of mode, which does not occur with any usual model. In addition, to solve the problem of the restricted parametric space of λ parameter, we proposed a new family of distribution called e-transmuted or exponential transmuted family of models. Goodness-of-fit by using two quadratic distance measures were presented to validate the results to the e-transmuted regression Weibull model. In order to continue this research, we are interested in studying the mixture transmuted models and its behavior in the presence of a latent variable. In parallel, we are interested in generalizing the ranking of transmutation maps proposing, for example, the kth order of a ranking transmutation map. Appendix A Hessian matrix A.1 The hessian matrix of transmuted log-logistic model of order 1 The Hessian matrix is given by A =  A11, A12, A13A21, A22, A23 A31, A32, A33  where  Vˆ11, Vˆ12, Vˆ13Vˆ21, Vˆ22, Vˆ23 Vˆ31, Vˆ32, Vˆ33  =  A11, A12, A13A21, A22, A23 A31, A32, A33  −1 and A11 = − ∂ 2l ∂µ2 = − n∑ i=1 ( eµxβi − λeµxβi 1 + eµxβi − λeµxβi + λ )( 1− e µxβi − λeµxβi 1 + eµxβi − λeµxβi + λ ) +3 n∑ i=1 ( eµxβi 1 + eµxβi )( 1− e µxβi 1 + eµxβi ) 88 89 HESSIAN MATRIX A12 = A21 = − ∂ 2l ∂µ∂β = − n∑ i=1 eµxβi lnxi − λeµxβi lnxi 1 + eµxβi − λeµxβi + λ ( eµxβi lnxi − λeµxβi lnxi )( eµxβi − λeµxβi ) ( 1 + eµxβi − λeµxβi + λ )2  +3 n∑ i=1 eµxβi lnxi 1 + eµxβi ( 1− e µxβi 1 + eµxβi ) A13 = A31 = − ∂ 2l ∂µ∂λ = n∑ i=1 ( eµxβi 1 + eµxβi − λeµxβi + λ )( 1 + (1− λ)(1− eµxβi ) 1 + eµxβi − λeµxβi + λ ) A22 = − ∂ 2l ∂β2 = n β2 − n∑ i=1 eµxβi lnx2i − λeµxβi lnx2i 1 + eµxβi − λeµxβi + λ − ( eµxβi lnxi − λeµxβi lnxi 1 + eµxβi − λeµxβi + λ )2 +3 n∑ i=1 eµxβi lnx 2 i 1 + eµxβi ( 1− e µxβi 1 + eµxβi ) A23 = A32 = − ∂ 2l ∂β∂λ = n∑ i=1 ( eµxβi lnxi 1 + eµxβi − λeµxβi + λ + (1− eµxβi )(eµxβi lnxi − λeµxβi lnxi) (1 + eµxβi − λeµxβi + λ)2 ) A33 = − ∂ 2l ∂λ2 = x n∑ i=1 ( 1− eµxβi 1 + eµxβi − λeµxβi + λ )2 A.2 The hessian matrix of transmuted Weibull of order 2 The Hessian matrix is given by A =  A11, A12, A13, A14 A21, A22, A23, A24 A31, A32, A33, A34 A41, A42, A43, A44  90 HESSIAN MATRIX where  Vˆ11, Vˆ12, Vˆ13, Vˆ14 Vˆ21, Vˆ22, Vˆ23, Vˆ24 Vˆ31, Vˆ32, Vˆ33, Vˆ34 Vˆ41, Vˆ42, Vˆ43, Vˆ44  =  A11, A12, A13, A14 A21, A22, A23, A24 A31, A32, A33, A34 A41, A42, A43, A44  −1 and A11 = − n µ2 + 2(λ1 + 2λ2 − 3) n∑ i=1 ki ωi [Sµi logSi − logSi − 2(λ1 + 2λ2 − 3)ki] + −6(1− λ2) n∑ i=1 ki [Z µ i logSi + 2ki] A12 = A21 = −n β + 1 β n∑ i=1 [µSµi logSi + S µ i ] + + 2(λ1 + 2λ2 − 3) β n∑ i=1 [ ki ωi (µ− µSµi + Sµi Zi) + Sµi Zi ] + 6(1− λ2) β n∑ i=1 Sµi Zi [1 + µZi logSi − 2µSµi Zi logSi] A13 = A31 = 2 n∑ i=1 ki ωi [ (λ1 + 2λ2 − 3)(2Zi − 1) ωi − 1 ] A14 = A41 = 2 n∑ i=1 ki Si [ (λ1 + 2λ2 − 3)(4Zi − 1) ωi − 2 ] + 6 n∑ i=1 kiZi A22 = nµ β2 + 2µβ2(λ1 + 2λ2 − 3)S µ i Zi ωi [ µSµi − 2(λ1 + 2λ2 − 3) Sµi Zi ωi − µ− 1 ] + +6(1− λ2)µ β n∑ i=1 Sµi Z 2 i [2µS µ i − µ− 1]− µ(µ+ 1) n∑ i=1 Sµi β2 A23 = A32 = 2 µ β n∑ i=1 Sµi Zi ωi [ 1− 2(λ1 + 2λ2 − 3)(2Zi − 1) ωi ] 91 HESSIAN MATRIX A24 = A42 = 2 µ β n∑ i=1 Sµi Zi ωi [ 1− 2(λ1 + 2λ2 − 3)(4Zi − 1) ωi ] − 6µ β Sµi Z 2 i A33 = − n∑ i=1 (2Zi − 1)2 ω2i A34 = A43 = − n∑ i=1 (2Zi − 1)(1 + 4Zi) ω2i A44 = − n∑ i=1 (4Zi − 1)2 ω2i where Zi = e −(xi/β)µ , Si = xiβ and ki = ( xi β )µ log ( xi β ) e−(xi/β)µ . A.3 The hessian matrix of transmuted log-logistic of order 2 The Hessian matrix is given by A =  A11, A12, A13, A14 A21, A22, A23, A24 A31, A32, A33, A34 A41, A42, A43, A44  where  Vˆ11, Vˆ12, Vˆ13, Vˆ14 Vˆ21, Vˆ22, Vˆ23, Vˆ24 Vˆ31, Vˆ32, Vˆ33, Vˆ34 Vˆ41, Vˆ42, Vˆ43, Vˆ44  =  A11, A12, A13, A14 A21, A22, A23, A24 A31, A32, A33, A34 A41, A42, A43, A44  −1 and 92 HESSIAN MATRIX A11 = 4 n∑ i=1 ηi (1 + ηi) [ ηi (1 + ηi)− 1 ] + 2 n∑ i=1 ηi { λ2 + ηi(6− 2λ2 − 2λ1) η2i (3− λ1 − λ2) + 2ηi + λ1 + − (λ2 + ηi(3− λ2 − λ1)) 2[ η2i (3− λ1 − λ2) + 2ηi + λ1 ]2 } A12 = A21 = 4 n∑ i=1 log xi ηi (1 + ηi) [ ηi (1 + ηi)− 1 ] + +2 n∑ i=1 ηi log xi [ λ2 + ηi(6− 2λ2 − 2λ1) η2i (3− λ1 − λ2) + 2ηi + λ1 ] + − n∑ i=1 ηi (λ2 + ηi(3− λ2 − λ1))(1 + log xi)[ η2i (3− λ1 − λ2) + 2ηi + λ1 ]2 A13 = A31 = n∑ i=1 [ ηi(λ2 + ηi(3− λ2 − λ1))(ηi − 1)[ η2i (3− λ1 − λ2) + 2ηi + λ1 ]2 − 2η2iη2i (3− λ1 − λ2) + 2ηi + λ1 ] A14 = A41 = n∑ i=1 [ η2i (λ2 + ηi(3− λ2 − λ1))(ηi − 2)[ η2i (3− λ1 − λ2) + 2ηi + λ1 ]2 − 2η2i (1− ηi)η2i (3− λ1 − λ2) + 2ηi + λ1 ] A22 = 4 n∑ i=1 (log xi) 2 ηi (1 + ηi) [ ηi (1 + ηi)− 1 ] − n β + +2 n∑ i=1 ηi { (log xi) 2 [λ2 + ηi(6− 2λ2 − 2λ1)] η2i (3− λ1 − λ2) + 2ηi + λ1 + − log xi(λ2 + ηi(3− λ2 − λ1)) 2[ η2i (3− λ1 − λ2) + 2ηi + λ1 ]2 } A23 = A32 = − n∑ i=1 { 2ηi log xi(λ2 − λ1ηi − λ2ηi + 3ηi)(1− η2i )[ η2i (3− λ1 − λ2) + 2ηi + λ1 ]2 + + 2η2i log xi η2i (3− λ1 − λ2) + 2ηi + λ1 } 93 HESSIAN MATRIX A24 = A42 = − n∑ i=1 { 2ηi log xi(λ2 − λ1ηi − λ2ηi + 3ηi)[ηi(2− ηi)]2[ η2i (3− λ1 − λ2) + 2ηi + λ1 ]2 + − 2ηi log xi(1− ηi) η2i (3− λ1 − λ2) + 2ηi + λ1 } A33 = − n∑ i=1 (1− η2i )2[ η2i (3− λ1 − λ2) + 2ηi + λ1 ]2 A34 = A43 = − n∑ i=1 (1− η2i )(ηi(2− ηi))2[ η2i (3− λ1 − λ2) + 2ηi + λ1 ]2 A44 = − n∑ i=1 (ηi(2− ηi))2[ η2i (3− λ1 − λ2) + 2ηi + λ1 ]2 where ηi = e µxβi . A.4 Hessian matrix of e-transmuted model In order to present the hessian matrix, we consider the model e-transmuted without covariates as follows:  I11 . . . I1(p+k) . . . ... I(p+k)(p+k)  (p+k)×(p+k) (A.1) with p the number of covariates and k the number of parameters of base distribution and the elements I as follow: I11 = ∂2l ∂ω2 = n [ ω2e−ω − e−2ω + 2e−ω − 1] ω2 (1− e−ω)2 I1j = ∂2l ∂ω∂θj = − n∑ i=1 G′θj (yi|θ) Ikj = ∂2l ∂θj∂θk = − n∑ i=1 g′′θjθk(yi|θ)g(yi|θ)− g′θj (yi|θ)g′θk(yi|θ) g2(yi|θ) − ω n∑ i=1 G′′θjθk(yi|θ) 94 HESSIAN MATRIX Ijj = ∂2l ∂θ2j = − n∑ i=1 [ g′θj (yi|θ) ]2 g2(yi|θ) + n∑ i=1 g′′θj (yi|θ) g(yi|θ) + ω n∑ i=1 G′′θj (yi|θ). Appendix B Tables Table B.1: MLEs for the original sample generated with different parameter values and sample sizes. Parameters Sample Estimated Estimated Generated Size Relative Difference(%) Standard Error (γ1, γ2, β, λ) n γ1 γ2 β λ γ1 γ2 β λ (−2,−1, 0.5,−0.8) 50 11.136 20.499 3.644 12.094 0.707 0.634 0.070 0.527 100 7.479 14.098 2.158 7.631 0.518 0.463 0.053 0.372 150 6.029 10.908 1.175 5.175 0.431 0.379 0.044 0.291 300 4.068 6.587 0.747 3.721 0.320 0.278 0.034 0.197 500 2.430 3.919 0.270 1.781 0.248 0.213 0.027 0.142 1000 1.244 2.397 0.228 0.987 0.173 0.149 0.019 0.096 (−2,−0.5, 1.5,−0.5) 50 6.667 29.552 0.051 15.737 0.727 0.643 0.210 0.598 100 4.237 21.604 0.860 2.789 0.535 0.466 0.155 0.437 150 4.479 19.244 0.609 4.435 0.465 0.400 0.134 0.386 300 3.155 10.868 0.466 1.911 0.347 0.294 0.098 0.291 500 2.710 8.546 0.202 4.239 0.282 0.237 0.079 0.240 1000 1.516 5.446 0.133 3.468 0.204 0.169 0.057 0.173 (−2, 0.5, 1,−0.5) 50 6.570 27.020 0.015 14.991 0.724 0.623 0.140 0.595 100 4.188 21.382 0.880 3.046 0.534 0.447 0.103 0.438 150 4.341 18.216 0.657 3.862 0.464 0.377 0.089 0.386 300 3.155 10.491 0.466 1.911 0.347 0.274 0.065 0.291 500 2.710 7.654 0.202 4.239 0.282 0.217 0.053 0.240 1000 1.516 4.722 0.133 3.468 0.204 0.152 0.038 0.173 (−3,−1, 0.5, 0.8) 50 2.204 0.718 3.115 20.596 0.677 0.592 0.069 0.491 100 1.009 0.442 1.984 14.506 0.471 0.416 0.052 0.364 150 1.031 0.397 1.570 10.668 0.378 0.333 0.046 0.294 300 0.829 0.079 0.646 6.547 0.263 0.231 0.033 0.197 500 0.730 0.094 0.284 4.466 0.201 0.176 0.026 0.147 1000 0.423 0.430 0.087 2.725 0.141 0.123 0.019 0.098 (−2,−0.5, 1.5, 0.5) 50 2.453 1.889 0.544 27.549 0.653 0.618 0.207 0.538 100 2.413 1.101 0.610 30.890 0.473 0.455 0.156 0.454 150 1.615 2.003 0.328 29.270 0.386 0.379 0.134 0.403 300 1.139 2.257 0.479 19.601 0.272 0.273 0.098 0.310 95 96 TABLES 500 0.300 1.557 0.330 13.072 0.207 0.213 0.080 0.248 1000 0.056 0.477 0.181 7.248 0.144 0.152 0.058 0.179 (−3, 0.5, 1, 0.5) 50 3.048 0.123 0.512 28.512 0.713 0.631 0.138 0.536 100 1.663 1.233 0.596 31.439 0.514 0.471 0.104 0.455 150 1.112 1.809 0.349 29.053 0.417 0.400 0.089 0.403 300 0.645 0.515 0.479 19.601 0.292 0.293 0.066 0.310 500 0.137 0.336 0.342 12.917 0.220 0.233 0.053 0.248 1000 0.015 0.247 0.181 7.248 0.152 0.168 0.039 0.179 Table B.2: Probability of coverage by considering 95% of confidence. Generated Sample Coverage Parameters Size Probability (γ1, γ2, β, λ) n γ1 γ2 β λ (−2,−1, 0.5,−0.8) 50 0.916 0.914 0.941 0.498 100 0.923 0.927 0.934 0.640 150 0.911 0.925 0.911 0.721 300 0.925 0.944 0.915 0.849 500 0.938 0.940 0.930 0.888 1000 0.951 0.947 0.954 0.930 (−2,−0.5, 1.5,−0.5) 50 0.914 0.909 0.937 0.634 100 0.920 0.908 0.919 0.716 150 0.900 0.913 0.922 0.792 300 0.899 0.931 0.924 0.838 500 0.916 0.937 0.932 0.889 1000 0.940 0.952 0.930 0.921 (−2, 0.5, 1,−0.5) 50 0.913 0.919 0.932 0.626 100 0.921 0.918 0.916 0.715 150 0.902 0.922 0.920 0.790 300 0.899 0.939 0.924 0.838 500 0.916 0.950 0.932 0.889 1000 0.940 0.955 0.930 0.921 (−3,−1, 0.5, 0.8) 50 0.933 0.934 0.933 0.494 100 0.952 0.942 0.924 0.674 150 0.945 0.959 0.925 0.763 300 0.952 0.956 0.923 0.841 500 0.940 0.950 0.925 0.894 1000 0.947 0.957 0.930 0.916 (−2,−0.5, 1.5, 0.5) 50 0.918 0.915 0.918 0.625 100 0.931 0.917 0.935 0.742 150 0.932 0.927 0.934 0.788 97 TABLES 300 0.955 0.937 0.925 0.843 500 0.952 0.941 0.929 0.894 1000 0.944 0.951 0.932 0.929 (−3, 0.5, 1, 0.5) 50 0.917 0.922 0.919 0.625 100 0.937 0.919 0.934 0.743 150 0.938 0.915 0.933 0.786 300 0.953 0.913 0.925 0.843 500 0.948 0.933 0.929 0.894 1000 0.947 0.955 0.932 0.929 References Akaike, H. (1973). Information theory and the maximum likelihood principle. International Simposium on Information Theory, eds. V. Petrov and F. 15, 59 Anaya-Izquierdo, K. and Marriott, P. (2007). Local mixture models of exponential families. Bernoulli , pages 623640. 72 Anderson, T. W. (1962). On the distribution of the two-sample cramer-von mises criterion. The Annals of Mathematical Statistics, pages 11481159. 82 Arnold, B. C., Balakrishnan, N., and Nagaraja, H. N. (1992). A first course in order statistics, volume 54. Siam. 6 Aryal, G. R. and Tsokos, C. P. (2009). On the transmuted extreme value distribution with application. Nonlinear Analysis: Theory, Methods & Applications, 71(12), e1401e1407. 1, 6 Aryal, G. R. and Tsokos, C. P. (2011). Transmuted weibull distribution: A generalization of theweibull probability distribution. European Journal of Pure and Applied Mathematics, 4(2), 89102. 1, 12, 14, 52 Barlow, R. E. and Campo, R. A. (1975). Total time on test processes and applications to failure data analysis. Technical report, DTIC Document. 17, 19 Barlow, W. E. and Prentice, R. L. (1988). Residuals for relative risk regression. Biometrika, 75(1), 6574. 41, 62 Barndorfe-Nielsen, O. E. and Cox, D. R. (1994). Inference and asymptotics. Chapman & Hall. 38 Barndorff-Nielsen, O. E. (1993). On a formula for the distribution of the maximum likeli- hood estimator. Biometrika. 41 Barndorff-Nielsen, O. E. and McCullagh, P. (1993). A note on the relation between modified profile likelihood and the cox-reid adjusted profile likelihood. Biometrika, 80(2), 321 328. 2, 40 Bennett, S. (1983). Log-logistic regression models for survival data. Applied Statistics, pages 165171. 2 Brickner, S. J. (1996). Oxazolidinone antibacterial agents. Current Pharmaceutical Design 2 , 2, 175194. 43 98 99 REFERENCES Burnham, K. P. and Anderson, D. R. (2002). Model selection and multimodel inference: a practical information- theoretic approach, Second Edition. Springer Series in Statistics. 20 Chechile, R. A. (2003). Mathematical tools for hazard function analysis. Journal of Math- ematical Psychology , 47, 478494. 11 Chen, M. h. and Ibrahim, J. G. (2006). The relationship between the power prior and hierarchical models. Bayesian Anal., 1(3), 551574. 24 Cox, D. and Reid, N. (1987). Parameter orthogonality and approximate conditional infer- ence. Journal of the Royal Statistical Society B . 40 Cox, D. R. and Snell, E. J. (1968). A general definition of residuals. Journal of the Royal Statistical Society , 30, 248275. 21 Dahiya, R. C. and Hossain, S. A. (1996). A modification of goel-okumoto model. pages 7784. 70 Daniels, M. J. (1999). A prior for the variance in hierarchical models. 24 Duchesne, T., Rioux, J., and Luong, A. (1997). Minimum cramér-von mises distance meth- ods for complete and grouped data. Communications in statistics-theory and methods, 26(2), 401420. 82 Fachini, B. J., Ortega, E. M. M., and Louzada, F. (2008). Influence diagnostics for poly- hazard models in the presence of covariates. Stat Meth Appl , 17, 413433. 2 Ferrari, S. L. P., Silva, M. F., and Cribari-Neto, F. (2007). Adjusted profile likelihoods for the Weibull shape parameter. Journal of Statistical Computation and Simulation. 2 Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models. Communications in Statistics. Theory and Methods, 1(3), 515533. 2, 24, 25 Gelman, A. and Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, (7), 457511. 26 Gentle, J. (2009). Computational Statistics. Statistics and Computing. Springer New York. 78 Ghitany, M. E. (2001). A compound rayleigh survival model and its application to randomly censored data. Statistical Paper, Springer , 42, 437450. 1 Ibrahim, J. G., Chen, M.-H., and Sinha, D. (2005). Bayesian survival analysis. Wiley Online Library. 1, 6 Johnson, N. L., Kotz, S., and Balakrishnan, N. (1996). Continuous univariate distributions, volume 2. Wiley New York. 52 Lawless, J. F. (2011). Statistical models and methods for lifetime data, volume 362. John Wiley & Sons. 2, 10, 18, 38, 52 100 REFERENCES Lee, S. Y., Lu, B., and Song, X. Y. (2006). Assessing local influence for nonlinear structural equation models with ignorable missing data. Computational Statistics & Data Analysis. 31, 48 Marshall, A. and Olkin, I. (2007a). Life distributions: structure of nonparametric, semi- parametric and parametric families. Springer, New York. 1 Marshall, A. and Olkin, I. (2007b). Life Distributions: Structure of Nonparametric, Semi- parametric, and Parametric Families. Springer Series in Statistics. Springer New York. 70 McCullagh, P. and Nelder, J. A. (1989). Generalized Linear Models. McGraw Hill, London. 41, 62 Migon, H. S., Gamerman, D., and Louzada, F. (2014). Statistical inference: an integrated approach. CRC press. 58 Mood, A. M., Graybill, F. A., and Boes, D. C. (1974). Introduction to the Theory of Statistics. McGraw Hill, 3 edition. 7 Nichols, M. and Padgett, W. (2006). A Bootstrap control chart for Weibull percentiles. Quality and Reliability Engineering International , 22, 141151. 52, 59, 87 Ortega, E. M., Bolfarine, H., and Paula, G. A. (2003). Influence diagnostics in generalized log-gamma regression models. Computational Statistics and Data Analysis, 42, 165186. 2 Paro, P. a. Z., Santos, A. L. Q., Maximiniano-Neto, A., Paro, J. L. N., Rodrigues, D. C., Cruz, G. C., Malta, T. S., Ribeiro, F. M., and Andrade, M. A. (2013). Anatomic study of the vascular casts of the testicular arteries in bovines of Tabapua race. Bioscience Journal , 19(1), 123132. 16, 17 Pereira, J. C. C. (2000). Contribuição genética do Zebu na pecuária bovina do Brasil. Informe Agropecuário, 21, 3038. 16, 17 Pham, H. and Lai, C.-D. (2007). On Recent Generalizations of the Weibull Distribution. IEEE Transactions on RELIABILITY , 56(3), 454458. 52 Polson, N. G. and Scott, J. G. (2012). On the half-Cauchy prior for a global scale parameter. Bayesian Anal., 7(4), 887902. 2, 24 Rayner, J. C., Thas, O., and Best, D. J. (2009). Smooth tests of goodness of fit: using R. John Wiley & Sons. 71 Ross, S. M. (2009). A First Course in Probability . Prentice Hall. 9 Severini, T. (1998). An approximation to the modified profile likelihood function. Biometrika. 2, 41 Shaw, W. T. and Buckley, I. R. C. (2009). The alchemy of probability distributions: beyond gram-charlier expansions, and a skew-kurtotic-normal distribution from a rank transmutation map. arXiv preprint arXiv:0901.0434 . 1, 5, 6, 53, 71 101 REFERENCES Therneau, T. M., Grambsch, P. M., and Fleming, T. R. (1990). Martingale-based residuals for survival models. Biometrika. 41, 42, 62 Tollenaar, N. and Mooijaart, A. (2003). Type i errors and power of the parametric bootstrap goodness-of-fit test: Full and limited information. British Journal of Mathematical and Statistical Psychology , 56(2), 271288. 80 Venzon, D. J. and Moolgavkar, S. H. (1988). A Method for Computing Profile-Likelihood Based Confidence Intervals. Applied Statistics. 38 Verdinelli, I., Wasserman, L., et al. (1998). Bayesian goodness-of-fit testing using infinite- dimensional exponential families. The Annals of Statistics, 26(4), 12151241. 69 Wilson, A. P. R., Cepeda, J. A., Hayman, S., Whitehouse, T., Singer, M., and Bellingan, G. (2006). In vitro susceptibility of Gram-positive pathogens to linezolid and teicoplanin and effect on outcome in critically ill patients. Journal of Antimicrobial Chemotherapy . 43