STANDARDIKÄYRÄN SOVITUS AIKAEROTTEISESSA FLUORESENSSI-IMMUNOMÄÄRITYKSESSÄ JA PASSINGBABLOK -MENETELMÄVERTAILU Valto Kuusisto Pro gradu -tutkielma 28. kesäkuuta 2021 MATEMATIIKAN JA TILASTOTIETEEN LAITOS Turun yliopiston laatujärjestelmän mukaisesti tämän julkaisun alkuperäisyys on tar- kastettu Turnitin OriginalityCheck-järjestelmällä TURUN YLIOPISTO Matematiikan ja tilastotieteen laitos KUUSISTO, VALTO: STANDARDIKÄYRÄN SOVITUS AIKAEROTTEISESSA FLUORESENSSI-IMMUNOMÄÄRITYKSESSÄ JA PASSINGBABLOK -MENETELMÄVERTAILU Pro gradu -tutkielma Tilastotiede 28. kesäkuuta 2021 Aikaerotteisessa fluoresenssi-immunomäärityksessä on usein tavoitteena selvittää tietyn proteiinin pitoisuus testiliuoksesta. Pitoisuuden selvittämiseen käytetään en- nalta määriteltyä standardikäyrää, joka muuntaa mitatun vasteen arvon tietyn pro- teiinin pitoisuudeksi. Standardikäyrän määrittämiseksi käytettävä asetelma pitää suunnitella huolellisesti, mukaan lukien aineiston kerääminen ja parametrien esti- mointiin käytettävät menetelmät. Tässä työssä standardikäyrän määrittämisessä käytetään simuloitua aineistoa, joka vastaa todellisia mittaustuloksia. Ennen mallien sovittamista sekä vasteet että pitoisuudet muunnetaan tasavälisemmiksi logaritmimuunnoksella. Standardikäyrä on usein epälineaarista muotoa ja etenkin sigmoidikäyrän muotoiset epälineaariset mallit sopivat standardikäyräksi. Lisäksi tarkastellaan silotetun splinimallin sovel- tamista standardikäyräksi. Muunnettujen vasteiden variaatio ei ole homogeenista, jolloin mallien parametrien estimointi suoritetaan painotetun pienimmän neliösum- man menetelmällä. Mallien välistä paremmuutta vertaillaan PassingBablok mene- telmävertailulla, jossa jokaista eri menetelmällä sovitettua standardikäyrää verra- taan referenssimenetelmään. Standardikäyrän sovittamisessa lineaarinen malli ei kykene kuvaamaan tarkasti proteiinin pitoisuutta, kun taas epälineaariset mallit pystyvät. Sigmoidikäyrän muo- toisista malleista symmetriset sigmoidimallit eivät myöskään kykene ennustamaan pitoisuuksia tarkasti, vaan epäsymmetriset sigmoidikäyrän muotoiset mallit tai si- lotettu splinimalli ennustavat tarkemmin proteiinin pitoisuutta. PassingBablok -menetelmävertailun sekä keskimääräistä ennustevirhettä tarkasteltaessa sigmoidi- käyrän muotoiset mallit suoriutuivat silotettua splinimallia heikommin. Silotettu splinimalli tuottaa pienemmän keskimääräisen ennustevirheen verrat- tuna sigmoidikäyrän muotoisiin malleihin. Lisäksi otoskoon kasvaessa keskimääräi- sen ennustevirheen suuruus pienenee ja sovitettu malli tulee tarkemmaksi. Passing Bablok -menetelmävertailun perusteella mallien ennustamien pitoisuuksien ja todel- listen pitoisuuksien välinen suhde ei ole lineaarista. Otoskoon ollessa suuri saattaa lineaarisuustarkastelun tulos olla virheellinen sekä hypoteesien testausmenetelmät olla liiankin kireitä. Silotettu splinimalli soveltuu parhaiten standardikäyräksi, mut- ta muita tarkasteluita on syytä tehdä ennen standardikäyrän käyttöönottoa. Asiasanat: Standardikäyrä, silotettu splini, logistinen regressiomalli, painotettu pie- nimmän neliösumman menetelmä, PassingBablok -menetelmävertailu Sisällys 1 Johdanto 1 2 Biologisesta määrityksestä 3 2.1 Aikaerotteinen fluoresenssi-immunomääritys . . . . . . . . . . . . . . 4 2.2 Standardikäyrä . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 Aineiston määrittäminen . . . . . . . . . . . . . . . . . . . . . . . . . 7 3 Lineaarinen malli biologisessa määrityksessä 9 3.1 Parametrien estimointi lineaarisessa mallissa . . . . . . . . . . . . . . 9 3.2 GaussNewton -menetelmä . . . . . . . . . . . . . . . . . . . . . . . . 10 4 Epälineaarinen mallintaminen biologisessa määrityksessä 13 4.1 Sigmoidikäyrämallit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 4.2 Painotettu pienimmän neliösumman menetelmä . . . . . . . . . . . . 17 5 Epäparametrinen mallintaminen biologisessa määrityksessä 19 5.1 Silotettu splini . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 5.2 Ristiinvalidointi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 6 PassingBablok -menetelmävertailu 24 6.1 Parametrien estimointi . . . . . . . . . . . . . . . . . . . . . . . . . . 26 6.2 Lineaarisuustarkastelu . . . . . . . . . . . . . . . . . . . . . . . . . . 28 6.3 Hypoteesintestaus PassingBablok -regressiosuoralle . . . . . . . . . . 29 7 Simulointiesimerkki 31 8 Päätelmät 36 Viitteet 38 1 Johdanto Tässä työssä tarkastellaan ns. biologisessa määrityksessä käytettäviä tilastollisia menetelmiä, erityisesti proteiinien pitoisuuksien arviointiin käytetyn standardikäy- rän sovittamista aikaerotteisessa fluoresenssi-immunomäärityksessä. Aikaerotteisen fluorerenssi-immunomäärityksen avulla tutkija kykenee määrittämään testiliuokses- ta tietyn proteiinin pitoisuuden tai lääkäri pystyy määrityttämään potilaan veri- näytteestä halutun antigeenin pitoisuuden, kuten sydäninfarktin-, raskaushormonin- tai koronaviruksen antigeenin pitoisuuden. Standardikäyrän avulla saadaan tarkas- ti määritettyä pitoisuus, mikäli standardikäyrä on luotu huolellisesti. Epätarkasti määritetty standardikäyrä voi johtaa virheellisesti arvioituihin antigeenin pitoisuu- lukemiin, jolloin tutkija saattaa ajautua vääriin johtopäätöksiin tai lääkäri saattaa tehdä virheellisiä tulkintoja potilaan terveydentilasta. Luvussa 2 esitellään aikaerotteisen fluoresenssi-immunomäärityksen perusidea, jolloin selvitetään kyseiseen määritykseen liittyviä oletuksia. Lisäksi luvussa esitel- lään biologisen määrityksen sellaisia erityispiirteitä, joita tulee ottaa mallintaessa huomioon. Luvussa kerrotaan standardikäyrän tarkoitus sekä selvitetään, millainen on standardikäyrän luomisessa käytetty otosaineisto. Luvussa 3 tarkastellaan biologisen määrityksen mallintamista. Aluksi selitetään yksinkertaisimman lineaarisen mallin sovittaminen aineistoon. Samalla esitellään pienimmän neliösumman menetelmä lineaarisen mallin tapauksessa, kun voidaan olettaa havaintojen olevan normaalisti jakautuneita. Pienimmän neliösumman me- netelmän ratkaisemiseen esitellään myös numeerinen GaussNewton -menetelmä. Neljännessä luvussa lineaarinen malli yleistetään epälineaarisiin parametrisiin malleihin. Epälineaarisista malleista esitellään sigmoidikäyrän muotoiset mallit, neli- ja viisiparametrinen logistinen regressiomalli sekä keskustellaan mallien erityispiir- teistä. Samalla tarkastellaan pienimmän neliösumman menetelmän toimivuutta epä- lineaarisen mallin tapauksessa. Luvussa tutkitaan heterogeenisen aineiston paramet- rien estimointia painotetulla pienimmän neliösumman menetelmällä. Viidennessä luvussa tarkastellaan epäparametrista mallintamista biologisessa määrityksessä. Luvussa tutkitaan silotetun splinimallin sovittamista sekä splinimal- lin solmukohtien asettamista havaintoaineistoon. Silotetun splinimallin sakkopara- metrin sopiva arvo optimoidaan ristiinvalidointimenetelmällä, josta esitellään kaksi tapaa. Kuudennessa luvussa esitellään PassingBablok -regressiomenetelmä sekä sovel- letaan tätä mentelmää biologisen aineiston pitoisuuksien vertailuun. PassingBablok -regressiomenetelmässä estimoidaan parametrit, tarkastellaan menetelmien lineaari- suutta sekä tutkitaan kahden menetelmän samankaltaisuus. Seitsemännessä luvussa tarkastellaan ja sovelletaan luvuissa 26 esiteltyjä mene- telmiä simulointiaineistoon. Simulointiaineiston avulla saadaan myös vertailtua eri mallien soveltuvuutta generoituun aineistoon sekä suoritettua menetelmävertailua PassingBablok -regressiomenetelmällä Työssä on käytetty pääosin D.J. Finneyn lähdeteosta [2], johon myös monet muut tämän työn lähdeteoksista viittaavat. Kappaleen 3 alkupuolella taas on seu- rattu T.V. Mzolon artikkelia [1] ja täydennetty esitystä sopivilta osin H. Huangin [4] ja D.J. Finneyn toisella teoksella [3]. Epälineaaristen menetelmien sovittamisessa on 1 seurattu P. Gotschalk et. al. [5] sekä H. Huangin [4] artikkeleita. Epäparametristen menetelmien yhteydessä on käytetty lähdemateriaaleina T. Hastien ja R. Tibshira- nin teoksia [8] sekä [9], joita on täydennetty enemmän biologiseen määritykseen so- veltavalla lähteellä [10]. PassingBablok -regressiomenetelmän lähteinä on käytetty H. Passingin ja W. Bablokin teosta [11] sekä L. Bilic-Zullen menetelmävertailuar- tikkelia [12]. 2 2 Biologisesta määrityksestä Tieteellistä koetta, jossa pyritään arvioimaan proteiinien biologista aktiivisuutta kutsutaan biologiseksi määritykseksi [1]. Biologisessa määrityksessä ja etenkin ke- mian tuotannossa on usein tavoitteena määrittää tiettyjen proteiinien pitoisuus liu- oksesta. Liuoksena toimii tyypillisesti verinäyte, johon liukenee eri sairauksien tai lääkkeiden vaikutuksesta proteiineja. Näiden proteiinien määrittämiseen on kehitet- ty useita menetelmiä, joista tässä työssä perehdytään aikaerotteisen fluoresenssi- immunomäärityksen (time resolved fluorescence immunoassay, TR-FIA) mallinta- miseen. Liuoksesta selvitetyn pitoisuuden perusteella tutkija kykenee tekemään päätel- miä kemiallisesta reaktiosta tai vastaavasti lääkäri kykenee tekemään arvioita poti- laan terveydentilasta. Väärin määritetty pitoisuus saattaa johtaa tutkijaa harhaan tai vastaavasti lääkäriä ohjaamaan potilaan vääriin jatkotutkimuksiin. Tässä työs- sä tarkastellaan ensin yleisellä tasolla ns. standardikäyrän sovittamista biologisessa määrityksessä. Erityisesti selvitetään standardikäyrän sovittamista aikaerotteiseen fluoresenssi-immunomääritykseen simulointiaineiston avulla. Kuva 1: Aikaerotteisen fluoresenssi-immunomäärityksen reaktiokuvaus. Tilastollisesta näkökulmasta tarkasteltuna biologinen määritys on ennalta suun- niteltu koeasetelma, jossa tutkitaan yhden tai useamman ärsykkeen biologista ak- tiivisuutta näytteessä [3]. Tässä ärsyke on yleinen nimitys lääkkeelle, myrkylle tai näytteelle, jota voidaan lisätä liuokseen tutkijan määrittämillä pitoisuuksilla. Täs- sä työssä liuos on verinäyte ja standardierä on verta vastaava liuos, johon on li- sätty tunnettu määrä antigeenejä. Kun antigeenien pitoisuus tunnetaan stadardie- rässä, saadaan määritettyä ns. standardikäyrä kyseiselle standardiliuoserälle tarkas- ti. Standardikäyrä on siis ärsykkeen vaste kuvattuna pitoisuuden funktiona, jonka avulla saadaan arvioitua myöhemmin saadun ärsykkeen vasteesta pitoisuus. Stan- dardikäyrä luodaan laitteen valmistajalla, joka toimitetaan tutkijalle tai lääkärille laitteen yhteydessä. Tällöin tutkija tai lääkäri ei itse sovita ja luo standardikäyrää, vaan tämän tekee asiaan perehtynyt asiantuntija. 3 2.1 Aikaerotteinen fluoresenssi-immunomääritys Tässä tutkimuksessa tarkasteltu proteiini on jokin antigeeni, jonka pitoisuutta mää- ritetään aikaerotteisella fluoresenssi-immunomäärityksellä. Tähän käytetään testia- lustaan kiinnitettyä vasta-ainetta sekä fluoresoivaa leimattua vasta-ainetta. Aikae- rotteisen fluoresenssi-immunomäärityksen prosessi havainnollistetaan kuvassa 1. Tes- tialustaan lisätään näyte, joka pestään hauduttamisen jälkeen pesuliuoksella. Hau- dutuksen aikana vasta-aine muodostaa näytteessä kiinnostuksen kohteena olevien antigeeni- tai proteiinimolekyylien kanssa sidoksia, johon liittyy vielä leimattu vasta- aine. Nämä vasta-aine-antigeenikompleksit jäävät testialustaan kiinni, kun muu si- toutumaton materiaali pestään pois testialustasta. Testialustaa valaistaan vielä fluo- resoivalla valolla, jolloin valon avulla saadaan määritettyä leimattujen vasta-aine- antigeenisidosten signaalien lukumäärä. Mitä enemmän antigeeniä on näytteessä, sitä enemmän myös sidoksia syntyy ja näin ollen määrityksessä saadaan suurempi signaalivasteen lukumäärä. Tämä oletus on voimassa kohtuullisilla antigeenipitoi- suuksilla, koska liuoksen saturoitumisen jälkeen antigeeni-vasta-ainesidoksia ei enää kykene syntymään, kun kaikki vasta-ainemolekyylien sitoutumispaikat on käytetty. Kuva 2: Aikaerotteisen määrityksen taustaidea. Kun LED-valo on päällä, kas- vaa fluoresenssin signaali. LED-valon sammuttua flueresenssin signaali laskee, kun antigeeni-vasta-ainemolekyylit fluoresoivat valoa. Luotettavan testituloksen saami- seksi viivästetään mittausta LED-valon sammuttamisesta. Näytettä valaistaan LED-valolla, jolloin leimatut vasta-aineet keräävät LED- valon energiaa itseensä ja heijastavat tätä valoenergiaa ulospäin, eli vasta-aineet fluoresoivat valoa. LED-valoa pidetään päällä, kunnes vasta-aineet ovat keränneet tarpeeksi energiaa itseensä, minkä jälkeen valo sammutetaan. Mittausta ei voida kui- tenkaan aloittaa heti valon sammuttamisesta, koska tällöin saataisiin mahdollisesti virheellinen mittaustulos. Valo myös fluoresoi ympärillä olevia materiaaleja, joten 4 mittauksen aloitusta viivästetään hieman ympäristön vaikutuksen tasaantumiseksi. Tällöin fluoresoivien vasta-aineiden lukumäärä saadaan laskettua tarkemmin ilman ympäristön tuomaa mittausvirhettä. Tätä prosessia on havainnollistettu kuvassa 2, joka esittää fluoresenssin määrää ajan funktiona. Mittauslaitteen laskeman signaalivasteen lukumäärän lukemisen jälkeen tutkijan tai lääkärin ongelmaksi tulisi se, että signaalien lukumäärää ei itsessään voida pi- tää pitoisuutena. Tästä syystä antigeenin pitoisuus näytteessä arvioidaan laitteelle syötetyllä standardikäyrällä. Standardikäyrä määritetään aina yhdelle standardiliu- oserälle kerrallaan, koska eri standardiliuoserien välillä saattaa olla eroavaisuuksia ja näin saadaan tarkempi määritys aikaiseksi. Standardikäyrää voidaan kutsua myös kalibrointikäyräksi tai vaihtoehtoisesti voidaan puhua myös standardisuorasta, mikä- li standardikäyrä on lineaarinen. Tältä standardikäyrältä laite kykenee kääntämään signaalivasteesta vastaavan antigeenin pitoisuuden arvion. Standardikäyrän on tar- koitus kuvata mahdollisimman tarkasti vasta-aine-antigeenikompleksien lukumäärää eri antigeenien pitoisuuksilla. Tämän vuoksi epätarkasti määritetty standardikäyrä tuottaa epätarkan pitoisuuden arvon todelliseen pitoisuuteen verrattuna. 2.2 Standardikäyrä Kuva 3: Illustraatio lineaarisen mallin sovittamisesta signaalivasteisiin. Standardikäyrän estimoinnin taustalla on käsite annosvasteen regressiofunktios- ta. Oletetaan, että yi on mitattu realisoitunut signaalivaste satunnaismuuttujasta Yi standardierän S mittausta vastaavalla pitoisuudella xi. Tämä voidaan kuvata regressiofunktiolla F (): yi = F (xi) + "i; "j  N(0; 2): (1) 5 Standardikäyrän määritys käsittää useita erillisiä mittauksia kullakin standardierän S pitoisuudella xi. Samalla pitoisuudella tehtyjä mittauksia kutsutaan toistomit- tauksiksi. Tavallisesti tällaisia toistomittauksia tehdään vähintään kolme [3] jokai- sesta erän S pitoisuudesta xi luotettavamman tuloksen saamiseksi. Yksinkertaisim- paan ja vähiten oletuksia vaativaan standardikäyrän määritykseen mitataan toisto- mittauksia vain kahdelta eri pitoisuudelta x1 ja x2. Tällöin valmistetun tuote-erän S standardikäyrä on lineaarinen suora kahden mittauspisteen välillä, kun taas to- dellisuudessa käyrän muoto saattaisikin olla epälineaarista (kuva 3). Kuvasta 4 sen sijaan voidaan nähdä, ettei tämä kuitenkaan aina kuvaa ilmiön todellista luonnetta, ja regressiosuoran käyttäminen standardikäyränä saattaa johtaa harhaiseen antigee- nin pitoisuuden estimointiin. Kuva 4: Annosvaste -ilmiön kuvaaja.Havaintoaineistoon on sovitettu neliparamet- rinen logistinen regressiofunktio sekä lineaarinen regressiomalli. Neliparametrinen logistinen regressio kykenee selittämään hyvin ilmiön todellisen luonteen. Sen sijaan lineaarinen malli ei kykene selittämään ilmiön S-mallista käyttäytymistä. Lineaari- sella mallilla voitaisiin selittää ilmiön käyttäytymistä välillä (-0.5, 1.0), jossa muutos on hyvinkin lineaarista. Kuvista 5 ja 4 voidaan myös havaita, että kun sekä pitoisuus että signaalivaste muutetaan logaritmiseen asteikkoon, saadaan ilmiö kuvattua selkeämmin verrattuna lineaariseen asteikkoon. Lineaarisessa asteikossa huomataan, että pienillä pitoisuuk- silla olevat arvot näkyvät kuvaajassa yhtenä havaintopisteenä. Tämän lisäksi tasavä- lisellä asteikolla suuremmilla pitoisuuksilla on suuri vipuvoima. Log-log -asteikkoon muuntamalla saadaan havaintopisteiden horisontaaliset välit tasaisemmiksi ja näin ollen myös kolmen pienen pitoisuuden pisteet tulevat näkyviin ja paremmin mallin- nettaviksi. Lisäksi sopivan mallin löytäminen lineaarisen asteikon havainnoille voisi olla ongelmallista tai pahimmassa tapauksessa saatettaisiin sovittaa kokonaan vää- 6 ränlainen mallityyppi aineistoon. (a) Lineaarisilla asteikoilla (b) Log-Log asteikoilla Kuva 5: Kuvassa (a) on kuvattu annosvaste -ilmö lineaarisilla asteikoilla, jolloin kolme alinta pitoisuuspisteen mittausta näkyvät kuvassa yhtenä pisteenä. Muun- tamalla molemmat muuttujat logaritmiseen asteikkoon, saadaan tasavälisempi ja helpommin mallinnettava annosvaste -ilmiö. 2.3 Aineiston määrittäminen Standardikäyrän määrityksessä havaintoaineiston pitoisuuspisteiden paikat tulee määrittää huolellisesti estimointivirheen minimoimiseksi ja tarkan standardikäyrän luomiseksi. Havaintoaineiston pitoisuuspisteiden paikat otetaan tässä työssä annet- tuina. Ennen havaintoaineiston mittaamista selvitetään, millä pitoisuusvälillä stan- dardikäyrän tulee olla määritetty ja tarkka. Tältä pitoisuusväliltä mitataan havain- toaineiston signaalivasteet yi eri standardipitoisuuksilla xi, jotta standardikäyrä saa- daan määritettyä tälle välille sopivaksi. Välin ulkopuolelle jäävällä alueella käyrän arvot joudutaan ekstrapoloimaan. Pitoisuusvälin ulkopuolella pitoisuudet ovat kui- tenkin joko todella matalia tai korkeita, jolloin standardikäyrän estimoiman pitoi- suuden harha ei ole joko tutkimuksen lopputuloksen tai potilaan terveyden tilan kannalta merkittävää. Mittausvälin määrittämisen jälkeen pitoisuusvälille sijoitetaan nc kappaletta pi- toisuuspisteitä, joilta kultakin määritetään k kappaletta toistomittauksia luotetta- vamman mittaustuloksen takaamiseksi. Pisteet on jaoteltu välille niin, että ne ovat tasavälisesti logaritmisella asteikolla, jolloin suurien pitoisuuksien signaalivasteilla ei ole vipuvoimaa pienien pitoisuuksien signaalivasteisiin verrattuna. Vipuvoimalla tarkoitetaan, että suurien pitoisuuksien mahdolliset virheet vaikuttaisivat enemmän standardikäyrän sovitukseen ja estimointitulokseen kuin pienempien pitoisuuksien signaalivasteiden arvot. Otoskooksi saadaan n = nck mittausta ja tällöin saadaan mitattua signaalivasteet yi, i = 1; : : : ; n. Standardikäyrän määrityksessä käytetään d mittauslaitetta mittausepävarmuu- den minimoimiseksi ja suuremman otoskoon tuottamiseksi. Tällöin saadaan otosko- 7 koa kasvatettua d-kertaiseksi yhden laitteen mittauksiin verrattuna. Tällä on myös käytännön ajankäytöllinen etu, koska useammalla laitteella kyetään mittaamaan d mittausta rinnakkain, kun yhdellä laitteella saadaan samassa ajassa vain yksi mit- taus. Lisäksi saadaan useammalta laitteelta mittauksia, jolloin yhden laitteen mah- dollisesti väärin mittaamat arvot saadaan suljettua pois määrityksestä. Tässä työssä valitaan laitteiden lukumääräksi d = 6. Laitteen lukumäärän ollessa d = 6 saadaan siis jokaiselle pitoisuudelle toistomittausten lukumääräksi 6k. Aineiston määrittämi- sen ja signaalivasteiden mittauksien jälkeen sovitetaan standardikäyrä aineistoon. 8 3 Lineaarinen malli biologisessa määrityksessä Varhaisissa eläinkokeissa biologiset vasteet olivat usein binäärisiä, kuten kuolema vs. selviytyminen tai lääke auttoi vs. lääke ei auttanut hoidossa. Näiden binääristen vasteiden mallintamiseen kehitettiin 1900-luvun alkupuolella erilaisia normaalisia sigmoidikäyrän muotoisia logistisia malleja log-log asteikolle (probit-mallit). Täs- sä työssä vasteet ovat kuitenkin kvantitatiivisia eli jatkuvia ja määrällisiä. Niiden mallintamiseen käytettiin 1900-luvun alkupuolella lineaarista regressiomallia. Line- aarinen malli saatiin sovitettua, kun ilmiö linearisoitiin muuntamalla pitoisuudet ja vasteet logaritmisille asteikoille. Oletetaan koeasetelma, jossa tehdään n:n suuruinen otos toisistaan riippumat- tomia mittauksia. Erän mittausvälille on asetettu ennalta määrätty määrä pitoi- suuspisteitä, joista mitataan toistomittauksia. Tällöin saadaan signaalivasteet yi pi- toisuudella x = (c1; : : : ; c1; c2; : : : ; c2; : : : ; cnc), missä kukin pitoisuustaso toistuu 6k kertaa. Tällöin voidaan havaintoaineistoon sovittaa lineaarinen regressiomalli, yi = F (xi; ; ) + "i = + log(xi)+ "i; "i  N(0; 2); (2) jossa on vakiotermi, on kulmakerroin ja "i virhetermi, jonka oletetaan olevan toisista virhetermeistä riippumaton ja normaalisti jakautunut nollaodotusarvolla ja vakiovarianssilla 2. Mallin parametrien estimoimiseksi lineaarinen malli (2) sovitetaan aineistoon pienimmän neliösumman menetelmällä (ordinary least squares). Pienimmän neliö- summan menetelmää voidaan käyttää, koska pienimmän neliösumman menetelmä on ekvivalentti suurimman uskottavuuden menetelmän kanssa, kun vasteiden yi ja- kaumat jokaisella pitoisuuspisteen xi arvolla ovat normaalisia [7]. Pienimmän neliö- summan menetelmässä pyritään löytämään parametrien ^ ja ^ arvot, jotka minimoi- vat neliösummalausekkeen todellisten vasteiden ja mallin F () antamien ennusteiden välillä: nX i=1 (yi F (xi; ; ))2 = nX i=1 (yi y^i)2; i = 1; : : : ; n; (3) jossa n on havaintojen lukumäärä ja F (xi; ; ) = y^i mallin mukainen vasteen odo- tusarvo pisteessä xi. 3.1 Parametrien estimointi lineaarisessa mallissa Lineaarisen mallin tapauksessa estimointi tapahtuu pienimmän neliösuman menetel- mällä logaritmisille signaalivaste-pitoisuuspareille (xi; yi), i = 1; : : : ; n. Tällöin em- piirinen malli saadaan kirjoitettua aiemmin mainittuun muotoon (2). Merkitsemällä F (X;) = (F (x1; ); F (x2; ); : : : F (xn; )) saadaan malli kirjoitettua matriisimuo- dossa: Y = F (X;) + "; (4) jossa  viittaa yleisesti p-ulotteiseen parametrivektoriin  = (1;    ; p). Para- metrivektorin  parametriavaruus  on osajoukko p-ulotteisesta reaaliavaruudesta 9   Rp ja X = 264x1... xn 375 0 ; jossa merkintä x0 viittaa matriisin transponointiin. Lisäksi Y = 264y1... yn 375 ; "1 = 264"1... "n 375 : Merkitään lisäksi F () = F (X;) merkintöjen selventämiseksi. Pienimmän neliösumman menetelmän estimaattoria parametrivektorille  mer- kitään ^, joka on piste parametriavaruudessa. Tällä pisteellä funktion F (^) arvo on lähimpänä havaittua arvoa y. Pienimmän neliösumman estimaattori saadaan joh- dettua minimoimalla residuaalineliösumma (RSS): S() = nX i=1 (yi F (xi;))2;  2   Rp: (5) Oletetaan, että funktio F () on differentioituva parametrin  suhteen. Etsitään pa- rametrien pienimmän neliösumman estimaatit ^ seuraavan yhtälöryhmän ratkaisu- na: @S() @v = 0; v = 1;    ; p: Tästä yhtälöryhmästä saadaan johdettua normaaliyhtälöt nX i=1 @F (xi;) @v (yi F (xi;)) = 0; v = 1;    ; p: (6) Kirjoittamalla Vi;v() = @F (xi;)=@v ja " = y F () saadaan normaaliyhtälöt kirjoitettua matriisimuodossa V()0" = 0; (7) jossa matriisia V kutsutaan vauhtimatriisiksi (velocity matrix). Lineaarisessa regres- siomallissa matriisia V kutsutaan mallimatriisiksi (design matrix), jonka arvot ovat riippumattomia parametrien arvoista . Normaaliyhtälöt kyetään lineaarisen mallin tapauksessa ratkaisemaan, mutta epälineaaristen regressiofunktioiden kohdalla tar- vitaan numeerisia ratkaisumenetelmiä. Yksi näistä on GaussNewton -menetelmä, joka etsii askeltavasti parhaimman estimaatin ^ arvon. 3.2 GaussNewton -menetelmä Normaaliyhtälöitä ratkaistaessa ajaudutaan usein tilanteeseen, jossa ei analyyttis- tä ratkaisua löydetä vaan täytyy turvautua numeerisiin menetelmiin. Näihin tilan- teisiin ajaudutaan, kun regressiofunktio on epälineaarista muotoa, jolloin tavalli- nen pienimmän neliösumman estimointi on harhaista. GaussNewton -menetelmä 10 ratkaisee epälineaariset regressiofunktiot askeltavasti hyödyntämällä Taylorin sarja- kehitelmää. Sarjakehitelmässä tarkastellaan pientä lähialuetta  parametrin  lä- heisyydessä, jolloin voidaan lineaarinen Taylorin sarjakehitelmä voidaan kirjoittaa muotoon F (X;)  F (X;) +rF (X;)0( ); (8) jossa rF (X;)0 =  @F (X;) @1 ;    ; @F (X; ) @p  : Approksimatiivinen RSS saadaan funktion F (xi;) lineaarisesta approksimaatiosta lähialueen  ympäristössä, kun S()  nX i=1 (yi F (xi;))2 pX v=1 Vi;v( )(v v)2: Approksimaatio on kuitenkin pätevä vain lokaalisti valitun lähialueen t sisällä ky- seiselle parametriestimaatille askeleella t. Lähialueen koko riippuu regressiofunktion kaarevuudesta; mitä kaarevampi (lineaarisempi) funktio, sitä pienemmässä (suurem- massa) lähialueessa lineaarinen approksimaatio pätee. Käyttämällä lineaarista Tay- lorin sarjakehitelmää askeltavaan päivittämiseen, saadaan vastaavat normaaliyhtä- löt (6) kuin lineaarisen pienimmän neliösumman tapauksessa t:nnen iteraatiokier- roksen kohdalla: V(t) 0V(t)( t) = V(t)0(y F (t)): (9) GaussNewton -menetelmän askeltava lähestymistapa saadaan yhtälöstä t+1 = t + tt; (10) jossa t on askeleen koko ja t on päivityksen uuden arvon suunta:  = (V(t) 0V(t))1V(t)0(y F (t)): Sopivan askelluksen koon ja suunnan löydyttyä menetelmä siirtyy seuraavaan askel- lukseen. Pienimmän neliösumman estimaattorin ^ asymptoottinen kovarianssi saadaan kirjoitettua Cov(^) = 2(V(^)0V(^))1: (11) Lineaarisessa mallissa vauhtimatriisi vastaa mallimatriisia, V(^) = X, eikä riipu pa- rametrien estimaateista. Lisäksi lineaarisen mallin kovarianssiyhtälö (11) ei tarvitse asymptoottista estimointia, jolloin se on luotettava otoskoosta riippumatta. Lineaarinen malli kykenee monissa tapauksissa selittämään haluttua ilmiötä tar- kasti ja usein logaritmisella muunnoksella saadaan myös epälineaariset tapaukset kuvattua lineaarisella mallilla. Kuitenkin aidosti epälineaarisia tapauksia ilmenee, jolloin lineaarisen mallin antamat estimointitulokset ovat harhaisia. Tällöin tulee pohtia epälineaaristen menetelmien käyttöä ilmiön kuvaamiseen. Usein vasteen ja pitoisuuden välinen suhde ei ole lineaarista vaan vasta logaritmiseen asteikkoon 11 muuntamalla saadaan suhde lineaariseksi ja näin voidaan käyttää pitoisuuden ker- toimen estimointiin samaa mallia kuin kuvattiin yhtälössä (2). Kuten aiemmin jo mainittiinkin luvussa 2.2, vasteen ja pitoisuuden välinen suhde on usein sigmoidi- käyrän muotoinen. Tällöin vasteen riippuvuutta konsentraatiosta ei saada logartimi- sellakaan muunnoksella lineaariseksi, vaan tällöin täytyy käyttää sigmoidifunktiota suhteen kuvaamiseen. Ilmiön lineaarisuuden lisäksi lineaarinen malli vaatii useita oletuksia annosvaste -ilmiöstä, jotta standardierän kertoimen estimointitulos olisi luotettava ja harhaton. Näitä oletuksia ovat ainakin vasteen ja pitoisuuden välinen lineaarisuus sekä resi- duaalien normaalisuus ja homogeenisuus yli pitoisuuksien [4]. Yhdenkin oletuksen rikkoutuessa tulee tehdä korjaavia toimenpiteitä, ja etenkin varianssin homoskedas- tisuusoletus rikkoutuu usein, jolloin tämä tulee huomioida estimoinnissa. Tämän lisäksi myös epälineaarisen mallin tapauksessa estimointitulos voi riippua Gauss Newton -menetelmälle annetuista alkuarvoista. Mikäli menetelmän alkuarvot on ase- tettu kauas todellisista arvoista, voi menetelmä pysähtyä mielestään parhaimpaan tulokseen, joka ei kuitenkaan vastaa todellista parasta arvoa. Virhetermien " ollessa normaalisti jakautuneita ja riippumattomia nollaodo- tusarvolla, pienimmän neliösumman estimaattori on harhaton, normaalisti jakautu- nut ja sillä on pienin varianssi kaikista harhattomista lineaarisista estimaattoreista (Best Linear Unbiased Estimator, BLUE). Biologisessa määrityksessä epälineaaris- ten regressiofunktioiden tapauksessa estimaattorilla ei kuitenkaan ole näitä samo- ja ominaisuuksia, koska epälineaarisen regressiofunktion tapauksessa estimaattori on BLUE vain asymptoottisesti. Kasvattamalla otoskokoa suureksi, tulee Gauss Newton estimaattori aymptoottisesti numeerisesti vakaaksi. Regressiofunktion ol- lessa hyvinkin epälineaarinen, jolloin vauhtimatriisin V() ensimmäisten derivaat- tojen arvot vaihtelevat rajusti, saattaa pienimmän neliösumman estimaattori olla numeerisesti epävakaa. 12 4 Epälineaarinen mallintaminen biologisessa määri- tyksessä Lineaarisen mallin sovituksessa oletuksena on, että vaste kasvaa aina samalla kul- makertoimella pitoisuudesta riippumatta. Biologisessa määrityksessä liuoksen satu- roituminen tai testin rakenne saattavat kuitenkin johtaa siihen, että vasteen arvon kasvu ei ole lineaarista, edes logaritmisella asteikolla. Saturoitumista kuvaa massa- vaikutuksen laki (law of mass action), jonka mukaan kemiallisen reaktion tapahtu- minen (vasta-aineiden ja antigeenin yhdistyminen) on lähes lineaarista ja nopeaa, kunnes reaktio alkaa lähestymään tasapainotilaansa, jolloin reaktio hidastuu. Tä- tä havainnollistetaan myös kuvassa 4, jossa havaitaan vasteen kasvavan kiihtyvästi, kunnes reaktio lähestyy tasapainotilaansa. Lähellä reaktion tasapainotilaa pitoisuu- den kasvattaminen ei enää kasvata vasteen arvoa. Tällaista epälineaarista kasvamis- ta voidaan kuvata erilaisilla sigmoidikäyrän muotoisilla malleilla, joiden yhtälöt ovat yleisesti samaa muotoa kuin yhtälössä (2). Yleisessä tapauksessa yhtälössä (2) F () on jokin monotoninen funktio,  pa- rametrivektori ja "i mallin virhetermit pitoisuudella xi. Virhetermien "i oletetaan olevan toisistaan riippumattomia, normaalisti jakautuneita nollaodotusarvolla ja he- teroskedastisella varianssilla 2i , N(0;  2 i ). Homoskedastisuutta ei voida olettaa, sillä logaritmisen vasteen hajonta pienenee pitoisuuden xi mukana ja homoskedastisuuso- letus saattaisi tuottaa harhaista estimointia. Tähän kuitenkin on löydetty ratkaisu, joka ottaa huomioon heteroskedastisuuden, ja tähän palataan luvussa 4.2. 4.1 Sigmoidikäyrämallit Kaksi perinteistä sigmoidifunktiomallia ovat symmetrinen neliparametrinen logis- tinen regressiomalli sekä laajennettu versio viisiparametrinen logistinen regressio- malli. Viisiparametrinen logistinen regressiomalli ottaa huomioon myös mahdollisen epäsymmetrian vasteen ja pitoisuuden välillä. Viisiparametrinen logistinen funktio on muotoa F (xi; b; c; d; e; f) = c+ d c (1 + exp(b(log(xi) log(e)))f ; (12) jossa xi pitoisuus, e on ns. EC50 (half maximal effective concentration) on se pitoi- suuden xi arvo, joka antaa puolet vasteen ylärajan d ja alarajan c väliltä (d c)=2. Alaraja c on funktion alempi asymptootti, joka vastaa keskimääräisen vasteen suu- ruutta nollapitoisuudella, kun funktio on monotonisesti kasvava. Vastaavasti mon- otonisesti kasvavalle funktiolle yläraja d kuvaa keskimääräisen vasteen ylempää asymptoottia. Lisäksi parametri b kuvaa kulmakertoimen pisteessä e. Parametri f säätelee funktion epäsymmetristä käyttäytymistä parametrin e ympärillä. Funktiot, jotka ovat muotoa 12 ovat aina joko kasvavasti tai laskevasti monotonisia paramet- rien b, c ja d arvoista riippuen. Taulukkoon 1 on kerätty eri parametrien vaihtoehto- jen vaikutukset siihen, onko funktio kasvava vai vähenevä. Epäsymmetriaparametrin ollessa f = 1, viisiparametrinen logistinen regressiomalli supistuu neliparametriseksi logistiseksi malliksi F (xi; b; c; d; e) = c+ d c 1 + exp(b(log(xi) log(e))) ; (13) 13 jossa muut parametrit ovat samat kuin viisiparametrisessa logistisessa regressiomal- lissa. Kuten aiemmin jo mainittiin, on neliparametrinen logistinen regressiomalli symmetrinen pisteen e suhteen, jolloin funktio kasvaa samalla nopeudella kohti ylä- rajaa d kuin se laskee kohti alarajaa c pisteen e suhteen. Mikäli alaraja d = 0, on neliparametrinen logistinen regressiomalli ekvivalentti kolmiparametrisen logisti- sen regressiomallin kanssa, jota ei käsitellä mainintaa enempää tässä työssä. Neli- Kuva 6: Viisi- ja neliparametristen logististen regressiomallien sovitukset simuloi- tuun aineistoon, joka on epäsymmetrinen pisteen (e; a+d 2 ) suhteen. Neliparametrinen logistinen regressiomalli ei kykene kuvaamaan ilmiötä yhtä hyvin kuin epäsymmet- rinen viisiparametrinen logistinen regressiomalli. tai viisiparametrisen regressiomallin sovituksen jälkeen halutaan arvioida havaituis- ta signaalivasteista yi tuntemattomat pitoisuudet xi, mikä onnistuu ratkaisemalla yhtälö (12) pitoisuuden xi suhteen log(xi) = F 1(yi; b; c; d; e; f) = log(e) ln  ( dc log(yi)c) 1=f 1  b ; i = 1; : : : ; n; (14) ja vastaavasti neliparametriselle logistiselle regressiomallille, kun f = 1 log(xi) = F 1(yi; b; c; d; e; f) = log(e) ln  dc log(yi)c 1  b : (15) Neliparametrinen logistinen malli on symmetrinen funktion keskipisteen (e; c+d 2 ) suhteen, jolloin havaitaan taulukon 1 tapauksien 1 ja 4 olevan ekvivalentit keske- nään, kun vain merkitään alarajan c arvo ylärajan d arvoksi ja toisinpäin sekä kun parametri b! b. Tämän välttämiseksi voidaan joko kiinnittää merkinnöissä c > d, 14 Taulukko 1: Parametrien c, b ja d vaikutukset monotonisten neli- ja viisiparametris- ten logististen regressiomallien derivaattoihin, eli kulmakertoimiin. Taulukossa (+) tarkoittaa kasvavaa derivaatta- ja () vähenevää derivaattafunktiota. Tapaus c:n ja d:n järjestys b:n merkki Kulmakerroin 1 c > d b > 0 - 2 c > d b < 0 + 3 c < d b > 0 + 4 c < d b < 0 - jolloin parametrin b merkki määrää kulmakertoimen neliparametriselle logistiselle regressiomallille tai voidaan kiinnittää parametri b > 0, jolloin alarajan c ja ylära- jan d järjestys määrää kulmakertoimen merkin. Kumpikaan näistä merkinnöistä ei vaikuta tai rajoita neliparametrisen logistisen regressiomallin toimintaa. Sen sijaan viisiparametrinen logistinen regressiomalli (f 6= 1) ei ole symmetri- nen pisteen e suhteen, jolloin taulukon 1 jokainen tapaus kuvaa eri tapausta. Ta- paukset 1 ja 4 sopivat kuvaamaan laskevaa annosvastetta, kun taas tapaukset 2 ja 3 sopivat kasvavien annosvastefunktioiden kuvaamiseen. Kuten kuvasta 6 havai- taan, neliparametrinen logistinen regressiomalli ei kykene kuvaamaan epäsymmet- risiä annosvasteilmiöitä yhtä tarkasti kuin viisiparametrinen logistinen regressio- malli. Hyvin pienillä pitoisuuksilla malli näyttäisi yliarvioivan riippuvuutta ja en- nustavan liian suuria vasteita näillä pitoisuuksilla. Samanlainen havainto tehdään myös suurilla pitoisuuksilla, joissa ensin pitoisuudella LOGPitoisuus = 1:0 nelipa- rametrinen logistinen regressiomalli aliarvioi vasteen arvoa, kun taas pitoisuudesta LOGPitoisuus = 2:0 eteenpäin malli yliarvioi todellisen ilmiön vasteen arvoja. 5PL sen sijaan näyttäisi taipuvan hyvin ensin loivemmin kaartuvaan alapäähän ja jyr- kemmin kaartuvaan yläpäähän epäsymmetrisen funktiotyypin vuoksi. Kuvassa 7 esitetään parametrien b f arvojen vaihtelun vaikutukset viisipara- metriseen logistiseen regressiomalliin, kun mallin muut parametrit ovat kiinnitet- tyinä. Siniset käyrät ovat keskenään ekvivalentit eri kuvaajien välillä. Kuvassa 7 (a) kasvatetaan ylärajan d arvoa, vastaavasti kuvassa 7 (b) kasvatetaan alarajan c arvoa. Kuvassa 7 (c) kasvatetaan kulmakerroinparametrin b arvoa, jolloin käy- rän kasvu hidastuu. Kuvassa 7 (d) kasvatetaan keskipisteparametrin e arvoa, jolloin käyrän sijainti vain muuttuu x-akselilla mitattuna. Viimeisessä kuvassa 7 (e) kasva- tetaan taas epäsymmetrisyysparametrin f arvoa, jolloin havaitaan, että malli tulee symmetrisemmäksi, kun f ! 1. Myös neli- ja viisiparametristen logististen regressiomallien parametrit saadaan estimoitua käyttämällä pienimmän neliösumman menetelmää. Yksinkertaistettuna parametrien estimointi on residuaalineliösumman minimointia, kuten yhtälössä (3). Tässä tapauksessa mallit eivät ole lineaarisia, jolloin tavanomainen pienimmän ne- liösumman estimointi saattaa olla harhaista. Tästä syystä pienimmän neliösumman estimointi täytyy suorittaa numeerisia menetelmiä käyttäen, kuten luvun 3.2 Gauss Newton -menetelmää. Tätä menetelmää täytyy kuitenkin laajentaa yleisemmäksi, jotta se toimisi myös epälineaaristen mallien tapauksessa. 15 (a) Parametrin d vaikutus (b) Parametrin c vaikutus (c) Parametrin b vaikutus (d) Parametrin e vaikutus (e) Parametrin f vaikutus Kuva 7: (a) - (e) parametrien b f eri arvojen vaikutus viisiparametriseen logis- tiseen regressiomalliin. Siniset käyrät ovat keskenään yhtäpitävät eli ne on tehty samoilla parametrien arvoilla, jolloin kuvat ovat keskenään vertailukelpoisia. Näissä kuvaajissa c > d ja b < 0, kuten kuvasta (c) havaitaan. 16 4.2 Painotettu pienimmän neliösumman menetelmä Pienimmän neliösumman menetelmä epälineaarisissa regressiomenetelmissä on har- haton, jos havaintojen variaatio on homogeenistä eri pitoisuuksilla. Näin ei kui- tenkaan usein ole, vaan logaritmisen pitoisuuden kasvaessa logaritmisten havainto- jen variaatio pienenee. Tällöin epälineaarisen mallin pienimmän neliösumman esti- moinnissa tulee huomioida myös heteroskedastisuus. Tämä huomioidaan painotetul- la pienimmän neliösumman estimoinnilla (weighted least squares), jota laskennan helppoudesta johtuen käytetään suurimman uskottavuuden menetelmän sijasta [6]. Tilastollisen regressioteorian mukaan [7] painotettu pienimmän neliösumman me- netelmä tuottaa suurimman uskottavuuden menetelmän kanssa täsmälleen samat, harhattomat parametrien estimaattorit, kun painotus on valittu oikein. Painotettu pienimmän neliösumman estimointi olettaa jokaisen pitoisuuspisteen xi havaintojen jakauman olevan suunnilleen normaalista [5]. Painotetussa pienim- män neliösumman menetelmässä tavoitteena on löytää parametrit, jotka minimoi- vat painotetun neliösumman havaittujen signaalivasteiden yi ja sovitusten F (xi;) välillä painokertoimella wi kerrottuna WRSS = ncX j=1 wi(yi F (xi;)): (16) Painokertoimella wi saadaan näin ollen tarkempi estimointi siihen osaan mallikäy- rää, johon halutaan tarkempaa estimointia. Painokertoimen valintaan on esitetty eri vaihtoehtoja, mutta yleisesti käytetty painokerroin wi = Var(yi) = 2j on havain- tojen zi otosvarianssi kyseisessä pitoisuuspisteessä xi. Tällä menetelmällä saadaan painotettua sovitusmenetelmä niin, että käyrää sovitetaan tiukemmin pienille pi- toisuuksille, joilla on suurempi varianssi ja löyhemmin suurille havaintojen arvoille, joilla on pienempi varianssi. Näin ollen saadaan estimointituloksena tarkempi tulos verrattuna tavalliseen pienimmän neliösumman estimointiin. Painotetussa estimointimenetelmässä korvataan residuaalineliösumma S() pai- notetulla residuaalineliösummalla WRSS (16), jolloin saadaan Taylorin laajennetulla sarjakehitelmällä (9) kirjoitettua normaaliyhtälöt muotoon: V(t) 0WV(t)( t) = V(t)0W(y F (t)); (17) jossa V() on mallin parametrien  vauhtimatriisi tai toiselta tulkinnaltaan Jacobin matriisi ja W on painokertoimien wi diagonaalimatriisi, jossa diagonaalialkioina ovat painokertoimet wi. Painokerroinmatriisin W diagonaalialkoioiden wi arvojen ollessa yksi on kyseessä tavallinen pienimmän neliösumman estimointi. Painotettu varianssi-kovarianssimatriisi saadaan, kun sijoitetaan painomatriisiW ja ratkaistaan varianssi-kovarianssimatsiisi kuten yhtälössä (11), jolloin Cov(^) = 2(V(^)0WV(^))1: Painotettu pienimmän neliösumman estimointi tuottaa tarkemman estimointitulok- sen, kun ei voida olettaa variaation olevan homoskedastista. Painotettu pienimmän neliösumman estimointikin vaatii muutamia oletuksia, kuten havaintojen approk- simatiivisen normaalisuuden jokaisella pitoisuuspisteellä xi. Lisäksi mallinnettavat 17 epälineaariset neli- ja viisiparametriset logistiset regressiomallit vaativat ilmiön mon- otonisuuden. Näitä oletuksia ei kuitenkaan aina voida todentaa, jolloin yhtenä vaih- toehtona on käyttää epäparametrisia mallinnusmenetelmiä, kuten splinejä. 18 5 Epäparametrinen mallintaminen biologisessa mää- rityksessä Epäparametriset regressiomenetelmät sopivat useimpiin mallinnusongelmiin, kuten myös annosvasteilmiön mallintamiseen. Näistä käydään tässä työssä läpi splinifunk- tioiden silotettu splini (smoothing spline), joka on kolmannen asteen paloittain mää- ritelty polynomifunktio. Muita mahdollisia vaihtoehtoja ovat yksinkertaisin, en- simmäisen asteen splini ja hieman mutkikkaampi toisen asteen splini. Ensimmäi- sen asteen splini on paloittain määritelty lineaarinen funktio pitoisuuspisteiden tj (t = (t0; t1; : : : ; tnc), tj  ti+1  : : : tnc) välillä. Näitä kutsutaan solmuiksi (knots), kun taas toisen asteen splinifunktio on paloittain kvadraattinen funktio välillä T = [t0; tnc ]. Tässä työssä käytetään kolmannen asteen silotettua spliniä. 5.1 Silotettu splini Silotetun splinin tarkoituksena on löytää funktio g(z), joka kuvaa annosvasteen il- miötä tarkasti muttei ylisovita mallia otosaineistoon. Tässä työssä rajoitetaan kaik- kien mahdollisten splinifunktioiden joukkoa sisältämään ainakin kahdesti jatkuvasti derivoituvat splinifunktiot välillä T . Halutaan siis löytää funktio, joka minimoi harhan ja variaation summaa. Har- han minimoiminen tarkoittaa residuaalineliösumman RSS minimointia, kun RSS =P i P j(yij g(zij))2. Harhan minimoiminen yksinään johtaisi havaintojen suoraan interpolointiin, jolloin saataisiin RSS = 0. Havaintojen interpolointi johtaisi siis aineiston ylisovittamiseen, joka kasvattaa funktiosovituksen variaatiota. Asetetaan välin T pisteet vastaamaan standardien pitoisuuksia (t0 = x0; t1 = x1; : : : ; tnc = xnc). Tällöin tarkastellaan tilannetta, jossa kaikista kahdesti jatkuvas- ti derivoituvista funktioista g(xi) pyritään löytämään se, joka minimoi sakotetun residuaalineliösumman (penalized residual sum of squares PRSS): PRSS(g(xi); ) = nX i=1 (yi g(xi))2 +  Z g00(v)2dv;  2 (0;1): (18) PRSS minimointi suoritetaan ennalta määrätyn sakkoparametrin  > 0 avulla. Yhtälössä (18) ensimmäinen summalauseketermi tappiofunktio mittaa sopivuutta aineistoon, kun taas toinen termi sakkofunktio toisen derivaatan arvon suuruutta ja sakottaa sitä enemmän, mitä suurempi sakkoparametri  on. Sakkoparametri  on ennalta valittu, ja sen suuruus vaikuttaa ratkaisuna saatavan kolmannen asteen splinifunktion sileyteen, jota havainnoidaan kuvassa 8. Kun sakkoparametri  ! 0 saadaan (rajalla) interpoloiva kolmannen asteen splinifunktio g(). Kun taas  ! 1, täytyy R g00(v)2dv ! 0, jotta minimointilauseke saisi pienen arvon. Näin ollen suurella sakkoparametrin  arvolla sovitetaan lineaarinen malli aineistoon, koska sen toinen derivaatta on nolla. Sakkoparametrin  avulla säädellään harhan ja variaation summan suuruutta ja vaikutusta. Funktiolla g(), joka minimoi yhtälön (18) on tiettyjä erityispiirteitä. Funktio g() on paloittain määritelty kolmannen asteen polynomifunktio yksikäsitteisillä sol- mujen tj arvoilla t1; t2; : : : ; tnc . Lisäksi funktiolla g() on jatkuvasti derivoituvat en- simmäisen ja toisen asteen derivaatat g0() ja g00() jokaisessa solmussa tj. Tarpeeksi 19 Kuva 8: Annosvaste -ilmiön kuvaajaan sovitettu kaksi Silotettua splinifunktiota eri sakkoparametrin  = 0:001 ja  = 0:04 arvoilla. Pienemmällä sakkoparametrin ar- volla splinifunktio on hyvin joustava ja käy jokaisen havaintopisteen läpi. Suurem- malla  arvolla splinifunktio on jäykempi ja muistuttaa enemmän sigmoidikäyrää. tarkan ja sileän splinin löytämisen ongelma on siis käytännössä löytää sopiva sakko- parametrin  arvo. Tämän löytämiseksi jaetaan osa otosaineistosta opetusaineistoon ja loput havainnot validointiaineistoon. Opetusaineiston avulla sovitetaan haluttu silotettu splinifunktio ennalta määrätyllä sakkoparametrin  arvolla ja tällä mal- lilla ennustetaan validointiaineiston havaintojen arvoja ja lasketaan ennusteiden ja todellisten havaintojen etäisyyksien neliösumma. Tätä menetelmää kutsutaan ris- tiinvalidoinniksi (Cross-Validation) ja tässä työssä käsitellään kahta eri ristiinvali- dointimenetelmää, jotka eroavat opetus- ja validointiaineistoon jaon suhteen. 5.2 Ristiinvalidointi Ristiinvalidoinnissa pyritään jakamaan otosaineisto f(x1; y1); : : : ; (xn; yn)g satunnai- sesti H kpl samansuuruiseen osajoukkoon. Tätä kutsutaan H-kertaiseksi ristiinvali- doinniksi (h-fold Cross-Validation). Ensin valitaan osajoukko h = 1 validointiaineis- toksi ja loput f2; 3; : : : ; Hg osajoukkoa opetusaineistoon. Tällöin validointiaineiston koko on n=H havaintoa ja opetusaineiston koko taas nn=H. Opetusaineistoon sovi- tetaan paras silotettu splini ennalta määrätyllä sakkoparametrin p (p = 1; : : : ; P ) arvolla tehtävän (18) mukaisesti. Tästä lasketaan sovitetun silotetun splinimallin avulla ennusteet validointiaineiston (h = 1) havainnoille ja lasketaan ennusteneliö- virheiden summa MSEh=1 = n=HX l=1 (yl yl^)2: (19) 20 Tämän jälkeen asetetaan osajoukko h = 2 validointiaineistoksi ja loput f1; 3; : : : ; Hg osajoukkoa opetusaineistoksi ja lasketaan ennusteneliövirheiden summa MSEh=2. Tätä jatketaan asettamalla jokainen osajoukko vuorollaan validointiaineistoksi, jol- loin saadaan ennusteneliövirheiden summat jokaiselle osajoukolleMSEh=1;MSEh=2; : : : ;MSEh=H samalla sakkoparametrin p arvolla. Näistä ennusteneliövirheiden sum- mista lasketaan keskiennusteneliövirhe CV(p) = 1 h hX l=1 MSEh=l; p = 1; 2; : : : ; P; (20) jolloin saadaan keskimääräinen ennusteneliövirheen estimaatti sakkoparametrin p arvolle. Tätä toistetaan käyttämällä useita eri sakkoparametrin p arvoja, joista valitaan keskiennusteneliövirheen mielessä paras sakkoparametrin p arvo. Keskien- nusteneliövirheen mielessä paras p arvo on se, joka tuottaa pienimmän keskiennus- teneliövirheen Min(CV(1); CV(2); : : : ; CV(P )). Kuva 9: Illustraatio H-kertaisen ristiinvalidoinnin toimintaperiaatteesta. Otosaineis- to jaetaan satunnaisesti kuvan tapauksessa seitsemään osajoukkoon (H = 7), jolloin aina jokaista osajoukkoa käytetään vuorotellen validointiaineistona (keltainen laa- tikko) ja loppuja opetusaineistona (vihreät laatikot). Opetusaineiston avulla sovite- taan mahdollisimman sopiva silotettu splinimalli, jonka avulla ennustetaan validoin- tiaineiston havaintojen arvoja. Jokaisen rivin kohdalla lasketaan ennusteneliövirheet MSEh=1 validointiaineiston ja sovitetun mallin välille. Seuraavalla kierroksella vali- taan seuraava osajoukkolaatikko validointiaineistoksi ja loput taas opetusaineistoksi, jolloin saadaan laskettua MSEh=2. Tätä jatketaan, kunnes on kaikki osajoukot H käyty lävitse ja saatu laskettua MSEh=1; : : : ;MSEh=H . Yllä on kuvattu H-kertainen ristiinvalidointimenetelmä, joka on varsin suosittu sen nopeuden ansiosta, osajoukkojen lukumäärän ollessa kohtuullinen esim. H = 10 tai H = 5. Tällöin voidaan käydä läpi tehokkaasti useitakin eri sakkoparametrin p (1; : : : ; P ) arvoja eli voidaan asettaa P suureksikin. Ristiinvalidoinnin etuja on myös sen sisäänrakennettu testaus validointiaineistoa hyödyntämällä. Voidaan todeta, et- tä opetusaineiston neliövirheen arvo aliarvioi todellista neliövirheen arvoa, jolloin siis opetusaineiston neliövirhe < validointiaineiston neliövirhe. Vain opetusaineis- ton perusteella tehdyt päätelmät voisivat olla harhaisempia ja optimistisempia kuin 21 validointiaineiston kanssa tehdyt päätelmät. Kuitenkin opetus- ja validointiaineis- toon jaossa saattaa tapahtua joitain virheitä, jolloin tämän jaon vuoksi tapahtuvaa harhaa saattaa syntyä. Tällöin voidaan todeta, että mitä useampaan osajoukkoon otosaineisto jaetaan, sitä useampia sovituksia tehdään ja voidaan olla varmempia valitun sakkoparametrin arvosta. Äärimmäisinä erikoistapauksina ristiinvalidoinnista ovat otosaineiston jako vain kahteen osajoukkoon (H = 2), jolloin suoritetaan vain kaksi sovitusta jokaiselle sakkoparametrille p. Toisena ääripäänä sen sijaan on otosaineiston jako niin useaan osajoukkoon kuin havaintojen lukumäärä on, jolloin siis osajoukkojen lukumäärä vastaa otosaineiston kokoa H = n. Näin ollen jokaiseen osajoukkoon h tulee vain yksi havaintopari (xh; yh), h = 1; : : : ; H. Tätä erikoistapausta kutsutaan yksi-pois -ristiinvalidoinniksi (Leave-one-Out Cross-Validation LOOCV) ja menetelmä toimii samoin kuin H-kertainen ristiinvalidointikin ja lasketaan yhtälön (20) mukaisesti. LOOCV-tapauksessa (H = n) yhtälö (20) voidaan kirjoittaa muotoon CV(p) = 1 nc ncX j=1 MSEh=j; p = 1; 2; : : : ; P: (21) Kuva 10: Illustraatio yksi-pois -ristiinvalidoinnin (LOOCV) toimintaperiaatteesta. Otosaineisto jaetaan H = n osajoukkoon, joihin kuuluu vain yksi havainto. Jokaista osajoukkoa käytetään vuorotellen validointiaineistona (keltainen laatikko) ja loppuja opetusaineistona (vihreät laatikot). LOOCV-estimaatilla on huomattavia etuja verrattuna suurempikokoisiin osajouk- koihin eli tapauksiin, jossa H = 2 tai H = 5. Ensinnäkin LOOCV-menetelmällä lasketulla ennustekeskineliövirheen estimaatti on huomattavasti vähemmän harhai- nen, koska lähes koko otosaineisto käytetään mallin opettamiseen. Muissa tapauk- sissa käytetään vähemmän havaintoja mallin opettamiseen ja enemmän testauk- seen, jolloin saatetaan yliarvioida validointiaineiston ennustevirhettä. Toisena etuna LOOCV-menetelmällä on sen toistettavuus. Koska jokainen havainto on vuorollaan validointiaineistossa, ei osajoukkoihin jaolla ole merkitystä ja saadaan aina sama keskineliöennustevirheen arvo sakkoparametrilla p. LOOCV-menetelmän haitta- puolena sen sijaan on sen hitaus, kun verrataan H-kertaiseen ristiinvalidointiin ja esimerkiksi H = 5. Otosaineiston ollessa hyvinkin suuri, esimerkiksi n = 100000, täytyy LOOCV-menetelmän tehdä silotetun splinin sovitus n = 100000 kertaa, kun 22 taas H-kertaisessa ristiinvalidoinnissa suoritetaan sovitus vain esimerkiksi H = 5 kertaa. Mikäli sakkoparametrin p arvoja on vain muutamia ei tämä välttämättä hidasta merkittävästi, mutta jos P = 100, LOOCV suorittaa sovituksen n  P = 100000100 = 10000000 kertaa ja H-kertainen ristiinvalidointi vain viidesosan tästä. 23 6 PassingBablok -menetelmävertailu Biologisessa määrityksessä on tavoitteena löytää malli, joka kuvaisi mahdollisimman tarkasti halutun ilmiön luonnetta. Tässä työssä halutaan, että sovitettu malli ky- kenee tuottamaan signaalivasteesta käännetyn pitoisuuden arvon mahdollisimman tarkasti. Empiirisessä esimerkissä luvussa 7 tarkastellaan epälineaaristen mallien ja silotetun splinin sopivuutta vertaamalla simulointiaineiston kontrollipisteiden pitoi- suuksien arvoja keskenään. Tätä kutsutaan menetelmävertailuksi (method compa- rison), jossa vertaillaan kahta menetelmää keskenään (käyrä, jolla simuloitu aineisto on luotu vs. sovitettu malli). Menetelmävertailu voidaan suorittaa mm. pääkom- ponenttianalyysillä, lineaarisella regressiomallilla tai vaihtoehtoisesti Blandt-Alman -kontrollikuvaajalla. Biologisessa määrityksessä on kuitenkin vakiintunut menetel- mävertailuksi PassingBablok -regressiomenetelmä [11], [12], joka olettaa kahden menetelmän välille lineaarisen rakenteellisen suhteen. Merkitään koeasetelmaa, jossa mitataan n = nc  k suuruinen otos toisistaan riippumattomia mittauksia populaatiosta, kuten luvussa 2.3 on kuvattu. Merkitään menetelmän 1 pitoisuuksia satunnaismuuttujalla X ja vastaavasti menetelmän 2 pi- toisuuksia satunnaismuuttujalla Z. Näin saadaan i:nnelle mittaukselle arvot xi ja zi, jotka ovat realisaatioita satunnaismuuttujista X ja Z, jossa X on tunnetun käyrän arvot ja Z taas mallin sovitteiden arvot. Tällöin voidaan kuvata jokainen satunnai- nen muuttuja kahden komponentin summana, jossa ensimmäinen muuttuja edus- taa populaation kaikkien mahdollisten näytteiden odotusarvon vaihtelua ja toinen muuttuja kuvaa mittausvirheen variaatiota kyseiselle otokselle. Näin ollen voidaan kuvata mittauksen i arvot: xi = x  i + i ja zi = z  i + i; (22) jossa xi ja z  i kuvaavat kyseisen mittauksen odotusarvoa ja i ja i mittausvirhettä. Kun kahden menetelmän välillä on lineaarinen rakenteellinen suhde, voidaan tämä kuvata lineaarisella mallilla: zi = + x  i ; jossa on mallin vakiotermi ja kulmakerroin. Tavoitteena on tällöin suorittaa n:stä mittausparista (xi; zi) menetelmävertailu seuraavasti: 1. estimoidaan parametrit ja ; 2. testataan mallin lineaarisuus; 3. testataan hypoteesi = 1; 4. testataan hypoteesi = 0. Mikäli kaikki testit hyväksytään, voidaan sanoa kahden menetelmän olevan saman- kaltaisia keskenään eli zi = x  i ja näin ollen molemmat menetelmät tuottavat samat pitoisuudet tarkastellulla välillä. Kuten yllä todettiin, mallinnetaan menetelmävertailututkimuksissa kahta mene- telmää lineaarisella mallilla. Biologisessa määrityksessä menetelmävertailua suori- tettaessa havaitaan, että: 24 Kuva 11: Identiteettisuora, PassingBablok -regressiosuora sekä lineaarinen malli sovitettuna kahden menetelmän pitoisuuksiin. ˆ Menetelmät 1 ja 2 sisältävät molemmat satunnaisvaihtelua. ˆ Mittausvirheiden jakaumat ovat harvoin normaalisti jakautuneita. ˆ Mittauksien odotusarvot xi ja z  i eivät ole normaalisti jakautuneita, koska mittaukset mitataan monelta eri pitoisuuspisteeltä. ˆ Äärimmäiset arvot (oudokit) eivät välttämättä ole mittausvirheitä, vaan ne saattavat olla eri menetelmien ominaisuuksia. Tästä syystä näitä arvoja ei saa poistaa ilman painavaa kokeellista syytä. ˆ Mittausvirheiden variaatio ei ole vakio yli pitoisuuspisteiden. Variaatio kasvaa suuremmaksi, kun pitoisuus kasvaa, kuten luvussa 2.2 todettiin. Tähän H. Passing ja W. Bablok kehittivät menetelmän, joka pystyy suorittamaan menetelmävertailun tarkasti, vaikka yllä mainitut ominaisuudet ovat läsnä. Vaikka menetelmä on lineaarinen regressiomalli, se eroaa parametrien ja estimoinnin suhteen tavallisiin lineaarisiin menetelmiin verrattuna. Tätä eroa havainnollistetaan kuvassa 11, jossa viimeisimmät havainnot vääntävät lineaarisen mallin kauemmas 25 identiteettisuorasta, kun taas PassingBablok -regressiosuora ei juurikaan eroa iden- titeettisuorasta. 6.1 Parametrien estimointi Ennen PassingBablok regressiomenetelmän sovittamista oletetaan, että xi ja z  i ovat satunnaismuuttujien X ja Z odotusarvoja mielivaltaisesta jatkuvasta jakau- masta. Oletetaan myös, että i ja i ovat mittausvirheiden realisaatioita, joiden jakaumat ovat samantyyppisiä. Lisäksi mittausvirheiden i ja i varianssien 2 ja 2 ei tarvitse olla vakioita pitoisuusvälillä, mutta niiden tulisi pysyä suhteellisina keskenään: 2 2 = 2: Nämä oletukset takaavat tarkan ja luotettavan parametrien estimoinnin ja hypo- teesintestauksen, kun  1 [11]. Tällöin saadaan estimoitua parametrien ja estimaatit. Aluksi tulee estimoida vakiotermin ja kulmakertoimen arvot. Kahden pis- teen välisen suoran kulmakerrointa voidaan hyödyntää, kun estimoidaan kulmaker- roinparametrin arvoa. Tällöin siis kahden pisteen välisen suoran kulmakerroin saadaan: Wij = zi zj xi xj ; kun 1  i < j  n; (23) jossa i ja j kuvaavat kahden eri pisteen arvoja. Tällöin, kun tarkastellaan koko ai- neiston läpi kaikki mahdolliset kulmakertoimet saadaan yhteensä lukumääräksi n 2  . Näistä kaikista kulmakertoimista Wij lasketaan korjattu mediaani, jonka laskemi- sessa on muutamia rajoitteita: ˆ Jos zi = zj ja xi = xj, saa Wij = 00 , joka on määrittelemätön. Nämä parit poistetaan estimoinnista. ˆ Jos zi > zj ja xi = xj, annetaanWij arvoksi jokin hyvin suuri positiivinen luku. Koska kulmakertoimista lasketaan mediaani, ei suuri arvo ole ongelmallinen mediaanin määrityksessä. ˆ Jos zi < zj ja xi = xj, annetaan Wij arvoksi jokin hyvin suuri negatiivinen luku. Koska kulmakertoimista lasketaan mediaani, ei suuri negatiivinen arvo ole ongelmallinen mediaanin määrityksessä. ˆ Jos Wij = 1, poistetaan tämä pari kulmakertoimen laskusta symmetriasyi- den nojalla [11]. Tämä korjaa harhaa, joka johtuu riippumattomuusoletuksen puuttumisesta arvojen Wij välillä. ˆ Jos Wij < 1, käytetään tätä havaintoparia estimointiin, mutta lisätään yksi laskuriin L. Tätä laskurin L arvoa käytetään mediaanin korjaamiseen. Kun kaikki yllä olevat askeleet on huomioitu järjestetään arvot Wij kasvavaan suu- ruusjärjestykseen W(1)  W(2)  : : :  W(N). Arvot Wij:n laskemisessa eivät ole 26 keskenään riippumattomia, jolloin niiden mediaanin käyttäminen suoraan kulma- kerroinparametrin estimaattorina tuottaa harhaisen tuloksen. Tästä syystä käy- tetään laskurin L arvoa korjaamaan tätä harhaa, jolloin kulmakertoimen arvo saadaan laskettua käyttämällä siirrettyä mediaania arvoista W(i): ^ = ( WN+1 2 +L; kun N pariton 1 2  (WN 2 +L +WN 2 +1+L); kun N parillinen: Kaksisuuntainen luottamusväli kulmakerroinestimaattorille ^ saadaan luotta- mustasolla koottua, kun merkitään standardinormaalijakauman (1 =2):ttä kvan- tiilia q =2. Tällöin saadaan laskettua luottamusvälin indeksit C: C = q =2 r n(n 1)(2n+ 5) 18 ; jota siirretään, M1 = N C 2 ; M2 = N M1 + 1: Tässä M1 on pyöristetty lähimpään kokonaislukuun, jolloin myös M2 on kokonaislu- ku ja N on havaintoparien kulmakertoimien todellinen lukumäärä. Tällöin saadaan luottamusväli laskettua parametrille ^ käyttämällä järjestettyjä kulmakertoimien Wij arvoja: bL = W(M1+L)  ^  W(M2+L) = bU : (24) Siirtolaskurin L arvo on siis havaintojen Wij < 1 lukumäärä, joka vastaa nol- lahypoteesia H0 : = 1 [11]. Tällöin tässä tapauksessa ^ on hyvä ja luotettava estimaattori kulmakerroinparametrista . Kulmakerroinparametrin estimaattorin ^ avulla lasketaan vakiotermin estimaat- tori ^ lineaarisesta mallista, kuten yhtälössä (2). Ratkaisemalla yhtälö vakiotermin suhteen ja asettamalla kulmakertoimen tilalle estimoitu parametri ^, saadaan jo- kaiselle havainnolle erikseen oma vakiotermin arvo: i = zi ^xi: Vakiotermit järjestellään myös suuruusjärjestykseen (1)  (2)  : : :  (N), joista lasketaan mediaani. Tämä mediaani on tällöin vakiotermin ^ estimaatti. Kun ollaan saatu estimoitua parametreille ^ ja ^ arvot, voidaan laskea näille luottamusvälit, joita käytetään menetelmien sopivuuden tarkasteluun. Vastaavasti taas saadaan laskettua luottamusväli vakiotermin estimaattorille ^ samalla luottamustasolla käyttämällä kulmakerroinparametrin luottamusvälin ala- ja ylärajan arvoja (bL; bU). Vakiotermin ^ luottamusvälin (aL; aU) estimointi vaatii, että vähintään puo- let havaintopisteistä sijaitsee regressiosuoran yläpuolella tai suoralla. Tällöin, koska (X;Z) on jatkuva kaksiarvoinen muuttuja, niin myös sama lukumäärä havaintoja sijaitsee regressiosuoran alapuolella tai suoralla todennäköisyydellä 1. Havaintopiste (xi; zi) sijaitsee suoran yläpuolella vain, kun a < zi ^xi, jolloin tulee todistettua, että: ^ = mediaani(zi ^xi) 27 on vakiotermiparametrin estimaattori. Tätä tietoa hyödyntämällä saadaan vakio- termiparametrin luottamusväli johdettua: aL = mediaani(zi bUxi) aU = mediaani(zi bLxi); (25) jossa aL kuvaa luottamusvälin alarajaa ja vastaavasti aU ylärajaa. Näitä luottamus- välejä (bL; bU) ja (aL; aU) käyttämällä voidaan testata nollahypoteeseja H0 : = 1 ja H0 : = 0, kun ollaan ensin varmistuttu, että menetelmien X ja Z välinen ilmiö on todellisuudessa lineaarista. Lisäksi voidaan havaita, että estimointimenetelmä ei ole riippuvainen menetelmien X ja Z järjestyksestä, koska estimaatit voidaan laskea kaikille n havaintoparille (xi; zi), kun merkitään: z = ^ + ^x x = A+Bz: Yllä on laskettu estimaattorit parametreille ja , joilla on seuraava ominaisuus: B = 1 ^ A = ^ ^ : Tästä voidaan päätellä, että on estimoinnin kannalta samantekevää kumpaa me- netelmää merkitään X:llä ja kumpaa Z:llä. Kun parametrit on estimoitu, täytyy varmistua, että menetelmien X ja Z välinen suhde on todellisuudessa lineaarista. Tämä tarkastelu on tarpeellista, koska menetelmiä vertaillaan lineaarisella mallilla, jolloin epälineaaristen suhteiden tarkastelu on ongelmallista. Lineaarisuutta testat- taessa tulee tarkastella regressiosuoran sopivuutta aineistoon tai tarkastella, kuinka satunnaisesti aineisto on hajautunut sovituksen z = ^ + ^x ympärillä. 6.2 Lineaarisuustarkastelu Lineaarisuutta voidaan tarkastella muokatulla Cusum-testillä, jossa tarkastellaan regressiosuoran kummallakin puolella peräkkäisten havaintojen etäisyyden kumu- latiivista summaa. Jos menetelmien X ja Z välillä on epälineaarista käyttäytymistä, voidaan havaita liian monta peräkkäistä mittausta joko sovitetun suoran ylä- tai alapuolella. Merkitään npos:llä havaintojen lukumäärää, jotka ovat sovitetun mal- lin yläpuolella (zi > ^ + ^xi) ja nneg:llä merkitään vastaavasti mallin alapuolelle jääviä havaintoja (zi < ^ + ^xi). Jokaiselle havaintopisteparille (xi; yi) lasketaan pistemäärä ri: ri = r nneg npos ; kun yi > ^ + ^xi; ri = r npos nneg ; kun yi < ^ + ^xi; ri = 0; kun yi = ^ + ^xi: Kuitenkin pistemäärien ri avulla laskettu tulos riippuu havaintopisteiden (xi; yi) jär- jestyksestä. Jos havaintopisteet järjestetään kasvavaan järjestykseen menetelmän X 28 mukaan saadaan eri tulos kuin menetelmän Z järjestyksellä. Tästä syystä järjestä- minen on parempi suorittaa havaintopisteiden (xi; yi) ja sovitetun suoran z = ^+ ^x välisen etäisyyden mukaan kasvavaan järjestykseen. Etäisyys Di saadaan laskettua projisoimalla: Di = yi + 1 ^ xi ^q 1 + 1 ^ 2 : Tällöin voidaan pistemäärät ri järjestää etäisyyksien mukaan kasvavaan järjestyk- seen D(1)  D(2)  : : :  D(N). Tämä menetelmä takaa, että lineaarisuustarkastelu ei ole riippuvainen kummankaan menetelmän järjestyksestä, vaan ainoastaan etäi- syydestä. Lineaarisuustarkastelun cusum-testisuure saadaan laskettua kumulatiivi- sen summan mukaan. Tarkastellaan koordinaatistojärjestelmää, jossa x-akseli kuvaa etäisyyksien Di järjestystä 1 N ja y-akseli taas kuvaa pistemäärien ri kumulatii- vista summaa. Kumulatiivinen summa saadaan laskettua indeksille i: cusum(i) = iX l=1 r(l); joka kertoo kyseisen indeksin positiivista tai negatiivista pistemäärän summaa pis- teiden 1 i välillä. Satunnainen pistemäärien järjestäminen johtaisi kohtuullisia jcusum(i)j arvoja, kun taas peräkkäiset negatiiviset tai positiiviset pistemäärien ri arvot johtavat suurempiin jcusum(i)j arvoihin. Cusum-testisuureen laskemisen jälkeen etsitään maksimaalinen peräkkäinen poik- keaminen regressiosuoralta max(jcusum(i)j). Mikäli valittu testisuureen arvo cusum(i) < 0, verrataan tätä niiden pistemäärien ri:n osajoukon jakaumaan, jos- sa ri < 0 ja vastaavasti verrataan cusum(i) > 0 osajoukon ri > 0 jakaumaan. Tätä testisuureen arvoa verrataan Kolmogorov-Smirnov -jakauman kriittiseen pisteeseen h luottamustasolla , taulukko 2. Mikäli testisuureen ja kriittisen pisteen välillä, jollekin indeksille i, (i = 1; : : : ; n), jcusum(i)j h p nneg + 1; voidaan todeta ettei satunnaista hajontaa regressiosuoran ympärillä ole ja lineaari- suusoletus hylätään havaittujen menetelmien x ja z välillä kyseisen otoksen perus- teella. Jos taas lineaarisuusoletus on jää voimaan, voidaan testata parametrien esti- maattorien ^ ja ^ arvoja. Suurilla otoksilla saattaa kuitenkin olla ongelmia cusum- lineaarisuustestin kanssa, koska suuressa otoksessa kaikki minimaalisetkin erot tu- levat tilastollisesti merkitseviksi ja johtavat hypoteesin hylkäämiseen.. Suuret testi- suureen arvot menevät yli Kolmogorov-Smirnov -jakauman kriittisen pisteen ja näin ollen lineaarisuusoletus tulee hylätyksi, vaikka menetelmien välillä olisikin todelli- suudessa lineaarinen suhde. 6.3 Hypoteesintestaus PassingBablok -regressiosuoralle PassingBablok regressiosuoran nollahypoteesina kulmakerroinparametrin estimaat- torin ^ oletetaan saavan arvon yksi, H0 : ^ = 1, jolloin kahden menetelmän väliset arvot (xi; yi) ovat keskenään samat. Tätä saadaan testattua tarkastelemalla lasket- tua luottamusväliä 24 halutulla luottamustasolla. Luottamusvälin sisältäessä arvon 29 Taulukko 2: Kolmogorov-Smirnov -jakauman kriittisten pisteiden arvot tietyllä luot- tamustasolla . h 1 % 1:63 5 % 1:36 10 % 1:22 yksi, bL  1  bU , voidaan olettaa menetelmien välisen suhteen olevan samankaltais- ta. Jos taas luottamusväli ei sisällä arvoa yksi, hylätään nollahypoteesi ja voidaan todeta, että menetelmien välillä on ainakin suhteellinen ero. Tämä testaus ei kuiten- kaan ole täysin riippumaton alla piilevien jakaumien suhteen, mutta testaus tuottaa yleisesti tarpeeksi luotettavia tuloksia [11]. Kulmakerroinparametrin luottamusvälin sisältäessä arvon yksi, voidaan testata vakiotermiparametrin arvoa. Nollahypoteesina vakiotermiparametrin estimaattorille H0 : ^ = 0, joka varmistaa, ettei kahden menetelmän välillä ole vakiotermin suu- ruista harhaa. Tästä voidaan varmistua, mikäli vakiotermiparametrin luottamusväli aL  0  aU sisältää arvon nolla. Jos taas luottamusväliin ei kuulu arvoa nolla, hylätään nollahypoteesi ja todetaan menetelmien eroavan ainakin vakion verran ja tuottavan harhaa menetelmien välille. Mikäli hyväksytään molemmat hypoteesit ^ = 1 ja ^ = 0, päätellään z = x ja menetelmien olevan identtiset. 30 7 Simulointiesimerkki Tässä työssä suoritettiin esimerkki simuloimalla oikeasta mittausdatasta johdettu aineisto, johon oli lisätty kohinaa ja harhaa. Aineiston luonti ja analysointi suoritet- tiin SAS JMP ohjelmistolla. Simulointiesimerkin tavoitteena oli, illustroida kuinka sigmoidikäyrän muotoiset mallit sekä silotettu splinimalli sopivat standardikäyräksi. Standardikäyrän sovitukseen käytettiin standardipitoisuuspisteitä, joista jokaisesta mitattiin haluttu määrä toistomittauksia. Tämä toimi opetusaineistona, kuten lu- vussa 5.2. Tämän lisäksi simuloitiin validointiaineistoksi 50 kappaletta ns. kontrol- lipisteitä tai -mittauksia tasavälisesti logaritmisella asteikolla sijoitettuja pitoisuus- pisteitä. Nämä kontrollipisteet oli sijoitettu siten, että pienimmän kontrollipisteen pitoisuus oli sama kuin pienimmän standardimittauksen pitoisuuden, sekä vastaa- vasti suurimman kontrollipisteen pitoisuus oli sama kuin suurimman standardimit- tauksen pitoisuus. Taulukko 3: Simulointiaineiston 14 toistomittauksen otosvarianssit yhteensä kahdek- sassa standardipitoisuuspisteestä. Log-standardipitoisuus Log-signaalivasteen otosvarianssi 2:046 0:25455 1:538 0:01888 1:102 0:00662 0:770 0:00085 0:215 0:00548 0:461 0:00006 1:204 0:00010 1:830 0:00008 Parhaimman mallin löytämisen lisäksi tarkasteltiin toistomittauksien lukumää- rien vaikutusta ja valittiin sopiva toistomittauksien lukumäärä. Toistomittauksien sopiva lukumäärä tuottaa riittävän pienen harhan, mutta ei ole liian suuri kustan- nuksien kannalta. Sopivan lukumäärän löydyttyä simuloitiin vielä yksi aineisto käyttäen valittua toistomittauksien lukumäärää. Tämän jälkeen sovitettiin sigmoidikäyrämallit ja si- lotettu splini sekä arvioitiin kontrollipisteiden pitoisuudet. Näitä kontrollipisteiden arvioituja pitoisuuksia vertailtiin PassingBablok menetelmävertailulla referenssipi- toisuuksiin verrattuna, jotka saatiin suoraan simulointiaineiston simulointiin käyte- tystä käyrästä. Aineisto simuloitiin käyttämällä kahdeksaa standardipitoisuuspistettä, 50 kont- rollipistettä ja kuutta mittauslaitetta. Toistomittauksien lukumäärä vaihteli kolmes- ta mittauksesta kahteenkymmeneen, jolloin yhden simuloinnin tuottama testiaineis- ton koko vaihteli 8 3 6 = 144 ja 8 20 6 = 960 välillä ja vastaavasti validoin- tiaineiston koko 900 6000 välillä. Standardi- ja kontrollipisteiden signaalit saatiin toisen asteen polynomifunktiosta, johon oli lisätty normaalijakautunutta kohinaa ja jonka variaatio kasvoi pitoisuuden kasvaessa. Tällöin saatiin heterogeenisia nor- maalisti jakautuneita mittaustuloksia jokaiselta pitoisuuspisteeltä. Lisäksi jokaisella toistomittauksen lukumäärällä simuloitiin 150 aineistoa. 31 Kuva 12: Keskimääräisen ennustevirheen muutos toistomittauksien funktiona, kun validointiaineistona on 50 kontrollipitoisuuspistettä. Logaritmisten signaalivasteiden otosvarianssi Var(yi) laskettiin jokaiselta stan- dardipitoisuuspisteeltä. Otosvarianssien avulla laskettiin painokertoimet pienimmän neliösumman estimointia varten. Logaritmisten signaalivasteiden otosvarianssit on kuvattu taulukossa 3. Logaritmisten signaalivasteiden otosvariansseista havaitaan, että ne pienenevät logaritmisen pitoisuuden kasvaessa. Tällöin, kun otosvarianssit si- joitetaan painotetun pienimmän neliösumman menetelmän kaavaan 16, niin pienien pitoisuuksien logaritmisilla signaalivasteilla on suurempi painokerroin kuin suurem- milla arvoilla. Jokaiselta toistomittauksen lukumäärältä (3 20) simuloitiin 150 eri aineistoa ja näihin jokaisen aineiston opetusaineistoksi valittiin standardipisteiden toistomit- taukset. Opetusaineistoon sovitettiin painotetulla pienimmän neliösumman mene- telmällä neliparametrinen-, viisiparametrinen- sekä silotettu splinimalli. Silotetulle splinimallille etsittiin jokaiselle aineistolle oma sakkoparametrin  arvo käyttäen H- kertaista ristiinvalidointia. Ristiinvalidoinnin osajoukkoihin jako tapahtui standar- dipitoisuuspisteiden mukaisesti, jolloin jokaisen standardipitoisuuspisteen toistomit- taukset olivat vuorotellen ristiinvalidoinnin validointiaineistona ja loput pisteet ope- tusaineistona. Sovitettujen mallien avulla ennustettiin viidenkymmenen kontrolli- 32 Kuva 13: silotettu splini, neliparametrinen logistinen malli sekä viisiparametrinen lo- gistinen malli sovitettuna simulointiaineiston opetusaineistoon, kun toistomittauk- sien lukumäärällä 14. pisteen pitoisuudet, joita verrattiin simuloituihin referenssipitoisuuksiin. Näin ollen saatiin laskettua jokaiselle toistomittaukselle keskimääräiset ennustevirheet, jotka esitetään kuvassa 12. Kuvasta 12 havaitaan, että silotetulle splinille sekä neliparametriselle logistisel- le mallille saadaan pienin keskimääräinen ennusteneliövirhe 18 toistomittauksella, kun taas viisiparametriselle logistiselle mallille minimi saatiin 20 toistomittauksella. Havaittiin myös, että silotettu splinimalli tuottaa jatkuvasti pienimmän ennuste- neliövirheen, muutamaa poikkeusta lukuun ottamatta. Kuvasta 12 nähdään myös, ettei keskimääräinen ennustevirhe laske juurikaan 14 toistomittauksen jälkeen, joten toistomittauksien lukumääräksi valittiin 14. Simuloitiin vielä yksi aineisto asettamalla toistomittausten lukumääräksi 14. Aineistosta laskettiin otosvarianssit jokaiselle standardipitoisuuspisteelle toistomit- tauksien logaritmisista signaalivasteista. Nämä lasketut otosvarianssit asetettiin pai- notetun pienimmän neliösumman painokertoimiksi. Aineistolle sovitettiin neli- ja viisiparametrinen logistinen regressio malli sekä siloitettu splinimalli, joita esitetään kuvassa 13. Mallit sovitettiin kuten yllä käyttämällä standardimittauspisteiden lo- 33 (a) Silotetun splinimallin pitoisuuksien vertailu (b) Neliparametrisen mallin pitoisuuksien vertailu(c) Viisiparametrisen mallin pitoisuuksien vertailu Kuva 14: Menetelmävertailu simulointiaineiston referenssipitoisuuden, (a) silotetus- ta splinistä, (b) neliparametrisesta- ja (c) viisiparametrisesta mallista arvioitujen pitoisuuksien välillä. Menetelmävertailussa kuvattuna identiteettisuora, Passing Bablok regressiosuora sekä lineaarinen malli. garitmisia signaalivaste ja -pitoisuuspareja opetusaineistona. Opetusaineistoon so- vitetuilla malleilla arvioitiin kontrollimittauksille pitoisuudet. Kuvasta 13 nähdään, että jokainen malli sopii hyvin opetusaineiston standardipitoisuuspisteiden toisto- mittauksiin. Validointiaineiston kontrollimittausten ennustevirheen perusteella silo- tettu splinimalli tuotti pienimmän ennustevirheen 7:3. Suurin keskimääräinen en- nustevirhe saatiin neliparametrinella logistinella mallilla 10:4 ja viisiparametrinen malli lähes saman kuin silotettu splinimalli 7:8. Jokaisesta sovitetusta mallista arvioiduille pitoisuuksille suoritettiin Passing Bablok menetelmävertailu (kuviat 14 (a) - (c)). Lisäksi kuviin 14 (a) - (c) on lisätty PassingBablok regressiomalli sekä identiteettisuora, jossa vakiotermi on 0 ja kul- makerroin 1. Identiteettisuora kuvaa täydellistä riippuvuutta menetelmien välillä, jolloin voidaan menetelmien sanoa olevan identtisiä keskenään. 34 PassingBablok -regressiomenetelmällä estimoitujen vakiotermien ja kulmaker- toimien arvot sekä 95 % luottamusvälit on annettu taulukossa 4. Näiden paramet- rien avulla varmistettiin mallin lineaarisuusoletus laskemalla Cusum-testisuureen arvot 5. Jokaiselle mallille havaittiin, että Cusum-testisuureen arvo oli huomatta- vasti kriittistä pistettä suurempi. Tällöin Cusum-testisuureen mukaan hylätään li- neaarisuusoletus referenssipitoisuuksien ja arvioitujen pitoisuuksien välillä. Otosko- ko n = 4158 oli hyvinkin suuri, jolloin Cusum-testisuureen lineaarisuustarkastelu saattoi olla harhaanjohtavaa. Taulukko 4: PassingBablok -regressiomallin parametriestimaatit referenssipitoisuu- den ja arvioitujen pitoisuuksien välillä sekä parametrien ja 95 % luottamusvälit. Sovitus Luottamusväli Luottamusväli Splini 0:0009 (0:0007; 0:00011) 1:0334 (1:0321; 1:0349) Neliparametrinen malli 0 (0:0002; 0:0003) 1:0289 (1:0282; 1:0297) Viisiparametrinen malli 0:0014 (0:0016;0:0013) 1:0399 (1:0389; 1:0409) PassingBablok -regressiomenetelmän hypoteesintestaus suoritettiin, kuten lu- vussa 6.3 on esitetty. Mikäli kulmakerroinparametrin luottamusväli sisältää arvon yksi ja vakiotermiparametrin luottamusväli sisältää arvon nolla, voidaan hyväk- syä nollahypoteesi menetelmien samankaltaisuudesta. Taulukosta 4 havaitaan, että yhdenkään mallin kulmakerroinparametrin 95 % luottamusväli eivät sisällä arvoa yksi. Silotetun splinimallin ja viisiparametrisen logistisen mallin vakiotermiparamet- rien luottamusvälit eivät sisältäneet arvoa nolla, mutta neliparametrisen logistisen mallin luottamusväli sisälsi. Tällöin hylättiin nollahypoteesi menetelmien samankal- taisuudesta 5 % luottamustasolla jokaisella sovitusmenetelmällä. Kulmakerroinpa- rametrien arvot olivat kuitenkin hyvinkin lähellä arvoa yksi ja vakiotermien ar- vot olivat käytännössä nolla. Molempien parametrien luottamusvälit olivat erittäin tiukkoja suuren otoskoon takia. Taulukko 5: Cusum-testisuureen arvo ja Kolmogorov-Smirnov -jakauman 5 % luot- tamustason kriittinen piste. Sovitus Cusum-testisuure Kriittinen piste Splini 880:44 51:21 Neliparametrinen malli 1457:45 40:18 Viisiparametrinen malli 1347:37 41:61 35 8 Päätelmät Biologisessa määrityksessä standardikäyrän tarkka sovittaminen on olennainen osa päätelmien tekoa ja pitoisuuksien arviointia liuoksesta. Biologisia määrityksiä on useita erilaisia. Niistä etenkin aikaerotteinen fluoresenssi-immunomääritys hyödyn- tää ns. standardikäyrää halutun aineen pitoisuuden arviointiin mittalaitteen anta- man signaalivasteen avulla. Standardikäyrä tulee määrittää mahdollisimman tarkas- ti, jotta arvioidut pitoisuudet olisivat tarkkoja. Standardikäyrä määritetään erilli- sellä määritysaineistolla etukäteen, ennen määritettävän näytteen pitoisuuden ar- viointia. Määritysaineisto on kuitenkin harvoin homogeenista ja täysin satunnaista, mikä täytyy huomioida mallinnuksessa. Heterogeenisen variaation lisäksi pitoisuuden ja vastesignaalin riippuvuus on harvoin lineaarista. Aineisto saatetaan saada logaritmi- muunnoksilla lineaariseksi tai ainakin tasavälisemmäksi. Muuntamattomalla määri- tysaineistolla erot pienimpien mittausten ja suurimpien mittausten välillä ovat suu- ria. Tällöin suurilla mittauksilla on huomattavan suuri vipuvoima mallintamisessa, mutta tämä saadaan logaritmimuunnoksella pienemmäksi. Määritysaineisto on harvoin kuitenkaan täysin lineaarista, kun toistomittauk- sia mitataan usealta eri pitoisuuspisteeltä. Tällöin täytyy turvautua epälineaarisiin menetelmiin. Epälineaarisista menetelmistä etenkin sigmoidikäyrän muotoiset mal- lit, kuten neli- ja viisiparametriset logistiset regressiomallit ovat yleisesti käytettyjä biologisessa määrityksessä. Mallien parametrit saadaan estimoitua painotetulla pie- nimmän neliösumman menetelmällä, joka huomioi määritysaineiston heterogenian. Lisäksi painottamalla saadaan painotettua pienien pitoisuuksien arvoja mallin so- vituksessa, jossa halutaan standardikäyrän olevan mahdollisimman tarkka. Epäline- aaristen mallien tapauksessa parametriestimointi tapahtuu myös numeeristen mene- telmien avulla iteratiivisesti ja tarkemmin sanottuna GaussNewton mentelmällä. Sigmoidikäyrän muotoisten mallien lisäksi silotettuja splinimalleja voidaan käyt- tää standardikäyrän luomiseen. Tällöin ongelmaksi käytännössä tulee sopivan sak- koparametrin arvon löytäminen, jonka etsimiseen käytetään ristiinvalidointimene- telmiä. Tässä työssä käytettiin 8-kertaista ristiinvalidointia, jossa määritysaineisto jaetaan pienempiin osajoukkoihin pitoisuuspisteittäin. Osajoukkojen toistomittaus- ten arvoja pyritään vuorotellen ennustamaan, kun silotettu splinimalli opetetaan käyttämällä loppujen pitoisuuspisteiden toistomittauksien arvoja. Sopivan mallin löydyttyä halutaan varmistua mallin sopivuudesta. Tämä tapah- tuu käytännössä vertailemalla sovitetusta standardikäyrästä arvioituja pitoisuuk- sia referenssipitoisuuksiin, jotka on saatu esimerkiksi vanhemmasta standardikäy- rästä arvioimalla. Tähän menetelmävertailuun on kehitetty useita tapoja, joista PassingBablok -regressiomenetelmä on vakiinnuttanut paikkansa biologisessa mää- rityksessä. PassingBablok -regressiomenetelmä on lineaarista mallia luotettavampi regressiosuora, koska se ei ole herkkä mahdollisesti suuren vipuvoiman omaaville havainnoille. PassingBablok -regressiomenetelmä kuitenkin olettaa lineaarisen suh- teen menetelmien välille, jota saadaan testattua Cusum-lineaarisuustestillä. Cusum- lineaarisuustesti kuitenkin hylkää herkästi suurella otoskoolla lineaarisuusoletuksen, vaikka menetelmien välinen suhde olisikin todellisuudessa lineaarista. Simulointiaineiston perusteella havaittiin, että toistomittauksien lukumäärän kas- 36 vaessa saadaan tarkempi sovitus keskimääräisellä ennustevirheellä mitattuna. Kui- tenkin kymmenen toistomittauksen jälkeen ennustevirheen ero alkaa pienenemään ja 14 toistomittauksen jälkeen ero on minimaalista. Lisäksi havaittiin, että silotet- tu splinimalli tuottaa pienimmän keskimääräisen ennustevirheen ja sigmoidikäyrän muotoisista malleista neliparametrinen logistinen malli vaikuttaisi olevan yleisempää viisiparametrisempaa logistista mallia keskimäärin hieman parempi. Sovitetuista malleista arvioituja pitoisuuksia vertailtaessa referenssipitoisuuksiin PassingBablok -regressiomenetelmällä havaitaan, ettei yksikään sovitus ole Cusum- lineaarisuustestin mukaan lineaarinen. Kuvaa katsellessa kuitenkin suhde vaikut- taisi olevan lineaarista. Kun tiedetään Cusum-lineaarisuustestin hylkäävän herkäs- ti lineaarisuusoletuksen suurella otoskoolla, voidaan suhteen olettaa olevan lineaa- rista. PassingBablok -regressiomenetelmän hypoteesintestaus mallien samankaltai- suudesta hylätään myös, koska parametrien luottamusvälit eivät sisällä arvoa yksi ja nolla. Parametrien estimaatteja lähemmin tarkasteltaessa tarkemmin voitiin kui- tenkin vakiotermien olevan käytännössä nolla jokaisella sovituksella sekä kulmaker- toimien olevan lähellä yhtä. Luottamusvälit eivät sisällä nollahypoteesin mukaisia arvoja, koska otoskoon ollessa suuri ovat luottamusvälit erittäin kapeat. Simulointiaineiston perusteella todettiin, että jokainen sovitettu malli sopii ai- neistoon, mutta silotettu splinimalli tuottaa huomattavasti pienimmän keskimää- räisen ennusteneliövirheen. Kuitenkin yksinkertaisemmat sigmoidikäyrän muotoiset mallit sopivat aineistoon hyvin, mutta kumpikaan näistä malleista ei vaikuttaisi so- pivan toista mallia paremmin. PassingBablok -regressiomenetelmän mukaan neli- parametrinen logistinen malli toimii kuitenkin parhaiten, kun taas silotetusta splini- mallista arvioitujen pitoisuuksien suhde referenssipitoisuuksiin on eniten lineaarista. Kaiken kaikkiaan silotettu splinimalli näyttäisi sopivan parhaiten standardikäyrän sovittamiseen, mutta muita menetelmiä ja malliperheitä on syytä tarkastella ennen standardikäyrän käyttöönottoa. 37 Viitteet [1] Mzolo T. V. Statistical methods for the analysis of bioassay data, Technische Universiteit Eindhoven, 2016 [2] Finney D. J., Statistical Method in Biological Assay, 2nd ed, 2nd impression, Charles Griffin & Company, 1971 [3] Finney D.J. Bioassay and the Practice of Statistical Inference, International Statistical Institute (ISI), Vol 47, No 1, 1979 [4] Huang H. Nonlinear Regression Analysis, International Encyclopedia of Educa- tion, 3rd ed., Oxford: Elsevier, 2010 [5] Paul G. Gottschalk, John R. Dunn The five-parameter logistic: A characteriza- tion and comparison with the four-parameter logistic, Analytical Biochemistry 343 (2005) 5465, Brendan Technologies, Inc., 2236 Rutherford Road, Suite 107, Carlsbad, CA 92008, USA, 2005 [6] R. J. Carrol , D. Ruppert Transformation and weighting in regression, New York, NY: Chapman and Hall, 1988 [7] J. Berkson, Estimation by least squares and by maximum likelihood, Berkeley Symposium on Mathematical Statistics and Probability, pp. 111, 1956 [8] T. Hastie, R. Tibshirani, D. Witten, G. James An Introduction to Statistical Learning with Applications in R, 7th printing, Springer Texts in Statistics, 2013 [9] T. Hastie, R. Tibshirani, J. Friedman The Elements of Statistical Learning: Data Mining, Inference, and Prediction., 2nd ed.,Springer Series in Statistics, 2009 [10] V. Guardabasso, D. Rodbard, P. J. Munson A model-free approach to estimation of relative potency dose-response curve analysis, the American Physiological Society, Vol 252 Issue 3, 357-364, 1987 [11] H. Passing, W. Bablok A new biometrical procedure for testing the equality of measurements from two different analytical methods, J. Clin. Chem. Clin. Biochem. Vol. 21, pp. 709-720, 1983 [12] L. Bilic-Zulle, Comparison of methods: Passing and Bablok regression, Bioche- mia Medica 21(1):49-52, 2011 38