Regresión lineal con datos censurados

Traballo Fin de Grao REGRESIÓN LINEAL CON DATOS CENSURADOS datos DATOS CENSURADOS María Barreira Miranda 2019/2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Regresión lineal con datos censurados María Barreira Miranda Xullo 2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Traballo proposto Área de Coñecemento: Estatística e Investigación Operativa Título: Regresión lineal con datos censurados Breve descrición do contido Os datos censurados son moi habituais na Análise de Supervivencia, que é a parte da Estatística que estuda os tempos de vida. E que os tempos de vida, que poden ser duracións dunha enfermidade, dun artigo de consumo (coches, teléfonos, ordenadores, etc.) ou calquera outro tempo entre dous eventos, normalmente requiren de certo seguimento. Se ese seguimento se interrompe, so coñeceremos que o tempo durou polo menos ata o momento da perda do seguimento. Nestas condicións pode seguir interesando considerar o efecto dalgunha variable sobre o tempo de vida. Por exemplo, pode interesar saber se a idade do ou da doente inúe no tempo de curación dunha lesión. Este traballo consiste en revisar as técnicas de estimación da regresión lineal cando a variable resposta está censurada. Exporanse os métodos xa existentes, estudaranse as súas propiedades mediante simulacións, e ilustraranse con datos reais. Recomendacións Ter un coñecemento básico do programa estatístico . iii Índice xeral Resumo viii ix 1. Introdución 1 1.1. Hipóteses do modelo de regresión lineal simple . . . . . . . . . . . . . . . . . 2 1.2. Estimación dos parámetros . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3. Propiedades dos estimadores . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4. Regresión lineal múltiple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2. Datos censurados 7 2.1. Introdución á Análise de Supervivencia . . . . . . . . . . . . . . . . . . . . . 9 2.2. Tiposdecensura ................................. 11 2.3. Funcións que caracterizan unha variable censurada . . . . . . . . . . . . . . 12 2.3.1. Función de Supervivencia . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.2. Funciónderisco ............................. 13 2.4. Medidas características dunha variable censurada . . . . . . . . . . . . . . . 13 2.5. Estimador de Kaplan-Meier . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3. Regresión censurada 17 3.1. O estimador de mínimos cadrados . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2. O estimador proposto por Miller . . . . . . . . . . . . . . . . . . . . . . . . 19 v vi ÍNDICE XERAL 3.3. O estimador proposto por Buckley e James . . . . . . . . . . . . . . . . . . 23 3.4. O estimador proposto por Jin, Lin e Ying . . . . . . . . . . . . . . . . . . . 25 4. Estudo de simulación 29 4.1. Introdución .................................... 29 4.2. Modeloconintercepto .............................. 31 4.2.1. Erro con distribución normal . . . . . . . . . . . . . . . . . . . . . . 32 4.2.2. Erro con distribución chi-cadrado . . . . . . . . . . . . . . . . . . . . 37 4.3. Modelosenintercepto .............................. 40 4.3.1. Erro con distribución normal . . . . . . . . . . . . . . . . . . . . . . 40 4.3.2. Erro con distribución chi-cadrado . . . . . . . . . . . . . . . . . . . . 45 5. Aplicación a datos reais 49 5.1. BasededatosUIS ................................ 49 5.2. Análise descritiva previa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.3. Estimación dun modelo de regresión . . . . . . . . . . . . . . . . . . . . . . 52 6. Conclusións 63 Anexo A: Comandos de R 65 Bibliografía 83 2 CAPÍTULO 1. INTRODUCIÓN 1.1. Hipóteses do modelo de regresión lineal simple Para poder estimar o modelo (1.1) necesitamos as seguintes hipóteses: Linealidade. Debido a que a función de regresión é unha liña recta podemos expresar este modelo como aparece na ecuación (1.1). Homocedasticidade. Para cada valor da variable explicativa x , a varianza do erro debe ser constante, é dicir: Var (ε|X=x) = σ2 para todo x . Independencia. Os erros teñen que ser independentes entre si. Normalidade. O erro ten distribución normal de media cero e varianza σ2 , é dicir: ε∈N0, σ2 1 . Para poder estimar β0 e β1 necesitamos unha mostra de datos que se obteñen ou ben dun deseño xo ou dun deseño aleatorio. No deseño xo fíxanse os valores da variable explicativa de forma que temos unha mostra do tipo {(x1, Y1), ... (xn, Yn)} . En canto ao deseño aleatorio, tanto a variable explicativa como a variable resposta son aleatorias, o que nos da unha mostra do tipo {(X1, Y1), ... (Xn, Yn)} . Para poder facer inferencia sobre os parámetros desexados usaremos un deseño de tipo xo. Neste traballo consideraremos que {W1, ... , Wn} son os datos e W(1), ... , W (p) son as correspondentes variables aleatorias. Resumindo, traballaremos cun modelo de regresión lineal simple, homocedástico, con errores normais e independentes e baixo deseño xo. 1.2. Estimación dos parámetros Para estimar os parámetros β0 , β1 e σ2 asociados a un modelo de regresión lineal simple supoñeremos as hipóteses antes explicadas. Para a predición do valor da variable resposta Y a partir do valor da variable explicativa X , temos os seguintes erros, chamados residuos da regresión: bεi=Yi−b β0−b β1xi sendo i∈ {1, ... , n}. 1 Cando escribimos Nµ, σ2 referímonos a unha normal de media µ e varianza σ2 , sendo a normal unha das distribucións de probabilidade máis frecuentes en Estatística. 1.3. PROPIEDADES DOS ESTIMADORES 3 Para efectuar a estimación faremos uso do método de mínimos cadrados cuxo obxectivo é minimizar a suma de residuos ao cadrado. Así, temos que buscar b β0 e b β1 de forma que fagan mínima esta suma, que vén dada por: n X i=1 Yi−b β0−b β1xi2= m´ın β0,β1 n X i=1 (Yi−β0−β1xi)2. Se derivamos a expresión anterior con respecto a β0 e β1 e logo igualamos a cero, obtemos as seguintes expresións: b β0=Y−SxY S2 x xb β1=SxY S2 x Por último estimaremos a varianza do erro σ2 do seguinte xeito: bσ2=1 n−2 n P i=1 bε2 i=1 n−2 n P i=1 Yi−b β0−b β1xi2 1.3. Propiedades dos estimadores Unha vez coñecidas as expresións dos estimadores, vainos ser interesante estudar as súas propiedades. Propiedades de c β1 Para calcular a esperanza de forma máis sinxela, expresaremos b β1 da seguinte forma: b β1=SxY S2 x = n P i=1 (xi−x)Yi−Y nS2 x = n P i=1 (xi−x) nS2 x (Yi−Y) = n P i=1 ωiYi−Y sendo ωi=(xi−x) nS2 x os pesos que so dependen da variable explicativa e, como estamos supoñendo que traballamos cun deseño xo, estes pesos non son aleatorios. Así, usando propiedades da media, podemos facer: Eb β1=En P i=1 ωiYi−Yi= n P i=1 ωiEYi−Y= n P i=1 (xi−x) nS2 x | {z } ωi β1(xi−x) | {z } E(Yi−Y) =β1 Para demostrar as últimas igualdades temos que ter en conta que: E(Yi) = β0+β1xi. 4 CAPÍTULO 1. INTRODUCIÓN EY=E1 n n P i=1 Yi=1 n n P i=1 E(Yi) = 1 n n P i=1 (β0+β1xi) = β0+β1x. En consecuencia temos EYi−Y=β0+β1xi−β0−β1x=β1(xi−x). Na última igualdade úsase que S2 x=1 n n P i=1 (xi−x)2 . Para poder calcular a varianza, imos expresar b β1 da seguinte forma: b β1= n P i=1 ωiYi−Y= n P i=1 ωiYi. debido a que n P i=1 ωi= 0 . Entón, Var b β1= Var n P i=1 ωiYi= n P i=1 ω2 i Var (Yi) = n P i=1 (xi−¯x)2 n2S4 x σ2=σ2 nS2 x . Para explicar estas igualdades recordamos que estamos traballando coa hipótese de independencia dos erros e homoceasticidade e usamos propiedades básicas do operador varianza. Finalmente, como β1 é combinación lineal de Y1, ... , Yn que son variables independentes e normais, este estimador ten distribución normal. En resumo: b β1∈Nβ1,σ2 nS2 x Propiedades de b β0 Calcularemos a media e a varianza de forma análoga ao caso anterior. Debido a que b β0=Y−b β1x , temos que: Eb β0=EY−xEb β1=β0+β1x−xβ1=β0. Para calcular agora a varianza, usaremos que Y=1 n n P i=1 Yi , xb β1= n P i=1 xωiYi e así temos: b β0=1 n n P i=1 Yi− n P i=1 xωiYi= n P i=1 1 n−xωiYi. Así, tendo en conta as hipóteses básicas do modelo de regresión lineal simple, verifícase que Var b β0= n X i=1 1 n−xωi2 Var (Yi) = σ2 n X i=1 1 n2+x2ω2 i−2¯xωi n=σ21 n+¯x2 nS2 x, 1.4. REGRESIÓN LINEAL MÚLTIPLE 5 onde para a última igualdade usamos que n P i=1 ωi= 0 e n P i=1 ω2 i=1 nS2 x . Finalmente, debido a que b β0 é combinación lineal de Y1, ... , Yn , tal como sucedía no caso de b β1 , temos que b β0 ten unha distribución normal. Polo tanto, podemos concluír: b β0∈Nβ0, σ21 n+x2 nS2 x Propiedades de bσ2 Aínda que neste caso non entraremos en detalles, o estimador da varianza do erro segue unha distribución do tipo chi-cadrado: (n−2) bσ2 σ2∈χ2 n−2 Nótese que cando estimamos a varianza, dividimos entre n−2 en lugar de facelo entre n para que o estimador σ2 sexa insesgado. 1.4. Regresión lineal múltiple Unha vez visto o modelo de regresión lineal simple, podemos estendelo a situacións máis complexas onde hai máis dunha variable explicativa, o que se coñece como regresión lineal múltiple. Neste modelo temos unha variable resposta Y e unha colección de variables explicativas X(1), ... , X(p−1) . Esta clase de modelos podemos escribilos en forma matricial do seguinte xeito:     Y1 . . . Yn    =    1x1,1··· x1,p−1 . . .. . .. . . 1xn,1··· xn,p−1        β0 . . . βp−1    +    ε1 . . . εn    . A expresión simplicada deste modelo é: Y=Xβ+ε, onde Y é o vector das variables respostas, X é unha matriz de n×p elementos, β∈Rp é o vector de coecientes e ε é o vector dos erros que verica ε∈Nn0, σ2In , sendo In a matriz identidade. Para estimar β , aplicaremos o método de mínimos cadrados que ao igual que no caso do modelo lineal simple ten como obxectivo minimizar a suma de residuos ao cadrado, é dicir: b β= arg m´ın β n X i=1 (Yi−xiβ)2, 6 CAPÍTULO 1. INTRODUCIÓN sendo xi a la i-ésima da matriz X . Este problema tamén se pode escribir en notación matricial de maneira equivalente como: b β= arg m´ın β(Y−Xβ)0(Y−Xβ). Se derivamos con respecto β e igualamos a cero, obtemos: X0Xβ=X0Y, e a súa solución é o estimador de β , que ven dado por: b β= (X0X)−1X0Y. Non escribiremos a demostración detallada pero é importante destacar as seguintes propiedades deste estimador: Media: Eb β=β , é dicir, trátase dun estimador inesgado. Covarianza: Cov b β, b β=σ2(X0X)−1. Distribución límite: b β∈Npβ, σ2(X0X)−1 Capítulo 2 Datos censurados Neste capítulo explicaremos que son os datos censurados moi empregados no ámbito da Análise de Supervivencia . Comezaremos dando unha idea sobre que é a censura e os datos censurados grazas a un exemplo motivador para logo explicar estes termos con máis profundidade. No ámbito da Bioestatística 1 podemos atopar moitos exemplos de datos censurados. Imaxinemos que traballamos cun estudo sobre o cancro. Neste caso os/as doentes poden morrer ou abandonar o estudo por diferentes causas e se isto ocorre non imos ter coñecemen- to dos seus datos completos. Isto é o que se coñece como datos censurados. Empregaremos como exemplo un ensaio clínico, que é un procedemento experimental dun medicamento ou tratamento en persoas para avaliar a súa seguridade ou a súa ecacia. Vexamos un exemplo concreto dun ensaio clínico que podemos atopar en [12]. Na Figura 2.1 representamos o tempo de vida de seis doentes e observamos que entran ao estudo durante un período de 2.5 anos que vai dende comezos do ano 2000 ata mediados do 2002 . Os/As individuos/as son seguidos durante 4.5 anos ata nais do 2007 , cando remata o ensaio clínico. Na Figura 2.1, as rectas verticais representan o comezo do ensaio, o peche do período para poder entrar ao estudo e o seu remate, respectivamente. Cada liña horizontal representa a un/unha individuo/a onde os puntos negros marcan a súa entrada ao estudo, os círculos representan os eventos censurados e as aspas denotan a morte do/da doente. Neste exemplo observamos un caso de censura xa que no momento que remata o estudo tres doentes (D1, D3 e D4) seguían vivos/as. Polo tanto, para estes/as tres doentes 1 A Bioestatística é unha rama da Estadística aplicada ás Ciencias da Vida como son a Bioloxía e a Medicina. 7 8 CAPÍTULO 2. DATOS CENSURADOS Figura 2.1: Ensaio clínico sobre seis doentes onde as rectas verticais representan o comezo do ensaio, o peche do período para poder entrar ao estudo e o seu remate, respectivamente. Cada liña horizontal representa a un/unha individuo/a onde os puntos negros marcan a súa entrada ao estudo, os círculos representan os eventos censurados e as aspas denotan a morte do/da doente. non temos información completa sobre o seu tempo de vida e polo tanto estes datos serán censurados. No caso do/da doente D1, por exemplo, sabemos que sobreviviu polo menos 7 anos pero non sabemos canto tempo máis vai vivir. Na Figura 2.1 temos os datos da súa morte (representada cunha aspa) pero este valor non sería coñecido no momento do estudo. En resumo, temos información completa de tres doentes dende que comeza o ensaio ata que morren e temos outros/as tres doentes que serán censurados. Nun estudo nestas circunstancias, no que non se coñecen todos os datos, poderíamos pensar en traballar so cos datos que temos e obviar o resto. Imos ver gracamente que ocorrería nesta situación empregando a función de distribución empírica. Dada unha mostra {Y1, ... , Yn} , recordemos que a función de distribución empírica dunha variable aleatoria Y é da seguinte forma: b Fn(y) = 1 n n X i=1 I(Yi≤y). (2.1) Na Figura 2.2 observamos dúas grácas. A parte (a) representa a función de distribución dunha variable normal N(0,1) para unha mostra de tamaño n= 100 , mentres que a parte (b) é para tamaño n= 1000 . O trazo en azul correspóndese coa función de distribución empírica para os datos cunha censura para valores menores que 1 (quedamos cos datos que cumpran y≤ 1 ) mentres que o vermello representa a función de distribución 2.1. INTRODUCIÓN Á ANÁLISE DE SUPERVIVENCIA 9 empírica dos datos completos (é dicir, sen presenza de censura). A curva negra representa a función de distribución teórica dunha N(0,1) e observamos como a curva vermella aproxima ben ao modelo teórico mentres que a azul non. Desta forma podemos ver que no caso dos datos censurados non podemos quedarnos cos datos que temos ata un certo tempo t xa que eses datos non se corresponde coa realidade. O problema tampouco é do tamaño da mostra xa que observamos a mesma situación no gráco da parte (a) e no da parte (b). Polo tanto, xorde a necesidade de buscar outras opcións para traballar con estes datos. Por exemplo, en lugar de empregar a función de distribución empírica, no caso de datos censurados usaremos a función Kaplan-Meier da que falaremos neste mesmo capítulo na Sección 2.5. Agora deniremos a censura e os datos censurados formalmente. Denición 2.1. Denomínase censura ao fenómeno que afecta ás variables de interese nun estudo cando existe unha limitación na información que temos delas. Os datos censurados son observacións que non poden ser cuanticadas, dado que so se coñece que o seu valor se atopa por debaixo ou por arriba dunha cota determinada ou ben que está incluído nun intervalo. 2.1. Introdución á Análise de Supervivencia Supoñamos que estamos interesados en estudar unha determinada variable Y que representa o tempo que pasa dende o comezo do experimento ata que ocorre un determinado suceso de interese que chamaremos morte ou fracaso. Así, cando nos reramos á variable no tempo t , estamos falando da variable Y . Este suceso tamén se pode entender como algo positivo, como por exemplo o tempo que pasa dende que un/unha doente entra nun ensaio clínico ata que responde favorablemente a un tratamento. O conxunto de técnicas estatísticas que se empregan para analizar este tipo de datos coñécense como a Análise de Supervivencia. Debido á presenza de censura, a Análise de Supervivencia tamén se coñece como análise de datos censurados. A principal característica da Análise de Supervivencia é que a variable Y é unha variable discreta ou continua non negativa e representa o tempo dende un inicio ata unha n denidos. Outra característica importante xurde cando o comezo ou a n dun evento non se observan completamente. Estes fenómenos coñécense como censura pola dereita e censura pola esquerda. Censura pola dereita: aparece cando o extremo nal é so coñecido por exceder un 10 CAPÍTULO 2. DATOS CENSURADOS (a) Tamaño: n= 100 (b) Tamaño: n= 1000 Figura 2.2: Representación gráca da estimación da función de distribución dunha variable normal no caso de datos censurados a partir de y= 1 e no caso de datos completos. 2.2. TIPOS DE CENSURA 11 valor particular. Formalmente, sexa Y unha variable que representa o tempo ata o fracaso e C unha variable que representa o tempo para un evento censurado. Entón denimos Z= m´ın (Y, C) e δ=I[Y⩽C] , que vai ser o indicador da censura. Desta forma, δ= 0 se Z é un tempo censurado e δ= 1 se Z é o tempo non censurado que se observa completamente. A Figura 2.1 mostra un exemplo de censura pola dereita. Censura pola esquerda: é menos frecuente e dáse cando os eventos teñen lugar antes do punto de inicio do estudo. Desta forma, o tempo de censura será o tempo de inicio do período de seguimento e non se coñece con exactitude. Ao longo deste traballo empregaremos datos censurados pola dereita, polo cal asumiremos que Z= m´ın (Y, C) denota a variable de interese observada e δ=I[Y⩽C] será a indicadora da censura. O obxectivo da Análise de Supervivencia é estimar a función de distribución (que veremos na Sección 2.3), comparar dúas ou máis distribucións de supervivencia e avaliar os efectos de certos factores sobre a variable Y . As técnicas que estudaremos teñen importantes similitudes coa clásica regresión lineal en media, coa importante diferenza de que a variable resposta é unha variable censurada. 2.2. Tipos de censura A partir de agora, como xa dixemos antes, cando falemos de censura estarémonos referindo a censura pola dereita. Esta censura divídese nos seguintes tipos, que explicaremos a partir de exemplos: Tipo I: Censura por tempo O tempo de censura está pre-denido. Nun estudo para deixar de fumar, séguese a cada doente dende que comeza o ensaio ata que sofre unha recaída, volve fumar, ou ben se despois de 180 días non se produce unha recaída. Os/As doentes que aos 180 días seguen sen fumar son censurados. Este exemplo podémolo ver en [12]. Tipo II: Censura por número de fallos Este caso ocorre cando os obxectos experimentais son seguidos ata que unha proporción deles fracasan. Dáse nas áreas da Biomedicina ou no ámbito da industria, onde o tempo de fallo dun dispositivo é o que máis interesa. Un exemplo deste tipo de censura pode darse cando no ámbito da industria realizamos un estudo que remata 18 CAPÍTULO 3. REGRESIÓN CENSURADA sendo εi os erros para cada i∈ {1, ... , n} . Desta forma, estimamos βdc =βdc 0, βdc 1 1 do seguinte xeito: b βdc = arg m´ın β 1 n n X i=1 bεi2= arg m´ın β 1 n n X i=1 Yi−βdc 0−βdc 1xi2. (3.1) O noso obxectivo é ser capaces de estimar o parámetro β no caso de que a variable resposta sexa censurada de forma similar ao que xa xemos con datos completos. No contexto de datos censurados, como vimos no capítulo anterior, para estimar a función de distribución (en lugar da función de distribución empírica), temos o estimador Kaplan- Meier b FKM (y) = n X i=1 WiI(Zi≤y). (3.2) onde Zi denota a variable resposta observada (posto que xa non coñecemos os valores de Yi ) e δi permítenos identicar a presenza de datos censurados. Tendo en mente o estimador de Kapla-Meier xunto co estimador (3.1), podemos pensar en estimar o parámetro b β asociado a un modelo de regresión con variable resposta Y censurada como: b β= arg m´ın n X i=1 Wibεi2, (3.3) sendo bεi=Zi−b β0−b β1xi os residuos do modelo de regresión. Esta técnica para atopar o estimador usando os pesos de Kaplan-Meier foi proposta por Stute en [16] (e por iso no seguinte capítulo lle chamaremos método de Stute) e estudada por Sánchez-Sellero na súa tese, que podemos ver en [14]. En concreto, demostraron que o estimador b β denido en (3.3) é asintoticamente normal, coa seguinte distribución límite para do modelo de regresión lineal simple que estamos tratando: √nb β−βd −−−−→ N0,Ω−1Π Ω−1, considerando a seguinte notación: Ω = E∂m (x) ∂βr ,∂m (x) ∂βsr,s∈{0,1} , sendo m(x) = E(Y|X=x) = β0+β1x, x ∈R , e Π=( Cov (ηϕr, ηϕs))r,s∈{0,1}, 1 O superíndice dc fai referencia a que estamos traballando con datos completos. 3.2. O ESTIMADOR PROPOSTO POR MILLER 19 onde ϕr(x, y) = [y−m(x)] ∂m (β) ∂βr , e ηϕ i=ϕ(xi, Zi)γ0(Zi)δi+γϕ 1(Zi) (1 −δi)−γϕ 2(Zi), γ0(y) = exp (Zy− −∞ e H0(dZ) 1−H(Z)), sendo F(u) e D(u) as funcións de distribución das variables Y e C respectivamente, entón temos: e H0(y) = P(Z⩽y, δ = 0) = Zy −∞ (1 −F(u)) D(du), e e H11 (x, y) = P(X⩽x, z ⩽y, δ = 1) , H(y) = P(Z⩽y). Ademais, para ϕ unha función real medible denida en R2 , denimos: γϕ 1(y) = 1 1−H(y)ZI{y<w}ϕ(x, w)γ0(w)e H11 (dx, dw), e γϕ 2(y) = Z Z I{v<y,v<w}ϕ(x, w)γ0(w) [1 −H(v)]2e H0(dv)e H11 (dx, dw). A demostración deste resultado non é obxectivo deste traballo pero podemos atopar todos os detalles explicados en [14]. 3.2. O estimador proposto por Miller Recordemos primeiro brevemente en que consiste o método de máxima verosimilitude que imos empregar para datos completos. Sexa {V1, ... , Vn} unha mostra aleatoria simple dunha variable V con función de distribución Fθ ou función de densidade fθ . O estimador de máxima verosimilitude (en adiante, EMV) é aquel valor que maximiza a masa de probabilidade (ou densidade) da mostra. Se denimos a función de verosimilitude como L(θ) = L(V1, ... , Vn;θ) = fθ(V1, ... , Vn) = n Y i=1 fθ(Yi). (3.4) Así, para cada mostra particular {Y1, ... , Yn} a estimación de máxima verosimilitude de β é o valor b βMV que maximiza a verosimilitude denida en (3.4), é dicir: LV1, ... , Vn;b θMV = m´ax θL(V1, ... , Vn;θ). Polo tanto, sexa θ o vector dos parámetros, o procedemento que hai que seguir para obter o EMV de θj , dada unha mostra {V1, ... , Vn} é: 20 CAPÍTULO 3. REGRESIÓN CENSURADA 1. Escribir a función de verosimilitude L(θ) = L(V1, ... , Vn;θ) . 2. Escribir o logaritmo da verosimilitude l(θ) = ln L(θ) , xa que o máximo non se ve alterado a través desta transformación e así os cálculos resultan máis sinxelos. 3. Obter o θj que cumpra ∂ ∂θj l(θ)=0, e denotámolo por b θj . Unha vez que xa falamos en termos xerais do método de máxima verosimilitude, vexamos agora que ocorre se estamos na situación dos datos censurados. Nesta sección seguiremos o proceso que aparece no libro de Miller que podemos atopar en [11] nas páxinas 11-14, onde aparece explicado para o caso dun modelo de regresión lineal múltiple e neste traballo adaptámolo para o noso contexto dun modelo lineal simple. Consideramos o par (Zi, δi) , recordando que Zi= m´ın (Yi, Ci) , sendo Ci os valores da variable de censura e δ=I[Y⩽C] é o indicador da censura. Asociado a este par temos a súa correspondente función de densidade: L(Zi, δi) =        f(Zi)se δi= 1, S(Zi)se δi= 0, que tamén se pode expresar da seguinte maneira: L(Zi, δi) = f(Zi)δiS(Zi)1−δi, (3.5) sendo f(Z) a función de densidade para os datos sen censura e S(Z) a función de supervivencia para os datos censurados. Consideramos agora a función de verosimilitude de toda a mostra L(β) = L(Z1, ... , Zn;δ1, ... , δn) = n Y i=1 L(Zi, δi) = Y U f(Zi)! Y C S(Zi)!, (3.6) denotando por QU e QC o produto sobre os datos sen censura 2 e os datos censurados respectivamente. Nesta igualdade estamos empregando a independencia dos Zi e usamos (3.5). Ademais, baixo a suposición de que o tempo censurado e o tempo de supervivencia son independentes, cabe destacar que a función de verosimilitude estaría multiplicada por 2 Empregaremos a notación que vén do inglés: uncensored e censored . 3.2. O ESTIMADOR PROPOSTO POR MILLER 21 unha constante que non depende do parámetro β e como o noso obxectivo é maximizar a función, podemos prescindir desta constante. Recordemos que estamos considerando o vector dos parámetros β= (β0, β1) . Para atopar o m´axβL(β) faremos unha transformación mediante o logaritmo e logo derivamos e igualamos a cero para obter o máximo tal e como detallamos no caso de datos completos. Desta forma, usando propiedades da función logaritmo quédanos: 0 = ∂ ∂βj log L(β) = n X i=1 ∂ ∂βj log Lβ(Zi, δi), con j∈ {0,1}. (3.7) Se substituímos agora o valor de L(β) da ecuación (3.6), temos o seguinte: 0 = X U ∂ ∂βj log fβ(Zi) + X C ∂ ∂βj log Sβ(Zi), con j∈ {0,1}. Debido á súa complexidade, para poder resolver esta igualdade e obter o valor de b β teremos que programar un método iterativo coa axuda dun ordenador. Neste traballo, por exemplo, detallaremos o método de Newton Rapson que é un algoritmo iterativo que se emprega con frecuencia para atopar aproximacións dos ceros dunha función G(x) necesitando un punto inicial X0 . A partir dunha iteración pasamos a seguinte e así sucesivamente mediante a seguinte ecuación: Xn+1 =Xn−G(Xn) G0(Xn). Para simplicar a notación escribiremos Li(β) = Lβ(Zi, δi) con i= 1, ... , n . Así, podemos reescribir (3.7) da seguinte maneira: 0 = n X i=1 ∂ ∂βj log Li(β), j = 1, ... , p, ou ben 0 = ∂ ∂β log L(β), sendo ∂ ∂β log L(β) = ∂ ∂β0 log L(β),∂ ∂β1 log L(β), ∂2 ∂β2log L(β) =        ∂2 ∂β2 0 log L(β)∂2 ∂β0∂β1 log L(β) ∂2 ∂β1∂β0 log L(β)∂2 ∂β2 1 log L(β)        . 22 CAPÍTULO 3. REGRESIÓN CENSURADA Supoñamos que a solución inicial é b β0=b β0 0,b β0 1 . Empregaremos o desenvolvemento de Taylor 3 para aproximar esa función e poder escribila da seguinte forma: 0 = n X i=1 ∂ ∂βj log Lib β= n X i=1 ∂ ∂βj log Lib β0+ 1 X k=0 b βk−b β0 kn X i=1 ∂2 ∂βk∂βj log Lib β0+... , con j= 0,1 , e tamén o podemos escribir como 0 = ∂ ∂β log Lb β=∂ ∂β log Lb β0+b β−b β0∂2 ∂β2log Lb β0+... . Usamos agora o método de Newton-Rapson e temos que a solución será b β1=b β0− ∂ ∂β log Lb β0 ∂2 ∂β2log Lb β0. (3.8) Este resultado proporciónanos as bases dunha aproximación, de maneira iterativa, para calcular o EMV. Así, dado un valor inicial b β0 , usamos a ecuación (3.8) para obter unha mellor estimación e repetimos este proceso ata xerar unha sucesión de estimadores que converxen ao EMV b β . Consideraremos a seguinte notación: ∂ ∂β log Lb β0 será o vector derivada en b β0 . −∂2 ∂β2log Lb β0 será a matriz de información en b β0 e denotarémola por ib β0 . Cabe destacar que E(i(β)) = −E∂2 ∂βk∂βj log L(β)=I(β), sendo I(β) a matriz de información de Fisher . Se agora substituímos a expresión do vector derivada e a matriz de información en (3.8) obtemos o seguinte: b β1=b β0+I−1b β0∂ ∂β log Lb β0. (3.9) Debido á dicultade para estimar a función de densidade ou a función de Supervivencia, podemos asumir que a distribución do erro é normal para realizar a estimación de máxima verosimilitude. Debemos destacar tamén que describimos o procedemento para o caso do modelo lineal simple pero de maneira intuitiva poderiamos estendelo ao contexto de múltiples variables explicativas. 3 Recordemos que o desenvolvemento de Taylor dunha función g(x) consiste en escribir g(x) = g(x0) + (x−x0)g0(x0) + (x−x0)2 2! g00 (x0) + ... 3.3. O ESTIMADOR PROPOSTO POR BUCKLEY E JAMES 23 3.3. O estimador proposto por Buckley e James Na Introdución deste traballo, no Capítulo 1, xa empregamos as ecuacións normais para estimar os parámetros no caso do modelo de regresión lineal con datos completos. Imos modicar estas ecuacións para o caso de observacións censuradas pero primeiro recordemos que no caso de datos completos temos que escoller como estimadores de β0 e β1 aqueles valores b β0 e b β1 que satisfagan: n X i=1 Yi−b β0−b β1xi= 0, (3.10) n X i=1 (xi−x)Yi−b β1xi= 0, (3.11) sendo x=1 nPn i=1 xi . Para o caso de datos censurados, como non podemos observar todos os datos {Y1, ... , Yn} , consideraremos Y∗ i=Yiδi+E(Yi|Yi> Ci) (1 −δi), i = 1, ... , n, sendo Ci os valores da variable de censura, δi=I[Yi⩽Ci] o indicador da censura, onde δi= 0 para os datos censurados e δi= 1 para os datos non censurados. Entón, E(Y∗ i) = β0+β1xi e polo tanto: E n X i=1 (xi−x) (Y∗ i−β1xi)!= 0. Por analoxía con (3.11), o ideal sería considerar un estimador b β1 para o cal se cumpra: n X i=1 (xi−x)Y∗ i−b β1xi= 0. Debido a que E(Yi|Yi> Ci) é descoñecida, imos considerar unha aproximación consistente empregando a función de Kaplan-Meier b FKM , onde b FKM (ε)=1−Y i;bεi≤bεn−i n−i+ 1δi , sendo bεib β0,b β1=Zi−b β0−b β1xi , Zi= m´ın (Yi, Ci) e Ci os valores da variable de censura. Imos substituír as observacións censuradas por Yib β1=b β1xi+X U Wik b β1Yk−b β1xk, (3.12) 24 CAPÍTULO 3. REGRESIÓN CENSURADA onde neste caso PU é un sumatorio dos datos non censurados sobre k e Wik b β1=             vkb β1 1−b FKM Ci−b β1xise bεi(0, b)<bεk0,b β1, 0 outro caso. onde vkb β1 é a masa de probabilidade asignada a función de distribución b FKM (ε) e Temos que escoller entón un estimador b β1 que satisfaga: b β1=PUYi(xi−x) + PCYi(β) (xi−x) Pn i=1 (xi−x)2. (3.13) Cabe destacar que isto é equivalente a substituír cada punto censurado (Ci, xi) polos puntos {b β1xi+Yk−b β1xk, xi}, k 6=i , darlle o peso Wik b β1 e nalmente substituílos nas ecuacións normais. Denotamos agora a parte dereita da ecuación (3.13) por γb β1 . Desta forma, en vez de querer minimizar unha función, o noso obxectivo será buscar un b que cumpra b=γ(b) , é dicir, queremos resolver a ecuación γ(b)−b= 0, (3.14) da que b β1 é raíz. Este problema pódese resolver empregando métodos de optimización e así (3.13) reescribímola como b β1=PU kYknxk−x+PC jWjk b β1(xj−x)o ηb β1 sendo ηb β1= n X i=1 (xi−x)2− C X j (xj−x)nxj−exjb β1o, onde exjb β1= U X k Wjk b β1xk. Este proceso lévase a cabo no artigo de Buckley e James que podemos ver con máis detalle en [1]. Unha vez que obtemos b β1 , podemos conseguir de forma análoga b β0 de forma que b β0=nPUYi+PCYib β1o n−b β1x. 3.4. O ESTIMADOR PROPOSTO POR JIN, LIN E YING 25 3.4. O estimador proposto por Jin, Lin e Ying O estimador obtido por Buckely e James é unha raíz da función de estimación que atopamos na ecuación (3.14) que non é nin continua nin monótona e as súas raíces poden non existir. O algoritmo iterativo de Buckley e James presenta algún problema, entre os que destacan que a converxencia do algoritmo non está garantida, é dicir, en ocasións non se atopa a solución ou esta solución oscila entre dous puntos non óptimos. Ademais, aínda que o algoritmo converxa, non está claro que nos leve a un estimador consistente xa que os resultados teóricos foron establecidos baseándose na hipótese de linearidade local. Debido a estes problemas, Jin, Lin e Ying intentaron mellorar este estimador no artigo que podemos ver con detalle en [8]. Para mellorar o método presentado por Buckley e James, será fundamental o papel que xoga o estimador inicial. Se o estimador inicial é consistente, entón para cada paso m , o estimador obtido na m -éstima iteración tamén será consistente. Ademais, se o estimador é asintoticamente normal, entón o estimador obtido na m -éstima iteración tamén o será. Isto que expomos está demostrado por Ritov ou Lai e Ying en [13] e [10], respectivamente. Así, este novo procedemento lévanos a unha clase de estimadores consistentes e asintoticamente normais. Desta forma, a idea que temos que seguir é a de darlle un bo valor inicial ao problema e aplicar o método iterativo. Este tipo de procedementos empréganse frecuentemente no ámbito da Matemática Aplicada como pode ser no método de Newton-Rapson que xa empregamos neste traballo. Jin, Lin e Ying propoñen escoller este valor inicial a partir da función de peso de Gehan, que se pode calcular aplicando técnicas de programación linear. Consideraremos a mesma notación que empregamos neste Capítulo 3, exceptuando que neste caso consideraremos a transformación logaritmo para a variable resposta. Por simplicidade de notación imos escribir Yi aínda que estamos considerando log Yi . Así, traballaremos co seguinte modelo de regresión linear: Yi=βxi+εi, é dicir, non imos ter en conta o intercepto. Na Sección 3.3 calculamos o estimador da pendente de Buckley e James, que lle chamaremos b βBJ , que é a raíz de U(β, β) = 0 , sendo U(β, b) = γ(b)−β, e recordemos que a función γ é a parte dereita da ecuación (3.13), é dicir, γ(b) = PUYi(xi−x) + PCYi(b) (xi−x) Pn i=1 (xi−x)2. 26 CAPÍTULO 3. REGRESIÓN CENSURADA É fácil ver que U(β, β) non é nin continua nin monótona en β e polo tanto é difícil calcular o estimador deste xeito, especialmente cando β é multidimensional. Podemos linealizar a función de estimación primeiramente dando un valor inicial b e logo resolvendo U(β, b)=0 para β . Esta operación lévanos a realizar β=γ(b) . Continúase este proceso co seguinte algoritmo iterativo: b β(m)=γb β(m−1), m ≥1. (3.15) Un estimador inicial consistente e asintoticamente normal de β0 pódese conseguir polo método rank-based de Jin, Lin, Wei e Ying, que podemos ver en [7]. Establecemos o estimador inicial b β(0) como o estimador tipo Gehan, b βG , descrito en [4] e que podemos calcular minimizando a seguinte función convexa: n X i=1 n X j=1 δi{εi(β)−εj(β)}−, onde a−=I{a < 0}|a| , e recordemos que εi(β) = Zi−βxi , Zi= m´ın (Yi, Ci) , Ci os valores da variable de censura, δi= 0 para os datos censurados e δi= 1 para os datos non censurados. Este problema de minimización é en realidade un problema de programación lineal simple e para ver máis detalles pódese consultar [7]. Para cada m , en [8] próbase que b β(m) é consistente e asintoticamente normal. Ademais, b β(m) é unha combinación linear do estimador tipo Gehan b βG e do estimador proposto por Buckley e James b βBJ , de forma que b β(m)=I−D−1Amb βG+I−I−D−1Amb βBJ +Opn−1 2, (3.16) onde I é a matriz identidade. D:= l´ımn→∞ n−1Pn i=1 (xi−x)2 é a matriz de pendentes da función de estimación de mínimos cadrados para os datos completos (datos non censurados). A é a matriz de pendentes da función estimada de Buckley e James que está denida e explicada en [8] pero neste TFG non imos profundizar sobre esta denición xa que non a empregaremos. Cando a porcentaxe de censura tende a cero, a matriz A aproxímase a D . Así, o primeiro termo da ecuación (3.16) vai ser cero e cada b β(m) aproxima ao estimador usual de mínimos cadrados. Se o algoritmo iterativo descrito en (3.15) converxe, entón b β(m) será o estimador 3.4. O ESTIMADOR PROPOSTO POR JIN, LIN E YING 27 de Buckley e James. Aínda que a secuencia iterativa non converxa, os estimadores seguen sendo consistentes e asintoticamente normais. Recordemos que unha hipótese do modelo de regresión lineal simple é a normalidade dos erros, polo que os erros seguen unha distribución normal de media cero e varianza σ2 , é dicir, ε∈N0, σ2 , sendo a normal unha función non decrecente. Pódese demostrar (ver [8]) que se a función de distribución do erro é non decrecente (como é o noso caso), entón cando D−A é denida positiva isto implica que I−D−1Am se aproxima a cero ou que b β(m) se aproxima a b βBJ (estimador proposto por Buckley e James) cando m tende a ∞ (para tamaños de mostra grandes). 34 CAPÍTULO 4. ESTUDO DE SIMULACIÓN b β0b β1 Sesgo Var ECM Sesgo Var ECM σ= 0.5 n= 100 M1 −675.90 149.58 195.26 −2135.01 620.24 1076.07 M2 94.68 202.96 203.86 −399.08 696.21 712.14 M3 1294.65 282.98 450.59 −5630.95 766.03 3936.79 M4 23.63 135.25 135.30 −162.03 525.92 528.55 n= 500 M1 −671.83 27.79 72.93 −2302.93 119.01 649.35 M2 14.19 42.04 42.06 −109.78 154.20 155.41 M3 1116.61 59.71 184.39 −5256.08 168.37 2931.01 M4 8.42 29.27 29.28 −76.61 129.43 130.03 n= 1000 M1 −627.14 14.57 53.90 −2384.88 63.80 632.57 M2 9.75 22.36 22.37 −39.34 82.10 82.26 M3 1077.09 33.95 149.96 −5156.54 95.99 2754.98 M4 38.41 18.11 18.25 −106.93 84.22 85.36 σ= 1 n= 100 M1 −2816.86 474.36 1267.83 −5428.37 1860.84 4807.56 M2 75.07 1361.29 1361.86 −1324.28 4525.82 4701.19 M3 1272.32 1455.49 1617.37 −5450.49 3923.49 6894.27 M4 −164.31 511.31 514.01 −194.59 1838.27 1842.06 n= 500 M1 −2806.58 97.49 885.19 −5902.56 367.56 3851.58 M2 19.31 462.55 462.59 −550.60 1544.08 1574.40 M3 901.49 567.38 648.65 −4156.82 1604.62 3332.53 M4 −101.53 104.52 105.55 −7.81 377.54 377.55 n= 1000 M1 −2720.51 48.09 788.20 −6123.53 182.60 3932.37 M2 −28.59 292.14 292.22 −260.67 936.91 943.70 M3 729.16 361.41 414.58 −3668.39 1012.26 2357.96 M4 −37.43 50.74 50.88 −56.47 185.97 186.29 Táboa 4.2: Sesgo, varianza e ECM dos estimadores obtidos (multiplicados por 10000 ) para o Modelo 1 a partir dos diferentes métodos M1 (estimador de mínimos cadrados ordinario), M2 (estimador proposto por Stute), M3 (estimador proposto por Stute con pesos Kaplan- Meier presuavizados) e M4 (estimador proposto por Buckley e James), sendo a porcentaxe de censura do 50 % para diferentes tamaños de mostra (denotado por n ) e desviacións do erro (denotado por σ ). 4.2. MODELO CON INTERCEPTO 35 Método M1: o estimador de mínimos cadrados clásico (detallada na Sección 1.2) aplicado só sobre os datos que observamos completamente, é dicir, cando δ= 1 . Para aplicar este método empregaremos a función lm de . Método M2: o estimador proposto por Stute, é dicir, un estimador de mínimos cadrados ponderado con pesos Kaplan-Meier (detallado na Sección 3.1). Para aplicar este método empregaremos a función lm de con argumento weight os pesos Kaplan- Meier calculados coa función KMW do paquete condSURV . Método M3: unha pequena modicación do estimador proposto por Stute onde se empregan uns pesos Kaplan-Meier presuavidazados (comentado na Sección 2.5). Para aplicar este método empregaremos a función lm de con argumento weight os pesos Kaplan-Meier presuavizados calculados coa función PKMW do paquete condSURV . Método M4: o estimador proposto por Buckley e James (detallado na Sección 3.3). Para aplicar este método empregaremos a función bf do paquete rms de . Para comparar os métodos anteriores calcularemos o sesgo, a varianza e o ECM de cada estimador. Os resultados pódense ver nas Táboas 4.1 e 4.2 que teñen asociadas unha porcentaxe do 25 % e 50 % de censura, respectivamente. Nestas táboas atopamos os resultados multiplicados por 10000 para poder comparar ben os ECM. Así evitamos a aparición de valores 0.000 ao aproximar os resultados. En cada táboa calculamos as medidas resumo para diferentes desviacións típicas do erro e diferentes tamaños de mostra onde en cada escenario se realizaron 1000 réplicas Monte Carlo. O código empregado para levar a cabo este estudo de simulación atópase no Anexo A deste traballo. En canto á programación, destacamos que hai que redenir o maior dato da variable resposta observada como non censurado para que o estimador de Stute, M2, sexa consistente (problema derivado da consistencia do estimador de Kaplan-Meier). Ademais, o método de Buckley e James, M4, non sempre converxe, polo que temos que eliminar as iteracións que non converxan para obter uns resultados que poidamos comparar co resto de métodos. En primeiro lugar, se observamos as Táboas 4.1 e 4.5 en conxunto, podemos sacar as seguintes conclusións xerais: O ECM dos estimadores aumenta canto máis grande sexa a varianza do erro, o cal é lóxico ao haber máis dispersión nos datos como se pode ver na Figura 4.1. 36 CAPÍTULO 4. ESTUDO DE SIMULACIÓN (a) Desviación típica de 0.5 . (b) Desviación típica de 1. Figura 4.1: Representación gráca dunha mostra de tamaño n= 100 do Modelo 1 xunto coa recta de regresión teórica para diferentes desviacións típicas da distribución do erro. O ECM diminúe conforme aumentamos o tamaño de mostra, o cal tamén é de esperar xa que ao ter un tamaño de mostra maior, temos máis información. O ECM aumenta cando aumenta a porcentaxe de censura. O ECM é menor en todos os métodos para o caso de censura do 25 % se os comparamos co obtido nos casos de censura de 50 % . Este feito débese a que ao aumentar a censura imos ter menos información e as estimacións serán menos precisas. Para todos os métodos considerados resulta mellor a estimación do intercepto que a da pendente, pois sempre ten un ECM máis baixo. A modo de exemplo, imos observar os datos obtidos para o intercepto e a pendente cando a desviación típica do erro é 0.5 e ímonos xar nas Táboas 4.1 e 4.2. En ambos ca- 4.2. MODELO CON INTERCEPTO 37 sos, o método M3 é o que peor resultados nos proporciona, pois os seus ECM son os máis elevados. Este feito resulta curioso posto que os pesos Kaplan-Meier presuavizados proporcionan mellores resultados que os clásicos pesos Kaplan-Meier na estimación da función de distribución baixo censura. Porén isto non se observa no contexto da regresión censurada onde os mellores resultados son os asociados aos clásicos pesos Kaplan-Meier. Ademais, pese a que a primeira vista poderíamos pensar que o método M1 non ía proporcionar bos estimadores dado que estamos tendo en conta só os datos que observamos completamente, tamén temos que ter en conta que estamos considerando unha censura dun 25 % ou dun 50 % e unha desviación do erro de 0.5 . Así, para un tamaño de mostra de n= 1000 , contamos con en torno a 750 ou 500 observacións respectivamente e é lóxico que este método sexa capaz de aproximar ben o modelo tendo en conta a distribución do erro. Imos ver que ocorre cando a desviación típica do erro é 1 . No caso de contar cunha porcentaxe de censura do 25 % , na Táboa 4.1, observamos que a pendente estímase peor sempre mediante M1. Destacamos que para un tamaño de 100 , M3 estima mellor tanto β0 como β1 que M2. En canto ao intercepto, para este escenario, non obtemos un resultado unánime con respecto a cal é o peor método para estimalo. Se temos en conta agora unha censura do 50 % , na Táboa 4.2 podemos ver que neste caso o peor método é M1. Simplemente facemos unha excepción para o tamaño de mostra de 100 , onde o peor volvería ser M3. Estes resultados poñen de manifesto a utilidade de presentar métodos de estimación especícos para o contexto de datos censurados. Como conclusión xeral desta simulación, observamos que, en termos de ECM, o que da mellores resultados sempre é o asociado ao método M4 xa que sempre se observan menores ECM para todos os casos. Este resultado coincide co que tiñamos pensado atoparnos antes de iniciar a simulación, pois o estimador proposto por Buckley e James é o máis able destes catro métodos. O seguinte método que estima mellor os parámetros, no caso da pendente, sería M2 e nalmente teríamos M1 e M3. Para o caso do intercepto, o segundo mellor método non está moi ben denido xa que para cada escenario contamos con diferentes conclusións. 4.2.2. Erro con distribución chi-cadrado Na Sección 4.2.1 observamos o caso no que o erro segue unha distribución normal, que é unha das hipóteses do modelo de regresión linear simple. Que ocorrería se consideramos outra distribución para o erro que non sexa normal? Nesta sección imos estudar este caso para salientar a importancia de comprobar as hipóteses para realizar unha boa Inferencia 38 CAPÍTULO 4. ESTUDO DE SIMULACIÓN b β0b β1 Sesgo Var ECM Sesgo Var ECM Censura: 25 % n= 100 M1 2.225 0.091 5.043 −0.494 0.271 0.516 M2 2.777 0.795 8.508 −0.20 2 2.728 2.769 M3 2.784 0.698 8.448 −0.219 2.362 2.409 M4 2.760 0.159 7.775 −0.002 0.480 0.480 n= 500 M1 2.212 0.017 4.912 −0.531 0.048 0.330 M2 2.815 0.548 8.471 −0.116 1.913 1.926 M3 2.813 0.754 8.665 −0.109 2.653 2.665 M4 2.825 0.034 8.016 −0.006 0.098 0.099 n= 1000 M1 2.207 0.009 4.879 −0.531 0.024 0.306 M2 2.845 0.441 8.533 −0.120 1.558 1.573 M3 2.844 0.731 8.821 −0.101 2.601 2.612 M4 2.845 0.018 8.112 −0.005 0.048 0.048 Censura: 50 % n= 100 M1 1.718 0.074 3.025 −0.548 0.242 0.542 M2 2.511 1.104 7.410 −0.323 3.816 3.920 M3 2.678 0.895 8.069 −0.673 2.961 3.413 M4 2.543 0.131 6.598 −0.006 0.366 0.367 n= 500 M1 1.718 0.014 2.964 −0.612 0.040 0.414 M2 2.613 0.795 7.620 −0.221 2.721 2.770 M3 2.740 0.933 8.443 −0.511 3.115 3.376 M4 2.654 0.034 7.077 −0.013 0.081 0.081 n= 1000 M1 1.713 0.006 2.942 −0.617 0.018 0.399 M2 2.684 0.719 7.924 −0.262 2.496 2.564 M3 2.804 0.926 8.786 −0.528 3.166 3.445 M4 2.685 0.018 7.225 −0.006 0.037 0.037 Táboa 4.3: Sesgo, varianza e ECM dos estimadores obtidos para o Modelo 1B a partir dos diferentes métodos M1 (estimador de mínimos cadrados ordinario), M2 (estimador proposto por Stute), M3 (estimador proposto por Stute con pesos Kaplan-Meier presuavizados) e M4 (estimador proposto por Buckley e James), sendo unha porcentaxe de censura do 25 % e 50 % para diferentes tamaños de mostra (denotado por n ). 4.2. MODELO CON INTERCEPTO 39 Estatística. Comezaremos xerando en valores do modelo Modelo 1B: Y= 1 + 2X+ε, onde X representa a variable explicativa que segue unha distribución uniforme no intervalo [0,1] , é dicir, X∈U[0,1] , e ε representa o erro do modelo que segue unha distribución chi-cadrado 1 con tres grados de liberdade, é dicir, ε∈χ2 3 . Ademais, a variable de censura C seguirá unha distribución normal de varianza 1 e con diferentes medias de cara a controlar a porcentaxe de censura nos distintos escenarios considerados. Aclaramos as medias que ten que ter a variable C para os distintos casos que imos tratar nesta sección: Censura do 25 % (aproximadamente):  A variable C ten que ter unha media de 6.34 para poder acadar esta porcentaxe de censura. Censura do 50 % (aproximadamente):  A variable C ten que ter unha media de 4.52 para poder acadar esta porcentaxe de censura. Calcularemos o sesgo, a varianza e o ECM de cada estimador para os diferentes métodos M1 (estimador de mínimos cadrados ordinario), M2 (estimador proposto por Stute), M3 (estimador proposto por Stute con pesos Kaplan-Meier presuavizados) e M4 (estimador proposto por Buckley e James) explicados na Sección 4.2.1. Os resultados pódense ver na Táboa 4.3 onde consideramos diferentes tamaños de mostra (denotado por n ) e en cada escenario se realizaron 1000 réplicas Monte Carlo. A única diferenza salientable en canto a programación é a distribución do erro, que neste caso é unha chi-cadrado de tres grados de liberdade. Polo tanto, para xerar a variable do erro, empregaremos o comando rchisq . Na Táboa 4.3 podemos observar os datos obtidos e, se nos xamos nos ECM, podemos ver que son valores máis altos do habitual. Isto quere dicir que as estimacións non son boas. 1 Sexan Z1, ... , Zm variables aleatorias normais estándar independentes. Diremos que a variable aleatoria X=Z2 1+··· +Z2 m segue unha distribución chi-cadrado con m grados de liberdade, onde χ2 m é a notación para a distribución chi-cadrado e o subíndice m representa os grados de liberdade. 40 CAPÍTULO 4. ESTUDO DE SIMULACIÓN Figura 4.2: Representación gráca dunha mostra de tamaño n= 100 do Modelo 1B xunto coa recta de regresión teórica. Sen embargo, o método M4 é capaz de obter outra vez as mellores estimacións, aínda que, ao non cumprir a hipótese de normalidade dos erros, estes resultados non os poderemos empregar para facer inferencia, pois non obteremos resultados ables. Se observamos a Figura 4.2, podemos ver como a nube de puntos se distribúe en torno á recta de regresión poñendo de manifesto a asimetría da distribución do erro. Con este exemplo ilustramos a importancia de validar as hipóteses dun modelo de regresión lineal. 4.3. Modelo sen intercepto 4.3.1. Erro con distribución normal Consideremos agora o modelo Modelo 2: Y= 2X+ε, onde X representa a variable explicativa que segue unha distribución uniforme no intervalo [0,1], é dicir, X∈U[0,1] e ε representa o erro do modelo que segue unha distribución normal de media 0 e varianza σ2 , que se denota por ε∈N0, σ2 . Ademais, como no caso anterior, a variable de censura C seguirá unha distribución normal de varianza 1 e con diferentes medias de cara a manter constante a porcentaxe de censura nos distintos 4.3. MODELO SEN INTERCEPTO 41 σ= 0.5σ= 1 Métodos Sesgo Varianza ECM Sesgo Varianza ECM n= 50 M1 −0.167 0.022 0.050 −0.523 0.069 0.343 M2 −0.011 0.023 0.023 −0.044 0.080 0.082 M3 −0.200 0.019 0.059 −0.149 0.062 0.084 M4 −0.017 0.080 0.080 −0.031 0.286 0.287 M5 −0.004 0.021 0.021 −0.012 0.076 0.076 M6 −0.015 0.081 0.081 −0.025 0.292 0.293 n= 100 M1 −0.171 0.010 0.040 −0.539 0.032 0.323 M2 −0.002 0.011 0.011 −0.019 0.041 0.041 M3 −0.191 0.009 0.046 −0.116 0.032 0.046 M4 0.010 0.039 0.039 0.021 0.144 0.144 M5 0.002 0.010 0.010 0.003 0.036 0.036 M6 0.011 0.039 0.039 0.022 0.144 0.144 n= 200 M1 −0.176 0.006 0.037 −0.553 0.017 0.324 M2 0.002 0.005 0.005 −0.012 0.021 0.021 M3 −0.192 0.005 0.042 −0.104 0.017 0.028 M4 0.003 0.019 0.019 0.002 0.072 0.072 M5 0.004 0.005 0.005 0.004 0.019 0.019 M6 0.005 0.019 0.019 0.003 0.072 0.072 Táboa 4.4: Sesgo, varianza e ECM dos estimadores da pendente obtidos para o Modelo 2 a partir dos diferentes métodos M1 (estimador de mínimos cadrados ordinario), M2 (estimador proposto por Stute), M3 (estimador proposto por Stute con pesos Kaplan-Meier presuavizados) e M4 (estimador proposto por Buckley e James), M5 (estimador proposto por Miller) e M6 (estimador proposto por Jin, Lin e Ying) sendo unha porcentaxe de censura do 25 % para diferentes tamaños de mostra (denotado por n ) e desviacións do erro (denotado por σ ). 42 CAPÍTULO 4. ESTUDO DE SIMULACIÓN σ= 0.5σ= 1 Métodos Sesgo Varianza ECM Sesgo Varianza ECM n= 50 M1 −0.315 0.037 0.136 −0.967 0.104 1.039 M2 −0.034 0.041 0.042 −0.135 0.127 0.145 M3 −0.390 0.026 0.178 −0.386 0.078 0.227 M4 −0.024 0.105 0.106 −0.035 0.376 0.377 M5 −0.010 0.034 0.034 −0.025 0.115 0.116 M6 −0.012 0.109 0.109 −0.019 0.380 0.380 n= 100 M1 −0.334 0.019 0.131 −1.012 0.052 1.076 M2 −0.014 0.023 0.023 −0.084 0.068 0.075 M3 −0.373 0.014 0.153 −0.339 0.043 0.157 M4 0.002 0.057 0.057 0.018 0.193 0.193 M5 −0.002 0.018 0.018 −0.006 0.059 0.059 M6 0.010 0.057 0.057 0.026 0.193 0.194 n= 200 M1 −0.335 0.010 0.123 −1.031 0.027 1.090 M2 −0.005 0.011 0.011 −0.056 0.035 0.038 M3 −0.369 0.007 0.143 −0.299 0.025 0.114 M4 −0.006 0.028 0.028 −0.005 0.093 0.093 M5 0.005 0.009 0.009 0.005 0.029 0.029 M6 0.001 0.028 0.028 0.003 0.092 0.092 Táboa 4.5: Sesgo, varianza e ECM dos estimadores da pendente obtidos para o Modelo 2 a partir dos diferentes métodos M1 (estimador de mínimos cadrados ordinario), M2 (estimador proposto por Stute), M3 (estimador proposto por Stute con pesos Kaplan-Meier presuavizados) e M4 (estimador proposto por Buckley e James), M5 (estimador proposto por Miller) e M6 (estimador proposto por Jin, Lin e Ying) sendo unha porcentaxe de censura do 50 % para diferentes tamaños de mostra (denotado por n ) e desviacións do erro (denotado por σ ). 4.3. MODELO SEN INTERCEPTO 43 escenarios considerados. Aclaramos as medias que ten que ter a variable C para os distintos casos que imos tratar: Censura do 25 % (aproximadamente):  Se ε∈N0,1 2 , entón a variable C ten que ter unha media de 1.85 .  Se ε∈N(0,1) , entón a variable C ten que ter unha media de 2 . Censura do 50 % (aproximadamente):  Se ε∈N0,1 2 , entón a variable C ten que ter unha media de 1 .  Se ε∈N(0,1) , entón a variable C ten que ter unha media de 1 . Para obter estes resultados, realizamos unha simulación para un tamaño de mostra grande (n= 10000) e comprobamos as censuras para os datos indicados. Como podemos observar, se queremos aumentar a porcentaxe de censura, teremos que diminuír a media da variable de censura C . Imos comparar seis estimadores diferentes con parámetro β1= 2 . Os métodos empregados serán os catro métodos que empregamos na Sección 4.2.1 no caso do modelo con intercepto, M1 (estimador de mínimos cadrados ordinario), M2 (estimador proposto por Stute), M3 (estimador proposto por Stute con pesos Kaplan-Meier presuavizados) e M4 (estimador proposto por Buckley e James), e os dous seguintes: Método M5: o estimador proposto por Miller (detallado na Sección 3.2). Para aplicar este método empregaremos a función CensReg.SMN , pero o argumento status sería 1−δ (o contrario da nosa notación). Para nós, recordemos, que δ= 0 é un dato censurado e δ= 1 é un dato non censurado. Método M6: o estimador proposto por Jin, Lin e Ying (detallado na Sección 3.4). Para aplicar este método empregaremos a función lss . A razón pola que separamos o modelo con intercepto do modelo sen intercepto é pola forma na que están programadas as funcións CensReg.SMN e lss en . Estas dúas funcións, que se corresponden cos métodos M5 e M6, so estiman a pendente e consideran so o caso dun modelo sen intercepto. Para levar a cabo isto, como na Sección 4.2.1, empregaremos o programa estatístico e calcularemos o sesgo, a varianza e o ECM de cada estimador. Os resultados pódense 50 CAPÍTULO 5. APLICACIÓN A DATOS REAIS No ensaio A formaron parte 444 persoas e consistiu nunha comparación de terapias que se modicaron en comunidades onde se incorporou un programa de educación sobre a saúde e sobre a prevención de recaídas. Este estudo tivo unha duración de tres a seis meses. Ademais, ás/aos participantes ensinóuselle a recoñecer as situacións que supoñen un alto risco que poden ocasionar unha posible recaída e tamén se lles ensinaron habilidades que lle permitían afrontar estas situacións complicadas sen ter que facer uso das drogas. No ensaio B participaron 184 persoas que recibían un programa terapéutico ou ben de seis meses ou de doce meses de duración. Este programa involucraba un estilo de vida perfectamente estruturado nun entorno de vida comunal onde os membros compartían diferentes aspectos da súa vida mediante un vínculo. A base de datos UIS facilítanos información sobre ambos ensaios. Máis concretamente, na Táboa 5.1 atópase o conxunto de variables recollidas ao longo do estudo. A variable TIME considerarase como a variable resposta na análise estatística que imos realizar e defínese como o número de días dende que se realiza a admisión do/a individuo/a en calquera dos dous centros posibles ata que esa persoa teña unha recaída no uso das drogas. Nótese que se trata dunha variable censurada onde a variable CENSOR será a variable indicadora da censura que se produce como consecuencia dunha perda do seguimento dun/- dunha doente. Neste estudo, que conta cun 19.3 % de censura, tamén se considera que se unha persoa sae do ensaio e, polo tanto, se perde o seguimento, é debido a unha recaída no uso das drogas. Poderíamos realizar diferentes estudos tendo en conta as posibles variables explicativas que poden ser signicativas para explicar o comportamento da variable resposta. Ademais, contamos cunha serie de variables categóricas que nos servirán para facer subgrupos e observar o comportamento dependendo, por exemplo, da raza do/a doente ou do tipo de droga consumida antes de entrar ao ensaio. As posibles variables explicativas consideradas serán: AGE, BECK, NDT, LEN.T e as variables categóricas que teremos en conta para os nosos subgrupos serán: HERCOC, IV, RACE, TREAT e SITE. 5.1. BASE DE DATOS UIS 51 Variable Descrición Códigos/Valores ID Código de identicación 1-628 AGE Idade de inscrición Anos BECKTOTA Puntuación de depresión na admisión 0.000 −54.000 HERCOC Uso de heroína ou cocaína durante 1=Heroína e cocaína meses antes da admisión 2=So heroína 3= So cocaína 4=Nin cocaína nin heroína IVHX Historia do uso da droga 1=Nunca 2=Previamente 3=Recentemente NDRUGTX Número de tratamentos anteriores 0−40 de drogas RACE Raza do/a individuo/a 0=Branca 1=Non branca TREAT Tratamento asignado 0=Curto 1=Largo SITE Sitio do tratamento 0=A 1=B LEN.T Período de estancia no tratamento Días TIME Tempo ata a recaída Días CENSOR Evento de perda do seguimento ou 1=Volver ás drogas ou recaída no uso das drogas perda do seguimento 0=Outro caso Y Logaritmo da variable TIME FRAC Lonxitude do tratamento fraccionado LEN.T/90, tratamento curto LEN.T/180, tratamento longo IV3 Uso recente de drogas 1=Si 0=Non Táboa 5.1: Explicación das variables procedentes da base de datos UIS. 52 CAPÍTULO 5. APLICACIÓN A DATOS REAIS 5.2. Análise descritiva previa Pode ser interesante realizar un estudo inicial básico para observar as características dos datos que imos analizar máis adiante. Na Figura 5.1 podemos observar os histogramas das posibles variables explicativas que poderíamos incluír nun modelo de regresión para explicar a variable TIME . Ademais destes grácos, tamén podemos calcular medidas características como a media, o mínimo ou o máximo de cada unha destas variables que obtemos mediante a función summary de . Así, observamos que a idade media dos/as individuos/as que participaron no estudo é de 32.38 anos mentres que as idades están comprendidas entre os 20 e os 56 anos. No caso da variable BECK , obtemos que o índice de depresión entre os e as participantes ten unha media de 17.37 e os seus valores atópanse entre 0 e 54 . En canto ao número de tratamentos anteriores de droga, NDT , observamos que o número medio é de 4.543 tratamentos e atopamos individuos/as con 0 tratamentos e tamén con 40 tratamentos, sendo este o máximo desta variable. Por último, destacamos que o número medio de días de estancia no tratamento, variable LEN.T , é de 100.8 días, sendo o mínimo 3 e o máximo de estadía 400 días. En canto ás variables categóricas, podemos denir subgrupos con certas características a partir da base de datos orixinal. En primeiro lugar, se comezamos estudando a variable HC , podemos ver como 104 persoas consumían heroína e cocaína, 107 so heroína, 172 so cocaína e 192 ningunha destas dúas drogas. Por outra banda, no estudo participaron 430 individuos/as de raza branca e 145 de raza non branca. Ademais, no tratamento curto formaron parte 289 persoas e no tratamento longo 286 . Finalmente, tal como se explicou ao comezo, 400 persoas participaron no ensaio no lugar A e 175 no lugar B. 5.3. Estimación dun modelo de regresión Realizaremos un estudo para estimar o intercepto e a pendente de diferentes modelos propostos a partir da base de datos que acabamos de introducir. Parece interesante estudar como afecta a idade ao tempo que pasa ata que se recae no consumo de drogas e comprobar como é a relación entre ambas variables. Para isto empregaremos un modelo de regresión lineal simple entre a variable explicativa AGE e a variable resposta TIME descritas na Táboa 5.1 que podemos escribir como: TIME =β0+β1 AGE +ε, (5.1) 5.3. ESTIMACIÓN DUN MODELO DE REGRESIÓN 53 (a) Histograma da variable AGE . (b) Histograma da variable BECK . (c) Histograma da variable NDT . (d) Histograma da variable LEN.T . Figura 5.1: Histogramas das posibles variables explicativas consideradas para explicar o comportamento da variable TIME . sendo ε o erro que segue unha distribución normal e neste caso a variable resposta é o tempo que pasa ata a recaída nas drogas, é dicir, TIME , e a variable explicativa será a idade, é dicir, AGE . Ademais, cos coñecementos obtidos durante a realización do estudo de simulación, para estimar o intercepto e a pendente sabemos que o método que nos proporciona mellores resultados é M4, é dicir, empregaremos o estimador proposto por Buckley e James. Por outra banda, se obtemos resultados que nos indiquen que o intercepto non é signicativo, entón sabemos tamén do estudo de simulación que o mellor método que podemos empregar é M5, é dicir, o estimador proposto por Miller. Estes métodos e a súa programación está explicada con máis detalle nas Seccións 4.2 e 4.3, respectivamente. Mostramos a continuación o código empregado en así como as saídas resultantes para describir a relación entre a idade e o tempo ata a recaída: >bj(Surv(TIME,CENSOR)~AGE,link="identity",x=TRUE, y=TRUE) Buckley-James Censored Data Regression 54 CAPÍTULO 5. APLICACIÓN A DATOS REAIS bj(formula = Surv(TIME, CENSOR) ~ AGE, link = "identity", x = TRUE, y = TRUE) Discrimination Indexes Obs 575 Regression d.f.1 g 33.396 Events 464 sigma124.2258 d.f. 462 Coef S.E. Wald Z Pr(>|Z|) Intercept 183.5572 31.2655 5.87 <0.0001 AGE 4.7533 0.9541 4.98 <0.0001 Observación 5.1 . Para comprobar que os valores estimados son signicativos, recordemos que temos que xarnos na columna de Pr(>|Z|) , onde se obtén o nivel crítico para o contraste de que o coeciente é cero. Se o nivel crítico é menor que os niveis de signicación habituais (10 %,5 % ou 1 %) entón diremos que dito coeciente é signicativamente distinto de cero. Por outra banda, no caso de que non sexa signicativo, entón poderemos considerar que dito parámetro é cero, e polo tanto, non aporta nada ao modelo de regresión. Se obtemos que a pendente non é signicativa, estamos dicindo que as variables non estarían relacionadas. Nesta saída observamos que tanto o intercepto como a pendente son signicativamente distintas de cero e teñen valores de 188.557 días e 4.753 días/anos respectivamente. Na Figura 5.2 mostramos o diagrama de dispersión xunto co axuste obtido e observamos unha relación lineal crecente entre ambas variables. Ademais, tamén podemos pensar en axustar o modelo (5.1) para os diferentes subgrupos que determinan as variables categóricas proporcionadas pola base de datos e, como xa introducimos antes, empregaremos a raza, variable RACE , e o tipo de droga consumida, variable HC . 5.3. ESTIMACIÓN DUN MODELO DE REGRESIÓN 55 Figura 5.2: Representación da idade dos/as participantes no ensaio clínico fronte ao tempo ata a recaída xunto co modelo (5.1) axustado. Raza Realizaremos este estudo separando as razas para comprobar se hai algunha diferencia entre os resultados obtidos para a raza branca e a non branca. Comezaremos analizando con detalle o caso dos/as individuos/as que forman parte do estudo que son de raza branca. Iremos escribindo as funcións empregadas en e as súas saídas resultantes. Para o caso das persoas de raza non branca podemos observar os resultados obtidos na Táboa 5.2 e unicamente comentaremos que tanto o intercepto como a pendente son signicativos e polo tanto son distintos de cero. Realizamos a estimación mediante o método M4, o proposto por Buckley e James, e usaremos a función bj : >bj(Surv(TIME,CENSOR)[RACE==0]~AGE[RACE==0],link="identity",x=TRUE, y=TRUE) Buckley-James Censored Data Regression bj(formula = Surv(TIME, CENSOR)[RACE == 0] ~ AGE[RACE == 0], link = "identity", x = TRUE, y = TRUE) Discrimination Indexes Obs 430 Regression d.f.1 g 26.442 Events 357 sigma122.9842 56 CAPÍTULO 5. APLICACIÓN A DATOS REAIS d.f. 355 Coef S.E. Wald Z Pr(>|Z|) Intercept 202.6077 34.8628 5.81 <0.0001 AGE 3.7079 1.0625 3.49 0.0005 Nesta saída observamos como tanto o intercepto coma a pendente son signicativos pois o nivel crítico de ambos parámetros é menor que os niveis de signicación habituais. Obtemos por tanto estimacións de 202.608 días e 3.708 días/anos respectivamente para cada parámetro. (a) Raza branca. (b) Raza non branca. Figura 5.3: Representacións da idade dos/as individuos/as de raza branca e non branca fronte ao tempo ata a recaída en drogas xunto cos axustes obtidos en cada caso. Se nos xamos na Táboa 5.2, podemos comprobar que as estimacións teñen o mesmo signo para ambas razas. Polo tanto, podemos observar efectos similares para os dous grupos, onde se obtén unha relación lineal crecente, ao ser a pendente positiva. Entón, podemos concluír que a medida que aumenta a idade dos/as doentes, tanto as persoas de raza branca como de raza non branca van a permanecer máis tempo sen tomar drogas, tal como se observa na Figura 5.3. Ademais, o efecto da idade é máis notable para as persoas de raza non branca xa que a pendente estimada ten un valor máis alto. Tipo de droga Realizaremos agora o mesmo procedemento feito para o caso das razas pero separando aos/ás individuos/as dependendo do tipo de droga consumida antes de entrar no ensaio clí- 5.3. ESTIMACIÓN DUN MODELO DE REGRESIÓN 57 Raza Raza 0: Branca Raza 1: Non branca b β0202.608 146.281 b β13.708 5.724 Táboa 5.2: Táboa das estimacións mediante o método Buckley e James do intercepto e a pendente para o modelo (5.1), considerando a variable explicativa AGE con respecto a variable resposta TIME diferenciando os resultados obtidos en función da variable RACE . nico: heroína e cocaína, so heroína, so cocaína e ningunha destas dúas drogas. Observamos todas as estimacións obtidas para os diferentes casos na Táboa 5.3. Veremos con detalle a obtención dos parámetros cando o tipo de droga consumida é a cocaína e cando non se consumiu ningunha destas posibles drogas. No caso da cocaína, atopamos que o intercepto non é signicativo, polo que teremos que empregar o método de Miller. Destacamos tamén que para o caso de consumir heroína e para ningunha destas drogas, as pendentes non son signicativas e polo tanto as variables resposta e explicativa non están relacionadas. Veremos que podemos facer cando nos atopamos esta situación. Tipo de droga Heroína e Cocaína Heroína Cocaína Nin cocaína nin heroína b β0−200.285 195.194 −3.790 (NS) 284.161 b β113.406 2.440 (NS) 10.602 −0.169 (NS) Táboa 5.3: Táboa das estimacións obtidas empregando o método M4 de Buckley e James do intercepto e a pendente no modelo (5.1) cando consideramos a variable explicativa AGE con respecto a variable resposta TIME diferenciando os resultados obtidos en función da variable HERCOC . Nótese que NS signica que dita estimación non resulta estatísticamente signicativa. A continuación, mostramos a saída de do correspondente estudo cando a droga consumida é a cocaína. Volvemos empregar o método M4 de Buckley e James e así, usamos a función bj : >bj(Surv(TIME,CENSOR)[HC==3]~AGE[HC==3],link="identity",x=TRUE, y=TRUE) Buckley-James Censored Data Regression 58 CAPÍTULO 5. APLICACIÓN A DATOS REAIS bj(formula = Surv(TIME, CENSOR)[HC == 3] ~ AGE[HC == 3], link = "identity", x = TRUE, y = TRUE) Discrimination Indexes Obs 172 Regression d.f.1 g 61.561 Events 131 sigma130.2886 d.f. 129 Coef S.E. Wald Z Pr(>|Z|) Intercept -3.7902 70.4801 -0.05 0.9571 AGE 10.6016 2.3038 4.60 <0.0001 Se nos xamos no nivel crítico do intercepto, que ten un valor de 0.9571 , podemos armar que este parámetro non vai ser signicativo e polo tanto temos evidencias de que o intercepto é cero. Unha vez concluído isto e usando os coñecementos adquiridos ao longo do Capítulo 4 deste traballo, sabemos que para un modelo sen intercepto o método que mellor estima a pendente é o método M5 proposto por Miller. Polo tanto, usaremos agora a función CensReg.SMN para obter unha mellor estimación da pendente: >CensReg.SMN(1-CENSOR[HC==3],AGE[HC==3],TIME[HC==3],cens="right", dist="Normal")$betas ------------------------------------------- EM estimates and SE ------------------------------------------- Estimates SE x1 9.69816 0.89852 sigma^2 67560.48197 13478.74770 ------------------------------------------ Model selection criteria ------------------------------------------- Loglik AIC BIC EDC Value -958.924 1921.849 1928.144 1923.095 5.3. ESTIMACIÓN DUN MODELO DE REGRESIÓN 59 ------------------------------------------- [,1] [1,] 9.698159 Obtemos un valor de 9.698 días/anos para a pendente, que é un resultado máis able que o obtido mediante o método de Buckley-James que observamos na Táboa 5.3. Analicemos agora a saída obtida cando as persoas que participaron no ensaio non consumiron nin cocaína nin heroína meses antes de entrar ao tratamento. Empregaremos en primeiro lugar a función bj : > bj(Surv(TIME,CENSOR)[HC==4]~AGE[HC==4],link="identity",x=TRUE, y=TRUE) Buckley-James Censored Data Regression bj(formula = Surv(TIME, CENSOR)[HC == 4] ~ AGE[HC == 4], link = "identity", x = TRUE, y = TRUE) Discrimination Indexes Obs 192 Regression d.f.1 g 1.235 Events 152 sigma127.2726 d.f. 150 Coef S.E. Wald Z Pr(>|Z|) Intercept 284.1605 52.2488 5.44 <0.0001 AGE -0.1687 1.5950 -0.11 0.9158 Se observamos o valor do nivel crítico para o intercepto vemos que é signicativo. Sen embargo, para o caso da pendente temos un valor de 0.9158 , claramente superior aos niveis de signicación habituais. Polo tanto podemos concluír que a variable resposta non está relacionada coa variable explicativa cando o/a individuo/a non consume nin cocaína nin heroína meses antes de entrar ao ensaio clínico. Despois de realizar este estudo, sabemos que a idade das persoas que consumen heroína e ningunha das dúas drogas consideradas non vai estar relacionada co tempo que tardan en recaer no consumo das mesmas. Ademais, tamén sabemos que o intercepto no caso da 66 ANEXO A: COMANDOS DE R points(2006,6,pch=4) points(2000,1, pch=16) points(2000,2, pch=16) points(2001,3, pch=16) points(2002,4, pch=16) points(2002,5, pch=16) points(2002,6, pch=16) abline(v=2000,lty=3,col="black") abline(v=2002.5,lty=3,col="black") abline(v=2007,lty=3,col="black") Gráca motivación datos censurados, Figura 2.2 #Tamaño n=100 n=100 xseq=seq(-3,3,by=0.01) plot(xseq,pnorm(xseq),type="l",lwd=2,xlab="y",ylab="Distribución") legend(x = "bottomright", legend = c("Datos censurados", "Datos completos"), fill = c("blue", "red")) x=rnorm(n) lines(ecdf(x),col="red",lwd=2) xc=x[x<=1] lines(ecdf(xc),col="blue",lwd=2) #Para tamaño n=10000 n=1000 xseq=seq(-3,3,by=0.01) plot(xseq,pnorm(xseq),type="l",lwd=2,xlab="y",ylab="Distribución") legend(x = "bottomright", legend = c("Datos censurados", "Datos completos"),fill = c("blue", "red")) 67 x=rnorm(n) lines(ecdf(x),col="red",lwd=6) xc=x[x<=1] lines(ecdf(xc),col="blue",lwd=2) Gráca Kaplan-Meier, Figura 2.3 n=10 x=runif(n,min=0,max=1) c=rnorm(n)+4 error=rnorm(n,sd=0.5) y=theta[1]+theta[2]*x+error z=pmin(y,c) delta=as.numeric(y<=c) datcens <- data.frame(t = x, cen = delta) library(survival) fit <- survfit(Surv(t, cen)~1, data = datcens) summary(fit) plot(fit, main = "Metodo de Kaplan-Meier") legend("bottomleft", c("F. Supervivencia", "Intervalo de confianza"), lty = 1:2) A.2. Comandos do Capítulo 4 Gráco comparativo de varianzas, Figura 4.1 #Varianza=0.5 n=100 theta=c(1,2) x=runif(n,min=0,max=1) 68 ANEXO A: COMANDOS DE R c=rnorm(n)+2.8# variable censurada error=rnorm(n,sd=0.5) # Error epsilon y=theta[1]+theta[2]*x+error # variable resposta z=pmin(y,c) delta=as.numeric(y<=c) plot(x,y,col=8,type="p",pch=19, xlab="x", ylab="y") abline(1,2) #valores reales #varianza=1 n=100 x=runif(n,min=0,max=1) c=rnorm(n)+3# variable censurada error=rnorm(n,sd=1) # Error epsilon y=theta[1]+theta[2]*x+error # variable resposta z=pmin(y,c) delta=as.numeric(y<=c) plot(x,y,col=8,type="p",pch=19,xlab="x", ylab="y") abline(1,2) #valores reales Gráco da Figura 4.2 n=100 theta=c(1,2) x=runif(n,min=0,max=1) c=rnorm(n)+2.8# variable censurada error<-rchisq(n,df=3) # Error Xi cadrado y=theta[1]+theta[2]*x+error # variable resposta z=pmin(y,c) plot(x,y,col=8,type="p",pch=19, xlab="x", ylab="y") abline(1,2) #valores reales Simulación do Modelo 4.2.1 #MODELO CON INTEREPTO (MODELO 1) 69 #CENSURA DUN 25% #Cargamos as librerias que usaremos library(survival) library(condSURV) #Librería para os pesos Kaplan-Meier library(xtable) #Librería para as táboas dos datos library(rms) #Librería para Buclkey-James #############--------------------------------------------------- #CASO DE DESVIACIÓN TÍPICA DO ERRO sd=0.5 #O resto de casos son análogos #############--------------------------------------------------- #Calculamos a censura: n=10000 x=runif(n,min=0,max=1) c=rnorm(n)+2.84 # Variable censurada error=rnorm(n,sd=0.5) # Error epsilon y=1+2*x+error # variable resposta z=pmin(y,c) delta=as.numeric(y<=c) censura=100-sum(delta)/length(delta)*100 #Porcentaxe de censura censura ###########----------------------------------------------------- set.seed(123456) #Fixamos semilla #Tamaño de mostra n=100 #n=500 #n=1000 70 ANEXO A: COMANDOS DE R m=1000 #Número de simulacións theta=c(1,2) #O beta teórico # As matrices theta_lm, theta_KM, theta_BJ e theta_Jin teran m filas # e 2 columnas: a primeira dos alphas (beta_0) e a segunda dos beta_1 theta_lm=matrix(nrow = m,ncol=2) theta_Stute=matrix(nrow=m,ncol=2) theta_KMS=matrix(nrow=m,ncol=2) theta_BJ=matrix(nrow=m,ncol=2) fallou=numeric(m) #vector empregado en Buclkey-James(M4) #Gardaremos os datos resultantes nunha matriz que defino agora: M=matrix(NA,nrow=4,ncol=6) for (i in 1:m){ x=runif(n,min=0,max=1) c=rnorm(n)+2.84 # Variable censurada error=rnorm(n,sd=0.1) # Error epsilon y=theta[1]+theta[2]*x+error # variable resposta z=pmin(y,c) delta=as.numeric(y<=c) indices=sort(z,index.return=T)$ix #Ordeno os datos z=z[indices] x=x[indices] delta=delta[indices] delta[n]=1 #Defino o ultimo dato como non censurado #Facemos agora os diferentes metodos: #MÉTODO M1:mínimos cadrados usual para datos sen censura #Escollemos os z e x que cumpran delta=1 (datos sen censura) lm(z[delta==1]~x[delta==1]) #O vector coef contén os alpha e beta estimados coa función lm coef=coefficients(lm(z[delta==1]~x[delta==1])) #Gardo na primeira columna os alphas e na segunda os betas 71 theta_lm[i,1:2]=coef #MÉTODO M2: Stute peso_KM=KMW(z,delta) #Calculo os pesos KM coa función KMW #Calculo os coeficientes empregando a función lm cos pesos KM coefStute=coefficients((lm(z~x,weights = peso_KM))) theta_Stute[i,1:2]=coefStute #MÉTODO M3: Stute con pesos KMS peso_KMS=PKMW(z,delta) #Calculo os pesos KM Suavizados coa función PKMW peso_KMS[n]=peso_KMS[n]+1-sum(peso_KMS) #Calculo os coeficientes empregando a función lm cos pesos KMS coefKMS=coefficients((lm(z~x,weights = peso_KMS))) theta_KMS[i,1:2]=coefKMS #MÉTODO M4: Buclkey-James #Empregamos a función bj modeloBJ=bj(Surv(z,delta)~x,link="identity",x=TRUE, y=TRUE) if (length(names(modeloBJ)) !=20){ fallou[i]=1 } else{ coefBJ=modeloBJ$coefficients theta_BJ[i,1:2]=coefBJ } } #Gardamos os datos nun arquivo: datos=cbind(theta_lm,theta_Stute,theta_KMS,theta_BJ) colnames(datos)=c("M1_I","M1_P","M2_I","M2_P", "M3_I", "M3_P", "M4_I", "M4_P") #Para n=100 write.table(datos, file = "modelo1_c25_100A.txt", sep = " " ,quote=F, col.names = colnames(datos),row.names=FALSE) 72 ANEXO A: COMANDOS DE R #Para n=500 #write.table(datos, file = "modelo1_c25_500A.txt", sep = " " ,quote=F, col.names = colnames(datos),row.names=FALSE) #Para n=1000 #write.table(datos, file = "modelo1_c25_1000A.txt", sep = " " ,quote=F, col.names = colnames(datos),row.names=FALSE) #Comprobamos que os escribimos correctamente: proba=read.table("modelo1_c25_100A.txt",header=TRUE) head(proba) #Calculo sesgo, varianza e ECM para cada caso #MÉTODO M1 media_lm=colMeans(theta_lm) sesgo_lm=media_lm-theta varianza_lm=c(var(theta_lm[,1]),var(theta_lm[,2])) ECM_lm=sesgo_lm^2+varianza_lm l1=c(sesgo_lm[1],varianza_lm[1],ECM_lm[1]) #correspondente a beta_0 l2=c(sesgo_lm[2],varianza_lm[2],ECM_lm[2]) #correspondente a beta_1 M[1,1:6]=c(l1,l2) #MÉTODO M2 media_Stute=colMeans(theta_Stute) sesgo_Stute=(media_Stute-theta) varianza_Stute=c(var(theta_Stute[,1]),var(theta_Stute[,2])) ECM_Stute=sesgo_Stute^2 +varianza_Stute Stute1=c(sesgo_Stute[1],varianza_Stute[1],ECM_Stute[1]) #beta_0 Stute2=c(sesgo_Stute[2],varianza_Stute[2],ECM_Stute[2]) #beta_1 M[2,1:6]=c(Stute1,Stute2) #MÉTODO M3 media_KMS=colMeans(theta_KMS) sesgo_KMS=media_KMS-theta varianza_KMS=c(var(theta_KMS[,1]),var(theta_KMS[,2])) 73 ECM_KMS=sesgo_KMS^2 +varianza_KMS KMS1=c(sesgo_KMS[1],varianza_KMS[1],ECM_KMS[1]) #beta_0 KMS2=c(sesgo_KMS[2],varianza_KMS[2],ECM_KMS[2]) #beta_1 M[3,1:6]=c(KMS1,KMS2) #MÉTODO M4 media_BJ=colMeans(theta_BJ[which(fallou !=1),]) sesgo_BJ=media_BJ-theta varianza_BJ=c(var(theta_BJ[which(fallou !=1),1]), var(theta_BJ[which(fallou !=1),2])) ECM_BJ=sesgo_BJ^2 +varianza_BJ BJ1=c(sesgo_BJ[1],varianza_BJ[1],ECM_BJ[1]) #beta_0 BJ2=c(sesgo_BJ[2],varianza_BJ[2],ECM_BJ[2]) #beta_1 M[4,1:6]=c(BJ1,BJ2) #Gardaremos os datos resultantes nunha matriz rownames(M)<-c("M1","M2", "M3", "M4") colnames(M)<-c("Sesgo", "Var","ECM","Sesgo", "Var","ECM") M M*10000 #Obtemos o sesgo, varianza e ECM de cada un dos metodos multiplicados por 10000 xtable(M*10000,digits = 3) #Para obter os comandos de Latex Simulación do Modelo 4.3.1 #MODELO SEN INTERCEPTO (MODELO 2) #CENSURA DUN 50% #Cargamos as librerias que usaremos library(survival) library(condSURV) #Libreria para os pesos Kaplan-Meier library(SMNCensReg) #Librería para o método de Miller 74 ANEXO A: COMANDOS DE R library(xtable) #Libreria para as táboas dos datos library(rms) #Libreria para o método de Buckley-James library(lss2) #Libreria para o método de Jin, Lin e Ying #############--------------------------------------------------- #CASO DE DESVIACIÓN TÍPICA DO ERRO: sd=0.5 #O resto de casos son análogos #############--------------------------------------------------- #Calculamos a censura: n=10000 x=runif(n,min=0,max=1) c=rnorm(n)+1 # Variable censurada error=rnorm(n,sd=0.5) # Error epsilon y=2*x+error # variable resposta z=pmin(y,c) delta=as.numeric(y<=c) censura=100-sum(delta)/length(delta)*100 #Porcentaxe de censura censura ###########----------------------------------------------------- set.seed(123456) #Fixamos semilla #Tamaño da mostra #n=50 #n=100 n=200 m=1000 #Numero de simulacións theta=2 #O beta teórico 75 # Os vectores theta_lm, theta_Stute, theta_KMS, theta_Miller, # theta_BJ e theta_Jin son vectores de m compoñentes onde # iremos gardando os resultados da estimación theta_lm=c(rep(NA,m)) theta_Stute=c(rep(NA,m)) theta_KMS=c(rep(NA,m)) theta_Miller=c(rep(NA,m)) theta_BJ=c(rep(NA,m)) theta_Jin=c(rep(NA,m)) fallou=numeric(m) #vector empregado en Buckley-James #Gardaremos os datos resultantes nunha matriz que defino agora: M=matrix(NA,nrow=6,ncol=3) for (i in 1:m){ set.seed(i) x=runif(n,min=0,max=1) c=rnorm(n)+1 # Variable censurada (para unha censura dun 25%) error=rnorm(n,sd=0.5) #Error epsilon y=theta*x+error #variable resposta z=pmin(y,c) delta=as.numeric(y<=c) #Ordeno os datos indices=sort(z,index.return=T)$ix z=z[indices] x=x[indices] delta=delta[indices] #Defino o ultimo dato como non censurado delta[n]=1 #Fago agora os diferentes metodos: 82 ANEXO A: COMANDOS DE R # Intercepto non significativo, entón usaremos o método de Miller CensReg.SMN(1-CENSOR[HC==3],AGE[HC==3],TIME[HC==3],cens="right", dist="Normal")$betas #------------------HC=4:Ningunha destas drogas----------------- bj(Surv(TIME,CENSOR)[HC==4]~AGE[HC==4],link="identity",x=TRUE, y=TRUE) bj(Surv(TIME,CENSOR)[HC==4]~AGE[HC==4]+BECK[HC==4]+NDT[HC==4]+ LEN.T[HC==4], link="identity",x=TRUE, y=TRUE) #------------------------------------------------- #RACE #------------------------------------------------- #---------------------RACE 0: branca------------------- bj(Surv(TIME,CENSOR)[RACE==0]~AGE[RACE==0],link="identity",x=TRUE, y=TRUE) plot(AGE[RACE==0],TIME[RACE==0]) abline(bj(Surv(TIME,CENSOR)[RACE==0]~AGE[RACE==0],link="identity", x=TRUE, y=TRUE)) #---------------------RACE 1: non branca------------------- bj(Surv(TIME,CENSOR)[RACE==1]~AGE[RACE==1],link="identity",x=TRUE, y=TRUE) plot(AGE[RACE==1],TIME[RACE==1]) abline(bj(Surv(TIME,CENSOR)[RACE==1]~AGE[RACE==1],link="identity", x=TRUE, y=TRUE)) Bibliografía [1] Buckley, J. e James, I., Linear regression with censored data , Biometrika, 66 (1979), 429436. [2] Cao R., López-de-Ullibarri, I., Janssen, P. e Veraverbeke, N., Presmoothed Kaplan- Meier and Nelson-Aalen estimators , Journal of Nonparametric Statistics, 17 (2005), 3156. [3] Dikta, G., On semiparametric random censorship models , Journal of Statistical Planning and Inference, 66 (1998), 253279. [4] Gehan, E.A., A generalized Wilcoxon test for comparing arbitrarily singlecensored samples , Biometrika, 52 (1965), 203223. [5] Gómez Villegas, M.A., Inferencia Estadística , Díaz de Santos (2005). [6] Hosmer, D. W., Lemeshow, S. e May, S., Applied Survival Analysis, Regression Modeling of Time-to-Event Data , Wiley, 2nd ed., (2008). [7] Jin, Z., Lin, D. Y., Wei, L. J e Ying, Z., Rank-based inference for the accelerated failure time model , Biometrika, 90 (2003), 34153. [8] Jin, Z., Lin, D. Y. e Ying, Z., On least-squares regression with censored data , Biometrika, 93 (2006), 147161. [9] Kaplan, E. L. e Meier, P., Nonparametric estimation from incomplete observations , Journal of the American Statistical Association, 53 (1958), 457481 [10] Lai, T. L. e Ying, Z., Large sample theory of modied Buckley-James estimator for regression analysis with censored data , Annals of Statistics, 10 (1991), 1370402. [11] Miller, R., Gong, G. e Muñoz, A., Survival Analysis , Thechnical report, Universidade de California (1980). 83 84 BIBLIOGRAFÍA [12] Moore, D. F., Applied Survival Analysis Using R , Springer (2016). [13] Ritov, Y., Estimation in a linear regression model with censored data , Annals of Statistics, 18 (1990), 30328. [14] Sánchez Sellero C., Inferencia Estadística en datos con censura y/o truncamiento , Universidade de Santiago de Compostela (2001). Tese de doutoramento. [15] Sánchez Sellero, C., Apuntes da materia Inferencia Estatística , Grao en Matemáticas, Universidade de Santiago de Compostela, 2018-2019. [16] Stute, W., Nonlinear censored regression , Statistica Sinica, 9 (1999), 10891102. [17] Stute, W., González Manteiga, W. e Sánchez Sellero, C., Nonparametric Model Checks In Censored Regression , Communication in Statistics-Theory and Methods, 29 (2000), 16111629.