Contributions to flexible bivariate regression models. Applications in Medicine and Environment

INTERNATIONAL DOCTORAL SCHOOL OF THE USC Óscar Lado Baleato PhD Thesis Contributions to flexible bivariate regression models. Applications in Medicine and Environment Santiago de Compostela, 2022 Doctoral Programme in Statistics and Operations Research DOCTORAL THESIS CONTRIBUTIONS TO FLEXIBLE BIVARIATE REGRESSION MODELS. APPLICATIONS IN MEDICINE AND ENVIRONMENT Óscar Lado Baleato INTERNATIONAL PHD SCHOOL OF THE UNIVERSITY OF SANTIAGO DE COMPOSTELA PHD PROGRAMME IN STATISTICS AND OPERATIONS RESEARCH SANTIAGO DE COMPOSTELA 2022 DECLARACIÓN DO AUTOR/A DA TESE Presento a miña tese, seguindo o procedemento axeitado ao Regulamento, e declaro que: 1) A tese abarca os resultados da elaboración do meu traballo. 2) De ser o caso, na tese faise referencia ás colaboracións que Bvo este traballo. 3) Confirmo que a tese non incorre en ningún Bpo de plaxio doutros autores nin de traballos presentados por min para a obtención doutros Htulos. 4) A tese é a versión definiBva presentada para a súa defensa e coincide a versión impresa coa presentada en formato electrónico. E comprométome a presentar o Compromiso Documental de Supervisión no caso de que o orixinal non estea na Escola. En SanBago de Compostela, 12 de setembro de 2022 Sinatura electrónica D. Óscar Lado Baleato Título da tese: ContribuBons to flexible bivariate regression models. ApplicaBons in Medicine and Environment AUTORIZACIÓN DOS DIRECTORES DA TESE Contributions to flexible bivariate regression models. Applications in Medicine and Environment D. Francisco Gude Sampedro D. Javier Roca Pardiñas INFORMAN: Que a presente tese, correspóndese co traballo realizado por D. Óscar Lado Baleato, baixo a nosa dirección e titorización, e autorizamos a súa presentación, considerando que reúne os requisitos esixidos no Regulamento de Estudos de Doutoramento da USC, e que como directores, e tutor, desta non incorre nas causas de abstención establecidas na Lei 40/2015. De acordo co indicado no Regulamento de Estudos de Doutoramento, declara tamén que a presente tese de doutoramento é idónea para ser defendida en base á modalidade de monográfica con reproducción de publicacións, nos que a participación do doutorando foi decisiva para a súa elaboración e as publicacións axustanse ao Plan de Investigación. En Santiago de Compostela, 8 de Novembro de 2022 Asinado: Francisco Gude Sampedro Javier Roca Pardiñas Universidade de Santiago de Compostela DISSERTATION Contributions to flexible bivariate regression models. Applications in Medicine and Environment Author: ´ Oscar Lado Baleato Advisors: Francisco Gude Sampedro Javier Roca Pardi˜ nas DEPARTMENT OF STATISTICS, MATHEMATICAL ANALYSIS, AND OPTIMIZATION Santiago de Compostela, 2022 ter defining the predictive accuracy of different time points, satisfactory results were obtained. This result highlights the general applicability and robustness of the model. In a following chapter, the model is further developed in terms of the predictive variables to be selected for use in a medical setting. A method is proposed for formally testing the effects of covariates on the reference region shape, based on the use of a bivariate Euclidean distance. The distribution under the null hypothesis was approximated by bootstrap resampling. This statistical test respected the nominal level when covariates had no effect, and offered a highly satisfactory power curve when their effects became stronger. The same test was then used for determining the effect of sex and age on children’s height and body weight. To the best of our knowledge no other statistical test exists for examining the effects of covariates on bivariate reference regions. This thesis concludes by discussing new research lines. A new method is proposed for estimating percentile curves and conditional reference regions for multivariate responses based on multivariate conditional transformation models (Klein et al.,2019). Estimates are based on the transformation of data to a reference distribution (usually the standard Gaussian) using a transformation function estimated from the data. This function offers a correspondence between the statistical measures of the reference distribution (e.g., mean, quantiles, density functions), and those of the original data. By applying this transformation function to a standard Gaussian tolerance region, multivariate conditional reference regions are obtained that show non-parametric restrictions, allowing non-linear structures of dependence between responses to be estimated. iv Resumo A regresi´ on multivariada mide o grao de asociaci´ on entre un conxunto de variables preditoras e varias variables resposta. O m´ etodo recom´ endase naqueles estudos onde existan m´ ultiples variables de interese correlacionadas. No ´ ambito da investigaci´ on m´ edica ´ e habitual encontrarse con este tipo de variables. Por exemplo, o peso dunha persoa non pode ser interpretado como saudable ou non saudable se non sabemos a s´ ua estatura. De feito para estimar o efecto de variables cl´ ınicas sobre os puntos diagn´ osticos de sobrepeso ou obesidade debemos considerar a informaci´ on de ambas variables de forma conxunta. No diagn´ ostico, e seguimento, de moitas enfermidades cr´ onicas, por exemplo a diabetes, tam´ en se empregan varios marcadores continuos cuxos resultados est´ an altamente correlacionados entre s´ ı. Na diabetes ´ e o caso da concentraci´ on de glicosa plasm´ atica en xax´ un e a porcentaxe de hemoglobina glicosilada. Noutras enfermidades atop´ amonos coa mesma problem´ atica, por exemplo; para medir a funci´ on hep´ atica empreganse as concentraci´ on de diversas transaminasas ou no diagn´ ostico de problemas tiroideos precisamos medir as concentraci´ ons das hormonas TSH, T3 e T4. Nestes casos, o uso dun modelo de regresi´ on multivariado ´ e de utilidade dado que caracteriza o efecto das variables cl´ ınicas sobre cada unha das respostas de forma marxinal, e a maiores estima o efecto das mesmas sobre a correlaci´ on entre as respostas. Por exemplo, ser´ ıa de interese comprobar se a idade modula a correlaci´ on entre as concentraci´ ons de glicosa e as da hemoglobina glicosilada. Deste xeito podemos comprobar se nun grupo de suxeitos ´ e esperable unha maior tasa de discordancias entre ambos criterios. Ademais, esta modelizaci´ on permitir´ ıanos identificar casos dunha combinaci´ on pouco frecuente dos marcadores, a cal pode ser indicativa dun peor progn´ ostico. Polo tanto, estes modelos permiten estimar de forma m´ ais certeira o efecto das variables cl´ ınicas sobre diversos marcadores diagn´ osticos, dando lugar a unha interpretaci´ on multivariada dos seus resultados, e en ´ ultimo termo mellorando as regras diagn´ osticas e o seguimento de diversas patolox´ ıas. A pesar da utilidade dos modelos de regresi´ on multivariada non se adoitan empregrar por parte de investigadores cl´ ınicos, nin noutras ´ areas do co˜ necemento. As primeiras propostas metodol´ oxicas datan dos anos 60 e permiten estimar o efecto das covariables sobre o vector de medias dunha variable resposta de distribuci´ on gaussiana bivariada. Sen embargo, esta proposta considera as variv anzas e covarianza da resposta como valores constantes e por tanto independentes das variables preditoras. No ´ ambito cl´ ınico ´ e habitual atopar variables de interese que seguen unha distribuci´ on non gaussiana e cuxa varianza aumenta coa idade e outras variables. Na ´ ultima decada a investigaci´ on estat´ ıstica en regresi´ on multivariada recibiu un pulo debido a proposta de modelos baseados en funci´ ons c´ opula param´ etricas. Estes permiten estimar o efecto das covariables sobre o vector de medias da resposta, e a maiores sobre as s´ uas varianzas e sobre a correlaci´ on entre elas. Ademais, a resposta pode seguir unha distribuci´ on complexa e o efecto das covariables pode ser estimado mediante preditores flexibles. A pesar das suas vantaxes, a regresi´ on multivariada baseada en funci´ ons c´ opula presenta restricci´ ons param´ etricas tanto para as respostas marxinais, como para a estrutura de dependencia entre as mesmas. As marxinais modelizanse mediante distribuci´ ons param´ etricas univariantes, mentres que a correlaci´ on require dunha funci´ on c´ opula, tam´ en param´ etrica. O investigador debe escoller a combinaci´ on destes tres elementos que ofreza o mellor axuste para os seus datos, sendo en moitos casos un axuste non satisfactorio. Ademais no contexto cl´ ınico presentan unha interpretabilidade complexa. Estes modelos permiten estimar o efecto de variables cl´ ınicas sobre a totalidade dos par´ ametros da resposta. Sen embargo, no contexto aplicado ser´ ıa de maior utilidade unha interpretaci´ on visual e probabil´ ıstica do efecto das covariables sobre os valores conxuntos da mesma. Nesta tese de doutoramento propo˜ nense novas formulaci´ ons da regresi´ on multivariada para a s´ ua aplicaci´ on a problemas de investigaci´ on cl´ ınica e con gran aplicabilidade noutros campos. Os modelos propostos non te˜ nen restricci´ ons param´ etricas para a resposta, co fin de ofrecer m´ etodos o m´ ais xerais posibles. Ademais permiten estimar os efectos das variables preditoras sobre os valores medios, variabilidade e a estrutura de correlaci´ on das respostas mediante preditores aditivos flexibles. Deste xeito os efectos non-lineais das covariables continuas estimanse mediante suavizadores locais tipo n´ ucleo, ou mediante splines penalizados. A partir do axuste destes modelos, podemos obter unha rexi´ on probabil´ ıstica que cont´ en unha porcentaxe dada dos datos. Mediante esta rexi´ on ´ e posible caracterizar de xeito non param´ etrico o cambio na forma da distribuci´ on bivariada da resposta en funci´ on do valor das covariables. En termos pr´ acticos, esta proposta ofrece a mesma utilidade ´ os investigadores m´ edicos que a regresi´ on cuantil. ´ E dicir, permite caracterizar que valores son m´ ais probables na poboaci´ on xeral axustando por caracter´ ısticas dos pacientes como a idade ou o sexo. Concretamente, podemos obter unha rexi´ on que conte˜ na o 95% dos resultados dos pacientes sans para cada idade da mostra, tendo as´ ı unha regra diagn´ ostica multivariada. O primeiro modelo de estimaci´ on desta rexi´ on presentase de forma detallada no Cap´ ıtulo 3. Ademais, no Ap´ endice A presentanse detalles computacionais. Mentres que no Ap´ endice C presentamos a implementaci´ on deste modelo en vi forma de software libre no paquete de R refreg. A estimaci´ on do modelo basease nun algoritmo complexo; nun primeiro paso estimanse as medias condicionais da resposta mediante preditores aditivos flexibles; no segundo paso, estimanse as varianzas e covarianza condicionais usando os residuos dos modelos obtidos no paso un. Finalmente, a estimaci´ on da rexi´ on realizase mediante a funci´ on de densidad bivariada dos residuos estandarizados. A estimaci´ on da densidade bivariada realizase cun estimador tipo n´ ucleo, cuxo par´ ametro de suavizaci´ on se obt´ en por validaci´ on cruzada minimizando a diferencia entre a cobertura esperada e a observada. A partir dos cinco modelos estimados, e coa rexi´ on obtida na escala dos residuos, ´ e posible obter a rexi´ on de referencia condicional para cada valor das covariables. Ademais, a rexi´ on estimada na escala dos residuos permite avaliar de forma visual que suxeitos se atopan f´ ora, ou dentro da mesma, tras axustar polas s´ uas caracter´ ısticas. Estudos de simulaci´ on num´ erica mostraron o bo desempe˜ no da estimaci´ on da rexi´ on de referencia condicional, cun erro de estimaci´ on baixo e cobertura pr´ oxima ´ o nivel n´ ominal, resultados que melloran con tama˜ nos de mostra maiores. A´ ında sendo os erros de estimaci´ on lixeiramente superiores nos escenarios con datos non gaussianos, o modelo tam´ en mostrou resultados satisfactorios nestes supostos. Esta proposta ten unha gran aplicabilidade por parte de investigadores cl´ ınicos na definici´ on de regras diagn´ osticas a partir dos resultados de marcadores continuos. Actualmente, est´ ımase que o 70% das decisi´ ons cl´ ınicas se basean nos resultados destes biomarcadores. Para a interpretaci´ on diagn´ ostica dos mesmos os m´ edicos empregan de forma rutineira os co˜ necidos como intervalos de referencia. Estes intervalos def´ ınense como dous puntos de corte entre os cales se atopan o 95% dos resultados dos pacientes sans. Desta forma, se un paciente non se atopa dentro deste intervalo, considerase que presenta un resultado an´ omalo indicativo dunha patolox´ ıa sen diagnosticar. Ademais, no caso de que unha variable cl´ ınica modifique os valores do marcador con independencia da patolox´ ıa, est´ ımanse intervalos de referencia condicionais, que reciben o nome de curvas de referencia. Deste xeito, def´ ınense varios puntos de corte dependendo de, por exemplo, a idade ou sexo do paciente. No caso das enfermidades que requiren dous ou m´ ais marcadores para o seu diagn´ ostico, precisamos dunha rexi´ on de referencia. Esta rexi´ on caracteriza os valores conxuntos dos marcadores continuos que se obte˜ nen na poboaci´ on san — tendo en conta deste modo a s´ ua correlaci´ on. A rexi´ on de referencia pode verse como a extensi´ on do intervalo de referencia para unha resposta bivariada continua. Esta interpretaci´ on permite obter unha maior especificidade e tam´ en unha maior sensibilidade no diagn´ ostico. De feito, a tasa de falsos positivos red´ ucese dado que se usa un s´ o criterio diagn´ ostico en lugar de dous criterios independentes univariantes. ´ O mesmo tempo tam´ en ofrece unha maior sensibilidade, dado que detecta pacientes que mostran resultados cunha combinaci´ on at´ ıpica de marcadores, casos que son indetectables mediante unha interpretaci´ on univii variante. Estas vantaxes foron definidas en traballos realizados a principio dos anos 80, sen embargo, a aplicaci´ on destas rexi´ ons ´ e anecd´ otica na investigaci´ on e pr´ actica cl´ ınicas. No noso caso, o m´ etodo estat´ ıstico proposto permitiu, por vez primeira, interpretar de forma conxunta os valores de d´ uas probas de control glic´ emico (niveis da glicosa en xax´ un e porcentaxe de hemoglobina glicosilada). Esta interpretaci´ on permite identificar distintos perf´ ıs de desregulaci´ on glic´ emica, avaliando dunha forma visual como cambia a distribuci´ on conxunta de ambos marcadores coa idade do paciente. Especificamente definimos 4 perf´ ıs de pacientes segundo a s´ ua posici´ on na rexi´ on de referencia. Na tipolox´ ıa 1 atopamos pacientes con valores altos en ambos marcadores; na 2, pacientes con valores altos de hemoglobina glicosilada e valores normais de glicosa; na 3, pacientes con valores extremadamente baixos en ambos marcadores; e na 4 pacientes con valores baixos de hemoglobina glicosilada e altos de glicosa en xax´ un. Os suxeitos dentro da tipolox´ ıa 1 probablemente padezan unha diabetes sen diagnosticar. Mentres que os individuos da tipolox´ ıa 2 poden ser descritos como alto glicadores. Un grupo xa descrito na literatura como pacientes con maior risco cardiovascular, dado que te˜ nen unha tasa de glicaci´ on alta, feito que se relaciona cun maior envellecemento vascular . No caso contrario estar´ ıan os pacientes do grupo 4 baixo a etiqueta de baixo glicadores. Os pacientes do grupo 3, en principio, non sofren problemas de desregulaci´ on glic´ emica. Ademais, describimos como os valores medios e variabilidade de ambos marcadores aumentan coa idade. As´ ı mesmo, a correlaci´ on entre ambos ´ e pr´ oxima a cero en suxeitos novos, e sobe de forma linear ata os 40 anos, a partires desa idade a s´ ua correlaci´ on ´ e constante, en torno a 0.40. Deste modo ´ e esperable obter unha porcentaxe maior de discordancias entre ambos m´ etodos en suxeitos de menor idade. Unha rexi´ on que cont´ en unha porcentaxe espec´ ıfica dos datos dunha distribuci´ on bivariada ´ e´ util na pr´ actica cl´ ınica, pero tam´ en noutros campos de investigaci´ on. Para exemplificar esta xenerabilidade, no Cap´ ıtulo 4 presentamos unha aplicaci´ on do modelo anteriormente mencionado a un problema medio ambiental. En concreto, usamos a rexi´ on de referencia como unha rexi´ on de predici´ on para estimar de forma conxunta os valores futuros de dous contaminantes atmosf´ ericos. Como variables preditoras contamos coas observaci´ ons previas dos contaminantes en distintos horizontes temporais. De forma emp´ ırica puidemos constatar que modelos con distinto n´ umero de variables preditoras amosaban unha cobertura dos datos similar. Sen embargo, estes modelos non amosaban a mesma precisi´ on a hora de modelizar un episodio de contaminaci´ on. Co fin de identificar os valores previos m´ ais informativos a hora de modelizar as concentraci´ ons futuras dos contaminantes, introducimos unha funci´ on de perda bivariada para a selecci´ on do modelo. Esta funci´ on mide a distancia entre o punto observado e o l´ ımite da rexi´ on correspondente a ese punto, tendo en conta ademais a cobertura desexada. Este novo estat´ ıstico demostrou a s´ ua utilidade para seleccionar o modelo m´ ais parsimonioso tanto con datos simulados como viii reais. O control dos niveis de contaminantes atmosf´ ericos ´ e de gran interese na sociedade actual, tanto polos seus efectos nocivos sobre a sa´ ude p´ ublica como polos danos medioambientais asociados. Ademais, a lexislaci´ on espa˜ nola recolle a obrigatoriedade de rexistrar a concentraci´ on de contaminantes atmosf´ ericos nas inmediaci´ ons de fontes de emisi´ on co˜ necidos (p.ex., f´ abricas electro-intensivas). No caso de superar unha media bihoraria establecida, a empresa emisora pode ser sancionada. Por tanto, existe dende as propias compa˜ n´ ıas un esforzo na recollida e modelizaci´ on de datos da calidade do aire, co fin de predicir futuros episodios de contaminaci´ on e as´ ı poder evitalos. A maior´ ıa dos m´ etodos de predici´ on de contaminantes atmosf´ ericos basease en modelos univariantes de autoregresi´ on. Sen embargo, estes modelos non te˜ nen en conta a correlaci´ on existente entre os distintos contaminantes. O feito de considerar esta informaci´ on nun modelo preditivo pode ofrecer predici´ ons m´ ais precisas. Nesta aplicaci´ on usamos os rexistros hist´ oricos dos contaminantes medidos nas inmediaci´ ons dunha industria electro-intensiva do norte de Espa˜ na para axustar un modelo de predici´ on de resposta bivariada. En concreto modelizamos as concentraci´ ons futuras do di´ oxido de xofre (SO2) e dos ´ oxidos de nitr´ oxeno (NOx), dous contaminantes atmosf´ ericos que amosan unha correlaci´ on alta. Por otra banda usamos datos correspondentes a un episodio de contaminaci´ on para a validaci´ on do mesmo. Estes datos de contaminaci´ on mostran unha alta asimetr´ ıa, con valores pr´ oximos a cero a maior´ ıa do tempo e valores extremos durante episodios de poluci´ on curtos. Tras o proceso de selecci´ on de variables, a predici´ on conxunta ofreceu resultados satisfactorios. Por tanto, a aplicaci´ on do modelo nestes datos complexos e cunha mostra de validaci´ on externa, demostra a xenerabilidade do modelo estat´ ıstico proposto e a s´ ua robustez na aplicaci´ on pr´ actica. No cap´ ıtulo 5 afondamos no problema de selecci´ on de variables destes modelos de regresi´ on bivariada. En concreto, propo˜ nemos un m´ etodo para contrastar a signficaci´ on estat´ ıstica do efecto dunha covariable sobre a forma da rexi´ on de referencia. Ademais, a formulaci´ on do m´ etodo orixinal extendeuse para incorporar un efecto de interacci´ on factor-por-rexi´ on (os detalles de estimaci´ on presentanse no Ap´ endice B). O estat´ ıstico proposto basease nunha distancia euclidiana bivariada entre a rexi´ on obtida polo modelo nulo, e baixo o alternativo. Neste contexto o modelo nulo considera un efecto aditivo das variables preditoras, mentres que o modelo alternativo considera unha interacci´ on entre elas. Sen embargo, o m´ etodo proposto poder´ ıa empregarse para contrastar outros efectos, tales como a inclusi´ on dunha nova variable, ou estimar un efecto nonlinear en lugar dun linear para unha covariable continua. A distribuci´ on do estat´ ıstico baixo a hip´ otese nula aproximouse mediante t´ ecnicas de remostraxe. O contraste proposto respecta os n´ ıveis nominais de significaci´ on no caso de hip´ otese nula. Ademais, ofrece unha curva de potencia satisfactoria, ´ e dicir, a medida que a magnitude do efecto da variable preditora aumenta, tam´ en o fai a porcentaxe de rexeitamentos da hip´ otese nula. Deste modo constatamos que ix a aproximaci´ on por t´ ecnicas de remostraxe da distribuci´ on do estat´ ıstico ´ e eficaz, e ademais mostra potencia. Esta ´ e a primeira proba estat´ ıstica proposta para contrastar o efecto de covariables sobre a forma dunha rexi´ on de referencia. Este m´ etodo aplicouse para contrastar un efecto factor por curva, fronte a un efecto aditivo na modelizaci´ on da distribuci´ on bivariada da (talla, peso) dunha mostra de infantes en funci´ on da s´ ua idade e sexo. Actualmente, a definici´ on de sobrepeso e obesidade en pediatr´ ıa basease no ´ Indice de Masa Corporal (IMC = peso/talla2). Especificamente, def´ ınense os puntos de corte diagn´ osticos a partir de curvas de referencia dependentes da idade e x´ enero para este ´ ındice. Se un infante est´ a por debaixo do cuantil 0.05 ou supera os cuant´ ıs 0.85 ou 0.90, considerase que padece infrapeso, sobrepeso ou obesidade respectivamente. A interpretaci´ on bivariada dos valores da talla e do peso, require por unha banda da estimaci´ on dunha rexi´ on de referencia condicional, e por outra dun contraste do efecto da idade e x´ enero para discernir se ´ e aditiva ou mostra interacci´ on. Co contraste proposto puidemos comprobamos que o efecto da idade e x´ enero sobre os valores conxuntos de (talla, peso) mostra interacci´ on. Especificamente, o efecto da idade ´ e distinto entre nenas e nenos a partir do inicio da puberdade. Ata os 13 anos o efecto da idade sobre medias e varianzas da talla e peso ´ e similar para ambos sexos, logo a partir desa idade a talla e peso medios dos nenos supera o das nenas. Ademais, a correlaci´ on entre a talla e o peso non ´ e constante coa idade, pois reducese cos anos. Os resultados obtidos mediante a aplicaci´ on da rexi´ on de referencia condicional compararonse coa clasificaci´ on baseada no ´ ındice de masa corporal, mostrando un exemplo interesante e visual, onde un modelo de regresi´ on de resposta bivariada pode ofrecer mellores resultados que unha alternativa univariante. A interpretaci´ on bivariada ofreceu novas perspectivas do efecto da idade e x´ enero sobre o peso e talla dos infantes, as´ ı como da relaci´ on entre ambos. O IMC s´ o ten en conta se un peso ´ e proporcional a talla ou non. Mentres que a rexi´ on de referencia espec´ ıfica para cada idade e sexo ofrece unha caracterizaci´ on total, identificando a maiores tallas e pesos que sendo proporcionais entre s´ ı, son at´ ıpicamente altos ou baixos. Finalmente no Cap´ ıtulo 6, esta tese abre unha nova li˜ na de investigaci´ on propo˜ nendo un m´ etodo de estimaci´ on de rexi´ ons de referencia. Este novo m´ etodo non asume unha dependencia linear entre as variables resposta, por tanto, ´ e esperable obter a´ ında mellores resultados en casos de distribuci´ ons complexas. Tomando como punto de partida a proposta de Klein et al. (2019) desenvolveronse unha serie de melloras encami˜ nadas a facilitar a aplicaci´ on dos Modelos Multivariantes de Transformaci´ on Condicional (MCTMs) a problemas de investigaci´ on con datos cl´ ınicos reais. A principal novidade dos MCMTs ´ e o seu m´ etodo de inferencia, que se basea na transformaci´ on m´ ais probable. Estes modelos en lugar de asumir unha distribuci´ on espec´ ıfica para modelizar un conxunto de datos transforman estes a unha distribuci´ on de referencia co˜ necida — habitualmente a gaussiana est´ andar. Esta transformaci´ on ofrece unha correspondencia entre as medidas caracter´ ısticas da distribuci´ on de referencia, e a distribuci´ on dos datos x observados. Deste xeito as medidas caracter´ ısticas da variable de interese poden obterse mediante unha funci´ on de transformaci´ on que se estima a partir dos datos. Esta funci´ on pode ademais estimarse de forma condicional a covariables no caso dunha variable multivariada. A estimaci´ on desta funci´ on de transformaci´ on basease na expresi´ on dos valores da variable resposta, e das covariables, en forma de bases polin´ omicas de Bernstein. O cal confire flexibilidade tanto na resposta como na forma dos efectos das variables preditoras continuas. Os par´ ametros do modelo estimanse mediante un estimador de m´ axima verosimilitude que permite considerar o efecto das covariables sobre a funci´ on de distribuci´ on das respostas marxinais e sobre a correlaci´ on entre as mesmas. A partires desta estrutura novidosa dos MCTMs propuxemos o m´ etodo de estimaci´ on das rexi´ ons de referencia condicionais. Este m´ etodo basease na existencia para datos gaussianos dunha definici´ on anal´ ıtica destas rexi´ ons. Aplicando pois a funci´ on inversa da transformaci´ on a estas rexi´ ons co˜ necidas, podemos obter unha estimaci´ on das mesmas para calquera resposta multivariada continua. Ademais, ´ o ser a transformaci´ on condicional a covariables, podemos obter esta rexi´ on para calquera valor das mesmas. Pola propia natureza do modelo, solventamos unha limitaci´ on da proposta anterior, a cal asum´ ıa un correlaci´ on linear entre as respostas. Neste caso, esta correlaci´ on pode ser non-linear e expresada en forma de ´ ındices tales como a rho de Spearman ou a tau de Kendall. Nos estudos de simulaci´ on obtivemos mellores resultados con este modelo en comparaci´ on con propostas anteriores. xi xii Contents List of Figures xvii List of Tables xxiii 1 Introduction 1 1.1 Thesis objectives .............................. 3 1.2 Thesis structure and achievements ................... 3 2 Theoretical background on flexible and multivariate regression 7 2.1 Flexible additive regression models ................... 7 2.1.1 Spline regression smoothers ................... 8 2.1.2 Kernel smoothers ......................... 13 2.1.3 Bernstein basis smoothers .................... 15 2.2 Regression models distributional assumptions ............ 17 2.3 Multivariate distributions analysis ................... 19 2.4 Multivariate regression models ..................... 23 2.4.1 Seemingly unrelated regression ................. 24 2.4.2 Conditional copula regression models ............. 25 2.4.3 Multivariate quantile regression ................ 26 3 Modelling conditional reference regions 29 3.1 Introduction ................................ 29 3.2 Conditional bivariate reference region model ............. 32 3.2.1 Model formulation ........................ 32 3.2.2 Estimation algorithm ....................... 33 3.2.3 Bivariate kernel bandwidth estimation ............ 35 3.3 Simulation Study ............................. 36 3.3.1 Estimation algorithm performance ............... 36 3.3.2 Reference region performance ................. 38 3.3.3 Bivariate kernel bandwidth performance ........... 43 3.4 Age-specific bivariate reference region for glycemic tests ...... 46 3.4.1 Motivating database ....................... 46 3.4.2 Conditional bivariate reference region estimation ...... 48 3.4.3 Software implementation .................... 50 xiii 5.6 Standarized 95% bivariate reference region compared with body mass index age dependent reference curves for males, and females. 78 5.7 Bivariate reference region (height, weight) real values characterization, compared with the WHO age x gender BMI cutpoints, and LMS 95% age x gender BMI reference curves. .............. 79 6.1 Graphical illustration of the estimation process to obtain the MCTMs’ multivariate reference regions. ..................... 87 6.2 True percentile curves for ⌧=0.05,0.50,0.95 (red line) with estimations mean (grey dashed line) and the corresponding 95% simulation intervals (blue dashed lines) for sample size n=200..... 89 6.3 Percentile curves estimation error for different sample sizes, and ⌧svalues for both estimation methods (MCTMs and quantreg). ... 90 6.4 Bivariate reference region estimation error for three methods (tolerance,refreg, MCTMs), for different sample sizes, ⌧=0.95, and two simulation scenarios. ........................ 91 6.5 Theoretical regions along with 1000 reference regions estimated using MCTMs, for different sample sizes, ⌧=0.95, and the second simulation scenario. ............................ 92 6.6 Conditional reference region estimation error for ⌧=0.95, and several basis orders for the response, and covariate effects. ..... 94 6.7 Glycemic markers’ bivariate distribution change with age for both genders. .................................. 95 6.8 Fasting plasma glucose and glycated hemoglobin age-dependent conditional cumulative distribution functions for both genders. .. 96 6.9 Fasting plasma glucose, and glycated hemoglobin age adjusted percentile curves, for men (red) and women (blue). ......... 96 6.10 Fasting plasma glucose, and glycated hemoglobin, correlation change with age, for both genders along with 95% pointwise confidence intervals. .................................. 97 6.11 Fasting plasma glucose, and glycated hemoglobin, reference region change with age for both genders. ................ 98 C.1 Estimated effects of age with the 95% bootstrap pointwise confidence interval on the FPG and HbA1c mean, variance models, and on their correlation. .........................126 C.2 Estimated region in the bivariate residuals scale for healthy patients (left), and patients with diabetes (right). ............127 C.3 Predicted reference regions for different ages. .............130 C.4 Two different representations of a pollution incident for SO2and NOx.....................................131 C.5 Estimated bivariate prediction region for a pollution episode. ...132 C.6 Scatter plot for three glycemic markers, with colour scale depending on age value. .............................134 xx C.7 Trivariate standarized, and conditional trivariate region for FPG, HbA1c, and Fr. ...............................135 xxi xxii List of Tables 3.1 Estimation algorithm bias and root mean square error for covariates effects on model components. ................... 38 3.2 Percentage coverage of the estimated conditional bivariate region for different sample sizes (500, 1000 and 2000), and nominal levels (5%,50%,90% and, 95%) in the three simulation scenarios. ...... 41 3.3 Coverage probability of bivariate data points for different sample sizes (500,1000 and 2000), covariate X1values, with X2fixed at zero, and bivariate kernel bandwidths estimators. Best-Coverage represents our equation (3.10) proposal, CV is the least-square crossvalidation method, and 5Han arbitrary large bandwidth. ..... 44 3.4 Percentage of test results contained within the estimated bivariate reference regions for different nominal levels. Cross validation evaluation refers to a leave-one-out cross validation evaluation. Coverage probability is presented for the entire dataset (Global) and for three age groups. ......................... 50 4.1 Scheme of the fitted models in the second simulation scenario. All covariates entered the models in the location and variancecovariance parameters. ......................... 58 4.2 Estimated coverage (ˆ⌧%) and L⌧value of the prediction regions obtained with each selected model of size qfor different values of ⌧. The cross Xindicates the covariates included in each model. ... 63 5.1 Estimated type I error (in per cent) at different significance levels (5, 10, 15, and 20 per cent) and different sample sizes (n=300, 500, and 1000). ................................. 74 6.1 MCTMs reference region data coverage for different sample sizes, covariate values, and two simulation scenarios. Coverage evaluation was performed in an out-sample design. ............. 93 6.2 Leave-one-out cross validation coverage evaluation. ......... 98 xxiii xxiv Chapter 1 Introduction Problems involving correlated response variables are common in the medical setting. By way of a simple example, a person’s weight cannot be judged healthy or unhealthy if no information is available about that person’s height, and if we are to estimate the effect of age on body size, information on both these variables would need to be examined simultaneously. Analogous problems are encountered in diabetes, liver and thyroid disease, the diagnosis of which is based on a number of correlated continuous markers. In such scenarios, a regression model involving multiple response variables is useful since the correlations between responses can be taken into account. This is an improvement over the use of several independent univariate models since the effects of covariates on the response variables’ association structure and joint distribution are characterized. This enriches our understanding of the effects of different clinical markers, offers a means of arriving at a ’multivariate interpretation’ of the results of continuous diagnostic tests, and ultimately allows for the formulation of better diagnostic criteria. Multivariate regression could provide new insights into the onset and progression of diabetes, a chronic disease characterized by high blood glucose concentrations. Given that early intervention can prevent or delay its appearance, it is crucial to identify high-risk individuals. To do so, general practitioners measure several markers of glycemia that are elevated in patients with glucose homeostasis dysfunctions. Yet, while these markers are highly correlated, their concentrations are usually interpreted separately. Hence, few clues are captured on how clinical variables might affect the association structure and joint distribution of the measured markers. Take for instance the two most commonly used glycemic markers, namely fasting plasma glucose (FPG) and glycated hemoglobin percentage (HbA1c). Figure 1.1 shows the FPG and HbA1c values for three age groups (20, 40 and 70 years) in a healthy population. It can be appreciated graphically how the shape of the markers’ joint values changes with age, and how the association between FPG and HbA1c does not remain constant. From a medical viewpoint, it would be of great interest to estimate 1 2CHAPTER 1. INTRODUCTION the effect of age on the joint distribution of (FPG, HbA1c), characterizing which bivariate values are more likely to be observed for each age group. This would reveal which joint values are “normal” for each age group, and which are “abnormal”, i.e., patients with glycemic dysregulation could be identified. In other words, the reference intervals (those containing 95% of healthy patients’ results) for each marker are extended into a bivariate setting. 60 80 100 120 140 160 180 4 5 6 7 20 years FPG, mg/dL HbA1c, % 60 80 100 120 140 160 180 4 5 6 7 40 years FPG, mg/dL HbA1c, % 60 80 100 120 140 160 180 4 5 6 7 70 years FPG, mg/dL HbA1c, % Figure 1.1: Fasting Plasma Glucose (FPG) and glycated hemoglobin (HbA1c) joint values for different age groups. This data were provided by the A-Estrada Glycation and Inflammation Study (AEGIS) project (Gude et al.,2017). Figure 1.2: Bivariate reference region (solid line) compared to univariate reference intervals (dashed lines) for independent and correlated bivariate data. 1.1. THESIS OBJECTIVES 3 Indeed, when there are several correlated markers, a reference region - a convex hull containing 95% of healthy patients’ multivariate results - can be constructed. A reference region is desirable over the use of several reference intervals for two reasons. Firstly, it avoids making separate interpretations (each of which may induce further error into the overall interpretation), thus reducing the number of false positive, resulting in higher specificity. Secondly, higher sensitivity is achieved since univariate methods ignore atypical combinations of markers. Figure 1.2 shows, in a simulated data example, the advantages of using a reference region over two reference intervals. When these two measurements are independent, both univariate and bivariate methods of interpretation can suggest the same conclusion. However, for correlated measurements, reference intervals ignore the atypical combinations of both markers. In other words, a value can be univariately ’normal’ but bivariately ’abnormal’. Identifying these multivariate atypical values will have clinical implications for patients. The literature, however, contains few proposals regarding how to estimate conditional reference regions, and there is no software that makes use of them in clinical medicine or medical research. The resulting need motivated the proposals discussed in Chapters 3-6. The proposed models could be used with any disease for which a diagnosis is based on several, correlated, continuous variables, and indeed in any research field that deals with multivariate responses. 1.1 Thesis objectives This thesis proposes and evaluates new regression models for interpreting multivariate continuous responses. Quantile regression extensions for bivariate continuous response variables are suggested. Quantile regression is usually used to better characterize the effects of covariates on a response variable over its total range. The extension of this idea to the bivariate setting might offer a better understanding of the association between correlated clinical variables. Software implementing the proposed models, and their application in real clinical problems, medical research, and in other fields, is presented. 1.2 Thesis structure and achievements Chapter 2: Theoretical background on flexible and multivariate regression Chapter 2offers an introduction to the statistical concepts used in Chapters 3-6. It focuses on flexible regression models based on spline functions, polynomial kernel and Bernstein basis smoothers. Flexible regression formulations beyond conditional means are introduced in the form of location-scale models. Finally, multivariate data characteristics and current multivariate regression models are introduced (parametric and non-parametric). 4CHAPTER 1. INTRODUCTION Chapter 3: Modelling conditional reference regions In Chapter 3a conditional reference region estimation method is proposed based on a bivariate regression model with location-scale effects. The estimation of the region relies upon a bivariate kernel density estimator, allowing covariate effects to be taken into account via the use of flexible additive predictors. A major novelty of this proposal is the lack of any parametric assumption regarding the response variable. The proposed methodology is validated in a simulation study that examines; i) the estimation error of the flexible additive model, ii) the general performance of the reference region, and iii) the effect of kernel bandwidth on the model fit. Finally, the method is applied to a real research problem, offering for the first time a joint interpretation for two correlated glycemic tests, depending on patient age. The contents of this chapter have been published in Statistics in Medicine journal (Lado-Baleato et al.,2021). Chapter 4: Variable selection for prediction regions The proposed model is here used with a real data problem in the environmental setting. Here, the reference region is used as a measure of uncertainty regarding joint predictions for two atmospheric pollutants. Specifically, the future concentrations of two correlated air pollutants SO2and NOxare predicted using the previously recorded concentrations of both as predictor covariates. Since data for several previous concentrations taken at different time points are available, the most informative need to be selected when constructing the final predictive model. For this purpose, a bivariate loss function is proposed (the usefulness of this function was not improved by the further inclusion of non-informative predictor variables). Two time points (values for the previous 45 and 60 minutes) were seen to be enough for predicting joint (SO2,NOx) values one hour in advance. This chapter is an extension of a co-authored paper that appeared in Stochastic Environmental Research and Risk Assessment (Roca-Pardi˜ nas et al.,2021). Chapter 5: Testing covariate effects on bivariate reference regions Chapter 5extends the conditional reference region estimation framework by proposing a statistical test examining the effects of the covariates. The proposed test is based on a bivariate Euclidean distance between the region obtained in the null model (with no covariate), and a competing alternative model (including the covariates’ effects). The null hypothesis distribution is approximated by bootstrap resampling. As a secondary aim, the model formulation presented in Chapter 3is extended by incorporating factor-by-curve interactions using penalized spline smoothers. A simulation study was performed to assess the test’s statistical power and type I error. Finally, an illustrative application modeling the effect of age and gender on children’s body height and weight is presented. The bootstrap-tested reference 1.2. THESIS STRUCTURE AND ACHIEVEMENTS 5 region proves that age has a different effect on the bivariate distribution shape of height and weight in males and females. Specifically, the age x gender model shows how the bivariate distributions of both anthropometric measurements overlap until puberty. Moreover, this novel interpretation allows four atypical combinations for body height and weight to be identified (whereas body mass index detects only obesity and underweight). Chapter 6: Multivariate reference regions based on conditional transformation models Chapter 6explores a more flexible model for estimating conditional reference regions. Even though the previous model has no parametric restrictions, it assumes a linear correlation between the response variables. This is associated with a larger estimation error in the case of responses showing non-elliptical structures of dependence. A more general and elegant formulation for reference region estimation is here presented via the use of multivariate conditional transformation models (Klein et al.,2019). Simulation studies showed that multivariate conditional reference regions are associated with smaller errors than are those resulting from the use of existing methods for non-Gaussian heterocedastic scenarios. The new method proved to be feasible with real data, providing age-gender specific reference regions for two glycemic tests. Chapter 7: General discussion This chapter closes the thesis by offering a general discussion on the presented results and by outlining future lines of research. Appendices Specific details on the models discussed in Chapters 3,4and 5, and their software implementation, are given in three additional appendices. Appendix A contains additional details on additive predictor estimation by means of polynomial kernel smoothers. Appendix Bcontains a modification of the aforementioned algorithm, incorporating factor-by-curve effects by means of penalized splines smoothers. Finally, Appendix Ccontains a user guide for the refreg R package (an R package implementing the statistical models developed). 12 CHAPTER 2. THEORETICAL BACKGROUND data overfitting. The main advantage of this strategy is that smoothness does not depend on the number, or position, of knots but on the smoothing parameter choice. The second derivative of B-splines might be used to take into account the curvature of the estimated function. More smoothly functions might be obtained by incorporating in the estimator the following penalty term Z(f00(X))2dx Spline basis first derivative is given by first order differences of the corresponding coefficients vector (see Formula (2.4)). Similarly, higher order derivatives can be obtained by means of rorder differences. By including this penalty term a least square estimator for f(X)can be written as: PLS()= n X i=1 Yi d X j=1 jBj(Xi)!2 + d X j=r+1 (0j)2 where 0denotes rth-order differences between adjacent regression coefficients. 0.0 0.2 0.4 0.6 0.8 1.0 −2−1 0 1 2 B−spline basis x y 0.0 0.2 0.4 0.6 0.8 1.0 −2−1 0 1 2 P−spline basis x y Figure 2.3: Comparison between B-spline and P-spline estimated coefficients. Adapted from (Durb´ an,2009). The rorder differences penalties can be written in matrix notation recursively. Given the first order differences as: 2.1. FLEXIBLE ADDITIVE REGRESSION MODELS 13 D1=0 B B B @ 11 11 ...... 11 1 C C C A(d1)⇥d Higher order differences can be obtained from D1following Dr=D1Dr1. Under this matrix notation, we can express the penalty term as:  d X j=r+1 (rj)2=>D> rDr=>Kr where Kr=D> rDr. Thus, the PLS estimator is given by: PLS()=(YZ)>(YZ)+>K Smoothing parameter choice is a big topic on the estimation of unknown non-linear effects by means of P-splines (Ruppert et al.,2003;Lee,2003). An automatic choice of this parameter can be obtained through minimizing several statistical criteria (e.g., GCV, AIC, UBRE, REML). For instance, the generalized cross validation (GCV) estimator for can be defined as ˆ GCV =argmin  n X i=1 (Yiˆ Yi)2 ntr(H) where tr(H)represents the trace of the matrix H=B(B>B+D>D)B>. 2.1.2 Kernel smoothers Under a spline regression model, continuous covariates non-linear effects was estimated from a global regression formulation. In this section, we focus on a regression smoother that approximates f(X)locally. This local approximation is based on a Taylor expansion of the function fat an specific point Xias f(Xi)⇡f(X)+(XiX)f0(X)+(XiX)2f00(X) 2! +...+(XiX)qf(q)(X) q!(2.5) where qis the polynomial degree. The derivatives f(q)(X) q!might be seen as a set of unknown parameters 1,..., q. Then, the expression (2.5) turns into a linear regression problem f(X)⇡0+(XiX)1+(XiX)22+...+(XiX)qq 14 CHAPTER 2. THEORETICAL BACKGROUND 0.0 0.2 0.4 0.6 0.8 1.0 02468 x y h = 0.30 0.0 0.2 0.4 0.6 0.8 1.0 02468 x y 0.0 0.2 0.4 0.6 0.8 1.0 02468 x y 0.0 0.2 0.4 0.6 0.8 1.0 02468 x y h = 0.10 0.0 0.2 0.4 0.6 0.8 1.0 02468 x y 0.0 0.2 0.4 0.6 0.8 1.0 02468 x y 0.0 0.2 0.4 0.6 0.8 1.0 02468 x y h = 0.02 0.0 0.2 0.4 0.6 0.8 1.0 02468 x y 0.0 0.2 0.4 0.6 0.8 1.0 02468 x y True effect Estimated Effect Local regression Kernel at xi Figure 2.4: Polynomial kernel smoother estimation performance depending on bandwidth value. Adapted from Garc´ ıa-Portugu´ es (2021). This approximation is only valid locally, that is for covariate values close to an specific point (Yi,X i). A local estimator can be obtained by including a weighting function in least square estimation formula (Wand and Jones,1994), as: ˆ =min  n X i=1 Yi q X j=0 j(XiX)j!2 Wh(X,Xi)(2.6) where the weights Wh(X, Xi), are derived from the kernel function Kdefined by Wh(X,Xi)=K✓XiX h◆with K(u)= 1 p2exp(1 2u2)(2.7) The core idea is that the closer a point Xis to the specific data point Xithe higher influence it has in the model estimation. Note that the equation (2.6) with q=0corresponds to the well known Nadaraya-Watson regression model 2.1. FLEXIBLE ADDITIVE REGRESSION MODELS 15 (Nadaraya,1964;Watson,1964). The kernel bandwidth h, in equation (2.7), plays a key role in the final performance of local kernel smoothers. In Figure 2.4, we depict a graphical illustration of this parameter effect on the regression function shape. As can be seen, very low hcorrespond to more local estimations resulting in data overfitting. On the other hand, large hcorresponds to a global polynomial regression estimator, thus ignoring the true shape of the non-linear effect. Bandwidth selection problem can be solved by several methods (see e.g., (K¨ ohler et al.,2014)). For instance, a cross-validation estimator for hmight be written as: ˆ hCV =argmin h 1 n n X i=1 (Yiˆ fi(Xi;q;h))2 where Yiis an specific observed value, and ˆ fi(Xi;q;h)represents the estimation of the smooth unknown effect without the idata point. 2.1.3 Bernstein basis smoothers In most practical applications, penalized splines, or polynomial kernel smoothers, are the preferred choice for estimating a non-linear association between a continuous predictor, and the response variable. Nevertheless, in specific data problems there exist subject-matter information about the regression function shape that must be incorporated into the model fit. For instance, if it is previously known that the regression function must be monotonically increasing, or strictly convex or concave. In such cases, the non-parametric estimator must satisfy a set of shape restrictions on the predictor variable support. This estimation problem is denominated in the statistical literature as nonlinear restricted shaped, or constrained, regression. Bernstein basis are among the preferred methods for estimating restricted shaped regression, because it is relatively easy to incorporate shape constraints under this model specification. Stadtm¨ uller (1986) introduced Bernstein basis for estimating a continuous covariate non-linear effect. This method was then extended for multiple predictor variables by (Tenbusch,1997). Moreover, this author proved the Bernstein regression uniform consistency, asymptotic normality, and that its estimation error was similar to kernel methods. More recently, Bernstein regression was extended, and used, by many authors (e.g. (Brown and Chen,1999;McKay Curtis and Ghosh,2011;Wang and Ghosh,2012)). A Bernstein polynomial of degree Mis defined as a number-to-vector function where the m-th output is: bm(X)=✓M m◆Xm(1 X)Mmm=0,1,···,M This function defines a mapping of the form b(X):[0,1] ![0,1]M+1. Thus, 16 CHAPTER 2. THEORETICAL BACKGROUND the predictor variable must be transformed to lie in the unit [0,1]. In Figure 2.5, Bernstein basis are depicted for several polynomial degrees M. As can be seen, a Bernstein basis is composed by M+1polynomials showing symmetry around the covariate midpoint (i.e. bm k(X)and bm mk(X)are mirror images of each other about X=0.50). One of the many remarkable properties of Bernstein polynomials is that all of the derivatives possess same convergence properties (see Farouki (2012)). This can be useful in determining appropriate shape restrictions for a particular application. 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Bernstein Basis (m = 3) x y 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Bernstein Basis (m = 7) x y 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Bernstein Basis (m = 10) x y Figure 2.5: Bernstein basis for polynomial degrees (m=3,7,10). −1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 x y True effect P−spline Restricted−shaped Figure 2.6: Bernstein regression with monotonic increasing restriction, compared with a penalized spline regression. In the regression context, for the observed data pairs {Yi,X i}n i=1, and a nonlinear model y=f(X)+✏. The smooth function f(X)might be represented in a Bernstein polynomial basis as 2.2. REGRESSION MODELS DISTRIBUTIONAL ASSUMPTIONS 17 f(X)= m X j=0 bj(X)j Chang et al. (2007) proved that if 01... m, then f0(X)0(i.e. monotonically increasing) for every X2[0,1] (see Figure 2.6). This restriction is very convenient for Chapter 6, because the transformation function h(used for model estimation and inference) must be monotonically increasing. 2.2 Regression models distributional assumptions Most statistical methods dealing with continuous data are based on the variable mean; a single number describing the center of distribution. Similarly, regression models focuses on estimating the response conditional expectation E[Y|X], and treat higher moments of the response distribution, e.g. variance or skewness, as constant parameters. This feature is shared by many regression formulations such as the simple linear models (LMs), generalized linear models (GLMs), or generalized additive models (GAMs). These distributional assumptions implies that for a model like, for instance Y=f(X)+"with "2N(0,2), we suppose that the response conditional distribution is given by Y2N(f(X),2). Thus, for Gaussian homocedastic data it is possible to obtain covariates effects not only on the mean but on the entire response distribution by means of the so called prediction intervals. For the aforementioned model, a prediction interval, also known as Z-score, can be defined as: Q⌧(Yi|Xi)= ˆ f(Xi)+2z⌧ where z⌧is the quantile from the distribution function N(0,2). These prediction intervals for ⌧=0.025,0.975 are depicted in Figure 2.7. Nonetheless, while for Gaussian distributed data the mean is equidistant with both distribution extremes, skewed or heterocedastic distributions are poorly described by this parameter. As a consequence the model inference, but also the information that can provide about the response conditional distribution, are not entirely reliable in such cases. Specifically, the model’s confidence, and their prediction, intervals show a poor performance in any real data application showing a non-Gaussian or heterocedastic distribution. Additionally, mean regression paradigm might offer but an incomplete characterization of the real influence that an specific covariate have on the response. For instance, in the medical case, the “mean patient” is likely to be the “healthy patient”, while those placed on the distribution extremes are likely to have some undiagnosed medical condition. Thus, a conditional distribution model is more interesting, informative, and flexible in clinical practice. 18 CHAPTER 2. THEORETICAL BACKGROUND 0.0 0.2 0.4 0.6 0.8 1.0 −1.5 −0.5 0.5 1.5 Conditional expectation x y 0.0 0.2 0.4 0.6 0.8 1.0 −1.5 −0.5 0.5 1.5 Prediction intervals x y Figure 2.7: Conditional mean regression model fit for homocedastic Gaussian data. In the right plot the model prediction intervals are depicted in red dashed lines for ⌧=0.025,0.975. Conditional distribution regression models was firstly approached by Koenker and Bassett Jr (1978). In this work, the classical least square estimation was replaced by a loss function to be minimized. In turn, it offers the estimation of the response conditional quantiles, imposing non-parametric assumptions. More recently, Cole and Green (1992) and then Stasinopoulos and Rigby (2007) proposed a distributional regression procedure based on the estimation of covariates effects not only on the response mean, but also on their higher moments such as variance. For that reason, these models are usually denominated location-scale models, if conditional variance are estimated, or location-scale and shape models if covariates effects on kurtosis, and skewness are also considered. After estimating covariates effects on each response parameter, and under an specific parametric distribution, the response distribution change with covariates can be easily reconstructed. A major advantage of these methods is that offers congruent conditional quantiles, avoiding the quantile crossing problem, very common in classical quantile regression. Given Lcategorical, and Kcontinuous predictor variables, a location-scale model for a continuous response Y2D(µ, )following a distribution Ddefined by the parameters µ2Rand 0, can be written as: (µ(X)=↵µ+PL l=1 µ lXl+PK k=L+1 fµ k(Xk) 2(X)=G⇣2 0+PL l=1 2 lXl+PK k=L+1 f2 k(Xk)⌘(2.8) where land fkrepresent respectively linear and non-linear effects. While Gis a 2.3. MULTIVARIATE DISTRIBUTIONS ANALYSIS 19 link function (e.g. exp(·)) which ensures that the additive predictor will always offer predicted variances higher than zero. From equation (2.8) model fit, conditional quantiles for the response can be easily computed as: Q⌧(X)=µ(X)+(X)"⌧ where "⌧is the ⌧quantile of the response distribution D(Kneib,2013). The improvement of this model upon classical regression is shown in Figure 2.8. 0.0 0.2 0.4 0.6 0.8 1.0 −2−1 0 1 2 3 Conditional expectation x y 0.0 0.2 0.4 0.6 0.8 1.0 −2−1 0 1 2 3 Location−Scale Model x y Figure 2.8: Comparison between prediction intervals obtained for conditional mean and a location-scale model. The prediction intervals are represented by red dashed lines, conditional means as black solid lines. In Chapter 3we extend the location-scale regression idea to a continuous bivariate response, but imposing no parametric restrictions. In this model the response mean vector, and its variance-covariance matrix, are modeled through flexible additive predictors taking into account continuous covariates non-linear effects. Before presenting this model, it is necessary to introduce multivariate data characteristics, and the already existing statistical techniques for modeling such data. 2.3 Multivariate distributions analysis Multivariate data is common in practical research. Whenever the value of several variables is recorded simultaneously for the same individual, we are dealing with this type of data. In medical studies, multivariate variables are, in a sense, the natural state of affairs. Each time a patient is consulted, several serological, anthropometric and psycho-social variables are taken for each of them. 20 CHAPTER 2. THEORETICAL BACKGROUND However, these variables information are usually studied side-by-side, and their multivariate distributions are rarely estimated in practice. At the beginning of this thesis project we stated our interest in estimating covariate effects on the continuous bivariate distribution of two correlated glycemic tests. Thus, a relevant issue to be addressed in this section is how the joint distribution of continuous multivariate variables are described. In the univariate context, distribution and density functions are characterized by a location, and a variance parameter. In a multivariate representation, the means and variances of the separate measures (known as marginal parameters) are relevant, but it is also necessary to take into account the dependence between the different variables. This dependence structure is usually measured by covariance, or its scaled version known as correlation coefficient. In the following subsections, we will illustrate the main techniques for describing the multivariate data joint distributions. Parametric methods Similarly to the univariate case the well-known Gaussian distribution are commonly applied for multivariate data description and inference. A multivariate variable Yof dimension dis said to follow a d-variate Gaussian (or Normal) if its density function is given by the expression: f(Y|µ,⌃)=(2⇡)r/2|⌃|1/2e1 2(Yµ)>⌃1(Yµ),Y2Rd where µ=E(Y)is the mean d-vector, and ⌃=Cov(Y)=E⇥(Yµ)(Yµ)>⇤ a positive-definite, symmetric (d⇥d)variance-covariance matrix. This matrix contains each variable variability, and the pairwise dependence between them. 115 120 125 130 135 140 145 150 20 30 40 50 60 Height, cm Weight, Kg Observed values Mean vector Guassian contours Figure 2.9: Children body height and weight measurements with density contours according to a Gaussian distribution assumption. 2.3. MULTIVARIATE DISTRIBUTIONS ANALYSIS 21 As an illustrative example for a d=2dataset – comprising 150 children’s height and weight records – we have a bivariate variable with the following mean vector and variance-covariance matrix ˆ µ=✓129cm 30Kg ◆and ˆ ⌃=69 50 50 67  Density contours drawn from these estimated parameters are depicted in Figure 2.9. Gaussian contours are elliptical-shaped symetric around the (height, weight) mean vector. From Gaussian distribution expression is possible to estimate tolerance regions containing an specific percentage of multivariate data – a method proposed by some physicians to derive reference regions (Boyd and Lacher,1982). Nonetheless, we can see that the Gaussian distribution fails to recover our dataset shape in the bottom-left corner. This is very common in practice, where the elliptical shaped Gaussian contours do not capture asymmetric structures of dependence. 110 120 130 140 150 10 20 30 40 50 60 Height, cm Weight, Kg Figure 2.10: Children body height and weight measurements with density contours according to a parametric copula representation. More flexible models for our bivariate data can be obtained by applying parametric copula functions. Since its proposal by Sklar (1973) copulas have become very popular in the multivariate statistical modeling because they can estimate a multivariate distribution in a modular fashion. Specifically, copula functions allow to join any parametric marginal distribution, using different parametric copulas. By varying the parametric copula we can design a bivariate distribution allowing different association structures (e.g. elliptical, upper-tail, lower-tail). For a d-dimensional random variable, with marginal distributions F1,...,F ddefined by a set of parameters (µ1, 1),...,(µd, d). Then there exists a copula, C, such that 28 CHAPTER 2. THEORETICAL BACKGROUND Chapter 3 Modelling conditional reference regions In this chapter we introduce a new statistical method for estimating bivariate conditional reference regions. Such regions can define which values are likely to be located in the most inner part of the joint distribution of two continuous response variables depending on predictor covariates. Therefore, they allow to identify excentric values of the response variable which could be labeled as atypical observations. This idea is useful in clinical primary care in order to define new diagnostic rules based on the joint values of two continuous diagnostic tests taking into account patient characteristics. A conditional reference region might be seen as the extension of the percentile curve ideas to the bivariate setting. Our statistical model counts with no parametric restrictions for the response, and covariates effects on the responses’ means and on their variance-covariance matrix might be estimated through flexible additive predictors. The remainder of this chapter is organized as follows. Section 3.1 gives background on the role of conditional reference regions in diabetes research. Then, Section 3.2 introduces the model structure and its estimation details. Section 3.3 evaluates the conditional reference region estimation, contemplating scenarios in which the model could be partially miss-specified. Section 3.4 makes use of the model for estimating an age-specific bivariate reference region for the diabetes markers FPG and HbA1c. Moreover, it gives some clues about the software implementation of our model. Finally, full details about flexible additive models estimation can be found in Appendix A. 3.1 Introduction The diagnosis and treatment of disease commonly rests on the results of measureable biomarker-based clinical laboratory tests. Indeed, some 70% of the decisions made in clinical practice are taken based on such results (Hallworth,2011). 29 30 CHAPTER 3. CONDITIONAL REFERENCE REGION For every test result, the clinical laboratory provides comparator values to help the clinician understand in context the information provided. These comparator values are often referred to as the reference interval (Siest et al.,2013) they usually reflect the range of values within which 95% of the results of the normal healthy population falls. When a single biomarker is examined, the reference interval is classically obtained using quantile estimation techniques (Wright and Royston,1999), or by conditional quantile regression (Koenker and Bassett Jr,1978) if any variable modifies the distribution of the response variable (i.e., the reference curve (Cole and Green,1992)). However, there is often more than one test for diagnosing a disease. For instance, the diagnosis of diabetes may be based on plasma glucose criteria, such as fasting plasma glucose (FPG) or the 2 h plasma glucose value obtained during a 75 g oral glucose tolerance test, or the glycated haemoglobin (HbA1c) test (American Diabetes Association,2019). The same tests may be used to screen for, diagnose and monitor the effects of treatment for diabetes. However, measuring glycaemic control is not foolproof; its clinical usefulness is affected by a number of biological and analytical factors. Disagreement between glycaemic control measurements are common, and clinicians need to know what might explain them (Sacks,2011). Following the guidelines of the American Diabetes Association, the diagnosis of diabetes is defined as FPG levels of 126 mg/dL, and of HbA1c6.5%. The same cut-offs are used for children, adolescents and adults (American Diabetes Association,2019). FPG and HbA1c levels are reported to be strongly correlated (Aleyassine et al.,1980), both in members of the general population and in patients with diabetes (Van’t Riet et al.,2010;Ramachandran et al.,2012). If two tests are indifferently used for the diagnosis of a disease, the correlation between them should be strong. Therefore, when diagnosing a disease using two markers, it might be reasonable to estimate their combined multivariate reference region instead of the reference interval for each test. The idea of combining reference regions for two or more laboratory tests has been discussed in the biomedical (Boyd and Lacher,1982;Harris et al.,1982;Slotnick and Etzell,1990) and statistical literature (Dong and Mathew,2015). However, the proposals made have required responses that follow a multivariate Gaussian distribution, condition that are not fulfilled by many markers used in the clinical setting. For instance, FPG and HbA1c concentrations both show a skewed distribution. Hence, a more general method for estimating multivariate reference regions is needed. Moreover, since HbA1c increases with age even after adjusting for glucose levels (Davidson, 1979;Pani et al.,2008;Espasand´ ın-Dom´ ınguez et al.,2019), this variable should be taken into account when establishing cut-offs or reference values. The literature reports but a few attempts to define reference regions for nonGaussian multivariate responses. Non-parametric reference regions for multivariate responses can be estimated using multivariate quantiles. However, there is no single definition of what a multivariate quantile is (Serfling,2002). Most 3.1. INTRODUCTION 31 current definitions are based on a centre-outward ordering of the data points (Chaudhuri,1996;Chakraborty,2001), defining a convex hull in which a proportion of the more central data points falls. Halfspace depth bivariate quantiles based on directional projections have recently been extended to the regression setting (Hallin et al.,2010) for estimating bivariate contours conditioned by covariates. These conditional bivariate quantile models have been used in the study of anthropometric characteristics affected by age (Geraci et al.,2020). However, the clinical interpretation of the results is not clear. Easily interpretable conditional bivariate quantile contours can be estimated using the Wei model (Wei, 2008), which is based on the estimation of: i) the marginal stratified conditional quantile regression for each response, and ii) bivariate quantiles using simulated data points from the above marginal stratified conditional quantile regression. The main drawback of this proposal is the influence of the univariate quantile regression outcome on the performance of the final bivariate quantile contour. Another non-parametric conditional method for detecting extreme combinations of two variables has also been proposed (Petersen,2009). This does not define a reference region, but four bivariate reference curves for detecting four possible atypical combinations of two variables. However, this alternative returns a higher false discovery rate for the reference region. Non-Gaussian bivariate regions for detecting joint outliers may be estimated using conditional copula regression models (Patton,2006). This was explored by Stander et al. (2019) to identify children with abnormal vision. Copula regression models allow a bivariate distribution to be constructed from a modular perspective, combining two univariate parametric distributions and the parametric copula joining them (Sklar,1973). Finally, the response parameters are made dependent on the covariates using flexible additive models (Klein and Kneib, 2016). Given copula regression model structure, and the effect of the estimated covariate on the means, variances and correlation of the responses, a conditional bivariate region for exploring the atypical combinations of two measurements can be estimated. The parametric representation of these models means several choices have to be made in the model building process, which complicates the use of this methodology when dealing with real data problems. The present chapter proposes a regression model for estimating a conditional bivariate reference region. The reference region is defined as the density function contour level which contains the bivariate data points with a given probability depending on the covariates. Unlike the existing methods, our statistical model places no parametric restrictions on the response, and the non-linear effects of continuous covariates may be estimated using local polynomial regression smoothers. Our proposal is an extension to bivariate data of a previous work (Mart´ ınez-Silva et al.,2016), where the authors used a location-scale model to estimate univariate percentile curves. The final performance of our conditional reference regions depends heavily on a bivariate kernel density estimator. In this work, we propose a new method for choosing the kernel bandwidth matrix in or- 32 CHAPTER 3. CONDITIONAL REFERENCE REGION der to obtain a reference region with a coverage of the bivariate data points close to the nominal level. Reference regions performance was succesfully tested using gaussian and non-gaussian simulated data. Finally, the proposed statistical methodology was applied to study jointly the FPG and HbA1c results, depending on age, in a cohort of nominally normoglycaemic persons. 3.2 Conditional bivariate reference region model This section presents a regression model for estimating a bivariate reference region conditioned by a set of covariates. The proposed model has no parametric restriction, and is based on the estimate of a location-scale regression, taking into account the effect of the covariate on the correlation between the response variables. 3.2.1 Model formulation Let X=(X1,...,X p)be a vector of p covariates, and let Y=(Y1,Y 2)be a continuous bivariate response of interest. In this context, the aim is to obtain a bivariate region ⌧of Yconditioned by the covariates Xdenoted as R⌧(X), and containing the 100⌧%of the bivariate data points. The following structure is assumed: Y=✓µ1(X) µ2(X)◆+⌃1/2(X)✓"1 "2◆(3.1) where µ1and µ2represent the response means, and the variance-covariance matrix is given by: ⌃(X)=✓2 1(X)12(X) 12(X)2 2(X)◆ The bivariate residuals ("1," 2)are assumed to be independent of the covariates and to have a mean of zero, zero unit variance, zero correlation, and an unknown density function f("1," 2). Note that, ⌃1/2(X)represents the Cholesky decomposition of the variance-covariance matrix ⌃(X)so that, ⌃1/2(X)(⌃1/2(X))T= ⌃(X). Thus, for any given Xthe bivariate region for (Y1,Y 2)is given by: R⌧(X)=✓µ1(X) µ2(X)◆+⌃1/2(X)R⌧for ⌧2[0,1] (3.2) where the R⌧is the unconditional bivariate region containing the 100⌧%of the model residuals ("1," 2)defined as: R⌧={(u, v)2R2|f(u, v)k}for ⌧2[0,1] where kvalue is chosen so that P(("1," 2)2R⌧)=⌧for ⌧2[0,1]. 3.2. CONDITIONAL BIVARIATE REFERENCE REGION MODEL 33 In equation (3.1), the response parameters (µ1,µ 2,2 1,2 2, 12)are related to the covariates vector Xvia additive predictors and known link functions G, which ensure that the restrictions on the parameter spaces are maintained. The following additive predictors are considered: µr(X)=↵r+ p X j=1 fjr(Xj)and 2 r(X)=G r+ p X j=1 gjr(Xj)!for r=1,2 (3.3) where fjr and gjr are smooth unknown functions, ↵rand rare intercept coefficients and the link function G= exp(·)to ensure that 2 r(X)0. Finally, the following additive structure is assumed for the responses’ association with one another, expressed as a linear correlation coefficient (⇢): ⇢(X)=G⇢ + p X j=1 mj(Xj)!(3.4) where mjare smooth unknown functions, an intercept coefficient and in this case the link function G⇢= tanh(·)to ensure that ˆ⇢(X)2[1,1]. For the sake of mathematical notational simplicity, only a non-linear effect of the continuous covariates is contemplated in equations (3.3) and (B.3), but they could easily be adapted to incorporate factor effects. For instance, if the first p1covariates define categories, the expression of µr(X)in (3.3) can be replaced by the following semiparametric structure µr(X)=↵r+Pp1 j=1 ↵jrXj+Pp j=p1+1 fjr(Xj). The variances and correlation structures can be similarly treated. 3.2.2 Estimation algorithm This section discusses the procedure for estimating the conditional bivariate region presented in equation (3.2). The methodology is based on the estimate of the covariate effects on the means of the responses via the use of a flexible additive predictor, and then on the variance-covariance matrix using the squared residuals of the conditional means estimates. Finally, the bivariate region ⌧is obtained using a bivariate kernel estimate of the standardized bivariate residuals density. Specifically, given a sample of {(Yi1,Y i2),Xi}n i=1, where Xi=(Xi1,...,X ip), the proposed estimation algorithm is as follows: Step 1: For r=1,2additive predictor is fitted to the original sample {Yir,Xi}n i=1 to obtain the estimates: ˆµr(Xi)=ˆ↵r+ p X j=1 ˆ fjr(Xij)for i=1,...,n (3.5) 34 CHAPTER 3. CONDITIONAL REFERENCE REGION Step 2: For r=1,2the squared residuals of the previous models ˆµr(Xi)are obtained and an additive predictor fitted to the sample {(Yirˆµr(Xi))2,Xi}n i=1 and obtain the estimates: ˆ2 r(Xi)=G ˆ r+ p X j=1 ˆgjr (Xij)!for i=1,...,n (3.6) Step 3: For i=1,...,nthe standardized residuals are computed: ˆri=(Yi1ˆµ1(Xi)) (Yi2ˆµ2(Xi)) ˆ1(Xi)ˆ2(Xi) and the correlation model ⇢(X)using the sample {ˆri,Xi}n i=1 is fitted as: ˆ⇢(Xi)=G⇢ ˆ+ p X j=1 ˆmj(Xij)!(3.7) Step 4: For i=1,...,nthe estimated standardized bivariate residuals are then obtained: ✓ˆ"i1 ˆ"i2◆=ˆ ⌃1/2(X)✓Yi1ˆµ1(Xi) Yi2ˆµ2(Xi)◆ where ˆ ⌃1/2(X)is the inverse of the Cholesky decomposition of ˆ ⌃(X). Using the sample {(ˆ"i1,ˆ"i2)}n i=1 the kernel estimation of the bivariate density ˆ f("1," 2)is obtained. Then, the bivariate region on the residual scale (R⌧) can be obtained as: ˆ R⌧={(u, v)2R2|ˆ f(u, v)ˆ k}(3.8) where ˆ kis the quantile ⌧of ˆ f(ˆ"1,ˆ"2). Step 5: Finally, the conditional bivariate region for each Xivalue is given by: ˆ R⌧(Xi)=✓ˆµ1(Xi) ˆµ2(Xi)◆+ˆ ⌃1/2(Xi)ˆ R⌧ In addition to the R⌧(X)estimate, the proposed algorithm can be used to obtain the univariate conditional reference curves for each response variable, applying the following expression: ˆ Q⌧r(X)=ˆµr(X)+ˆr(X)ˆ"⌧rfor r=1,2(3.9) where ˆ"⌧ris the empirical ⌧-quantile of the univariate errors "1r,...," nr. 3.2. CONDITIONAL BIVARIATE REFERENCE REGION MODEL 35 3.2.3 Bivariate kernel bandwidth estimation The covariate dependent bivariate reference region R⌧(Xi)requires the estimation of a region containing the standardized bivariate residuals with a given probability ⌧derived from the bivariate residuals’ density (see equations (11) and (12)). Given the sample {(ˆ"i1,ˆ"i2)}n i=1, the fdensity estimator at a given point (u, v)is given by: ˆ f((u, v),H)=1 n n X i=1 KH✓uˆ"i1 vˆ"i2◆ where K(·)represents the kernel (a bivariate symmetric probability density function, usually the standard bivariate Gaussian distribution), and Ha diagonal matrix defining the kernel bandwidth, the selection of which is crucial for obtaining good estimate of R⌧and hence the final region R⌧(Xi). −2 0 2 4 6 −2 0 2 4 6 Y1|x=0.5 Y2|x=0.5 Theoretical Estimated −2 0 2 4 6 −2 0 2 4 6 Y1|x=0.5 Y2|x=0.5 Theoretical Estimated Figure 3.1: Example estimating the R⌧(Xi)for ⌧=0.10,0.50 and 0.90, using the plug-in bandwidth estimator (Sheather and Jones,1991) (left) and following the method proposed in equation (3.10) (right). The theoretical bivariate regions were obtained using the parametric density function of the response variable, and 100000 simulated data points. For the selection of the bandwidth H, a plug-in (Sheather and Jones,1991) or cross-validation estimator can be used (Bowman,1984), as in any density estimation problem. However, the optimal bandwidth for density estimation might not be optimal for the coverage properties of the conditional bivariate region ˆ R⌧(Xi) (see Figure 3.1 and Appendix 3.3.3). Thus, a bandwidth needs to be chosen such that it minimizes the difference between the estimated and nominal coverage of ˆ R⌧(Xi). This is achieved here by basing the selection on the expression ˆ H=ˆ h, 36 CHAPTER 3. CONDITIONAL REFERENCE REGION where ˆ his the former bandwidth estimate obtained using the plug-in estimator, and is a parameter modulating the final shape of the estimated region. This parameter is estimated as: ˆ =arg min  n1 n X i=1 I{(Yi1,Y i2)2R(i)(Xi)}!⌧ (3.10) where ⌧is the desired coverage and ˆ R(i) ⌧(Xi)is the estimated bivariate region without the i-th observation. Given the high computational cost of (3.10), a kfold cross-validation scheme could be used instead. 3.3 Simulation Study In this section we evaluate the general performance of the estimation algorithm presented in Section 3.2.2. Furthermore, the estimations of the conditional reference region, see equation (3.2), are evaluated for bivariate Gaussian and non Gaussian data. Finally, we study the effect of the bivariate kernel bandwidth on the estimation accuracy of the conditional reference region. 3.3.1 Estimation algorithm performance The performance of the bivariate location-scale model estimator is evaluated in this section. Specifically, we will assess if the proposed algorithm is able to estimate covariates effects on response means, and on their variance-covariance matrix. Two explanatory covariates X=(X1,X 2)were drawn from a uniform distribution U[0,2]. The bivariate response variable Y=(Y1,Y 2)was generated under the model ✓Y1 Y2◆=✓µ1(X1,X 2) µ2(X1,X 2)◆+✓2 1(X1,X 2)12(X1,X 2) 12(X1,X 2)2 2(X1,X 2)◆1/2✓"1 "2◆ One thousand independent samples {Xi,Yi}n i=1 were generated from this model, under the scenario given by µ1(X1,X 2)=10+5sin(X1),µ 2(X1,X 2)=10+5(X21)2, 1(X1,X 2)=0.15 ·µ1(X1,X 2), 2(X1,X 2)=0.15 ·µ2(X1,X 2) and 12(X1,X 2)=1(X1,X 2)2(X1,X 2)⇢(X1,X 2) with ⇢(X1,X 2)=tanh1.5+sinX2 1+X2 3.3. SIMULATION STUDY 37 Errors "1and "2were drawn from independent distributions N(0,1), and sample sizes fixed to n=500,1000,2000. 0.0 0.5 1.0 1.5 2.0 4 6 8 10 12 14 16 X1 µ1(X1) 0.0 0.5 1.0 1.5 2.0 01234 X1 σ1(X1) 0.0 0.5 1.0 1.5 2.0 10 11 12 13 14 15 16 X2 µ2(X2) 0.0 0.5 1.0 1.5 2.0 1.0 1.5 2.0 2.5 X2 σ2(X2) 0.0 0.5 1.0 1.5 2.0 −1.0 −0.8 −0.6 −0.4 −0.2 0.0 X1 ρ(X1) 0.0 0.5 1.0 1.5 2.0 −1.0 −0.5 0.0 0.5 1.0 X2 ρ(X2) True Estimated mean Simulation 95% CI Figure 3.2: True effects of the model (solid red line) versus the average of the estimated effects (black dashed line), and the corresponding 95% simulation intervals (blue dashed lines) for a sample of n= 1000. 44 CHAPTER 3. CONDITIONAL REFERENCE REGION density estimator, this result is not surprising. This also explains the good coverages obtained for a large kernel bandwidth (see Table 3.3). Hence, if the standardized bivariate residuals are Gaussian-distributed, a parametric expression would be better used in the estimation of R⌧. For non-Gaussian data (scenario 2), our proposed method obtains a smaller estimation error than that returned by the plug-in, and cross-validation estimators. Moreover, the use of a large bandwidth returns a worse estimation error. In addition, theoretical region perimeters are better approximated by our kernel bandwidth estimator. Table 3.3: Coverage probability of bivariate data points for different sample sizes (500,1000 and 2000), covariate X1values, with X2fixed at zero, and bivariate kernel bandwidths estimators. Best-Coverage represents our equation (3.10) proposal, CV is the least-square cross-validation method, and 5Han arbitrary large bandwidth. X1 Bandwidth 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Scenario 1 n = 500 Best-Coverage 93.0 93.9 93.9 94.0 94.0 94.0 93.8 93.7 93.6 93.4 91.7 Plug-in 91.5 92.5 92.6 92.6 92.7 92.6 92.5 92.3 92.3 92.1 90.1 CV 87.2 88.4 88.5 88.5 88.6 88.5 88.4 88.2 88.1 87.9 85.6 5H93.6 94.4 94.4 94.5 94.5 94.5 94.3 94.2 94.2 94.0 92.3 n = 1000 Best-Coverage 93.8 94.3 94.3 94.3 94.3 94.3 94.2 94.1 94.1 94.0 93.2 Plug-in 93.1 93.6 93.6 93.6 93.6 93.6 93.5 93.4 93.4 93.3 92.4 CV 89.8 90.5 90.5 90.5 90.5 90.5 90.4 90.3 90.3 90.2 89.0 5H94.3 94.7 94.7 94.7 94.7 94.7 94.6 94.6 94.5 94.4 93.7 n = 2000 Best-Coverage 94.2 94.5 94.5 94.5 94.5 94.4 94.4 94.4 94.3 94.3 93.8 Plug-in 93.8 94.1 94.1 94.1 94.1 94.1 94.0 94.0 93.9 93.9 93.4 CV 91.6 91.9 91.9 91.9 91.9 91.9 91.8 91.8 91.7 91.7 91.1 5H94.5 94.8 94.8 94.8 94.8 94.8 94.7 94.7 94.7 94.6 94.2 Scenario 2 n = 500 Best-Coverage 92.5 93.6 94.0 94.3 94.3 94.2 94.0 94.2 93.6 93.4 91.7 Plug-in 89.9 91.1 91.8 92.0 92.2 92.2 92.0 92.1 91.5 91.1 88.4 CV 86.1 87.5 88.3 88.6 88.8 89.0 88.7 88.7 88.2 87.6 84.2 5H94.0 94.8 95.2 95.4 95.3 95.2 94.9 95.1 94.5 94.4 93.5 n = 1000 Best Coverage 92.7 93.5 93.9 94.0 94.0 94.0 94.1 94.4 94.0 94.1 93.5 Plug-in 91.2 92.1 92.5 92.7 92.9 92.9 93.0 93.4 92.9 92.9 91.6 CV 87.4 88.5 89.1 89.3 89.5 89.6 89.7 89.9 89.6 89.5 87.6 5H94.2 95.0 95.2 95.2 95.1 95.0 95.0 95.4 95.0 95.1 95.3 n = 2000 Best Coverage 93.1 93.6 93.9 94.1 94.1 94.1 94.2 94.6 94.1 94.1 93.7 Plug-in 92.2 92.7 93.0 93.3 93.3 93.5 93.5 93.9 93.4 93.4 92.8 CV 90.0 90.5 91.0 91.2 91.3 91.5 91.5 91.9 91.4 91.4 90.6 5H94.4 94.9 95.1 95.2 95.1 95.1 95.0 95.4 94.9 95.0 95.1 Figure 3.6 shows the estimated bivariate reference region for every kernel bandwidth. Plug-in, and cross-validation estimators returned greater variabilities in their estimates, while the large bandwidth ignored the true shape of the region. Moreover, the best coverages were obtained using the bandwidth selector described in Section 2.3 (see Table 3.3). Specifically, the plug-in, and crossvalidation methods seems to overfit the training dataset, resulting in a coverage below 95% for the estimated reference region, while the large bandwidth (5H) offer a coverage of over 95%. Hence, when dealing with non-standard responses, the density estimator plays a key role in the performance of the bivariate reference region. The present method shows better performance than the plug-in, 3.3. SIMULATION STUDY 45 and cross-validation estimators, but the bandwidth selection problem requires further work. n = 500 n = 1000 n = 2000 Perimeter RMSE 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 155 205 255 305 355 405 0.0 0.1 0.2 0.3 0.4 0.5 Covariate Value Best−Coverage Plug−in CV 5H Scenario 1: Gaussian data n = 500 n = 1000 n = 2000 Perimeter RMSE 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 240 290 340 390 440 490 540 590 640 690 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Covariate Value Best−Coverage Plug−in CV 5H Scenario 2: non−Gaussian data Figure 3.5: Bivariate reference region performance depending on kernel bandwidth estimator for different X1predictor variable values, with X2fixed at zero, sample sizes (500, 1000 and 2000), for Gaussian and non-Gaussian data. Red line represents the theoretical region’s perimeter. Best-Coverage represents our equation (3.10) proposal, CV is the least-square cross-validation method, and 5Han arbitrary large bandwidth. 46 CHAPTER 3. CONDITIONAL REFERENCE REGION Best Coverage Plug−in Cross−validation 5H Figure 3.6: Estimation of the bivariate reference region for 50 replicates (grey) along with the theoretical region (red), for n=1000, X1=0.5and X2=0, and the bivariate kernel bandwidths used in the estimation of the model’s residual density function. Best Coverage represents our equation (3.10) proposal, and 5H an arbitrary large bandwidth. 3.4 Age-specific bivariate reference region for glycemic tests 3.4.1 Motivating database The A-Estrada Glycation and Inflammation Study (AEGIS) is a cross-sectional, population-based study that was performed in the municipality of A Estrada (Galicia, NW Spain). The study objective was to investigate the association between glycation, inflammation status, lifestyle and common diseases, and to investigate any discordance between glycaemic marker results (Gude et al.,2017). An age-stratified random sample of the population aged 18 years was drawn from Spain’s National Health System Registry. From November 2012 until March 2015, all subjects were successively convened for one day at the A Estrada Primary Care Centre for an evaluation which comprised fasting venous blood sampling, questionnaire interviews, and the description of subjects’ lifestyles. 3.4. AGE-SPECIFIC BIVARIATE REFERENCE REGION FOR GLYCEMIC TESTS47 FPG was determined in subjects’ plasma samples using the glucose oxidaseperoxidase method. HbA1c was determined by high performance liquid chromatography using a Menarini Diagnosticcs HA-8160 analyser; all HbA1c values were converted to DCCT-aligned values (Hoelzel et al.,2004). A total of 1516 subjects (55% female) agreed to participate in the study; their mean age was 52 years, (range 18-91). Among them, 187 (12%) had been diagnosed with diabetes, and among these 66.8% took oral antidiabetics, 3.7% took insulin alone, and 13.3% took insulin and oral drugs. The remaining 16.2% took none of these medications. 4.00 4.25 4.50 4.75 5.00 5.25 5.50 5.75 6.00 6.25 6.50 6.75 7.00 60 70 80 90 100 110 120 130 140 150 160 170 Fasting Plasma Glucose, mg/dL Glycated Haemoglobin, % 4.50 4.75 5.00 5.25 5.50 5.75 6.00 6.25 6.50 70 80 90 100 110 120 130 Fasting Plasma Glucose, mg/dL Glycated Haemoglobin, % 80 years 20 years Figure 3.7: In the left plot an scatter plot and univariate density estimations of the glycaemic test results for subjects not previously diagnosed with diabetes is given, the black lines represent the current diagnostic criteria. In the right plot a comparison of the test results between the younger (20 years) and older patients (80 years) is depicted. Figure 3.7 shows the glycaemic marker concentrations for the subjects along with the current diagnosis cut-off points. As can be seen, according to the criteria in current use, a physician may encounter discordant FPG and HbA1c results, i.e., i) an FPG value outside its reference interval but the HbA1c concentration within the normal range, ii) an HbA1c value outside its reference interval but the FPG level inside the normal range, or iii) both values inside their reference intervals but representing an unlikely combination. To the best of our knowledge, there is no criterion for interpreting such results. Moreover, the use of the same diagnosis cut-off points for the younger and older patients seem unreasonable since the mean values, the variability in the results, and the correlation between the glycaemic markers, increased in older, diabetes-free subjects. 48 CHAPTER 3. CONDITIONAL REFERENCE REGION 3.4.2 Conditional bivariate reference region estimation This section examines the 1329 AEGIS subjects with no diagnosis of diabetes. The age-specific bivariate region containing 95% of these subjects was estimated using the formula below: ✓FPG HbA1c ◆=✓µ1(age) µ2(age)◆+✓2 1(age)12(age) 12(age)2 2(age)◆✓"1 "2◆(3.14) 20 30 40 50 60 70 80 90 5.0 5.2 5.4 5.6 5.8 Age, years E[Glycated Haemoglobin, %] 20 30 40 50 60 70 80 90 0.2 0.4 0.6 Age, years SD[Glycated Haemoglobin, %] 20 30 40 50 60 70 80 90 80 85 90 95 Age, years E[Fasting Plasma Glucose, mg/dL] 20 30 40 50 60 70 80 90 6 8 12 16 Age, years SD[Fasting Plasma Glucose, mg/dL] 20 30 40 50 60 70 80 90 −0.2 0.2 0.6 1.0 Age, years Glycaemic markers correlation Figure 3.8: Estimated effect of age (black) and the corresponding 95% confidence interval (grey) on the mean (E) and standard deviation (SD) of both glycaemic markers and the correlation between them. The effect of age on expectations, the variance, and the correlation between the markers was estimated using polynomial kernel smoothers to account for 3.4. AGE-SPECIFIC BIVARIATE REFERENCE REGION FOR GLYCEMIC TESTS49 possible non-linear trends. As shown in Figure 3.8, the mean HbA1c and FPG levels, and their variability, increase with age, while the strength of their correlation appears not to change. Figure 3.9 shows the bivariate reference region displayed in the standarized residuals scale, after adjusting for age, including approximately 95% of the diseasefree subjects. From a clinical point of view, these subjects would have “normal” values for both glycaemic tests taking into account their age. The other 5% of the participants might be classified in four different groups: (I) first quadrant, individuals with high values for both tests; (II) second quadrant, discordant individuals with high HbA1c concentrations and low/medium FPG; (III) third quadrant, individuals with low values for both tests; and (IV) fourth quadrant, individuals with low/medium HbA1c concentrations and high FPG values. Previous results may have clinical implications, especially for subjects with both markers showing high values, but also for those with discordant results. Subjects returning high values for both (first quadrant) very likely have undiagnosed diabetes. Individuals who fall outside of the reference region in the second quadrant could be labelled as high glycators, i.e., people with normal glucose values but who are at higher risk of cardiovascular disease (McCarter et al.,2004;Cohen,2007) because of their glycation rate. In contrast, individuals in the fourth quadrant, who show high glucose levels but normal glycated haemoglobin levels, could be labelled as low glycators. −2 0 2 4 6 −4−2 0 2 4 6 8 Fasting Plasma Glucose, mg/dL Glycated hemoglobin, % 0.90 0.95 0.975 III III IV Figure 3.9: Bivariate reference region for the glycemic markers (FPG, HbA1c) using standarized bivariate residuals for ⌧=0.90,0.95,0.975. In Figure 3.10 the bivariate reference region is depicted for several ages. These regions shift towards the upper right corner and expand as age increase. This agrees with the non-linear effects of age on the expected means and variability of both markers (Figure 3.8). This may also have clinical implications. For instance, a subject older than 40 years of age with FPG=100mg/dL and HbA1c=6.0% should 50 CHAPTER 3. CONDITIONAL REFERENCE REGION 60 70 80 90 100 110 120 130 4.5 5.0 5.5 6.0 6.5 7.0 Fasting Plasma Glucose, mg/dL Glycated Haemoglobin, % 20 years 30 years 40 years 50 years 60 years 70 years 80 years Figure 3.10: Reference region (⌧=0.95) for age decades. be considered diabetes-free, while a younger subject with the same levels of both markers should be considered to have glycaemic dysregulation. Table 3.4: Percentage of test results contained within the estimated bivariate reference regions for different nominal levels. Cross validation evaluation refers to a leave-one-out cross validation evaluation. Coverage probability is presented for the entire dataset (Global) and for three age groups. Nominal Apparent Cross-validation evaluation Global (18,40] (40,60] (60,91] 54.9 4.9 5.7 5.1 4.0 50 50.0 49.7 50.9 49.6 48.5 90 90.0 89.7 91.6 87.5 90.2 95 95.0 94.1 95.0 93.1 94.3 97 97.1 96.7 97.1 96.2 96.9 98 97.8 97.6 97.6 97.7 97.6 99 99.0 98.2 98.1 98.1 98.5 The performance of the proposed reference region was evaluated in terms of coverage for different ⌧using a leave-one-out cross-validation scheme and three age groups. Table 3.4 shows the estimated region to have a coverage close to the nominal level for every ⌧and age group considered. 3.4.3 Software implementation This methodology was implemented into an R package already available on CRAN (Lado-Baleato et al.,2022). The refreg package contains a set of functions for estimating a conditional reference regions. Its working framework was designed so that people without a strong statistical background can use it. Indeed, only two functions need to be taken into account by the user: 1) the effects 3.4. AGE-SPECIFIC BIVARIATE REFERENCE REGION FOR GLYCEMIC TESTS51 of the covariates on responses need to be estimated using the bivRegr function, a step that requires the user choose which variables may influence the region; 2) bivRegion needs to be applied to a bivRegr object so that the reference region can be estimated. Numerical and graphical summaries of bivRegrand bivRegion-fitted objects for both objects’ classes can be obtained using e.g., summary boot,summary,predict, and plot routines. The bivRegr() function has the following structure: bivRegr(f = formulas, data = data) The fargument contains a list of five R formulae corresponding to the additive predictors for the means, variances and correlation models shown in equation (3.1). Because bivRegr() uses internally the function mgcv::gam(), the user can estimate covariates non-linear effects using the operator s(). For instance the code: mu1 <- y1 ˜ s(x1) mu2 <- y2 ˜ s(x1) var1 <- ˜ x2 var2 <- ˜ x2 rho <- ˜ s(x3) formula = list(mu1,mu2,var1,var2,rho) assumes a smooth effect of x1 on the response means, a parametric effect of x2 on their variances, and a smooth effect of x3 on the response correlation. The bivRegion() function is designed for non-parametrically estimating a bivariate reference region: bivRegion(object, tau = 0.95, bandwidth = "plug-in") The object may be a set of bivariate data points, or a bivRegr object, while tau defines the desired coverage(s) for the reference region, which might be a single value or a vector. Finally, “bandwidth” specifies the kernel bandwidth selection method. The user can chose between the plug-in, cross-validation, or the best coverage method (see equation (3.10)). Appendix Ccontains an R vignette implementing the case study of this chapter and also Chapter 4one. Furthermore, as a proof of concept, we exemplify how this methodology could be extended for a trivariate response variable. 52 CHAPTER 3. CONDITIONAL REFERENCE REGION Chapter 4 Variable selection for prediction regions In this chapter, the introduced bivariate regression model is applied to a real problem in environmental research in order to offer a joint forecasting of two correlated air pollutants. In this context, the conditional reference region is understood as a probabilistic region measuring the joint forecasting uncertainty. Specifically, the future concentrations of two air pollutants, namely SO2and NOx, are simultaneously predicted providing a prediction region which covers a given percentage of the data. The air pollution data comprises several predictor variables corresponding to several time lags. In order to know which covariates are more informative in predicting the course of a pollution episode, a bivariate loss function is introduced for model selection purposes. This loss function was validated with simulated data, showing a satisfactory performance. In practice, it allows to obtain a parsimonious model offering satisfactory results in the joint monitoring of both air pollutants. The remainder of the chapter is organized as follows; Section 4.1 offers a nonexhaustive review about multivariate methods applied in air quality modeling. Then, Section 4.2 revisits our bivariate regression model formulation introducing a loss function for variable selection purposes. Section 4.3 evaluates this loss function performance for model choice. Finally, Section 4.4 develops a model for the joint forecasting of future concentrations of SO2and NOxduring a pollution episode. 4.1 Introduction Air quality is a topic that has been and continues to be object of intensive debate concerning many people, especially those living in industrial areas. One of the main conclusions from this debate is the need to take corrective measures to reduce the emissions of several pollutants into the atmosphere. Predicting 53 60 CHAPTER 4. VARIABLE SELECTION bivariate loss function outperforms test coverage evaluation for every ⌧considered. tau = 0.10 tau = 0.50 tau = 0.95 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 25 50 75 100 a True model selection, (%) Cov. Test Cov. Train Biv. loss Figure 4.4: Percentage of true model selection for coverage and loss function criteria. 4.4 Case study The dataset used in this chapter belongs to a coal-fired power station located in Northern Spain. In Spain, a limit on the mean of 24 consecutive determinations (taken at 5 minutes intervals) of main pollutants is placed around industrial areas. Therefore, there exist a data collection effort from the electro-intensive industrial companies in order to control the concentration of such air pollutants. In our case we have a full record of SO2and NOxconcentrations over a year and several pollution episodes registered. 01:00 06:00 11:00 16:00 21:00 0 100 200 300 400 500 time 5 10 15 20 25 30 SO2 NOx SO2 NOX 0 10 20 30 40 0 100 200 300 400 500 600 NOx SO2 12:00 12:30 13:00 13:30 14:00 14:30 15:00 15:30 16:00 16:35 Figure 4.5: Two different representations of a pollution incident for SO2and Nox. 4.4. CASE STUDY 61 In this section we present the main characteristics of the air pollution data and its restructuring in order to defining predictor variables from time series observations. Let tbe the present time, and SO2tand NOxtthe mean concentrations obtained respectively by the series of bi-hourly SO2and NOxat instant t(5minute temporal instants). Being zthe prediction horizon, the interest is to predict Yt=SO2t+z,NO xt+zusing the vector of predictive covariates Xt=((SO2t,NO xt),(SO2t1,NO xt1),...,(SO2ti,NO xti)) (4.5) providing a prediction region R⌧for these estimations. We are interested in predicting an hour in advance, and therefore, we will consider z=12from now on. The performance of the proposed model will be evaluated in the monitoring of an air pollution incident. The course of this incident is shown in Figure 4.5. In the left plot of this figure, SO2and NOxconcentrations along time are represented with solid and open circles, respectively. Each point in the right plot shows the concentration of both polutants at a specific instant of time. Both plots show an evident correlation between SO2and NOxconcentrations. Figure 4.6 represents the scatter plot of the air pollutants historical matrix. As can be seen the nature of air pollutants time series is characteristic: the values are close to zero most of the time, but in certain points they rise to high levels and then fall back to zero gradually. In order to obtain a reasonably large number of incidents to test our algorithm, avoiding non informative values close to zero, we took as our sample 5000 rows of a historical matrix Mt,{(Xi,Yi)}5000 i=1 . This dataset was constructed removing concentrations below 20µgm3, then the observations time period was divided in 20 strata with a similar number of observations. From the whole number of observations in each stratum, 250 values were selected by random sampling. Finally, we randomly split the historical matrix into a training set MI t=(XI i,YI i) 3000 i=1 and a test set MII t=(XII i,YII i) 2000 i=1 . The prediction region ˆ R⌧(XI i)will be estimated from samples in the first matrix MI t, its performance will be evaluated to the second matrix MII t. 4.4.1 Joint forecasting of two air pollutants In this study, we have considered a maximum of 5 lags, so the predictor covariates are given by X=(SO2t,NOxt)with SO2t=(SO2t,SO 2t1,...,SO 2t4)and NOxt=(NOxt,NO xt1,...,NO xt4). Future joint values of both air pollutants are predicted by the model: ✓SO2t+12 NOxt+12 ◆=✓µ1(SO2t) µ2(NOxt)◆+✓2 1(SO2t)12(SO2t,NOxt) 12(SO2t,NOxt)2 2(NOxt)◆✓"1 "2◆ Table 4.2 shows the coverage and loss function value of the estimated prediction regions evaluated in the test dataset, corresponding to the best combinations 62 CHAPTER 4. VARIABLE SELECTION ●●●●●●●●●●●●● ● ●●●●●●●●●●●●●●●●●●●●●●●●●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ●●●● ● ● ●●● ● ● ● ● ● ●●●●●●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ●●●● ● ● ●● ● ● ● ●●● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●● ● ●● ● ● ● ● ●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ●●● ● ●●●●●●●●●● ● ● ●●●●●●● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ●● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ●●● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ●● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ●● ● ● ● ● ● ● ● ●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ●●●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●●●●●●●●●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●●●●● ●● ●●● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●● ● ● ● ● ● ● ● ● ●●●●●●●●●●● ● ● ● ●● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●● ● ● ● ●●●●●●● ● ● ● ● ● ● ●● ●● ● ●●●●●●●●●●●●●●●●●●●●●●●●●●● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●●● ● ●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●● ● ●● ● ● ● ● ● ●● ● ● ● ●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ●●●● ●● ● ● ● ● ● ● ● ● 0 200 400 600 0 200 400 600 So2(t) So2(t+12) ●●● ●● ●●●●●● ●● ●●●●●●●● ● ●●●●●●●●●●●●●● ● ●●●●●●●●●● ● ●●●●●●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ●● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ●● ●● ● ● ● ●●●● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●●●●●●●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●●● ● ●●●●●●●●●●●● ●● ● ● ● ● ●●●●● ● ● ● ● ● ● ●● ● ●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ●●●●●● ● ● ●● ● ● ●●●●●●●●●●●●●●●●●●●●●●● ● ● ● ●●●●●● ● ●●●●●●●●●●●●●●●●●●●●●●●● ● ●●●●● ● ● ●●● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ●●●● ● ● ●●● ● ● ● ●●●●●●●●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●●●●●● ● ●●●●●● ●●●●●●●● ● ● ●●●●●●●●● ● ● ● ●●●●●●●● ● ●●● ● ● ● ● ●●●●●● ● ● ● ● ● ● ● ● ●● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ●● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●●●●●● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●●●●●●●● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ●● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ●● ● ●● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●●●●● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ●● ●● ●● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ●● ●● ● ●● ● ●● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ●●● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ●● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ●● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●●●● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ●● ●● ● ● ● ●● ● ● ●● ● ● ● ● ● ●● ● ●●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●●●●●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●● ● ●● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ● ●●●●●●●●●●●●●●●●●●●●● ● ● ● ●●●●●●●●●●● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ●●●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ●● ● ●●● ● ●●●●●●●●● ● ●● ●● ● ●● ● ●● ●●●●●●●●●●●●●●● ● ●●●●●●●● ● ●●● ●●●●●●●● ● ● ● ● ● ● ● ● ● ● ●●● ● ●●●●●●●●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●●● ● ● ● ● ● ● ● ● ●● ● ●●●●●●●●●●●●●●●● ● ● ● ●●●●●●●●●●●●●●●●●●●●● ● ● ● ●●●●●●●●●●●●●●●●●●●●●●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●● ● ●●● ● ● ● ●●●● ●● ● ● ●●●●●●●●● ● ● ●●●●●●●●●●● ● ● ●●●●●● ●●●●● ●●●●●●●●●●●●●●●●●●●●●● ● ● ● ● ● ● ●●●●●●●●●●●●●●●●●● ● ● ● ● ● ● ●●●●●●●●●●●●●●●●●●●●● ● ●●●●●●●●●●●●●●●●● ● ●●●●●●●●●●●●●●●●●●●●●● ● ● ●●●●●●●●●●●● ● ● ● ● ●●● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ● ● ● ● ●● ● ●●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ●● ●●● ● ● ●● ● ● ● ● ● ● ● ● ● ●●●●●●● ●● ●●●● ● ● ● ● ● ●● ● ●● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ●●●●●● ● ● ● ● ● ● ● ● 0 10 30 50 0 20 40 60 Nox(t) Nox(t+12) ●●●●●●●●●●●●● ● ●●●●●●●●●●●●●●●●●●●●●●●●●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ●● ●●●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ●●●● ● ● ●●● ● ● ● ● ● ●●●●●●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ●●●● ● ● ●● ● ● ● ●●● ●● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ●● ● ● ● ● ●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ●● ● ●●●●●●●●●● ● ● ●●●●●●● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ●●●● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●● ● ● ●● ● ●● ● ● ● ●● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ●●● ● ● ●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●●● ● ●● ● ● ● ● ●● ● ● ●● ● ●●●● ● ● ● ● ● ●● ●● ● ● ●●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ●● ● ● ● ● ● ●●● ●●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●● ● ●● ●● ● ● ● ●● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ●● ●● ● ●●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ●●●● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●●●●●●●●●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ●● ●●●●● ●● ●●● ● ● ● ● ●●●● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●● ● ● ● ●●● ● ● ●●●●●●●●●●● ● ● ● ●● ● ●● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ●●●●●●● ● ● ● ● ● ● ●● ●● ● ●●●●●●●●●●●●●●●●●●●●●●●●●●● ● ●● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●●● ● ●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●● ● ● ● ●● ●● ● ●● ● ● ●●● ● ● ● ● ● ● ● ● ●●● ● ●● ● ● ●● ● ● ● ●●● ● ●● ● ● ● ●● ●● ● ● ● ●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ●●●● ●● ● ● ● ● ● ● ● ● 0 10 30 50 0 200 400 600 Nox(t) So2(t) ●●● ●● ●●●●●● ●● ●●●●●●●● ● ●●●●●●●●●●●●●● ● ●●●●●●●●●● ● ●●●●●●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ●● ●● ● ● ● ●●●● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●●●●●●●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●●● ● ●●●●●●●●●●●● ●● ● ● ● ● ●●●●● ● ● ● ● ● ● ●● ● ●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●● ● ● ● ● ●● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ●● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ●●●●●● ● ● ●● ● ● ●●●●●●●●●●●●●●●●●●●●●●● ● ● ● ●●●●●● ● ●●●●●●●●●●●●●●●●●●●●●●●● ● ●●●●● ● ● ●●● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ●●●● ● ● ●●● ● ● ● ●●●●●●●●● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●●●●●● ● ●●●●●●●●●●●●●● ● ● ●●●●●●●●● ● ● ● ●●●●●●●● ● ●●● ● ● ● ● ●●●●●● ● ● ● ● ● ● ● ● ●● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ●● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●●●●●●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ●● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●●●●● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ●●● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●●● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ●● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●●●● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●●●●●● ● ● ● ● ●● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●● ● ●● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ● ●●●●●●●●●●●●●●●●●●●●● ● ● ● ●●●●● ●●●●●● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ●●●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ●●● ● ●●●●●●●●●● ●● ●● ● ●● ● ● ●●●●●●●●●●●●●●●● ● ●●●●●●●● ● ●●● ●●●●●●●● ● ● ● ● ● ● ● ● ● ● ●●● ● ●●●●●●●●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ●●●● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ●●● ● ● ● ● ● ● ● ● ●● ● ●●●●●●●●●●●●●●●● ● ● ● ●●●●●●●●●●●●●●●●●●●●● ● ● ● ●●●●●●●●●●●●●●●●●●●●●●● ● ● ●● ● ● ● ● ● ● ● ●●● ●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●● ● ●●● ● ● ● ●●●● ●● ● ● ●●●●●●●●● ● ● ●●●●●●●●●●● ● ● ●●●●●● ●●●●● ●●●●●●●●●●●●●●●●●●●●●● ● ● ● ● ● ● ●●●●●●●●●●●●●●●●●● ● ● ● ● ● ● ●●●●●●●●●●●●●●●●●●●●● ● ●●●●●●●●●●●●●●●●● ● ●●●●●●●●●●●●●●●●●●●●●● ● ● ●●●●●●●●●●●● ● ● ● ● ●●● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ● ● ● ● ●● ● ●●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ●● ●●● ● ● ●● ● ● ● ● ● ● ● ● ● ●●●●●●● ●● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ●●●●●● ● ● ● ● ● ● ● ● 0 200 400 600 0 20 40 60 So2(t) Nox(t+12) Figure 4.6: Graphic representation of some of the data stored in the historical matrix. of the lags in X, and for different values of ⌧. The results obtained show an estimated coverage ˆ⌧close to the true one, ⌧, in all cases. It should be taken into account that the covariate selection may depend on ⌧. That is, the covariates of two models with the same number of terms can change with ⌧. Also note that variable Intappears in all models, while Int1is present in almost all models. As expected, the higher the ⌧, the lower the loss function value of the prediction region. Figure 4.7 shows the value of the loss function for models with a different number of covariates, given the ⌧value, calculated using the test data. As can be appreciated, there is a substantial reduction of the loss value from q=1to q=2, and then the error becomes almost stable. Therefore, in general the most recent measured concentrations have more influence on the prediction than the oldest. Moreover, two or three covariates are enough to predict adequately the concentration values of both pollutants one hour in advance. 4.4. CASE STUDY 63 Table 4.2: Estimated coverage (ˆ⌧%) and L⌧value of the prediction regions obtained with each selected model of size qfor different values of ⌧. The cross Xindicates the covariates included in each model. selected lags ⌧%qˆ⌧% loss tt1t2t3t4 50 3 49.9 8.8 XXX 2 48.6 9.6 XX 4 48.4 9.6 XX X X 4 48.6 9.7 XX X X 75 3 74.8 6.7 XX X 2 76.0 7.1 XX 2 74.0 7.2 XX 5 75.1 7.2 XXXXX 90 4 90.0 3.0 XXXX 3 90.7 3.1 XXX 4 90.3 3.1 XX X X 2 90.0 3.1 XX 95 2 94.9 1.8 XX 2 95.0 1.9 XX 3 94.5 2.0 XX X 2 95.3 2.1 XX ● ● ● ● ●●●● ●● ● ●● ● ● ●●● ● ●● ● ● ● ● ●● ● ● ● ● number of covariates loss value 12345 tau=0.5 ● ● ● ●● ●● ● ●●● ● ● ●● ●●●● ●● ● ●● ● ● ● ●● ● ● number of covariates loss value 12345 tau=0.75 ● ● ●● ● ●● ● ● ●● ●● ●●● ● ●● ●● ●● ●● ● ● ● ● ● ● number of covariates loss value 12345 tau=0.9 ● ● ● ● ●● ●● ●● ●●● ●● ●●●● ● ●● ● ●● ● ● ●●● ● number of covariates loss value 12345 tau=0.95 Figure 4.7: Loss function (L⌧) depending on the number of covariates qfor different values of ⌧. 64 CHAPTER 4. VARIABLE SELECTION Thus, we chosen the following model to monitoring (SO2,NOx) joint values one hour in advance: ✓SO2t+12 NOxt+12 ◆=✓µ1(SO2t,SO 2t2) µ2(NOxt,NO xt2)◆+✓2 1(SO2t,SO 2t2)12(SO2t,NO x2) 12(SO2t,NO x2)2 2(NOxt,NO xt2)◆ Figure 4.8 depicts the joint (SO2,NOx) concentration during one pollution episode along with their bivariate prediction region. Each plot corresponds to different time instants. For each of them, an observation (gray point), its prediction (black point) and the prediction region at 95% are also represented. Note that the shape and the size of these regions change along time. Some of the regions are quite big, which means that there is a significant uncertainty in the estimation of the pollutant concentrations. In fact, the largest prediction regions are located at the top of the episodes, that is, where the curves reach a maximum. This is somewhat expected, it seems the maximum corresponds to a transition between the increase and the decrease of the concentrations, whose estimate is more difficult. The prediction regions in third row are elongated along the X direction, so the estimation of NOxconcentrations has more uncertainty than the corresponding to SO2concentrations. It is also appreciated that most of the times the observed point is inside the region. 4.4. CASE STUDY 65 5 10 20 30 0 100 200 300 400 500 NOx SO2 time: 12:25 ● ● ● ● Observed Predicted 5 10 20 30 0 100 200 300 400 500 NOx SO2 time: 13:15 ● ● 5 10 20 30 0 100 200 300 400 500 NOx SO2 time: 13:40 ● ● 5 10 20 30 0 100 200 300 400 500 NOx SO2 time: 14:30 ● ● 5 10 20 30 0 100 200 300 400 500 NOx SO2 time: 14:55 ● ● 5 10 20 30 0 100 200 300 400 500 NOx SO2 time: 15:20 ● ● 5 10 20 30 0 100 200 300 400 500 NOx SO2 time: 15:45 ● ● 5 10 20 30 0 100 200 300 400 500 NOx SO2 time: 16:35 ● ● 5 10 20 30 0 100 200 300 400 500 NOx SO2 time: 17:25 ● ● Figure 4.8: Example of an SO2-NOxpollution incident (dotted gray line). Measured (gray point) and predicted (black point) concentrations for an individual observation at a specific time. The 95% prediction region for this observation is shown as a black solid line. 66 CHAPTER 4. VARIABLE SELECTION Chapter 5 Testing covariates effects on bivariate reference regions In this chapter, the methodology previously introduced in Chapter 3is extended by proposing a testing procedure for covariates effects on conditional bivariate reference regions. Moreover, the model formulation will be completed by proposing a factor-by-region interaction term. An estimation algorithm based on smoothing splines was used to construct the bivariate reference region for a pediatric anthropometric dataset, and a bootstrap-based hypothesis test proposed for testing age and gender effect on the shape of the reference region. (Height, weight)’s bivariate distribution was shown to depend on the interaction between age and gender. This bivariate growth chart was compared with univariate age-gender body mass index (BMI) percentile curves offering additional insights about children’s body shape. Whereas the well known BMI criterion detects only two atypical situations (i.e., underweight, overweight), the bootstraptested bivariate reference region detected also abnormally large or small body frames for different ages and gender. The layout of this chapter is as follows. Section 5.1 introduces the factor by curve interaction problem for a conditional reference region. Then, Section 5.2 provides a description of the proposed bivariate reference region regression model taking into account a factor-by-curve effect. This section also proposes an estimation algorithm based on penalized splines. Section 5.3 introduces a bootstrap-based procedure for testing covariate interactions in the context of reference regions. In Section 5.4, the adequacy of this bootstrapping method is assessed in a simulation study. Finally, Section 5.5 then sees the model used with real pediatric data, and derives a bivariate interpretation of body height and weight, testing the effect of age and gender on the joint distribution of both measurements. An exhaustive explanation of the model construction details is provided in Appendix B. 67 68 CHAPTER 5. TESTING FACTOR-BY-REGION EFFECT 5.1 Introduction Multivariate reference regions (MVRs) were proposed 50 years ago as a means of interpreting multiple laboratory test results (Winkel and Lyngbye,1972). Despite their theoretical advantages, namely their greater sensitivity and specificity (Boyd,2004), they are rarely used by physicians. Indeed, in practice, when several measurements are available for the same patient, they are commonly transformed into a single entity, or treated as mutually independent variables. These methods ignore the multivariate nature of many patient characteristics, and the information arising from their joint distribution. Moreover, as for percentile curves in the univariate setting, the estimation of covariate effects on the reference region shape can be desirable. For instance, the distribution of continuous test results are reported to be highly influenced by patient age and gender, independent of disease status. The literature contains few proposals for estimating conditional reference regions. A number of statistical proposals and clinical applications based on multivariate Gaussian distribution properties do exist (Dong and Mathew,2015; Mattsson et al.,2008;Selmeryd et al.,2018), but in most clinical settings Gaussiandistributed data are rarely observed. More general methods based on parametric conditional copula models (Klein and Kneib,2016;Marra and Radice,2017) have been used in medicine (Espasand´ ın-Dom´ ınguez et al.,2019;Stander et al.,2019) but selecting a copula regression model requires one choose i) suitable parametric distributions for margins, and ii) a copula that faithfully describes the response association structure. Thus, selecting a model for a real dataset can be time-consuming, and that which best fits the data might not always be entirely reliable. Multivariate quantile regression models have been tested by defining bivariate growth charts or estimating the effects of covariates on the multivariate distribution of anthropometric measurements (Wei,2008;Hallin et al.,2010; Geraci et al.,2020). A method for estimating reference regions based on data depth hyper-rectangular prediction and tolerance regions was recently proposed (Young and Mathew,2020), but while the authors thoroughly studied reference region performance depending on different data depth measurements, the significance of covariate effects on the region shape was not formally tested. In Chapter 3we described a location-scale model for modeling a bivariate continuous response in which the means and variance-covariance matrix were approximated by flexible additive models (Lado-Baleato et al.,2021). Our method provided not only an estimate of the effect of the covariates at both response moments, but also a non-parametric estimate of the region containing 100⌧%of the bivariate data points. The model was then used to forecast the concentrations of two air pollutants (Roca-Pardi˜ nas et al.,2021). In the present chapter, this statistical framework is extended by proposing a bootstrap resampling technique for evaluating the effect of covariates on the shape of the probabilistic region. The secondary aim was to develop and computationally implement an algorithm for 5.2. MATHEMATICAL MODEL 69 estimating the regression model of (Lado-Baleato et al.,2021) in the setting of a more complex formulation involving interaction terms. The proposed statistical test was used with real data from a pediatric anthropometric database. The characterization of body size in children is currently based on age-gender-specific body mass index (BMI =weight, kg/height, m2) cut-offs. In the present chapter we explore the feasibility of a bivariate interpretation of children’s joint (height, weight) values. This strategy requires a bivariate reference region containing most (height, weight) joint values, so that it characterizes which bivariate values are more commonly observed among healthy children. When defining a (height, weight) reference region, the effects of age and gender on the region’s shape must be characterized. The bivariate interpretation proved to be more adequate than BMI cut-offs for describing age and gender effects on child body size. BMI only takes into account whether weight is proportional to height, and only identifies problems of overor underweight. The present age-gender specific reference region, however, offers a total characterization of body size, identifying body frames that are atypically large or small. Identifying such cases in children would be helpful in the early diagnosis of growth hormone disorders (Clayton et al.,2000). Moreover, BMI assumes a specific and constant, association between body height and weight (i.e., body weight is given by the square of body height). This relationship, however, was estimated from the data collected for a homogeneous adult sample, and it might not realistically reflect the condition of children. Indeed, the proposed model describes a reduction in the correlation between both anthropometric measurements. 5.2 Mathematical model Let X=(X1,...,X p)be a vector of pcontinuous covariates and Ya bivariate response of interest Y=(Y1,Y 2). We are interested in obtaining an accurate ⌧th bivariate uncertainty region of Yconditional on X, denoted as R⌧(X), verifying that P(Y2R⌧(X)|X)=⌧, that is R⌧(X)contains the ⌧%of the bivariate response observations. In this context, in Chapter 3we assumed the following locationscale regression model ✓Y1 Y2◆=✓µ1(X) µ2(X)◆+⌃1/2(X)✓"1 "2◆(5.1) where ⌃1/2(X)represents the Cholesky decomposition of the variance-covariance matrix ⌃(X)=Var(Y|X), and the bivariate residuals ("1," 2)are assumed to be independent of the covariates, with zero mean, unit variance, and zero correlation. In practice, the relationship between Yand each of the continuous covariates Xj can vary among the subsets defined by the levels of a categorical factor Fwith M 76 CHAPTER 5. TESTING FACTOR-BY-REGION EFFECT and gender. As can be seen, the bivariate values shift to the upper tail corner as age increase, corresponding to larger bodies. Differences between the bivariate distribution for females and males can easily be seen. 5.5.1 Age-gender effect on the (height, weight) reference region This section reports the main results for the (height, weight) reference region, with a formal test of the effect of age and gender on its shape. Since 80% to 90% of the healthy child population is contained between the current BMI cut-offs, the present reference region R⌧(age, gender)was estimated as ⌧=0.80 and 0.90. Age = 5, years Age = 7, years Age = 9, years Age = 11, years Additive predictor Interaction predictor 90 110 130 150 170 190 90 110 130 150 170 190 90 110 130 150 170 190 90 110 130 150 170 190 0 20 40 60 80 100 120 0 20 40 60 80 100 120 Height, cm Weight, Kg Age = 13, years Age = 15, years Age = 17, years Age = 19, years Additive predictor Interaction predictor 90 110 130 150 170 190 90 110 130 150 170 190 90 110 130 150 170 190 90 110 130 150 170 190 0 20 40 60 80 100 120 0 20 40 60 80 100 120 Height, cm Weight, Kg females males Figure 5.4: Predicted regions for the additive predictor (age + gender), and interaction predictor model (age x gender). 5.5. APPLICATION IN THE CLINICAL SETTING 77 Age and gender entered the reference region model under two predictor specifications. The null model considered an additive effect of age and gender, while an interaction between them was considered in an alternative model. In mathematical terms, the null model is given by equation (2) in which X=(age, gender), while the alternative model is given by equation (3) in which X=age, and F=gender. Figure 5.4 shows the change in the (height, weight) bivariate reference distribution depending on age and gender for the additive and alternative interaction model. Clear differences can be seen between these modelings. The additive effect model (age + gender) shows the female reference region to be closer to the bottom left corner for all ages, while the bivariate distribution of both body measurements in the age x gender interaction model overlaps until puberty. From 13 years onward, the male reference region shifts towards the upper right corner, corresponding to larger body frames. The bootstrapping test shows that the change in the reference region is better described by the interaction term (p value <0.01). 110 120 130 140 150 160 170 180 5 6 7 8 9 10111213141516171819 Age, years Height, cm Height mean 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 5 6 7 8 9 1011 12 1314 15 1617 18 19 Age, years Weight, Kg Weight mean 1 2 3 4 5 6 7 8 9 5 6 7 8 9 10 11 12 13 14 15 16 17 18 Age, years Height, cm Height standard deviation 0 5 10 15 20 25 30 5 6 7 8 9 1011 12 1314 15 1617 18 19 Age, years Weight, Kg Weight standard deviation 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 6 7 8 9 10 11 12 13 14 15 16 17 18 Age, years (Weight, Height) correlation females males (Weight, Height) correlation Figure 5.5: Height and weight mean, standard deviation, and correlation change with age for males and females. The estimated effect is depicted along with 95% point-wise bootstrap confidence interval. Figure 5.5 shows the effect of age x gender on mean height and weight and their standard deviations, as well as on their correlation. Mean height increases linearly with age until puberty in both genders, but from 13 years onward the boys are taller than the girls. Height variability is greater in older subjects, and similar for both genders. Weight increases linearly with age in girls, but in boys a sharp increase is seen after 13 years of age. The standard deviation for weight shows a similar increase with age for both males and females. According to the interaction model, the interdependence between weight and height becomes 78 CHAPTER 5. TESTING FACTOR-BY-REGION EFFECT weaker with age. In 6 year-old boys the anthropometric measurements show a correlation of close to ⇢=0.75, while in 18 year-old young men it is ⇢=0.40. The weakening the correlation between height and weight with age is more pronounced in males than in females. The 95% bootstrap-based point-wise confidence interval confirms that the correlation between weight and height cannot be considered constant in either children or adolescents. 5.5.2 Comparison with BMI interpretations In this section, the conditional reference region results are compared with univariate age-gender BMI percentile curves, estimated using an LMS model Cole and Green (1992) for ⌧=0.05,0.85,and 0.95. These ⌧scorrespond to underweight, overweight, and obesity cut-offs. −2 0 2 4 −20246 Standarized reference region Height, cm Weight, Kg Big Body Obesity Small body Undernutrition 80% reference region 90% reference region III III IV BMI age−dependent cutpoints (boys) Age, years Body Mass Index, Kg/m2 5 6 7 8 9 11 13 15 17 19 5 10 20 30 40 50 Big body Overweight Small body Underweight 5% LMS ref. curve 85% LMS ref. curve 955% LMS ref. curve BMI age−dependent cutpoints (girls) Age, years Body Mass Index, Kg/m2 5 6 7 8 9 11 13 15 17 19 5 10 20 30 40 50 Big body Overweight Small body Underweight 5% LMS ref. curve 85% LMS ref. curve 955% LMS ref. curve Figure 5.6: Standarized 95% bivariate reference region compared with body mass index age dependent reference curves for males, and females. Figure 5.6 shows the difference between the bivariate interpretation and traditional BMI reference curves. Whereas the BMI cut-offs identify two clinical groups (underweight and overweight), the reference region recognizes four. In the standardized scale, the reference region places extreme values in four different quadrants, each one with a different clinical interpretation. The first quadrant identifies those children with extremely high values for both anthropometric measurements (large body frame); the second quadrant identifies children with obesity; the third identifies those with extremely low values for both height and weight (small body frames); and the fourth identifies children suffering from underweight. 5.5. APPLICATION IN THE CLINICAL SETTING 79 Weight, Kg 0 20 50 80 110 140 Males' bivariate growth chart Subject 1 6 years 8 years Subject 3 10 years 12 years Height, m Weight, Kg 90 130 170 0 30 60 90 120 14 years Height, m 90 110 140 170 Subject 7 Subject 5 Subject 6 15 years Height, m 90 110 140 170 16 years Height, m 90 110 140 170 18 years Weight, Kg 0 20 50 80 110 140 Females' bivariate growth chart 6 years Subject 2 8 years 10 years Subject 4 12 years Height, m Weight, Kg 90 130 170 0 30 60 90 120 14 years Height, m 90 110 140 170 15 years Height, m 90 110 140 170 16 years Height, m 90 110 140 170 18 years 5% LMS ref. curve 85% LMS ref. curve 95% LMS ref. curve Figure 5.7: Bivariate reference region (height, weight) real values characterization, compared with the WHO age x gender BMI cutpoints, and LMS 95% age x gender BMI reference curves. 80 CHAPTER 5. TESTING FACTOR-BY-REGION EFFECT The BMI age-gender percentile curves and the proposed model were highly coincident in terms of diagnosing obesity. Most patients in the reference region’s first quadrant are located on the upper part of the overweight or obesity percentile curves. Underweight was not very common according to the proposed model, and no agreement was seen with the BMI conclusions. Most children who, according to the bivariate reference region, had an atypically large or small body frame, were considered as being of normal weight by the BMI reference curves. Figure 5.7 shows the predicted reference region for several ages (each age includes children 6 months younger or older than the stated age in years) and both genders, along with the recorded data. The BMI cut-offs for height and weight values are overlain on the graph. It can be seen how the latter are very close to the reference region contour. The estimated region also takes into account extreme (height, weight) values. Subjects 2, 3, 4 and 7 correspond to boys and girls of different ages with a large body frame for their age; in terms of BMI, however, these children would be regarded as within the normal range. Subject 5 is a similar case, but for a small body size. It can be seen how the reference region properly describes the real shape of the data for each gender and age. Chapter 6 Multivariate reference regions based on conditional transformation models As discussed in previous chapters, our location-scale bivariate regression model impose no distributional restrictions for the response variable. Nevertheless, the model has the underlying assumption of a linear association between the response variables. Chapter 3simulation study showed how this model feature was associated with poorer performance in case of non-elliptical structures of dependence. In this chapter, we introduce a new formulation for estimating conditional reference regions in order to overcome this linear correlation limitation. To do so, we take advantage of Multivariate Conditional Transformation Models (MCTMs Klein et al. (2019)), which allows an elegant formulation for multivariate reference regions. Additionally, MCTMs marginal results are expressed for the first time as percentile curves, instead of using conditional cumulative distributions functions. Both proposals were validated by simulation studies, and applied to a real case in diabetes research. The remainder of this chapter is organized into five additional sections. Section 6.1 gives an overview about transformation models core ideas. Then, Section 6.2 provides an overview of MCTMs structure and inference. Section 6.3 is dedicated to MCTMs percentile curves estimation and the formulation of conditional multivariate reference regions. In Section 6.4, a qualitative analysis of the MCTMs percentile curves, and multivariate conditional regions was carried out. Finally, Section 6.5 presents a detailed analysis for our clinical study including guidelines for model building and results interpretations. 81 82 CHAPTER 6. MCTMS’ REFERENCE REGIONS 6.1 Introduction Recently, Klein et al. (2019) proposed the Multivariate Conditional Transformation Models (MCTMs) - a leap forward in multivariate regression. MCTMs are an attractive statistical tool in multivariate data modeling since they have no parametric restrictions, and are scalable to more than two dimensions. Further, they allow the effects of non-linear covariates on the response marginal moments, marginal and joint quantiles, and dependence structures, to be estimated. These novel features arise from an estimation strategy based on most likely transformation (MLT) statistical inference (Box and Cox,1964). Unlike classical inference procedures, MLT estimation is not based on parametric assumptions about the data, but on their transformation to a known reference distribution. This requires the estimation of a transformation function, hereinafter denoted as h. Once the data are transformed, a correspondence exists between the original data values and the reference distribution, mediated by the transformation function h, and its inverse function h1. The variable moments, density function and cumulative distribution function are obtained by applying h1to the reference distribution. In other words, the problem of estimating the characteristics of the response variables is reduced to the problem of estimating the transformation function h. Hothorn et al. (2014) proposed a conditional transformation function for taking into account covariate effects on a single response variable for continuous, binary or censored data. A maximum likelihood estimator for MLT was also proposed (Hothorn et al.,2018), and extended to MCTMs inference, which provides consistent and asymptotically normal estimators. Therefore, for a multivariate variable Ywith a distribution that depends on a covariate vector X, MCTM estimation is given by the conditional transformation function h(Y|X).Klein et al. (2019) expressed h(Y|X)using Bernstein bases for the multivariate response variables and covariate effects. The present chapter considers the role of MCTMs in the modeling of multivariate biomedical data, focusing on the effects of age and gender on the multivariate distribution of the glycemic markers FPG and HbA1c. Based on the structure of MCTMs, estimates were made for the reference regions covering different percentages of the bivariate distribution of the responses depending on covariate values. The literature contains some recent proposals in this respect (Young and Mathew,2020;Lado-Baleato et al.,2021), but the first of these is valid only for Gaussian-distributed data, and the latter assumes a linear correlation between responses. The conditional MVRs of MCTMs, however, aim to characterize which bivariate response values are more likely to be observed for each covariate value, with the advantage that they deal with non-linear structures of dependence between the responses. The conditional MVR estimator discussed in the present work proved to be reliable with simulated data, returning a smaller estimation error than existing methods involving non-Gaussian heterocedastic response variables. Moreover, these regions provide a way of evaluating the multivariate 6.2. MULTIVARIATE CONDITIONAL TRANSFORMATION MODELS 83 distribution shape for the examined glycemic markers with respect to age and gender. Thus, reference values that take into account the correlation between these markers can be obtained for use by clinicians. Finally, we proposed the use of MCTMs’ conditional percentile curves which shows a smaller estimation error than that associated with classical quantile regression (Koenker and Bassett Jr,1978). Although, the proposed methods were developed for use with the mentioned glycemic markers, they could be used with any disease that requires a diagnosis based on a number of continuous markers, or indeed for examining any problem that involves correlated responses. 6.2 Multivariate conditional transformation models The following describes the main characteristics of the conditional transformation function h(Y|X)for use in multivariate regression. For a multivariate response variable, his expressed as a linear combination of the marginal transformation functions, and a set of parameters measuring the dependence between responses. The estimation of the transformation function requires its expression in terms of basis functions. Thus, the shape of the distribution of the response variables, and the effects of the covariates, are approximated via a known basis function and unknown parameters that need to be estimated from the data. MCTMs estimates of these conditional distributions of the responses, and their inference, rely upon a maximum likelihood estimator and parametric bootstrap resampling. It is beyond the scope of the present work to give a full overview of transformation models, for further details (Hothorn et al.,2014,2018;Klein et al., 2019). 6.2.1 MCTMs structure Multivariate transformation models were developed for J-dimensional, absolutely continuous response vectors Y=(Y1,...,Y J)>2RJ. The most likely transformation estimation procedure then relies on a transformation function h:RJ!RJthat maps the vector Yto a set of Jindependent and identically distributed variables Zj⇠PZ,j =1,...,J where PZis a pre-defined distribution (usually, Zj⇠N(0,1)), i.e. h(Y)=(h1(Y1),...,h J(YJ))>=(Z1,...,Z J)=Z2RJ Multivariate conditional transformation models (MCTMs) are obtained by extending the previous transformation function to include a covariates vector X, resulting in h(Y|X)=(h1(Y1|X),...,h J(YJ|X)).(6.1) 84 CHAPTER 6. MCTMS’ REFERENCE REGIONS A triangular structure of the transformation function his imposed in order to simplify calculations, i.e., it is assumed that the jth component of the transformation function depends only on the first jelements of Y. Therefore, for each j21,...,J, the components of (6.1) can be written as: hj(Y|X)=hj((Y1,...,Y J)|X)=hj((Y1,...,Y j)|X).(6.2) Finally, the transformation functions hjare assumed to be linear combination of marginal transformation functions ˜ hj:RJ!RJ, such that hj(Y|X)= j1 X j=1 j|(X)˜ h|(Y||X)+˜ hj(Yj|X) where j,|controls the linear combination of marginal transformations. Hence, MCTMs are characterized by a set of marginal conditional transformations ˜ hj(Yj|X), j21,...,J, and by a triangular (J⇥J)matrix of transformation coefficients ⇤(X) ⇤(X)=0 B B B B B @ 10 21(X)1 31(X)32(X)1 . . .. . .... J1(X)J2(X)...  J,J1(X)1 1 C C C C C A , Under the standard normal reference distribution PZ= N(0,1), the coefficients of ⇤(X)characterize the dependence structure of the responses and the model specification can be shown to result in a Gaussian copula model with regression effects on the dependence structure and arbitrary marginals. This is so since ˜ Z(X)=( ˜ Z1(X),..., ˜ ZJ(X))>, defined by the random variables ˜ Zj(X)= ˜ hj(Yj|X), follows a multivariate Gaussian distribution ˜ Z(X)⇠N(0J,⌃(X)), with ⌃(X)=⇤(X)1(⇤(X)1)>. Even though the pairwise correlations of the transformed vector ˜ Zare restricted to linear dependencies, the original responses Yare allowed to have nonlinear dependence structure due to the inverse marginal transformation functions Yj=˜ h1 j(˜ Zj). Conditional dependence measures, such as the well known Spearman’s rho (in the following ⇢s), may be derived from the MCTMs’s fit as ⇢sj|(⇤)= 6 ⇡arcsin(R[j, |]/2) (6.3) where R[j, |]is defined as R=S⌃S, with S=diag(1 1,...,1 j), and ⌃1= ⇤1⇤>. 6.2. MULTIVARIATE CONDITIONAL TRANSFORMATION MODELS 85 6.2.2 Transformation function parametrization As explained above, for mathematical convenience h(Y|X)is expressed as a linear combination of marginal transformation functions ˜ hj(Y|X), mediated by parameters j|(X). Both quantities are approximated by basis function expansions. For the marginal (with respect to the response yj) the conditional transformation function may be approximated by ˜ h(Yj|X)=c(Yj,X)#j=a(Yj)>#j,1X>j(6.4) where a(Yj)is a basis representation of the response element Yj,#is a parameter vector that comprise the response’s basis coefficients #j,1and jcontains covariates parametric effects. Note that c(Yj,X)#j, and therefore the marginal transformation function, depends on both the response element Yjand covariates vector X. In equation (6.4)ajrepresents Pjdimensional basis functions aj:R!RPj with basis coefficients #j. For a continuous response ˜ hjshould be smooth in Yj, so any polynomial or spline basis is a suitable choice for a.Hothorn et al. (2018) proposed the use of Bernstein polynomials basis of order M (with M +1 parameters). Under this representation aj(Yj)>results from evaluating densities of beta distributions (characterized by two parameters; mand M), a choice that is computationally convenient because strict monotonicity can be formulated as a set of Mlinear constraints on the parameters #m<# m+1 for all m=0,...,M (McKay Curtis and Ghosh,2011;Hothorn et al.,2018). Apart from the marginal transformation function ˜ hthat allow us to estimate the covariate effects on each response variable conditional density fj(Yj|X)or conditional CDF Fj(Yj|X), MCTMs allows to estimate the covariates effects on the multivariate response dependence structure. This feature is implemented by covariate dependent coefficient of the ⇤matrix as j|(X)=↵j|+X>j|,1|<jJ(6.5) where ↵j|represents the model intercept and X>j|a linear effect of the covariate vector on the j|correlation. Even though, the covariates effects on (6.4) and (6.5) were expressed as a linear shift, non-linear effects may be included using a Bernstein basis representation for a continuous covariate bj(Xj), as a(yj)>#j,1b(X)>j. More complex effects as factor-by-curve may be included as cj=(a> j⌦(1,X>)>)>. 6.2.3 MCTMs inference The conditional transformation function ˜ hj(Yj|X)for each response element Yj, and therefore its conditional density and CDF, is given by the parameters vector 92 CHAPTER 6. MCTMS’ REFERENCE REGIONS improvement in such regions estimation is a consequence of MCTMs structure, which allows it to capture non-linear association structures between response variables. −20246 −6−4−2 0 2 4 6 y1 y2 x=0.10 n=200 −2 0 2 4 6 −6−4−2 0 2 4 6 y1 y2 x=0.50 −2 0 2 4 6 −6−4−2 0 2 4 6 y1 y2 x=0.90 −20246 −6−4−2 0 2 4 6 y1 y2 n=500 −2 0 2 4 6 −6−4−2 0 2 4 6 y1 y2 −2 0 2 4 6 −6−4−2 0 2 4 6 y1 y2 −20246 −6−4−2 0 2 4 6 y1 y2 n=1000 −2 0 2 4 6 −6−4−2 0 2 4 6 y1 y2 −2 0 2 4 6 −6−4−2 0 2 4 6 y1 y2 Theoretical region Estimated regions 1 Figure 6.5: Theoretical regions along with 1000 reference regions estimated using MCTMs, for different sample sizes, ⌧=0.95, and the second simulation scenario. Figure 6.5 shows how the obtained MCTMs reference region faithfully capture the theoretical region’s complex shape. Finally, the MCTMs reference region data coverage was assessed using a test sample (Table 6.1). In general, the coverage was close to that desired for each sample size, covariate value, and simulation scenario. As expected, better data coverage were obtained as the sample size increased. 6.4. SIMULATION STUDY 93 Table 6.1: MCTMs reference region data coverage for different sample sizes, covariate values, and two simulation scenarios. Coverage evaluation was performed in an out-sample design. x Nominal 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Scenario 1 n = 200 10 10.1(1.2) 9.7(1.1) 9.3(1.0) 9.5(1.1) 9.6(1.0) 9.7(1.1) 9.7(1.0) 10.1(1.1) 9.5(1.1) 50 50.2(3.4) 48.6(3.2) 49.4(3.2) 48.9(3.1) 49(3.1) 48.4(3.1) 49.1(3.2) 49.3(3.2) 48.5(3.5) 90 88.2(2.6) 88.4(2.1) 89.1(1.9) 88.8(1.8) 88.9(1.8) 88.9(1.8) 88.3(2.1) 88.5(2.2) 87.1(2.6) 95 93.2(2.1) 93.4(1.6) 94(1.4) 93.7(1.3) 93.8(1.2) 93.9(1.3) 93.5(1.5) 93.6(1.8) 92.1(2.2) 99 96.7(1.7) 97.3(1.2) 97.9(0.8) 98(0.7) 98(0.7) 98.2(0.7) 97.7(0.9) 97.3(1.3) 95.8(2.0) n = 500 10 10.1(0.8) 9.8(0.7) 9.4(0.6) 9.8(0.7) 10(0.6) 10.1(0.7) 10.2(0.7) 10.5(0.8) 9.9(0.8) 50 50.7(2.1) 49(2.0) 50(2.0) 49.7(1.8) 50(1.9) 49.5(1.9) 50.4(2.0) 50.5(2.1) 50(2.3) 90 89.5(1.4) 89.3(1.2) 89.9(1.1) 89.5(1.0) 89.6(1.0) 89.7(1.0) 89.3(1.1) 89.6(1.2) 88.7(1.4) 95 94.6(1.0) 94.2(0.8) 94.7(0.8) 94.4(0.7) 94.5(0.7) 94.6(0.7) 94.4(0.8) 94.7(0.9) 93.8(1.2) 99 98.3(0.8) 98.3(0.5) 98.6(0.4) 98.6(0.3) 98.6(0.4) 98.8(0.3) 98.6(0.4) 98.5(0.5) 97.7(0.9) n = 1000 10 10.1(0.8) 9.8(0.7) 9.4(0.6) 9.8(0.7) 10(0.6) 10.1(0.7) 10.2(0.7) 10.5(0.8) 9.9(0.8) 50 50.7(2.1) 49(2.0) 50(2.0) 49.7(1.8) 50(1.9) 49.5(1.9) 50.4(2.0) 50.5(2.1) 50(2.3) 90 89.5(1.4) 89.3(1.2) 89.9(1.1) 89.5(1.0) 89.6(1.0) 89.7(1.0) 89.3(1.1) 89.6(1.2) 88.7(1.4) 95 94.6(1.0) 94.2(0.8) 94.7(0.8) 94.4(0.7) 94.5(0.7) 94.6(0.7) 94.4(0.8) 94.7(0.9) 93.8(1.2) 99 98.3(0.8) 98.3(0.5) 98.6(0.4) 98.6(0.3) 98.6(0.4) 98.8(0.3) 98.6(0.4) 98.5(0.5) 97.7(0.9) Scenario 2 n = 200 10 9.2(1.7) 9.8(1.6) 10(1.6) 9.9(1.6) 9.8(1.6) 9.7(1.6) 9.6(1.6) 9.5(1.7) 9.3(2.0) 50 45.4(5.8) 47.5(5.1) 48.5(5.0) 48.4(5.0) 48.0(5.0) 47.5(4.9) 47.3(5.1) 46.9(5.4) 45.5(6.3) 90 82.8(5.1) 84.9(4.0) 86(3.6) 85.9(3.6) 85.7(3.7) 85.5(3.6) 85.1(3.8) 84.5(4.2) 82.3(5.4) 95 88.5(4.3) 90.3(3.2) 91.1(2.9) 91.1(2.9) 90.9(2.9) 90.8(2.9) 90.5(3.1) 89.9(3.4) 87.9(4.6) 99 94.6(2.9) 95.7(2.1) 96.2(1.8) 96.2(1.8) 96.1(1.9) 96(1.9) 95.7(2.0) 95.2(2.3) 93.8(3.2) n = 500 10 9.8(1.1) 10.0(1.1) 10.2(1) 10.0(1.0) 10.0(1.0) 9.9(1.0) 9.8(1.1) 9.9(1.2) 9.9(1.3) 50 48.3(3.6) 49.3(3.4) 49.7(3.2) 49.5(3.2) 49.2(3.2) 48.9(3.1) 48.8(3.4) 48.9(3.6) 48.6(4.1) 90 87.2(2.6) 88.1(2.4) 88.6(2.1) 88.3(2.1) 88.2(2.2) 88.1(2.1) 88.1(2.2) 88(2.3) 87.2(2.9) 95 92.5(2.0) 93.2(1.8) 93.6(1.6) 93.4(1.6) 93.3(1.6) 93.3(1.6) 93.2(1.6) 93(1.7) 92.4(2.2) 99 97.4(1.1) 97.8(0.9) 98.0(0.9) 97.9(0.9) 97.8(0.9) 97.8(0.9) 97.7(0.9) 97.5(1) 97.1(1.3) n = 1000 10 10.0(0.8) 10.1(0.8) 10.2(0.7) 10.1(0.8) 10.0(0.9) 10.0(0.8) 9.8(0.8) 10.0(0.9) 10.1(1.0) 50 49.4(2.5) 49.8(2.5) 50.0(2.3) 49.8(2.2) 49.5(2.3) 49.2(2.3) 49.3(2.5) 49.6(2.7) 49.6(2.9) 90 88.6(1.7) 89.1(1.6) 89.4(1.4) 89.1(1.4) 89.0(1.5) 89.0(1.4) 89.0(1.5) 89.1(1.6) 88.8(1.9) 95 93.7(1.3) 94.2(1.2) 94.4(1.0) 94.2(1.0) 94.1(1.1) 94.1(1.0) 94.1(1.1) 94.1(1.1) 93.9(1.4) 99 98.2(0.7) 98.4(0.5) 98.5(0.5) 98.4(0.5) 98.4(0.5) 98.4(0.5) 98.3(0.5) 98.3(0.6) 98.1(0.7) Bernstein basis order and MVR estimation In the MCTMs original work, the authors suggested a Bernstein basis of order 6 to be the best choice for estimating non-linear covariate effects, and for representing the joint distribution of response variables. To test this, the performance of the reference region was compared (in terms of RMSE) for bases of different order. These five models were contemplated: bY1order bY2order bXorder fit 1 6 6 6 fit 2 3 3 6 fit 3 9 9 6 fit 4 6 6 3 fit 5 6 6 9 MCTMs performance was estimated for sample sizes of n=200,500 and 1000 (⌧=0.95,1000 simulation replications). The data was generated for a non- 94 CHAPTER 6. MCTMS’ REFERENCE REGIONS Gaussian bivariate response. Specifically, a logistic and a reverse Gumbel distribution joined by a Gaussian copula. Marginal, and association, parameters were made dependent on a single continuous covariate X2U[1,1] as: 8 > < > : µ1(X) = sin(⇡X)2 1(X)=0.50 + 0.25X µ2(X)=cos(⇡X)2 2(X)=0.50 + 0.25X ✓(x)=tanh(1.5+X) n = 200 n = 500 n = 1000 0.1 0.3 0.5 0.7 0.9 0.1 0.3 0.5 0.7 0.9 0.1 0.3 0.5 0.7 0.9 −2.0 −1.5 −1.0 −0.5 0.0 x log(RMSE) fit 1 fit 2 fit 3 fit 4 fit 5 Figure 6.6: Conditional reference region estimation error for ⌧=0.95, and several basis orders for the response, and covariate effects. As shown in Figure 6.6, the MCTMs’ reference region estimate was similar for each proposed model, except for fit 4 which underestimated the complex shape of the covariate effect. Increasing covariate basis order (fit 5) did not improve the estimate of the region. Fit 3, which used a higher basis order for (Y1,Y 2) than 6 did not improve the estimation error of fit 1. Similarly, a reduced basis order (fit 2) showed similar results to those obtained using 6 order basis. Thus, a basis order of M = 6 for both the response and covariate effect would seem a reasonable compromise between model complexity and flexibility. 6.5 Glycemic markers joint modeling In this section, the joint modeling of the correlated glycemic markers FPG and HbA1c is presented. The distributions and correlation of these two markers may be influenced by age and gender. 6.5.1 Glycemic dataset and model specification MCTMs was used to study the relationship between FPG, and HbA1c, adjusting by age and gender. The data used came from the A-Estrada Glycation and 6.5. GLYCEMIC MARKERS JOINT MODELING 95 Inflammation Study (AEGIS); a population-based study conducted in the AEstrada municipality of northwestern Spain to assess the levels of inflammation and glycation markers in the general adult population, and to examine their association with patient clinical variables (Gude et al.,2017). Women Men 60 80 100 120 140 160 180 60 80 100 120 140 160 180 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 FPG, mg/dL HbA1c, % 20 30 40 50 60 70 80 90 Age Figure 6.7: Glycemic markers’ bivariate distribution change with age for both genders. Figure 6.7 shows the bivariate distribution of both markers for patients not previously diagnosed with diabetes. The joint distribution of both markers seems to be displaced from the bottom left towards the upper right corner as age increases in both genders. Based on a reference distribution Pz=N(0,1), the marginal conditional distributions of the glycemic markers are parameterized as follows: ˆ hj(Yj|Age)=a(Yj)>#j,1Gender ⌦b(Age)>j,j2{FPG,HbA1c} where aand bare Bernstein basis of order 6. The ⇤coefficients, which measures the glycemic markers correlations, were parameterized as: j|(Age, Gender)=Gender ⌦b(Age)j|,|<j2{FPG,HbA1c} Conditional correlation was expressed as Spearman’s rho following equation (6.3). 6.5.2 Marginal results Figure 6.8 depicts conditional marginal cumulative distribution (CDF), and density, functions for both glycemic markers. We may observe a change on the CDF for each marker. Specifically, FPG and HbA1c CDFs presents higher values for 96 CHAPTER 6. MCTMS’ REFERENCE REGIONS Men Women 60 80 100 120 140 160 180 60 80 100 120 140 160 180 0.00 0.25 0.50 0.75 1.00 FPG, mg/dL F(FPG|Age) FPG conditional marginal distribution Men Women 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 0.00 0.25 0.50 0.75 1.00 HbA1c, % F(HbA1c|Age) Age 20 40 60 70 80 90 HbA1c conditional marginal distribution Figure 6.8: Fasting plasma glucose and glycated hemoglobin age-dependent conditional cumulative distribution functions for both genders. FPG's percentiles of Men Age, years FPG, mg/dL 20 30 40 50 60 70 80 90 60 75 90 110 130 150 Men percentiles Women percentiles Glycemic markers values FPG's percentiles of Women Age, years FPG, mg/dL 20 30 40 50 60 70 80 90 60 75 90 110 130 150 HbA1c's percentiles of Men Age, years HbA1c, % 20 30 40 50 60 70 80 90 4 4.5 5 5.5 6 6.5 7 HbA1c's percentiles of Women Age, years HbA1c, % 20 30 40 50 60 70 80 90 4 4.5 5 5.5 6 6.5 7 1 Figure 6.9: Fasting plasma glucose, and glycated hemoglobin age adjusted percentile curves, for men (red) and women (blue). each quantile as age increase, but for FPG in males which show a decrease in subjects older than 60 years. 6.5. GLYCEMIC MARKERS JOINT MODELING 97 Figure 6.9 shows the estimated conditional percentile curves for both genders. FPG increased with age, more so for higher distribution quantiles. Moreover, men between 35 and 70 years old showed higher median blood sugar values and outer quantile than younger subjects. No differences were seen between men and women for HbA1c but were detected for the upper quantiles in middle age (between 45 and 60 years). 6.5.3 Multivariate results Figure 6.10 shows the correlation between FPG and HbA1c to change with age for men and women (Spearman correlation scale). The shaded area represents the 95% pointwise confidence interval for 1000 parametrically drawn bootstrap replicates. According to the confidence interval, no correlation exists between FPG and HbA1c in normoglycemic patients until 20 years of age in women, and 25 in men. Thus, until that age, both markers can be considered independent measures in healthy patients. The correlation then becomes stronger until 40 years of age (and thereafter) when it reaches a value close to ⇢s=0.40 for both genders. FPG − HbA1c) correlation in Men Age, years ρs(FPG −HbA1c) 20 30 40 50 60 70 80 90 −1−0.7 −0.4 −0.1 0.2 0.5 0.8 1 (FPG − HbA1c) correlation in Women Age, years ρs(FPG −HbA1c) 20 30 40 50 60 70 80 90 −1−0.7 −0.4 −0.1 0.2 0.5 0.8 1 Figure 6.10: Fasting plasma glucose, and glycated hemoglobin, correlation change with age, for both genders along with 95% pointwise confidence intervals. Figure 6.11 shows the change in the bivariate reference region with age for men and women. The data shape is faithfully described by MCTMs for every age. Table 6.2 shows the results of a cross validation evaluation for the reference region; data coverage was close to the desired level for several age groups, both for men and women. 98 CHAPTER 6. MCTMS’ REFERENCE REGIONS FPG, mg/dL HbA1c, % 60 80 100 120 140 160 3 3.75 4.75 5.75 6.75 7.75 8.75 Age = 20 years Men FPG, mg/dL HbA1c, % 60 80 100 120 140 160 3 3.75 4.75 5.75 6.75 7.75 8.75 Age = 30 years FPG, mg/dL HbA1c, % 60 80 100 120 140 160 3 3.75 4.75 5.75 6.75 7.75 8.75 Age = 40 years FPG, mg/dL HbA1c, % 60 80 100 120 140 160 3 3.75 4.75 5.75 6.75 7.75 8.75 Women FPG, mg/dL HbA1c, % 60 80 100 120 140 160 3 3.75 4.75 5.75 6.75 7.75 8.75 FPG, mg/dL HbA1c, % 60 80 100 120 140 160 3 3.75 4.75 5.75 6.75 7.75 8.75 FPG, mg/dL HbA1c, % 60 80 100 120 140 160 3 3.75 4.75 5.75 6.75 7.75 8.75 Age = 50 years Men FPG, mg/dL HbA1c, % 60 80 100 120 140 160 3 3.75 4.75 5.75 6.75 7.75 8.75 Age = 60 years FPG, mg/dL HbA1c, % 60 80 100 120 140 160 3 3.75 4.75 5.75 6.75 7.75 8.75 Age = 70 years FPG, mg/dL HbA1c, % 60 80 100 120 140 160 3 3.75 4.75 5.75 6.75 7.75 8.75 Women FPG, mg/dL HbA1c, % 60 80 100 120 140 160 3 3.75 4.75 5.75 6.75 7.75 8.75 FPG, mg/dL HbA1c, % 60 80 100 120 140 160 3 3.75 4.75 5.75 6.75 7.75 8.75 1 Figure 6.11: Fasting plasma glucose, and glycated hemoglobin, reference region change with age for both genders. Table 6.2: Leave-one-out cross validation coverage evaluation. Men Women ⌧20-40 40-60 60+ 20-40 40-60 60+ 0.10 12.6 9.7 10.9 11.0 8.2 9.1 0.50 48.4 53.2 48.8 51.0 42.7 49.5 0.90 88.9 90.3 90.8 93.1 89.3 92.4 0.95 92.6 95.1 93.7 96.5 92.7 93.9 0.99 97.9 98.4 95.9 97.9 97.5 96.9 Chapter 7 Discussion 7.1 Chapter 3: Modelling conditional reference regions This work proposes a new means of estimating reference regions under the assumption of non-parametric conditions, that can be used to help diagnose and treat patients with diseases for which the results of two different tests are considered. The proposed method overcomes the problem of the Gaussian distribution restriction of previously introduced multivariate reference regions (Boyd and Lacher,1982). Moreover, it can estimate non-linear effects of continuous covariates using polynomial kernel smoothers. In simulation studies, it was shown that the procedure for estimating the conditional bivariate reference region was efficient, even for datasets with complex bivariate response distributions. The use of the proposed model revealed that the two biomarkers routinely used in diabetes screening and control (FPG and HbA1c) are better interpreted jointly. Patient age should be taken into account in all interpretations. Disagreements between the measured concentrations of different biomarkers can hinder decision-making when they are so-examined, but they may also provide insight into the future progress of the disease (Rodr´ ıguez-Segade et al.,2011;Nayak et al.,2013;Kim et al.,2018;Nayak et al.,2019). Screening for, and the control of, diabetes mellitus is based largely on the results for two biomarkers, but in other diseases three or more biomarkers may be taken into consideration. For example, thyroid dysfunction is assessed by measuring blood concentrations of thyroid stimulating hormone, tri-iodothyronine (T3) and tetriodothyronine (T4). Currently, the results for each test are compared with their respective univariate reference intervals. Bivariate and trivariate reference regions for a thyroid-healthy control group and for a sample of patients was previously reported (Hoermann et al.,2016), comparing their diagnostic efficiency with the standard assessment method. However, these authors applied the current definition of reference regions that required multivariate Gaussian99 100 CHAPTER 7. DISCUSSION ity, and the variation of thyroid hormone concentrations caused by a number of biological factors(Jonklaas and Razvi,2019) was not taken into account. Other clinical studies have estimated trivariate reference regions based on the Gaussian distribution for cancer (Mattsson et al.,2008) and cardiovascular disease (Selmeryd et al.,2018). Given the need for reference regions for more than two tests, it would be of interest to extend the proposed model beyond the bivariate case. The lack of parametric restrictions in the proposed model, and the possibility of estimating the non-linear effects of the predictor variables would be advantageous in the clinical setting. The proposed model suffers the limitation that the correlation between the response variables is measured using a linear correlation coefficient. This is not optimal for non-Gaussian margins and non-elliptical dependence structures. Although in the simulation studies the proposed model was shown to be quite robust for this type of miss-specification, a slightly increased error in estimates might be expected. More general correlation measures, such as those derived from copula functions (Gijbels et al.,2011;Veraverbeke et al.,2011), might help overcome this problem. 7.2 Chapter 4: Variable selection in conditional prediction regions In this chapter we developed a method to forecast simultaneously concentrations of two pollutants. Given historical concentrations of both pollutants, the method predicts the concentrations in the future, including a bivariate probabilistic region that covers a specific percentage of the concentrations of both pollutants. The size, shape and orientation of these regions provide information regarding the relationships between the concentrations of both pollutants and the uncertainty associated to the estimation change with time. The relationship between the two covariates is taking into account through the variance-covariance matrix. In order to evaluate the performance of the developed model, we define a loss function that takes into account the distance of each observation to the contour of the bivariate uncertainty region, as well as if the observation is inside or outside this region. The method was first validated with simulated data, and the results were very satisfactory in terms of coverage and shape of the uncertainty region. Then, this method was applied to a real case aimed to forecast, one hour in advance, a pollution episode, and the corresponding probabilistic bivariate uncertainty regions, for SO2and NOxconcentrations emitted by a coal-fire power station. The results were also satisfactory, and the analysis of the importance of the predictor variables showed that the best results are reached from data within the three closest data in time. 7.3. CHAPTER 5: TESTING COVARIATES EFFECTS 101 7.3 Chapter 5: Testing covariates effects This work reports a bootstrap-based hypothesis test for examining the effect of covariates on bivariate reference regions. It thereby extends a previously reported statistical framework for estimating such regions. Simulation studies revealed a satisfactory power curve and good approximation of the statistic null hypothesis distribution. Its use with real data showed the bivariate distribution of children’s body height and weight to be different depending on age and gender. According to the present results, body size sexual dimorphism is observed only from puberty. The use of the children’s body size data highlights a clinical setting in which the bivariate regression strategy is preferable to the use of a univariate statistical model. Indeed, several cases were detected in which the bivariate model more reliably identified atypical body sizes. A result with likely clinical implications is that the correlation between body weight and body height becomes weaker with age in children. Since the BMI system was originally defined using the data for a homogeneous adult population, in children the correlation of this index with body composition measurements has remained uncertain (Vanderwall et al.,2017;Mei et al.,2002;L´ opez et al.,2005). The present results agree with those of a recent study (Johnson et al.,2020) in which the authors show a change in the correlation between weight and height when using an age-specific Benn Index (Benn,1971;Garn and Pesick,1982). The proposed testing procedure might be used to help define conditional reference regions in clinical practice. This is also of great importance in biomedical research in which general health controls are based on several continuous and correlated markers. Identifying how clinical variables (e.g., age, gender, body size) affect the multivariate distribution of those markers might offer useful clinical information - information that is currently ignored. For instance, diagnosing hyperor hypothyroidism is based on the blood concentration of three hormones (TSH, T3 and T4). Formally testing whether the multivariate distribution of the hormones differs between genders, or changes with age, would offer a better understanding of thyroid dysfunction (Hoermann et al.,2016). 7.4 Chapter 6: Multivariate reference regions based on conditional transformation models In the final chapter we extend the use of the MCTM regression framework by proposing a method for estimating conditional reference regions, which was shown reliable for describing Gaussian and non-Gaussian bivariate responses. A comparison with existing methods showed a smaller estimation error when re- 108 REFERENCES Klein, N. and Kneib, T. (2016). Simultaneous inference in structured additive conditional copula regression models: a unifying bayesian approach. Statistics and Computing, 26(4):841–860. Klein, N. and Kneib, T. (2019). Directional bivariate quantiles: a robust approach based on the cumulative distribution function. AStA Advances in Statistical Analysis, pages 1–36. Klein, N., Kneib, T., Klasen, S., and Lang, S. (2015). Bayesian structured additive distributional regression for multivariate responses. Journal of the Royal Statistical Society: Series C: Applied Statistics, pages 569–591. Kneib, T. (2013). Beyond mean regression. Statistical Modelling, 13(4):275–303. Koenker, R. (2021). quantreg: Quantile Regression. R package version 5.85. Koenker, R. and Bassett Jr, G. (1978). Regression quantiles. Econometrica., 46(1):33–50. K¨ ohler, M., Schindler, A., and Sperlich, S. (2014). A review and comparison of bandwidth selection methods for kernel regression. International Statistical Review, 82(2):243–274. Kojadinovic, I. and Yan, J. (2010). Modeling multivariate distributions with continuous margins using the copula r package. Journal of Statistical Software, 34(9):1–20. Kong, L. and Mizera, I. (2012). Quantile tomography: using quantiles with multivariate data. Statistica Sinica, 22:1589–1610. Kostova, S. P., Rumchev, K. V., Vlaev, T., and Popova, S. B. (2012). Using copulas to measure association between air pollution and respiratory diseases. In Proceedings of World Academy of Science, Engineering and Technology, number 71 in 1, page 749. World Academy of Science, Engineering and Technology (WASET). Kreuzer, A., Valle, L. D., and Czado, C. (2019). A bayesian non-linear state space copula model to predict air pollution in beijing. arXiv preprint arXiv:1903.08421. Lado-Baleato, O., Roca-Pardinas, J., Cadarso-Suarez, C., and Francisco, G. (2022). refreg: Conditional Multivariate Reference Regions. R package version 0.1.2. Lado-Baleato, ´ O., Roca-Pardi˜ nas, J., Cadarso-Su´ arez, C., and Gude, F. (2021). Modeling conditional reference regions: Application to glycemic markers. Statistics in Medicine, 40(26):5926–5946. Lang, S., Adebayo, S. B., Fahrmeir, L., and Steiner, W. J. (2003). Bayesian geoadditive seemingly unrelated regression. Computational Statistics, 18(2):263–292. REFERENCES 109 Lee, T. C. (2003). Smoothing parameter selection for smoothing splines: a simulation study. Computational statistics & Data analysis, 42(1-2):139–148. Liao, K., Huang, X., Dang, H., Ren, Y., Zuo, S., and Duan, C. (2021). Statistical approaches for forecasting primary air pollutants: A review. Atmosphere, 12(6):686. L´ opez, J. A. F., Remesar, X., and Alemany, M. (2005). Ventajas te´ oricas del ´ ındice de rohrer (p/a3) sobre el ´ ındice de masa corporal (p/a2) para la estimaci´ on de la adiposidad en humanos. Revista Espa˜nola de Obesidad, 3(1):47–55. Marra, G. and Radice, R. (2017). Bivariate copula additive models for location, scale and shape. Computational Statistics & Data Analysis, 112:99–113. Mart´ ınez-Silva, I., Roca-Pardi˜ nas, J., and Ord´ o˜ nez, C. (2016). Forecasting so2 pollution incidents by means of quantile curves based on additive models. Environmetrics, 27(3):147–157. Mattsson, A., Svensson, D., Schuett, B., Osterziel, K. J., and Ranke, M. B. (2008). Multidimensional reference regions for igf-i, igfbp-2 and igfbp-3 concentrations in serum of healthy adults. Growth Hormone & IGF Research, 18(6):506– 516. McCarter, R. J., Hempe, J. M., Gomez, R., and Chalew, S. A. (2004). Biological variation in hba1c predicts risk of retinopathy and nephropathy in type 1 diabetes. Diabetes Care, 27(6):1259–1264. McKay Curtis, S. and Ghosh, S. K. (2011). A variable selection approach to monotonic regression with bernstein polynomials. Journal of Applied Statistics, 38(5):961–976. Mei, Z., Grummer-Strawn, L. M., Pietrobelli, A., Goulding, A., Goran, M. I., and Dietz, W. H. (2002). Validity of body mass index compared with other bodycomposition screening indexes for the assessment of body fatness in children and adolescents. The American journal of clinical nutrition, 75(6):978–985. Microsoft and Weston, S. (2020). foreach: Provides Foreach Looping Construct.R package version 1.5.1. Mohd, Z. I., Roziah, Z., Marzuki, I., Muhd, S. L., et al. (2009). Forecasting and time series analysis of air pollutants in several area of malaysia. American Journal of Environmental Sciences, 5(5):625–632. Nadaraya, E. A. (1964). On estimating regression. Theory of Probability & Its Applications, 9(1):141–142. 110 REFERENCES Nayak, A. U., Nevill, A. M., Bassett, P., and Singh, B. M. (2013). Association of glycation gap with mortality and vascular complications in diabetes. Diabetes Care, 36(10):3247–3253. Nayak, A. U., Singh, B. M., and Dunmore, S. J. (2019). Potential clinical error arising from use of hba1c in diabetes: effects of the glycation gap. Endocrine Reviews, 40(4):988–999. Nelder, J. A. and Wedderburn, R. W. (1972). Generalized linear models. Journal of the Royal Statistical Society: Series A (General), 135(3):370–384. Nelsen, R. B. (2007). An introduction to copulas. Springer Science & Business Media. Nieto, P. G., Lasheras, F. S., Garc´ ıa-Gonzalo, E., and de Cos Juez, F. (2018). Estimation of pm 10 concentration from air quality data in the vicinity of a major steelworks site in the metropolitan area of avil´ es (northern spain) using machine learning techniques. Stochastic Environmental Research and Risk Assessment, 32(11):3287–3298. Nilforooshan, M. A. (2020). mbend: Matrix Bending. R package version 1.3.1. Onis, M. d., Onyango, A. W., Borghi, E., Siyam, A., Nishida, C., and Siekmann, J. (2007). Development of a who growth reference for school-aged children and adolescents. Bulletin of the World health Organization, 85:660–667. Pani, L. N., Korenda, L., Meigs, J. B., Driver, C., Chamany, S., Fox, C. S., Sullivan, L., D’Agostino, R. B., and Nathan, D. M. (2008). Effect of aging on a1c levels in individuals without diabetes: evidence from the framingham offspring study and the national health and nutrition examination survey 2001–2004. Diabetes Care, 31(10):1991–1996. Patton, A. J. (2006). Modelling asymmetric exchange rate dependence. International Economic Review, 47(2):527–556. Pearson, K. and Lee, A. (1900). Mathematical contributions to the theory of evolution. viii. on the inheritance of characters not capable of exact quantitative measurement. part i. introductory. part ii. on the inheritance of coat-colour in horses. part iii. on the inheritance of eye-colour in man. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 195:79–150. Perperoglou, A., Sauerbrei, W., Abrahamowicz, M., and Schmid, M. (2019). A review of spline function procedures in r. BMC medical research methodology, 19(1):1–16. REFERENCES 111 Petersen, J. H. (2009). A non-parametric conditional bivariate reference region with an application to height/weight measurements on normal girls. Biometrical Journal, 51(4):697–709. Ramachandran, A., Riddle, M. C., Kabali, C., Gerstein, H. C., Selmeryd, J., Henriksen, E., Dalen, H., and Hedberg, P. (2012). Relationship between a1c and fasting plasma glucose in dysglycemia or type 2 diabetes: an analysis of baseline data from the origin trial. Diabetes Care, 35(4):749–753. Roca-Pardi˜ nas, J., Gonz´ alez-Manteiga, W., Febrero-Bande, M., Prada-S´ anchez, J., and Cadarso-Su´ arez, C. (2004). Predicting binary time series of so2 using generalized additive models with unknown link function. Environmetrics, 15(7):729–742. Roca-Pardi˜ nas, J., Ord´ o˜ nez, C., and Lado-Baleato, O. (2021). Nonparametric location–scale model for the joint forecasting of so2 and nox pollution episodes. Stochastic Environmental Research and Risk Assessment, 35(2):1–14. Rodr´ ıguez-Segade, S., Rodr´ ıguez, J., Cabezas-Agricola, J. M., Casanueva, F. F., and Camina, F. (2011). Progression of nephropathy in type 2 diabetes: the glycation gap is a significant predictor after adjustment for glycohemoglobin (hb a1c). Clinical Chemistry, 57(2):264–271. Rousseeuw, P. J., Ruts, I., and Tukey, J. W. (1999). The bagplot: a bivariate boxplot. The American Statistician, 53(4):382–387. Ruppert, D., Wand, M. P., and Carroll, R. J. (2003). Semiparametric regression. Number 12. Cambridge university press. Sacks, D. B. (2011). A1c versus glucose testing: a comparison. Diabetes Care, 34(2):518–523. Santos, B. and Kneib, T. (2020). Noncrossing structured additive multiple-output bayesian quantile regression models. Statistics and Computing, 30(4):855–869. Selmeryd, J., Henriksen, E., Dalen, H., and Hedberg, P. (2018). Derivation and evaluation of age-specific multivariate reference regions to aid in identification of abnormal filling patterns: the hunt and vamis studies. JACC: Cardiovascular Imaging, 11(3):400–408. Serfling, R. (2002). Quantile functions for multivariate analysis: approaches and applications. Statistica Neerlandica, 56(2):214–232. Sheather, S. J. and Jones, M. C. (1991). A reliable data-based bandwidth selection method for kernel density estimation. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 53(3):683–690. 112 REFERENCES Siest, G., Henny, J., Gr¨ asbeck, R., Wilding, P., Petitclerc, C., Queralt´ o, J. M., and Petersen, P. H. (2013). The theory of reference values: an unfinished symphony. Clinical Chemistry and Laboratory Medicine, 51(1):47–64. Siman, M. and Bocek, P. (2019). modQR: Multiple-Output Directional Quantile Regression. R package version 0.1.2. Sklar, A. (1973). Random variables, joint distribution functions, and copulas. Kybernetika., 9(6):449–460. Slotnick, H. and Etzell, P. (1990). Multivariate interpretation of laboratory tests used in monitoring patients. Clinical Chemistry, 36(5):748–751. Stadtm¨ uller, U. (1986). Asymptotic properties of nonparametric curve estimates. Periodica Mathematica Hungarica, 17(2):83–108. Stander, J., Dalla Valle, L., Taglioni, C., Liseo, B., Wade, A., and Cortina-Borja, M. (2019). Analysis of paediatric visual acuity using bayesian copula models with sinh-arcsinh marginal densities. Statistics in Medicine, 38(18):3421–3443. Stasinopoulos, D. M. and Rigby, R. A. (2007). Generalized additive models for location scale and shape (gamlss) in r. Journal of Statistical Software, 23(7):1–46. Tang, L. L. and Zhou, X.-H. (2012). A semiparametric separation curve approach for comparing correlated roc data from multiple markers. Journal of Computational and Graphical Statistics, 21(3):662–676. Tenbusch, A. (1997). Nonparametric curve estimation with bernstein estimates. Metrika, 45(1):1–30. Trivedi, P. K. and Zimmer, D. M. (2007). Copula modeling: an introduction for practitioners. Now Publishers Inc. Tsay, R. S. (2013). Multivariate time series analysis: with R and financial applications. John Wiley & Sons. Tukey, J. W. (1975). Mathematics and the picturing of data. In Proceedings of the International Congress of Mathematicians, Vancouver, 1975, volume 2, pages 523–531. Vanderwall, C., Clark, R. R., Eickhoff, J., and Carrel, A. L. (2017). Bmi is a poor predictor of adiposity in young overweight and obese children. BMC pediatrics, 17(1):1–6. Van’t Riet, E., Alssema, M., Rijkelijkhuizen, J. M., Kostense, P. J., Nijpels, G., and Dekker, J. M. (2010). Relationship between a1c and glucose levels in the general dutch population: the new hoorn study. Diabetes Care, 33(1):61–66. REFERENCES 113 Vatter, T. and Chavez-Demoulin, V. (2015). Generalized additive models for conditional dependence structures. Journal of Multivariate Analysis, 141:147–167. Vatter, T. and Nagler, T. (2018). Generalized additive models for pair-copula constructions. Journal of Computational and Graphical Statistics, 27(4):715–727. Veraverbeke, N., Omelka, M., and Gijbels, I. (2011). Estimation of a conditional copula and association measures. Scandinavian Journal of Statistics, 38(4):766– 780. Wand, M. P. and Jones, M. C. (1994). Kernel smoothing. CRC press. Wang, J. and Ghosh, S. K. (2012). Shape restricted nonparametric regression with bernstein polynomials. Computational Statistics & Data Analysis, 56(9):2729– 2741. Watson, G. S. (1964). Smooth regression analysis. Sankhy¯a: The Indian Journal of Statistics, Series A, pages 359–372. Wei, Y. (2008). An approach to multivariate covariate-dependent quantile contours with application to bivariate conditional growth charts. Journal of the American Statistical Association, 103(481):397–409. Winkel, P. and Lyngbye, J. (1972). The normal region—a multivariate problem. Scandinavian journal of clinical and laboratory investigation, 30(3):339–344. Wood, S. (2017). Generalized Additive Models: An Introduction with R. London, UK: Chapman and Hall/CRC, 2 edition. Wood, S. (2021). mgcv: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation. R package version 1.8-38. Wood, S., N., Pya, and Safken, B. (2016). Smoothing parameter and model selection for general smooth models (with discussion). Journal of the American Statistical Association, 111:1548–1575. Wright, E. M. and Royston, P. (1999). Calculating reference intervals for laboratory measurements. Statistical Methods in Medical Research, 8(2):93–112. Yee, T. W. (2015). Vector generalized linear and additive models: with an implementation in R. Springer. Young, D. S. and Mathew, T. (2020). Nonparametric hyperrectangular tolerance and prediction regions for setting multivariate reference regions in laboratory medicine. Statistical Methods in Medical Research, 29(12):3569–3585. Yu, K. and Lu, Z. (2004). Local linear additive quantile regression. Scandinavian Journal of Statistics, 31(3):333–346. 114 REFERENCES Zellner, A. (1962). An efficient method of estimating seemingly unrelated regressions and tests for aggregation bias. Journal of the American statistical Association, 57(298):348–368. Zellner, A. (1963). Estimators for seemingly unrelated regression equations: Some exact finite sample results. Journal of the American Statistical Association, 58(304):977–992. Zhanqiong, H., Sriboonchitta, S., and Jing, D. (2013). Modeling dependence dynamics of air pollution: time series analysis using a copula based garch type model. In Uncertainty analysis in econometrics with applications, pages 215–226. Springer. Appendix A Supplementarial material to Chapter 3 A.1 Flexible additive models estimation In order to obtain the estimated additive models in equations (3.5), (3.6), and (3.7), we have used a backfitting algorithm based on local polynomial kernel smoothers. For mathematical notation simplicity we denote, in this section, Yas our response variable, and X=(X1,...,X p)the p vector of covariates. In this regression framework, we consider the transformed additive model: E[Y|X]=G ↵+ p X j=1 fj(Xj)! where G(·)is a known link function, ↵is a constant and fjunknown functions representing the effects of continuous covariates. Given a sample {(Yi,Xi)}n i=1, this model can be estimated using the following iterative process based on a Newton Raphson procedure which extends the ACE (Alternating Conditional Expectation) algorithm (Hastie and Tibshirani,1990). Initialize: compute the initial estimates, ˆ↵=G1(¯ Y)with ¯ Y=n1Pn i=1 Yi, ˆ f0 1,..., ˆ f0 p=0. Step 1: for i=1,...,n construct the linearised response ˜ Yand the weights Wso that: ˜ Yi=ˆ⌘0 i+YiG(ˆ⌘0 i) G0(ˆ⌘0 i)and Wi=G0(ˆ⌘0 i)2 ˆ2 i where ˆ⌘0 i=ˆ↵+Pp j=1 ˆ f0 j(Xj),G0(⌘)=G ⌘ , and ˆ2 iis an estimation of the variance 2(Yi|G(ˆ⌘0 i)). The estimated ˆ2 ican be obtained fitting an additive model to (Yi G(ˆ⌘0 i))2. 115 116 APPENDIX A. CHAPTER 3 APPENDIX Step 2: fit an additive model to ˜ Yweighted by Wand compute the updates ˆ↵ and ˆ fjfor j=1,...,p. At this step we have used an inner backfitting algorithm based on a local polynomial kernel smoother: Step 2.0: update the constant ˆ↵=( Pn i=1 Wi)1Pn i=1 Wi˜ Y Step 2.1: for j=1,...,pcalculate the partial residuals Rj i=˜ Yiˆ↵ j1 X k=1 ˆ fk(Xik) p X k=j+1 ˆ f0 k(Xik) and for i=1,...,n, compute the polynomial kernel estimator updates: ˆ fj(Xij)= ˆ ⇣Xij,(Xlj,R j l,W i) n l=1 ,h f j⌘ being hf jthe smoothing bandwidth associated with the estimation of fj. Step 2.2: Repeat Step 2.1 replacing ˆ f0 j(Xij)by ˆ fj(Xij)for j=1,...,pand i=1,...,n, until the convergence criterion Pn i=1 ⇣ˆ fj(Xij)ˆ f0 j(Xij)⌘2 Pn i=1 ⇣ˆ f0 j(Xij)⌘2+0.001 "for all j=1,...,p is reached. Step 3: repeat the Steps 1 and 2with ˆ⌘0 ibeing replaced by ˆ⌘i=ˆ↵+Pp j=1 ˆ fj(Xij) for i=1,...,nuntil: |MSE(ˆ⌘,Y )MSE(ˆ⌘0,Y)| MSE(ˆ⌘0,Y)✏ where ✏is a small and the mean squared error MSE(ˆ⌘,Y )is defined as MSE(ˆ⌘,Y )= n1Pn i=1 Wi(YiG(ˆ⌘i))2 The proposed algorithm use two loops: (i) an external loop, for adjusting the transformed response models (Step 1 and 3), (ii) an internal loop which estimates the non-linear effect of continuous covariates using a backfitting method (Step 2). This algorithm was used to estimate the additive models for means (equation (3.5)) with identity link , the additive models for variances (equation (3.6)) with exponential link and the correlation additive model (equation (3.7)) with tanh(·) link. Note that in the first case the algorithm is reduced to the internal loop. Finally, the polynomial kernel smotther used in Step 2 may be replaced by an alternative estimator (e.g. penalized splines) (Wood,2017). A.2. LOCAL POLYNOMIAL KERNEL SMOOTHERS 117 A.2 Local polynomial kernel smoothers Given a sample {(Xi,Y i)}n i=1 with a vector of weights {Wi}n i=1 the local linear kernel smoother at a location x,ˆ (x)= ˆ (x, {(Xi,Y i,W i)}n l=1 ,h)is defined as ˆ (x)= ˆ , where ˆ =(ˆ 0,...,ˆ q)is a vector which minimizes: n X i=1 Wi"Yi q X j=0 j(Xix)j#2 K✓Xix h◆ where K(·)denotes a kernel function (a symmetric density), qthe polynomial degree and h>0is the smoothing parameter, chosen using: CV =1 n n X i=1 Wi⇣Yiˆ (i)(Xi)⌘2 where ˆ (i)(Xi)indicates the fit at Xileaving out the ith data vector A.3 Bootstrap inference for flexible additive predictors In this section we present a bootstrap procedure to obtain punctual confidence intervals, given a specific vector of covariates X0, for the components (mean, deviation and correlation components ) of the model presented in equation (1). The steps for construction of the bootstrap confidence intervals are: Step 1. From the sample data {(Yi1,Y i2),Xi}n i=1 obtain the estimates ˆµr(X0), ˆr(X0)(r=1,2) and ˆ⇢(X0). Step 2. For b=1,...,Bgenerate bootstrap samples {(Y• i1,Y• i2),Xi}n i=1 with ✓Y• i1 Y• i2◆=✓ˆµ1(Xi) ˆµ2(Xi)◆+ˆ ⌃1/2(Xi)✓ˆ"• i1 ˆ"• i2◆ where {(ˆ"• i1,ˆ"• i2)}n i=1 is a sample of size nfrom the residuals {(ˆ"i1,ˆ"i2)}n i=1 with replacement, and compute ˆµ•b r(X0),ˆ•b r(X0)and ˆ⇢•b(X0)as in Step 1. The limits for the 100(1 ↵)% confidence intervals of the true components µr(X0),r(X0)and ⇢(X0)are given repectively by ˆµ↵/2 r(X0),ˆµ1↵/2 r(X0) ˆ↵/2 r(X0),ˆ1↵/2 r(X0) and ˆ⇢↵/2(X0),ˆ⇢1↵/2(X0), where ˆµp r(X0)represents the p-percentil of ˆµ•1 r(X0),...,ˆµ•B r(X0) ˆp r(X0)represents the p-percentil of ˆ•1 r(X0),...,ˆ•B r(X0), and ˆ⇢p(X0)is the ppercentil of ˆ⇢•1(X0),...,ˆ⇢•B(X0).