Parameterskattning av spatiala klusterprocesser med en inblick i nervdata Parameter estimation of spatial cluster processes with a view towards nerve data Examensarbete för kandidatexamen i matematik vid Göteborgs universitet Kandidatarbete inom civilingenjörsutbildningen vid Chalmers Hanna Ekelund Elsa Helle Lirim Kadriu Casper Opperud Hanna Skytt Samuel Thorén Institutionen för Matematiska vetenskaper CHALMERS TEKNISKA HÖGSKOLA GÖTEBORGS UNIVERSITET Göteborg, Sverige 2020 Parameterskattning av spatiala klusterprocesser med en in- blick i nervdata Examensarbete för kandidatexamen i matematik inom Matematikprogrammet vid Göteborgs universitet Elsa Helle Kandidatarbete i matematik inom civilingenjörsprogrammet Teknisk matematik vid Chalmers Hanna Ekelund Lirim Kadriu Casper Opperud Hanna Skytt Samuel Thorén Handledare: Aila Särkkä Institutionen för Matematiska vetenskaper CHALMERS TEKNISKA HÖGSKOLA GÖTEBORGS UNIVERSITET Göteborg, Sverige 2020 Förord Vi vill främst tacka vår handledare Aila Särkkä vars handledning och insikter har varit ovärderliga. Vi vill även tacka Esmée Berger, Christoffer Ekgrim, Julia Jansson och David Larsson för gott samarbete när det gäller feedback på varandras rapporter. Vi vill dessutom tacka varandra för gott samarbete. Gruppen har löpande fört loggbok där allas bidrag till projektet loggats och där alla problem som vi stött på under arbetets gång har dokumenterats. Nedan i tabellen specificeras vilka med- lemmar som tillskrivits som huvudförfattare till de olika avsnitten av projektet. Förutom det som syns som resultat i rapporten har det även spenderats många timmar åt litte- raturstudier, kodning i R, illustrerande av figurer, samt korrekturläsning och feedback till varandras text. Alla i gruppen har ägnat sig åt litteraturstudier för att få en förståelse för projektets teoretiska bakgrund. Samuel läste på lite extra om Ripleys funktioner samt Monte Carlo-simuleringnar, Elsa läste om parkorrelationsfunktionen och Hanna S läste om Minimum Contrast-metoden. Största delen av simuleringsstudierna samt övrig kod i R har skapats av Lirim och Casper med undantag för violindiagrammen i appendix D.3 som skapats av Samuel. Det är också de som genererat alla figurer i rapporten med undantag för figur 4, figur 5, 6 samt figur 20 som illustrerats av Elsa. Lirim har även deltagit aktivt i avsnitten teori och analys av nervdata samt engagerat sig i den generella utformningen av rapporten. Under de sista dagarna lästes rapporten genom noggrant av Hanna S och Elsa där språket var i fokus. Utöver arbete kring rapporten så har Hanna S och Hanna E även varit på flera föreläsningar om hur man skriver vetenskaplig text, de har använt sina kunskaper från dessa och delat med sig till övriga gruppmedlemmar. Avsnitt Skrivet av Övrig information Populärvetenskaplig presentation Hanna E Diskuterat med Elsa, Figur: Lirim Sammanfattning och abstract Lirim 1 Inledning Elsa 1.1 Hanna S 1.2 Casper 1.3 Samuel 2 Teori Lirim 2.1 Lirim, Samuel 2.1.1 Lirim 2.1.2 Elsa, Lirim Figur: Elsa 2.2 Samuel 2.2.1 Samuel, Elsa Figur: Elsa 2.2.2 Hanna S 2.3 Elsa 2.3.1 Elsa, Samuel Figur: Elsa 2.3.2 Hanna E 2.4 Casper 2.5 Hanna S 2.6 Samuel, Lirim 2.7 Lirim, Casper 2.7.1 Lirim 2.7.2 Casper 3 Simuleringsstudie Lirim 3.1 Lirim 3.2 Lirim 4 Analys av nervdata Lirim 4.1 Lirim 4.2 Lirim, Elsa 4.3 Lirim, Elsa 4.4 Lirim, Elsa 4.4.1 Elsa Figur och figurtext: Lirim 4.4.2 Lirim 4.5 Elsa Figur och figurtext: Lirim 5 Sammanfattning Samuel, Casper, Lirim Figur: Elsa A Härledningar A.1 Hanna S A.2 Samuel B Simulering B.1 Lirim, Elsa B.2 Lirim C Tabeller och resultat Lirim D Figurer D.1 Lirim, Samuel Text: Samuel D.2 Lirim, Samuel Text: Samuel D.3 Samuel E Kod E.1 Lirim E.2 Lirim E.3 Lirim E.4 Casper E.5 Lirim E.6 Lirim E.7 Casper E.8 Samuel E.9 Lirim Populärvetenskaplig presentation Spatiala punktprocesser beskriver hur punkter i ett punktmönster förhåller sig till varandra. För att kunna beskriva ett punktmönster matematiskt används modeller för punktprocesser. En mo- dell använder sig av två olika sorters punkter, föräldra- och dotterpunkter. Punktprocessmodellen placerar först ut föräldrapunkter i rummet och sedan ett antal dotterpunkter runt varje föräld- rapunkt. Både antalet föräldrapunkter och antalet dotterpunkter per föräldrapunkt beskrivs med funktioner som beskriver en fördelning, även dotterpunkternas position förhållande till föräldrarna beskrivs med dessa funktioner. Detta gör det möjligt att beskriva många olika punktmönster med en och samma modell. Av det följer att olika punktmönster kan se olika ut och därmed har även vissa mönster fått namn som beskriver dem. I figur 1(a) ser vi ett punktmönster som är helt slump- mässigt, det vill säga att dotterpunkterna är helt oberoende och likformigt fördelade i hela planet. Det är värt att notera att det endast är dotterpunkter som visas i figuren. Ett annat punktmönster är klustrat som vi kan se i figur 1(b). Punkterna i ett klustrat punktmönster är samlade i grupper. I vissa fall kan det vara svårt att avgöra om ett punktmönster är klustrat eller slumpmässigt med blotta ögat, men med statistik visas hur klustrat ett punktmönster är. (a) Slumpmässigt punktmönster. (b) Klustrat punktmönster. Figur 1: Olika punktmönster på ett begränsat område i planet (R2). Det är intressant att undersöka punktmönster då de kan modellera många fenomen i vår vardag. Till exempel går det att beskriva stjärnornas position i himlen eller trädens position i en skog. I detta arbete har vi använt punktprocesser för att beskriva positionen för nervändar i överhuden för friska individer och individer med sjukdomen diabetesneuropati. Diabetespatienter riskerar att utveckla en rad olika komplikationer som bland annat diabe- tesneuropati, vilket först ofta visar sig i fötter och lår i form av domningar och smärtor, men kan sprida sig till andra delar av kroppen och leda till nedsatt känsel. I nuläget upptäcks neuropati sent i förloppet, när symptom visat sig ett tag. Den sena upptäckten beror på att neuropati utvecklas långsamt vilket leder till att patienterna är sjuka långt innan de upptäcker symptomen. För att kunna behandla neuropati krävs att patienten diagnoseras så tidigt som möjligt. En skillnad som observerats i tidigare studier är att patienter med diabetesneuropati har färre nervtrådar i överhuden, samt att nerverna ej försvinner slumpmässigt. Nerverna försvinner så att resterande nerver formar kluster. Nervmönstrerna kan modelleras som en spatial punktprocess, där nerverna i överhuden representeras av dotterpunkter och deras rötter som föräldrapunkter. Detta punktmönster tenderar att vara klustrat för alla individer, oberoende av neuropati, men mönstren hos patienter med neuropati är mer klustrade än hos friska individer. Med hjälp av den så kallade Ripleys K-funktion och parkorrelationsfunktionen är det möjligt att mäta hur klustrat ett punktmönster är. Dessa funktioner kan användas för att anpassa en klustrad punktprocessmodell på given da- ta. Utifrån dessa kunde vi kvantifiera att punktmönstrena för både friska och sjuka individer är klustrade. Dessutom såg vi att patienterna med diabetesneuropati hade mer klustrade punktmöns- ter är de friska. Detta betyder att det är möjligt att använda denna metod för att diagnostisera diabetesneuropati, men metoden behöver testas på en större mängd data. Sammanfattning Ny teknik har lett till möjligheten att ta fram detaljerade bilder på nerverna i ytterhuden med hjälp av mikroskopi, vilket i sin tur avslöjat olika detaljer kring fördelningen av nervfibrer. Hos patienter med olika former av neuropati har det kunnat uttydas en större grad av klustring för nervtrådsändarna, jämfört med hos friska individer. Exempelvis är det intressant att se hur ändpunkterna av nervfibrernas spatiala mönster förändras i takt med att diabetesneuropatin avancerar. Eftersom mikroskopibilderna starkt tyder på att klustringen blir starkare ju mer utvecklad sjukdomen är, blir det viktigt att undersöka huruvida det är möjligt att i ett tidigare stadie diagnostisera patienter baserat på statistisk analys av ändpunktmönstren. Centralt i detta arbete var att modellera spridningen av nervfibrers ändpunkter med hjälp av så kallade thomasprocesser, som är modeller för klustrade punktmönster, och därefter se hur väl parametrarna kan skattas för modellen. Populära skattningsmetoder som maximum likelihood-metoden fungerar väldigt bra vid arbete med analytiska uttryck av likelihoodfunk- tionen. Vid spatiala punktprocessmodeller och klusterprocesser är det dock väldigt svårt eller rent av omöjligt att härleda ett analytiskt uttryck. Därför användes skattningsmetoden mi- nimum contrast-metoden, som är en minsta kvadrat-metod som bygger på sammanfattande statistikor. I detta arbete undersökte vi hur valet av den sammanfattande statistikan i minimum contrast-metoden påverkar parameterskattningarna av thomasprocessen. Vi valde att jämföra hur två statistikor, parkorrelationsfunktionen och Ripleys K-funktion, påverkar skattningarna för parametrarna i thomasprocessen. Från simuleringsstudien tog vi fram parametervärden för vilka minimum contrast-metoden kombinerat med en av de två sammanfattande statistikorna ger tillförlitliga skattningar. Vidare anpassades modellen efter den givna nervdatan som tagits fram genom biopsi av huden och vi undersökte möjligheten av att modellera punktmönstrens nervtrådsändar genom thomasprocessen. Från simuleringsstudien framgick det, om två av parametrarna fixeras och den sista va- rieras, att minimum contrast-algoritmen endast går att lita på för vissa värden av den vari- erande parametern. Dessa intervall skiljde sig åt beroende på om parkorrelationsfunktionen eller Ripleys K-funktion användes. Intervallet då Ripleys K-funktion användes var större än motsvarande för parkorrelationsfunktionen. Det kan konstateras från analysen av nervdatan att punktmönstret av nervfibrernas änd- punkter är klustrade för både friska och sjuka patienter och att thomasprocessen fungerar bra som modell för detta. De sjuka patienterna verkar ha färre basnerver samt mer täta och separerade kluster än de friska. Eftersom den givna datan är ett väldigt litet stickprov så är det svårt att ta fram intervall för parametrarna i thomasprocessen som kan hjälpa oss avgöra ifall ett givet punktmönster kommer från en frisk eller en sjuk patient. Nyckelord: diabetesneuropati, minimum contrast-metoden, Monte Carlo-simulering, parkor- relationsfunktionen, Ripleys K-funktion, spatiala punktmönster, thomasprocessen. Abstract Advances in image technology have given medical scientists the ability to take microscopic images of nerve fibers in the epidermis, which is the outermost part of the skin. This has subsequently unveiled some interesting characteristics of the spatial distribution of the nerve fibers. Patients suffering from different forms of neuropathy have exhibited a greater level of clustering of the nerve fiber end points as opposed to healthy individuals. Of particular interest is to study how the spatial pattern of the nerve fiber end points changes based on how advanced the diabetic neuropathy, in our case diabetic neuropathy, is. This leads to the question whether we can use spatial statistical analysis of the end point patterns in order to diagnose patients earlier. Of major importance in this project was to model the distribution of nerve fibers using the so-called Thomas process, which is a model for clustered point patterns. This is used to evaluate how well we can estimate the parameters of the model. The popular parameter estimation method maximum likelihood method works well when we have analytical expressions of the likelihood function in closed form available. However, regarding spatial point processes and clustered processes, analytical expressions are very rarely available. This means that we had to resort to other means of estimation such as the minimum contrast method, which is a least square estimation technique that uses different summary statistics as input. In this project we investigated how the choice of summary statistic in the minimum contrast method affects the estimations of the parameters of the Thomas process. We chose to compare the two summary statistics, the pair correlation function and Ripley’s K-function. From the simulation work we found a range of parameter values that can be estimated reliably by using the minimum contrast method combined with one of the two summary statistics. Then, the model was fitted to the given data that has been collected by skin biopsies and we investigated whether it is possible to model the nerve end points data using the Thomas process. The pair correlation function and Ripley’s K-function influence the accuracy of the esti- mations of the parameters in the Thomas process. The simulation study is essential in this project since it was used in order to conclude proper intervals of parameters which gives rea- sonably accurate estimates, after which the model is fitted to the data that has been generated from skin biopsies. Keeping two parameters fixed and only varying the third, the results of the simulation study indicated that we could only obtain reliable estimates of the parameters when the varied parameter was kept within a certain interval. This interval differed depending on whether we used the pair correlation function or Ripley’s K-function as summary statistic. The interval obtained using the K-function was larger than the corresponding interval for the pair correlation function. The analysis of the nerve end points data indicated that both point patterns for sick and healthy individuals alike were clustered and the Thomas process was indeed a good model for the data. However the sick patients had fewer base nerves and tighter and more separated clusters. Since the data is a very small sample, we could not calculate specific values of the parameters for which we could ascribe any point pattern to a healthy or sick individual. Keywords: diabetic neuropathy, minimum contrast method, Monte Carlo simulation, pair correlation function, Ripley’s K-function, spatial point processes, Thomas process. Förkortningslista Förkortning Betydelse CSR Complete spacial randomness MCM Minimum contrast-metoden MCS Monte Carlo-simulering PCF Parkorrelationsfunktionen Notationslista Notation Förklaring T (κ, µ, σ) Thomasprocessen med parametrar κ, µ, σ. κ Intensitet av föräldrapunkter för thomasprocessen. µ Genomsnittliga antalet dotterpunkter per kluster. σ Standardavvikelsen för avståndet av en punkt från sin föräldrapunkt. λ Förväntade antalet totala dotterpunkter i ett givet punktmönster. V[X] Varians av slumpvariabeln X. E[X] Väntevärde av slumpvariabeln X. S(t) Beteckningen för en godtycklig sammanfattande statistika. θ Godtycklig parametervektor, till exempel θ = (κ, µ, σ). D(θ) Funktionen som MCM bygger på. Minimerar felet mellan teoretisk och skattad parameter. K(r) Ripleys K-funktion, beror av radien r. L(r) Ripleys L-funktion, beror av radien r. g(r) Parkorrelationsfunktionen, beror av radien r. Kpoi(r) K-funktionen motsvarande den för CSR-punktmönster. Lpoi(r) L-funktionen motsvarande den för CSR-punktmönster. gpoi(r) PCF motsvarande den för CSR-punktmönster. Kiso(r) K-funktionen skattad från data. Liso(r) L-funktionen skattad från data. giso(r) PCF skattad från data. IK,σ Intervall över σ som ger tillförlitliga parametrar, skattat av Ripleys K-funktion. IP,σ Intervall över σ som ger tillförlitliga parametrar, skattat av PCF. Innehåll 1 Inledning 1 1.1 Syfte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Avgränsningar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Teori 2 2.1 Punktprocesser och spatial statistik . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.1.1 Statistiska fördelningar och Poissonprocessen . . . . . . . . . . . . . . . . . 3 2.1.2 Thomasprocessen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Ripleys K-funktion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2.1 Uppskattning av K-funktionen . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2.2 Teoretiska K-funktionen för thomasprocessen . . . . . . . . . . . . . . . . . 5 2.3 Parkorrelationsfunktionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.3.1 Uppskattning av parkorrelationsfunktionen . . . . . . . . . . . . . . . . . . 6 2.3.2 Teoretiska parkorrelationsfunktionen för thomasprocessen . . . . . . . . . . 7 2.4 Viktat medelvärde från data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.5 Minimum Contrast-metoden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.6 Monte Carlo-simulering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.7 Modellanpassning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.7.1 Ripleys K-funktion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.7.2 Parkorrelationsfunktionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3 Simuleringsstudie 11 3.1 Skattning av σ, κ och µ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.2 Sammanfattning av resultaten simuleringsstudien . . . . . . . . . . . . . . . . . . . 14 4 Analys av nervdata 14 4.1 Förundersökning av nervdata . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 4.2 Visualisering av nervdata . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4.3 Anpassning av thomasprocessen till nervdata . . . . . . . . . . . . . . . . . . . . . 16 4.4 Utvärdering av thomasprocessen som modell för nervdata . . . . . . . . . . . . . . 17 4.4.1 Grafisk utvärdering med hjälp av envelopes . . . . . . . . . . . . . . . . . . 17 4.4.2 Utvärdering via Monte Carlo-test . . . . . . . . . . . . . . . . . . . . . . . . 18 4.5 Sammanfattning av resultaten från analysen av nervdata . . . . . . . . . . . . . . . 19 5 Sammanfattning 19 Referenser 21 A Härledning av statistikor för thomasprocessen 22 A.1 Ripleys K-funktion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 A.2 Parkorrelationsfunktionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 B Simulering 25 B.1 Ändring av data och skalning av simuleringsfönster . . . . . . . . . . . . . . . . . . 25 B.2 Framtagning av lämpliga parametrar för simuleringsstudien . . . . . . . . . . . . . 26 C Tabeller och resultat 27 D Figurer 30 D.1 Punktmönster från nervdatan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 D.1.1 Grupp 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 D.1.2 Grupp 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 D.2 Grafisk undersökning för anpassningen av thomasprocessen . . . . . . . . . . . . . 33 D.2.1 Grupp 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 D.2.2 Grupp 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 D.3 Violindiagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 D.3.1 Ripleys K-funktion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 D.3.2 Parkorrelationsfunktionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 E Kod 39 E.1 Kod för histogram för skattningar då medelvärde används . . . . . . . . . . . . . . 39 E.2 Skillnaden i avvikande datapunkter, median kontra medelvärde . . . . . . . . . . . 40 E.3 Kod för parameterskattningar då medelvärde används . . . . . . . . . . . . . . . . 42 E.4 Kod för visualisering av statistikorna . . . . . . . . . . . . . . . . . . . . . . . . . . 44 E.5 Kod för parameterskattningar då median används . . . . . . . . . . . . . . . . . . . 47 E.6 Kod för K- och L-figurerna och deras envelopes . . . . . . . . . . . . . . . . . . . . 49 E.7 Kod för envelopes på nervdatan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 E.8 Kod för violindiagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 E.9 Kod för skattning/figurer av nervdata . . . . . . . . . . . . . . . . . . . . . . . . . 55 1 Inledning En komplikation för diabetespatienter är att de även kan utveckla nervsjukdomen diabetesneuro- pati. Sjukdomen bryter ofta ut i patientens fötter och kan ge upphov till smärtor, nedsatt känsel och i de värsta fall infekterade fotsår som kan leda till amputation. När sjukdomen väl har uppstått kan inte skadorna i nerverna gå tillbaka, det är därför av stor vikt att sjukdomen upptäcks i tid så att det finns en chans att förebygga skador i bästa mån [1]. Tidigare forskning har visat att personer med diabetesneuropati har färre nervändar i överhu- den. Dessutom verkar inte nerverna förtvina slumpmässigt, utan de försvinner så att resterande nervändar tenderar att bilda ett mer klustrat punktmönster än hos friska personer [16]. Ett punkt- mönster är data i form av punkter som är spridda över ett område och är ett resultat av en spatial punktprocess. Genom att analysera spridningen av punkterna är det möjligt att avgöra om punkt- mönstret är helt slumpmässigt, klustrat eller reguljärt. Vi behandlar så kallade klusterprocesser på grund av att spridningen av nervändar tenderar att vara mer klustrad än slumpmässig. Figur 2: Nervändarna i överhuden sett från sidan i ett tvärsnitt av huden. Sett uppifrån fås en bild av hur spridningen ser ut. Bild från Claes Andersson, 2017. (a) Slumpmässigt punktmönster. (b) Klustrat punktmönster. Figur 3: Punkterna representerar nervändarna i ytterhuden, sett uppifrån. För att skatta parametrar av modeller till data används traditionellt maximum likelihood- metoden, men för klusterprocesser är det inte möjligt att skriva ner likelihood-funktionen analy- tiskt och därför kommer vi istället använda oss av den så kallade minimum contrast-metoden. Det är en typ av minsta kvadrat-metod som har som avsikt att mäta avvikelsen mellan modell och data. Den gör detta genom att minimera skillnaden mellan en godtyckligt sammanfattad statistika som uppskattas från nervmönsterna och dess simulerade motsvarighet givet modellen [15]. Exempel på sammanfattande statistikor är parkorrelationsfunktionen och Ripleys K-funktion, som uppskattar antal punkter inom en given radie till en godtycklig punkt av processen. Genom att studera nerv- data i form av punktmönster, från såväl friska som diabetessjuka patienter som är i ett tidigt stadie av neuropati, är förhoppningen att kunna hitta en modell som beskriver dessa typer av klustrade punktmönster. Detta för att kunna diagnostisera diabetesnuropati innan sjukdomsförloppet har utvecklas avsevärt. 1 1.1 Syfte Syftet med detta projekt är att med hjälp av simuleringar undersöka hur parameterskattning för punktprocessmodeller påverkas av olika statistikor när minimum contrast-metoden används som skattningsmetod. Målen med projektet är att få en djupare förståelse för teorin om spatiala punkt- processer, undersöka klustrade punktmönster med hjälp av paketet spatstat i programspråket R samt kunna anpassa den modell som undersöks på given data. Projektet strävar efter att vara en del av huvudmålet, vilket är att kunna ligga ett steg före i processen av att upptäcka tidig diabetesneuropati genom att undersöka patienters nervändsfördelning. 1.2 Problem Projektet har för avsikt att i första hand utföra en simuleringsstudie på genererad data, detta genom att göra ett lämpligt val av modell och statistika i minimum contrast-metoden. I detta arbete kommer vi att använda oss av thomasprocessen som är en modell som lämpar sig väl när det rör sig om mönster som uppvisar högre grad av klustring jämfört med ett mönster som är helt slumpmässigt genererat. De sammanfattande statistikorna som ligger i fokus för denna studie är Ripleys K-funktion och parkorrelationsfunktionen. Målet med simuleringsstudien är att undersöka om valet av statistika påverkar parameterskatt- ningarna på den genererade datan. Utifrån detta kommer vi få en bild av vilken av statistikorna som ger bäst resultat för den genererade datan, som sedan kommer att användas när vi anpassar thomasprocessen på nervdatan. Förutsättningarna säger att både individer med och utan diabe- tesneuropati kommer att följa ett klustrat punktmönster. Däremot förväntas det att insjuknade individer kommer att ha ett mer klustrat punktmönster än de som inte är drabbade av diabe- tesneuropati. 1.3 Avgränsningar I detta projekt har vi valt att avgränsa valet till två olika statistikor. De metoder som kommer att undersökas är Ripleys K-funktion och parkorrelationsfunktionen. Anledningen till valet av dessa statistikor är att de ofta används på spatiala punktprocesser samt att båda de teoretiska funktionerna för thomasprocessen är kända. Dessutom är tidskravet för projektet ytterligare en bidragande faktor till varför vi har valt att avgränsa projektet till dessa två statistikor. Projektet kommer att använda sig av data från totalt 14 personer, vilket är ett ganska litet urval. Detta kan komma att påverka vilka slutsatser som kan dras och det är inte säkert om det går att generalisera resultatet i samma utsträckning som om stickproven varit fler. 2 Teori I detta avsnitt presenteras den matematiska teorin som ligger till grund för de metoder och sta- tistiska modeller som tillämpas. Thomasprocessen är en klusterprocess som bygger på ett antal fördelningar och processer som vi definierar nedan. Vi skattar parametrarna med hjälp av mi- nimum contrast-metoden (MCM) där vi kommer att använda oss av Ripleys K-funktion samt parkorrelationsfunktionen (PCF). En beskrivning av dessa samt redogörelse för hur de uppskat- tas hittas i detta avsnitt. Vi visualiserar sedan resultaten med hjälp av Monte Carlo-simuleringar (MCS). 2.1 Punktprocesser och spatial statistik Inom punktmönsteranalys studeras utfallen av punkter relativt med varandra samt deras position inom given yta. Eftersom vi i detta projekt enbart kommer att behandla två dimensioner så pratar vi om yta, men generellt sett är det möjligt att arbeta i rummet Rn, n = 1, 2, .... Ett vanligt exempel är att analysera hur en viss trädart växer inom ett slutet geografiskt område. Då marke- ras positionerna av träden på en karta och använder sig av modeller för att beskriva den spatiala fördelningen av träden. Om positionerna av träden är samlade i flera olika grupper, så kallade kluster, är trädens placering klustrad. I det fallet kan klusterprocesser användas, som är en viss typ 2 av punktprocesser, som modell för trädens positioner. Detta görs ofta genom att använda sig av så kallade föräldra- och dotterpunkter. I en klusterprocess genereras en del punkter som föräldrapunk- ter med ett antal dotterpunkter fördelade runt dem enligt en fördelning. Klusterprocesser består bara av dotterpunkterna och ofta finns ingen information om föräldrapunkterna. I exemplet med träden skulle föräldrapunkterna vara de ursprungliga träden och dotterpunkterna vara de träden som skapas från frön från föräldrapunkterna [9]. Spatiala punktmönster kan analyseras genom deskriptiv och inferentiell statistik. I deskriptiv statistik används icke-parametriska metoder för att belysa strukturen av det spatiala mönstret. I inferentiell statistik anpassas en parametrisk modell för datan och skattar dess parametrar från data. I detta projekt kommer vi att använda oss av deskriptiv statistik i form av sammanfattande statistikor och inferentiell statistik i form av anpassning av en klusterprocess. 2.1.1 Statistiska fördelningar och Poissonprocessen Den process som vi har valt att använda för simuleringar är thomasprocessen. För att närma- re kunna beskriva den måste vi presentera lite förkunskaper. Thomasprocessen använder sig av fördelningsfunktionerna normalfördelning och poissonfördelning. Normalfördelning är den vanligaste fördelningen som uppkommer naturligt i väldigt många sammanhang. Det är en kontinuerlig sannolikhetsfördelning för en reell stokastisk variabel. Definition 1 (Normalfördelning). För en kontinuerlig stokastisk variabel X ges dess täthetsfunk- tion av fX(x : µ, σ2) = 1 σ √ 2π exp ( −1 2 [ x− µ σ ]2 ) , x ∈ R. (1) Där parametrarna är väntevärdet E[X] = µ ∈ R och variansen V[X] = σ2 > 0. För att beskriva att en variabel X är normalfördelad med väntevärde µ och varians σ2 skriver vi X ∼ N (µ, σ2). Utfallen i en normalfördelning är symmetriskt fördelade kring medelvärdet och ju längre ifrån medelvärdet desto färre utfall förekommer. Poissonfördelningen är den vanligaste diskreta fördelningen. Denna fördelning används i pois- sonprocessen för att modellera antalet gånger en händelse har inträffat inom ett givet tidsintervall eller ett intervall i ett godtyckligt rum [13], [11], [10]. Vi definierar fördelningen enligt: Definition 2 (Poissonfördelning). En diskret stokastisk variabel X sägs vara poissonfördelad med parameter λ > 0, om för varje x = 0, 1, 2, ..., ges täthetsfunktionen av fX(x : λ) = P[X = x] = λxe−λ x! . (2) Väntevärdet och variansen av X är lika med λ, E[X] = V[X] = λ. För en stokastisk variabel som är poissonfördelad med parameter λ skriver vi X ∼ P(λ). Som nämnts tidigare används denna fördelning för att modellera händelser i tid och endimen- sionellt rum i poissonprocessen. Det är dock möjligt att generellt arebta i Rd, d ≥ 1. När d > 1 brukar man prata om den spatiala poissonprocessen. Definition 3 (Spatial poissonprocess). Låt A ⊂ R2 och låt NA vara antalet punkter i A. Beteckna ytan av A som |A|. För en poissonprocess med parameter λ, även kallat intensiteten, gäller att 1. NA har, för varje sluten mängd A, en poissonfördelning med parameter λ|A|, 2. om två mängder A och B är disjunkta, A ∩ B = ∅, så är NA och NB oberoende stokastiska variabler. Den likformiga fördelningen uppkommer naturligt i denna process eftersom punkterna får en likformig fördelning över given yta. Av denna anledning refereras denna process även som Complete Spatial Randomness (CSR). Vid simulering av denna process så genereras x- och y-koordinaterna likformigt och oberoende av varandra och antalet punkter genereras enligt poissonfördelningen. 3 2.1.2 Thomasprocessen Thomasprocessen är en klusterprocess som vi kommer att använda oss av i simuleringsstudien och i anpassning till data. Den bygger på alla ovanstående definitioner. Först genereras ett antal föräldrapunkter enligt den spatiala poissonprocessen som beskrivs ovan, med intensitet κ, ∼ P(κ). Kring dessa föräldrapunkter genereras, för varje förälder, dotterpunkter som i antal är poissonför- delade med µ, ∼ P(µ). Dotterpunkternas avstånd från sina respektive föräldrar är normalfördelade med väntevärde 0 och varians σ2, ∼ N (0, σ2), längs x- och y-axeln med föräldrapunkten i origo. Slutligen tas föräldrapunkterna bort så att endast dotterpunkterna kvarstår. Vi har valt att an- vända oss av thomasprocessen då den motsvarar nervdatan genom att låta basnerverna, se figur 3(b), betecknas av föräldrapunkter som inte observeras och deras förgreningar (ändpunkter) av dotterpunkter. Figur 4 visar visuellt hur denna process går till för att modellera ett punktmönster. För att uppskatta parametrarna av thomasprocessen krävs mer bakgrundsteori som presenteras i följande kapitel. Figur 4: Thomasprocessen i tre steg i område A. Modellering av föräldrapunkter (vänster), föräld- rapunkter med tillhörande dotterpunkter (mitten) och bara dotterpunkter (höger). En statistika är ett tal som beskriver stickprov, som medelvärde och varians, men det kan även beskriva karaktärsdrag hos ett stickprov som är betydligt mer komplicerade. Följande kapitel pre- senterar de statistikor vi undersöker och kommer att använda oss av för att beskriva punktmönster och skatta parametrarna av thomasprocessen. 2.2 Ripleys K-funktion Ripleys K-funktion är en fördelningsfunktion som används för att beskriva beteendet hos spatiala punktprocesser. Den beskriver, i det tvådimensionella fallet, antal punkter inom en radie, r, från en godtycklig punkt i punktprocess. Funktionen är definierad enligt: K(r) = λ−1E[N0(r)], (3) där N0(r) är antalet punkter hos processen inom ett avstånd r till en godtycklig punkt av processen och intensiteten λ är det genomsnittliga antalet punkter per enhetsarea [6]. Ripleys K-funktion ger information om hur klustrat ett spatialt punktmönster är. För ett punktmönster som uppfyller CSR, gäller det att K(r) = πr2,∀r > 0. Om värdet för K-funktionen är större än πr2 betyder det att punktmönstret är mer klustrat, och om det har ett lägre värde är punktmönstret mer utspritt än CSR [7]. En variant av K-funktionen är Ripleys L-funktion, vilket ger en fördelningsfunktion med sta- biliserad varians. Ripleys L-funktion definieras enligt L(r) = ( K(r) π )1/2 . (4) Precis som medK-funktionen går det att jämföra L-funktionen mot CSR. För ett punktmönster som uppfyller CSR är L(r) = r, ∀r > 0. Det är vanligt att rita upp en graf för L(r) − r. Om 4 L(r) − r > 0 uppvisar punktmönstret en högre grad av klustring och om L(r) − r < 0 så rör det sig om ett regelbundet mönster. Ju närmare värdet 0 desto bättre följer punktmönstret CSR. 2.2.1 Uppskattning av K-funktionen För att det ska vara möjligt att använda K-funktionen krävs det att vi definierar en lämplig uppskattning av funktionen. För att uppskattningen ska vara väntevärdesriktig krävs det även att vi introducerar en viktfunktion wij . Den är nödvändig då punkter som ligger nära områdets kant annars riskerar att bidra till att uppskattningen blir för låg, på grund av att inte hela området med radie r kring dessa punkter ligger i området A [7], [8]. Detta fenomen kallas kanteffekt och illustreras i figur 5. Då kan K-funktionen skattas av K̂(r) = λ̂−1 n∑ i=1 ∑ j 6=i w−1 ij I(rij ≤ r) n , (5) där λ̂ = n− 1 |A| (6) är en uppskattning av processens intensitet och n är antal punkter på området A. I ekvation (5) betecknar rij avståndet mellan punkt i och j. I är indikatorfunktionen, alltså är den lika med noll om rij > r och lika med ett om rij ≤ r [6]. Figur 5: Bilden ovan illustrerar fenomenet kanteffekt och hur det bidrar till felaktig uppskattning av K-funktionen. Viktfunktionen, wij , används i syfte att minska kanteffekten. Vi kommer använda Ripleys vikt- funktion som definieras genom att låta w(xi, r) beteckna andelen av omkretsen av cirkeln med centrum i xi och radie r som ligger i området A. Definiera sedan wij = w(wi, rij). Resultatet av detta blir att de punkter som ligger nära kanten av området vi undersöker kommer att viktas högre i summan av uppskattningen i ekvation (5). 2.2.2 Teoretiska K-funktionen för thomasprocessen I sektion 2.2 definierades Ripleys K-funktion enligt ekvation (3). Denna funktion kan explicit härledas för vissa processer som exempelvis thomasprocessen [12]. Den teoretiska K-funktionen för thomasprocessen med parametrar θ = (κ, µ, σ) blir K(r) = πr2 + κ−1 ( 1− exp ( − r2 4σ2 )) . (7) Härledningen hittas i appendix A.1. Det kan noteras från ekvation (7) att K-funktionen inte beror på µ. Därmed kan två olika thomasprocesser där det enda som skiljer är parametern µ fortfarande ha samma K-funktion. 5 2.3 Parkorrelationsfunktionen Parkorrelationsfunktionen, förkortad PCF, för ett punktmönster beskriver hur densiteten av res- terande utfall varierar i förhållande till avståndet från ett godtyckligt utfall i systemet. På så sätt beskriver PCF derivatan av Ripleys K-funktion. Funktionen definieras med hjälp av funktionen ρ(ξ)dξ, vilket är sannolikheten att vi kan hitta ett utfall i den oändligt lilla omgivningen med centrum i ξ, och funktionen ρ(2)(ξ, η)dξdη som är sannolikheten att hitta utfall i båda oändligt små omgivningar med centrum i ξ och η [14]. PCF ges sedan av g(ξ, η) = ρ(2)(ξ, η) ρ(ξ)ρ(η) . (8) Vi kan anta att vårt system är stationärt och isotropiskt, det vill säga att vi antar att det inte spelar någon roll var i systemet vi befinner oss eller hur punkterna ξ och η ligger i förhållande till varandra förutom avståndet, r = |ξ−η|, mellan dem. I det fallet blir g invariant under omvandling så att g(ξ, η) = g(r). För CSR är ρ(2)(ξ, η) = ρ(ξ)ρ(η) och därför är g(r) = 1. Om g(r) > 1 betyder det att det är mer sannolikt att hitta två punkter på avståndet r jämfört med CSR, det vill säga en homogen poissonprocess med samma intensitet. Vilket betyder att punktprocessen är klustrat, medan om g(r) < 1 betyder det att systemet är mer regelbundet än CSR. Notera att lim r→∞ g(r) = 1 oavsett om systemet är helt slumpmässigt, klustrat eller regelbundet. 2.3.1 Uppskattning av parkorrelationsfunktionen PCF uppskattas genom en algoritm som beräknar antal punkter som ligger inom ett avstånd från r och r+dr från en referenspunkt, vilket illustreras i Figur 6. Antag ett punktmönster medN punkter Figur 6: Bilden ovan illustrerar hur PCF uppskattas med hjälp av att införa avståndet dr och på så sätt räkna antal punkter som ligger på ett visst avstånd, r, från en godtycklig punkt. och area A, då är medeldensiteten ρ = N/A. Algoritmen för uppskattningen av g(r) innefattar upprepade beräkningar för de värden av r som är av intresse. För varje r beräknas medelantalet punkter som ligger inom området mellan r och r+ dr från en godtycklig referenspunkt i systemet. Därefter divideras det antalet med ytarean av området, ungefär 2πr dr. Slutligen divideras talet även med ρ, för att säkerställa att g(r) = 1 för punktmönster som inte har någon struktur. Detta beskrivs enligt g(r) = N(r + dr) Ω(r + dr) 1 ρ , (9) där N(r + dr) är antalet punkter som ligger innanför intervallet r + dr och Ω(r + dr) är arean av det området som undersöks [5]. 6 2.3.2 Teoretiska parkorrelationsfunktionen för thomasprocessen Tidigare i kapitlet definierades PCF enligt ekvation (8) och thomasprocessen i sektion 2.1.2. Den teoretiska PCF för thomasprocessen med parametrar θ = (κ, µ, σ) är g(r) = 1 + 1 4πκσ2 exp ( − r2 4σ2 ) . (10) En härledning för detta finns i appendix A.2. Som tidigare noterat så gäller att lim r→∞ g(r) = 1, (11) vilket följer direkt från ekvation (10). Notera även att PCF ej beror på µ. 2.4 Viktat medelvärde från data Anta att vi har data som består av m punktmönster. Om de underliggande punktprocesserna antas ha samma K-funktion, kan vi uppskatta den gemensamma K-funktionen genom att ta ett medelvärde av K-funktionerna som skattats från de m punktmönstrena. Det kan vara fördelaktigt att ta hänsyn till hur många punkter varje punktmönster har, för att få en bättre bild av hur processen beter sig. Vi kan då beräkna ett viktat medelvärde av K-funktionen, enligt ekvation (12). Här är ni antalet punkter för punktmönstret i, och Ki är K-funktionen till process i = 1, ...,m. K̂(r) = m∑ i=1 niK̂i(r)/ m∑ i=1 ni. (12) Vi behöver inte anta att intensiteten till de underliggande processerna är samma på grund av attK-funktionen inte påverkas punkter tas bort slumpmässigt. Vi kommer använda oss av detta för att kunna jämföra de två grupperna i nervdatan, genom att använda funktionen pool i R. Framöver kommer vi därför referera till de gemensamma K-funktionerna som poolade K-funktioner. 2.5 Minimum Contrast-metoden Maximum likelihood-metoden (MLM) är en parameterskattningsmetod som ofta används tack vare dess bra teoretiska egenskaper för att beskriva data under en viss sannolikhetsfördelning. Då maximeras likelihood-funktionen, L(θ), som fås från täthetsfunktionen eller sannolikhetsfunktionen med avseende på parametrarna θ. Men för klusterprocesser så är det sällan möjligt att skriva ner likelihood-funktionen analytiskt, varför det istället är möjligt att använda minimum contrast- metoden, förkortad MCM. MCM är en minsta kvadrat-metod som skattar parametrar som minimerar avvikelsen mellan det teoretiska värdet (eller simuleringsbaserade värdet), S(t; θ), och det uppskattade värdet från datan, Ŝ(t), för den givna sammanfattande statistikan. Parametrarna uppskattas genom att minimera integralen D(θ) = ∫ t0 ν w(t) ( Ŝ(t)c − S(t; θ)c )p dt, (13) med avseende på θ. Integralen beror på flera parametrar; ν, t0, w(t), c och p. R-paketet spatstat har en metod för MCM, mincontrast, som används för simuleringsstudien, där c = 0.25 och p = 2 är de förvalda värdena [3], och om inte annat specificeras är t0 = 0.25, ν = 0. Detta är lämpligt eftersom både vid Ripleys K-funktion och PCF så rör det sig om radiella avstånd. Parametervärdet p = 2 ges av likheten med minsta kvadrat-metoden. Enligt Diggle så är w(t) = 1 ett bra val till klustrade punktmönster [6]. Diggle rekommenderar även t0 = 0.25 min(a, b) där a och b är sidlängderna för området. I ett tidigare arbete undersöktes hur valet av integrationsgränserna påverkade parameterskattningarna [12]. Då det redan har gjorts en gång, undersöks det inte igen, utan de förvalda värdena för t0 och ν används. För vissa statistikor såsom Ripleys K-funktion och PCF är det möjligt att explicit definiera det teoretiska värdet. Exempelvis sätter vi S(t) i MCM som Ripleys K-funktion med de valda parametrarna och erhåller följande integral: 7 D(θ) = ∫ 0.25 0 ( K̂(t)0.25 −K(t; θ)0.25 )2 dt. (14) Analogt görs för PCF. 2.6 Monte Carlo-simulering Monte Carlo-metoden är en simuleringsmetod som används för att göra numeriska uppskattningar då andra metoder blir för komplexa. Den grundas på ett stort antal simuleringar, så kallade Monte Carlo-simuleringar (MCS), som sker genom att slumpmässigt generera en stor mängd variabler från en sannolikhetsfördelning över ett område. Trots att varje enskild simulering är slumpmässig så är idén att genomsnittet ska representera den verkliga fördelningen. Monte-Carlo metoden är möjlig att använda för att uppskatta integralen I = ∫ Ω f(x)µ̂(dx), (15) med summan Î = 1 m m∑ i=1 f(xi). (16) Parametern m är antalet genererade slumpvariabler, xi, med fördelningen µ̂. Området Ω är det som variablerna slumpas över och f är en godtycklig funktion som är intressant i sammanhanget. Med hjälp av MCS är det möjligt att bland annat göra hypotesprövningar och uppskatta konfi- densintervall [4]. I vårt fall kommer vi småningom göra ett p-test för anpassningen av thomasprocessen. Det enk- laste sättet att ta fram ett p-värde är att beräkna skattningen av den sammanfattande statistikan, Ŝ(t), för observerat ändpunktmönster och därefter skattningarna Ŝi(t), i = 1, ...,m från m antal simuleringar av den anpassade thomasprocessen med parametrar θ̂ = (κ̂, µ̂, σ̂). Dessa parametrar erhålls genom att använda Ŝ(t) och MCM, se ekvation (13). Dessa Ŝi(t) jämförs mot det analytiska uttrycken för Si(t : κ̂, µ̂) som i vårt fall antingen kommer vara det analytiska uttrycket för Ripleys K-funktion eller PCF för thomasprocessen vilka ges i avsnitt 2.2.2 och 2.3.2. Notera att de analy- tiska uttrycken Si(t : κ̂, µ̂) innehåller de olika skattningarna av κ, µ som erhålls av ekvation (13), därav indexeras dem med i nedan i ekvation (17). Statistikan skattas för alla t ∈ [0,min(a, b)]. Vi inför en ny teststatistika Gi som vi definierar som summan av den kvadratiska skillnaden mellan skattningarna Ŝi(t) och den teoretiska statistikan Si(t), alltså Gi = ∑ t∈[0,min(a,b)] [ Ŝi(t)− Si(t : κ̂i, µ̂i) ]2 . (17) Notera att dessa Gi enbart kommer från simuleringarna, men vi behöver även en referens av G-statistikan att jämföra mot. Denna referens är Gobs som erhålls genom att sätta in Ŝ(t) och S(t : κ̂, µ̂) i summanden i ekvation (17) enligt Gobs = ∑ t∈[0,min(a,b)] [Ŝ(t)− S(t : κ̂, µ̂)]. (18) För att beräkna p-värdet för varje patient räknar vi hur många simulerade statistikor Gi är större än eller lika med Gobs och beräknar andelen enligt p = 1 m m∑ i=1 I(Gi ≥ Gobs). (19) 2.7 Modellanpassning Som nämnts tidigare kan det vara komplicerat, om inte omöjligt, att ta fram slutna analytiska uttryck. Som följd av detta tillgripes ofta simulering av nollhypotesen som ett sätt att skapa en 8 referensfördelning. Analysen av punktmönster börjar ofta genom att testa nollhypotesen H0, då punktspridningen representeras av en CSR-modell. I detta avsnitt visar vi hur statistikorna används för att skapa överskådliga grafer, som enkelt kan läsas av för att få en uppskattning av de olika egenskaperna hos punktmönster. 2.7.1 Ripleys K-funktion Vid första anblick är det lockande att bara visuellt titta på ett punktmönster för att avgöra hur klustrat mönstret är. Det är dock svårt att bedöma egenskaper hos punktmönster genom att endast observera det. Därför använder vi oss av Kest i spatstat för att skatta K-funktionen för givna punktmönster. Vi jämför K-funktionen Kiso som är skattad från data med K-funktionen Kpoi som är skattad från CSR-mönster med samma intensitet som datan. Observera att Kiso motsvarar K(r) i integranden i ekvation (14). För de r som ger Kiso(r) > Kpoi(r) kan vi förvänta oss fler punkter på det avståndet än vid CSR med samma intensitet, dvs. det finns högre tendens för ett klustrat mönster. För de r som medför Kiso(r) < Kpoi(r) förväntas ett mer regelbundet punktmönster. Figur 7 visar hur Kiso(r) ochKpoi(r) förhåller sig till varandra beroende på om vi tittar på ett ett helt slumpmässigt mönster så som figur 3(a) eller ett mer klustrat mönster, som figur 3(b). (a) Slumpmässigt mönster. (b) Klustrat mönster. Figur 7: Uppskattade K-funktionen Kiso(r) jämfört med den teoretiska K-funktionen under CSR Kpoi(r). För att få en tydligare bild av situationen för små r använder vi istället L-funktionen. Detta beräknas i R via funktionen Lest i spatstat. Se resultaten av detta i figur 8. (a) Slumpmässigt punktmönster. (b) Klustrat punktmönster. Figur 8: Uppskattade L-funktionen Liso(r) jämfört med den teoretiska L-funktionen under CSR Lpoi(r). För att få en ännu tydligare bild gör vi ytterligare en modifikation av K-funktionen och ritar L(r)−r istället. Då är Lpoi(r)−r alltid lika med 0. Figur 9 nedan visar resultaten. Notera skillnaden 9 mellan dessa funktioner beroende på om vi använder ett klustrat mönster eller ett CSR-mönster. För det slumpmässiga mönstret är skillnaden mellan Lpoi(r) och Liso(r) betydligt mindre längs hela intervallet för r medan denna skillnad ökar väldigt snabbt för det klustrade mönstret. (a) Slumpmässigt punktmönster. (b) Klustrat punktmönster. Figur 9: Uppskattad Liso(r) − r funktion jämfört med den teoretiska Lpoi(r) − r funktion under CSR. Problemet med att enbart jämföra Liso(r) med Lpoi(r) är att det inte ger oss någon information om variationen. Genom att simulera ett stort antal CSR-mönster och använda envelope i R är det möjligt att konstruera ett band baserat på de K-funktionerna uppskattade från punktmönstrena. Bandet är räknat så att för varje r lämnas 2.5% av de lägsta och 2.5% av de högsta värdena utanför bandet. Figur 10(a) visar att L(r)−r det slumpmässiga mönstret ligger inom det gråa bandet vilket indikerar att det inte finns något som talar emot att detta faktiskt är ett CSR-mönster. (a) Slumpmässigt punktmönster. (b) Klustrat punktmönster. Figur 10: 95%-band på respektive L(r)− r grafer. I figur 10(b) avviker Liso(r)− r mycket ifrån Lpoi(r)− r vilket medför att den hamnar utanför det grå CSR-bandet. Det är därmed inte möjligt att med stor säkerhet säga att det rör sig om en CSR-process och vi måste förkasta nollhypotesen vilket betyder att punktmönstret som används i figur 10(b) verkar klart vara mer klustrat än CSR. 2.7.2 Parkorrelationsfunktionen Som nämns i 2.3 så beskriver PCF densiteten av förväntat utfall för ett kluster i ett punktmönster. Vi kommer använda oss av 3(a) och 3(b) som referensmönster. Även om det går att avgöra visuellt vilken typ av punktmönster som de följer, är detta inte alltid fallet. Vi använder oss av PCF i spatstat för att uppskatta g(r), r ∈ [0, 0.25], för samma punktmönster som i figurerna ovan. Vi jämför sedan giso(r) med gpoi(r). Här är giso(r) PCF som uppskattas från datapunkterna med Ripleys viktfunktion, se 2.2.1, och gpoi(r) betecknar PCF för CSR. För de r där giso(r) > gpoi(r) så observeras det flera par av punkter vid detta avstånd jämfört med CSR, alltså 10 är mönstret mer klustrat. När giso(r) < gpoi(r) är mönstret mer regelbundet än CSR. Figur 11 visar hur giso(r) och gpoi(r) ser ut för ett helt slumpmässigt mönster och ett klustrat mönster. För CSR-mönstret följer det observerade PCF, giso(r), gpoi(r). För det klustrade mönstret observeras att flera punkter ligger i närheten av föräldrapunkterna för intervallet r ∈ [0.00, 0.05]. Sedan avtar antal punkter i förhållandet till föräldrapunkten och börjar följa gpoi(r). Detta tyder på att enbart det slumpmässiga mönstret följer CSR, som det bör göra. (a) Slumpmässigt mönster. (b) Klustrat mönster. Figur 11: Uppskattad g-funktion giso(r) jämfört med den teoretiska g-funktionen under CSR gpoi(r). Vi använder oss nu av MCS och funktionen envelope i R. Vi simulerar 99 stycken CSR punkt- mönster, där PCF uppskattas. Detta visualiseras i figur 12(a) och 12(b). För det slumpmässiga mönstret kan vi konstatera att det verkar ha genererats från CSR, så nollhypotsen kan inte för- kastas. Detta representeras med att de observerade värdena giso(r) ligger mellan ghi(r) och glo(r). För det klustrade mönstret ligger giso(r) utanför bandet, så vi kan förkasta nollhypotesen och dra slutsatsen att det inte finns skäl till att anta att det rör sig om en CSR-process. (a) Slumpmässigt punktmönster. (b) Klustrat punktmönster. Figur 12: 95%-band på respektive g(r)-grafer. 3 Simuleringsstudie I simuleringsstudien undersöker vi hur de olika statistikorna påverkar skattningen av parametrar- na i thomasprocessen. Vi simulerar realisationer av thomasprocessen med olika parametrar som vi har valt utifrån observationer från förundersökningen av nervdata, se appendix B.2. Eftersom thomasprocessen karaktäriseras av parametrarna κ, µ och σ, betecknas processen med T (κ, µ, σ). Parametern κ motsvarar intensiteten av föräldrapunkter, µ är intensiteten av dotterpunkterna och σ är standardavvikelsen av avståndet av en punkt (längs varje koordinataxel) från sitt klustercent- rum. Den totala intensiteten för hela processen är λ = κµ som är antalet förväntade punkter. 11 3.1 Skattning av σ, κ och µ Vi vill få en inblick i hur väl vi kan skatta parametrarna med olika parametervärden och de sammanfattande statistikor som vi arbetar med. Vi har valt att simulera med fasta värden för κ och µ, nämligen κ = 22.9 och µ = 4, och låta σ variera. Värden på κ och µ har valts så att de motsvarar värdena i nervdatan. Hur vi kommit fram till dessa värden förklaras i appendix B.2 där vi gör en förundersökning av nervdatan. Alla simuleringar är gjorda i enhetskvadraten [0, 1]× [0, 1]. Vi börjar med att titta på hur punktmönstrena för tre olika värden på σ, se figur 13. (a) Punktprocess nära CSR. (b) Uppenbar klustring. (c) Extrem klustring. Figur 13: Punktmönster med varierande σ. När σ är relativt stort ser punktmönstret ut att vara slumpmässigt och minskande σ ökar klustringen. Skulle vi välja för stora värden på standardavvikelsen så närmar sig thomasprocessen CSR som är ointressant, se figur 13(a). Ett för litet värde på standardavvikelsen ger å andra sidan ett väldigt tätt klustrat mönster som kan ses i figur 13(c). Vi bör hålla oss till σ > 0.01 för att inte klustringen ska bli så extrem och icke-representativ av nervdatan. För att undersöka ett lite bredare spektrum av värden har vi valt att undersöka alla σ ∈ [0.01, 0.15]. De funktioner i R som användes var rThomas, thomas.estK och thomas.estpcf. Den funktion som genererar slumpmässiga realisationer av thomasprocessen är rThomas, som genererar algorit- men för föräldrapunkter och dotterpunkter. Funktionerna thomas.estK och thomas.estpcf skat- tar parametrarna κ och µ givet ett thomas-punktmönster genom att använda Ripleys K-funktion respektive PCF i MCM. Se detaljerna för hur skattningarna beräknas i avsnitt 2.2.1 och 2.3.1. Enligt ekvationerna (7) och (8) är statistikorna obereoende av µ, och därför så kommer inte funk- tionerna att skatta denna parameter. Detta gör vi istället genom att från punktmönstret ta fram totala antalet punkter och dividera dessa med κ̂, ty vi vet att κ · µ = λ är förväntade antalet punkter. För att gå vidare med frågeställningen måste vi se hur väl vi kan skatta parametrarna från den simulerade thomasprocessen med varierande σ ∈ [0.01, 0.15]. Detta görs enklast genom att partitio- nera intervallet [0.01, 0.15] i några mindre delintervall med likformig steglängd h = 0.15−0.01 14 = 0.01 och loopa igenom med rThomas och thomas.estK respektive thomas.estpcf för dessa värden. För varje värde på σ kördes 1000 simuleringar av thomasprocessen och lika många skattningar av varje parameter. Dessa parametrar läggs in i 3 olika 1× 1000-vektorer och medelvärdet av alla värden i respektive vektor beräknades för att slutligen ta fram en rimlig skattning. Detta gjordes en gång med Ripleys K-funktion och en gång med PCF, tabell C.7 och C.8 i appendix visar resultaten. Tabellerna visar en stor varians hos skattningarna och att generella trender saknas. Vi kan därför inte säga så mycket om intervallet Iσ där skattningsmetoden fungerar bra. För att förstå vad som har gått snett tittar vi på de empiriska fördelningarna av parameterskattningarna. Vi kan visualisera dessa avvikande värden genom att använda så kallade violindiagram, se figur 14 och appendix D.3. Observera att y-axeln är logaritmerad för de violindiagrammen i appendix. Vi ser att för både PCF och Ripleys K-funktion har alla empiriska fördelningar för κ och µ en stor spridning efter σ ≈ 0.05. Detta leder självklart till att även σ får en stor spridning eftersom σ direkt beror av κ och µ. Detta leder till att beräkningar av medelvärde påverkas väldigt mycket 12 Figur 14: Ett violindiagram av empiriska fördelningar av κ̂ då K-funktionen används. När σ är lågt koncentreras alla tusen utfallen från simuleringen kring ett värde, i detta fall kring 23. I takt med att σ ökar uppkommer fler avvikande värden upp till två tusen. ”Kroppen” på violinen innehåller de flesta utfall och de långa strecken är avvikande värden. och därmed ger missvisande resultat. Eftersom de flesta värden för µ och κ hamnar i intervallen [0, 10] respektive [0, 40] så måste vi ändra på hur vi får gemensamma skattningen från simuleringarna för att eliminera påverkan från avvikelserna. Det enklaste sättet är att använda medianvärdet istället för medelvärdet. Figurerna 15 och 16 visar hur värden på skattningarna beror på σ. För både K-funktionen och PCF erhålls betydligt bättre skattningar baserat på simuleringarna när median används. Notera skalan på y- axlarna för medelvärde kontra de för median. (a) K-funktion och medelvärde. (b) K-funktion och median. Figur 15: Skattningar av de tre parametrar vid användning av Ripleys K-funktion och medelvärde (vänster) respektive median (höger) av de 1000 skattningarna. Notera att skalan på y-axeln i den vänstra figuren skiljer sig åt från skalan i den högra. 13 (a) PCF och medelvärde. (b) PCF och median. Figur 16: Skattningar av de tre parametrar vid användning av PCF och medelvärde (vänster) respek- tive median (höger) av de 1000 skattningarna. Notera att skalan på y-axeln i den vänstra figuren skiljer sig åt från skalan i den högra. Från tabell C.7 och C.8 i appendix och figurerna 15(b) och 16(b) ser vi att alla skattningarna κ̂, µ̂, σ̂, för K-funktionen då median används, är väldigt stabila kring det verkliga värdet fram till och med σ = 0.09. När PCF används försämras värdena för κ̂ betydlig fortare, redan vid σ = 0.03. 3.2 Sammanfattning av resultaten simuleringsstudien Som slutsats för simuleringsstudien kan vi säga att ett bra intervall för σ då vi använder oss av K-funktionen är Iσ,K = [0.02, 0.09], ty vi tidigare infört restriktionen att σ > 0.01 för att begränsa klustringen. Då vi använder PCF fås det tillförlitliga intervallet till Iσ,P = [0.02, 0.03]. Spontant ser det ut som att K-funktionen fungerar bättre eftersom det erhålls högre stabilitet på skattningarna över ett större intervall för σ. Notera att dessa intervall baseras på att simuleringarna av thomasprocessen genererats av (κ, µ) = (22.9, 4). Hade dessa värden varit annorlunda hade även intervallen Iσ,K och Iσ,P skiljt sig från erhållna resultat. Detta säger oss att MCM-algoritmen i R konvergerar endast för σ som leder till att punktmönstret är klustrat. Ju längre ifrån CSR, desto bättre skattningar. 4 Analys av nervdata Datan är tillgänglig som en rds-fil, endpoints.rds, visualisering av datan finns i appendix D.1. Nervdatan och övriga punktmönster lagras som ett ppp-objekt i R och står för planar point pattern. Struktureringen av ett sådant objekt innehåller alla x- och y-koordinater för varje dotterpunkt, antal punkter n i varje mönster, samt fönsterstorleken i vilken dessa punkter ska genereras. Dessa punktmönster har samlats genom att undersöka prover av hudbiopsin, vilket är små fragment av hud som skurits ut från patienter. Datan består av två grupper, en grupp med 8 friska patienter, kallad grupp 1, samt grupp 2 som innehåller data från 6 patienter som lider av diabetesneuropati. Vi kommer i detta kapitel börja med att analysera den givna nervdatan för att se hur mycket punktmönstren skiljer sig från en homogen, slumpmässig poissonprocess, samt i vilken riktning; är de mer klustrade än CSR eller mindre? 4.1 Förundersökning av nervdata Vi vill undersöka hur punktmönstrerna för de två grupperna av data ser ut i jämförelse med CSR och vill därför plotta liknande figurer med envelopes som vi gjort i tidigare avsnitt. För detta måste vi veta intensiteten av ändpunktsprocessen i de två grupperna. Intensiteten av processen ges av λ = κµ, vilket ger det förväntade antalet punkter per enhetsyta. Vi använder enhetskvadrat i simuleringarna och därför kan λ uppskattas med antal punkter i 14 mönstret delat med fönstrets storlek. Vi sammanställer medelantalet punkter i alla grupper. Antalet punkter för patient j i grupp i betecknas av nij , där i ∈ {1, 2} och j ∈ {1, 2, 3, 4, 5, 6, 7, 8}. Grupp 1 Patient 1 n11 111 Patient 2 n12 144 Patient 3 n13 154 Patient 4 n14 143 Patient 5 n15 113 Patient 6 n16 175 Patient 7 n17 157 Patient 8 n18 72 Medelantal n̄1. 142.4 Tabell 1: Antal punkter i grupp 1, rödmarkerade exkluderade i beräkningen av medelantal. Grupp 2 Patient 1 n21 153 Patient 2 n22 113 Patient 3 n23 83 Patient 4 n24 94 Patient 5 n25 30 Patient 6 n26 101 Medelantal n̄2. 108.8 Tabell 2: Antal punkter i grupp 2, rödmarkerade exkluderade i beräkningen av medelantal. För att inte få missvisande medelvärden behöver vi ta bort avvikelser. Eftersom vi har ett väldigt litet dataset med endast 14 patienter så får vi vara försiktiga med att eliminera allt för många avvikande värden. Vi nöjer oss med att ta bort det mest extrema fallen i respektive grupp, patient 8 i grupp 1 och patient 5 i grupp 2. Vi har därmed inte tagit hänsyn till dessa avvikelser i beräkningarna av medelvärdena n̄1. och n̄2.. Nervdatan är skalad, se appendix B.1, och de skalade fönsterdimensionerna är [0, 1.31]× [0, 1]. Detta betyder att den uppskattade intensiteten för de två grupperna är λ1 = n̄1. 1.31 respektive λ2 = n̄2. 1.31 . 4.2 Visualisering av nervdata För att se hur varje enskilt punktmönster ser ut hänvisar vi till appendix D.1. Notera att datan där är skalad, enligt beskrivning i appendix B.1. Det kan ges en bättre uppfattning om hur klus- (a) Grupp 1. (b) Grupp 2. Figur 17: Envelopes för PCF för grupp 1 (vänster) och grupp 2 (höger). Notera att skalan på y-axeln skiljer sig åt avsevärt mellan grupp 1 och grupp 2. tertrenderna ser ut för de olika grupperna genom att generera och studera figurer som visar den uppskattade statistikan, antingen Ripleys K-funktion eller PCF, för ett punktmönster tillsammans med envelopes för CSR, precis som vi redogjort för i avsnitt 2.7. Målet är att kunna jämföra de två grupperna och få en tydligare bild av hur de förhåller sig till CSR. Vi visar figurer för de ge- nomsnittliga poolade uppskattningarna av PCF för den samlade datan i grupp 1 respektive grupp 2, se figur 17 och D.2 i appendix. 15 Kom ihåg att vad PCF visar är antalet par av punkter som befinner sig på radie r från varandra, i förhållande till vad som förväntas för CSR. Det vi kan se i figurerena är att det verkar som nervmönsterna för båda grupperna är generellt mer klustrade än CSR. Framförallt kan vi se att för grupp 1 finns det en klustertrend i intervallet r ∈ [0.01, 0.07] från en godtycklig punkt, medan för grupp 2 verkar klustertrenden vara något snävare inom intervall r ∈ [0.00, 0.05]. Dessutom är avvikelsen från CSR större i grupp 2 än grupp 1. Tolkningen av detta blir att vi kan förvänta oss ungefär lika stora kluster, det vill säga antal dotterpunkter, för de två grupperna, men att grupp 2 genererar något tätare och mer separerade kluster än grupp 1. 4.3 Anpassning av thomasprocessen till nervdata Förundersökningen visar att nervpunktmönstrerna är klustrade och därför kan thomasprocessen vara en bra modell för ändpunkterna. Vi skattar parametrarna av thomasprocessen genom att an- vända både Ripleys K-funktion och PCF för varje patient i varje grupp, resultatet är presenterat i tabellerna 3 och 4. Notera återigen att nervdatan är skalad, vilket betyder att dessa skattningar för σ och κ inte är detsamma som för datans originalstorlek, utan betydligt mindre. Se appen- dix B.1 för mer information om skalningen. Notera även att värdena i tabellen nedan är utifrån enhetsfönstret [0, 1] × [0, 1]. Vid första anblick av resultaten lägger vi märke till att den största Grupp 1 Ripleys PCF Patient µ̂ κ̂ σ̂ µ̂ κ̂ σ̂ 1 3.65061 30.40585 0.03974 2.72907 40.67326 0.02851 2 1.59493 90.28592 0.02061 1.71094 84.16438 0.01831 3 3.48571 44.18040 0.03328 2.49982 61.60454 0.02325 4 2.27386 62.88867 0.02209 1.79656 79.59639 0.01660 5 6.61844 17.07352 0.04996 6.58899 17.14983 0.04500 6 8.53601 20.50139 0.05050 5.15993 33.91520 0.03015 7 4.70917 33.33922 0.03523 2.54671 61.64821 0.01868 8 3.20799 22.44392 0.03298 3.17679 22.66436 0.03060 Medel 4.25959 40.13986 0.03555 3.27612 50.17702 0.02639 Tabell 3: Parameterskattningar för grupp 1. Grupp 2 Ripleys PCF Patient µ̂ κ̂ σ̂ µ̂ κ̂ σ̂ 1 6.81699 15.26133 0.03895 6.75069 21.73420 0.02686 2 5.03477 63.21135 0.01496 4.98580 76.92758 0.01166 3 3.69811 20.46026 0.02143 3.66214 19.22366 0.02020 4 4.18822 20.44933 0.02899 4.14748 20.14053 0.02675 5 1.33667 18.83654 0.02727 1.32366 15.97330 0.02628 6 4.50011 12.41746 0.03316 4.45634 17.23940 0.02292 Medel 4.26248 22.44392 0.03298 4.22102 22.66436 0.03060 Tabell 4: Parameterskattningar för grupp 2. skillnaden mellan grupperna verkar vara gällande parametern κ. Minns att κ är intensiteten av föräldrapunktprocessen, vilket i vår data motsvarar antal basnerver i enhetsfönstret. Vi kan se att antalet basnerver verkar var större i grupp 1, friska patienter, än i grupp 2, sjuka patienter, vilket troligen är grund av att för patienter drabbade av diabetesneuropati så dör en del av nerverna. Hela 16 nervträdet, dvs. föräldrapunkter och dotterpunkter, försvinner. Vi ser att parameterskattningarna för µ, dvs. det förväntade antalet dotterpunkter per kluster, inte verkar skilja sig anmärkningsvärt mellan de olika skattnigsmetoderna eller grupperna i datan. Detta är förenligt med den analys som gjordes i avsnittet ovan med stöd av figur 17. Slutligen, om vi vänder blickarna mot paramterskattningarna av parametern σ så verkar det även här inte vara någon avsevärd stor skillnad mellan grupperna. Dessutom fastställdes i avsnitt 3.1 att det pålitliga intervallet gällande parameterskattningar för PCF var Iσ,P = [0.02, 0.03] och då våra skattningar ligger precis på den övre gränsen kanske vi inte ska förlita oss helt på dem utan fokusera mer på skattningarna gjorda med Ripleys K-funktion. Där skiljer sig skattningarna av σ mellan grupp 1 och 2 minimalt. Skulle den skillnaden vara signifikant så stöder resultatet hypotesen om att grupp 2 har tätare kluster än grupp 1. 4.4 Utvärdering av thomasprocessen som modell för nervdata I detta avsnitt ska vi undersöka om thomasprocessen är en bra modell för vår data. Vi har valt att göra detta på två sätt. Den första metoden bygger på att göra en grafisk jämförelse genom upp- repade simuleringar av thomasprocessen där vi använder de skattade parametrarna för nervdatan enligt tabellerna 3 och 4. Den andra metoden bygger på konstruktion av Monte Carlo test baserat på PCF [4]. Den- na metod är bra för att kunna kvantifiera hur bra thomasprocessen anpassar punktmönstrena i nervdatan. Vi kommer i detta test att basera teststatistikan på PCF och bifoga de för Ripleys K- funktion i appendix. Metoden är i grund och botten ett hypotestest där vi beräknar p-värdet för att avgöra om nollhypotesen, som i detta fall är att ändpunktmönstret kommer från thomasprocessen med skattade parametervärden, bör förkastas. 4.4.1 Grafisk utvärdering med hjälp av envelopes Utvärderingen görs genom att utföra 9999 Monte Carlo simuleringar och uppskatta en teststatistika för varje patients skattade parametrar, enligt tabeller 3 och 4. Sedan skapar vi envelopes som representerar ett intervall för vart vi förväntar oss att kurvan ska ligga för ett punktmönster med nervdatans skattade parametrar. Till sist skapar vi en figur som visar dessa envelopes tillsammans med PCF uppskattad för patienten i fråga och undersöker hur de placerar sig i förhållande till varandra. Om kurvan faller inom området som illustrerar acceptansintervallet av vad vi förväntar oss kan vi anta att thomasprocessen är en god modell för datan. Vi undersöker dessa grafer för varje patient i de två grupperna. För att undvika ett cirkulärt resonemang väljer vi att använda skattningarna som gjordes med hjälp av Ripleys K-funktion och sedan använda PCF som teststatistika för att utvärdera modellen. (a) Kurvor motsvarande grupp 1. (b) Kurvor motsvarande grupp 2. Figur 18: De skattade PCF för två patienter (en frisk till vänster och en sjuk till höger) från nervda- tan tillsammans med envelopes för motsvarande simulerad data, samt en röd linje som representerar den genomsnittliga kurvan för thomasprocessen med motsvarande parametrar. 17 Vi valde ut ovanstående figurer från ovannämnda process för att illustrera trenden vi såg bland resultaten. Övriga figurer finns att hitta i Appendix D.2. Genom att observera figurerna ovan verkar det som vår data har liknande mönster som den data som är resultat från en thomasprocess, eftersom kurvan faller inom det förväntade området samt ligger nära den teoretiska röda linjen som ovan betecknas med gpoi(r). Detta stämde för alla patienter och vi kan säga att thomasprocessen är en bra modell för denna nervdata. 4.4.2 Utvärdering via Monte Carlo-test Vi formulerar nollhypotesen enligt H0 = Punktmönstren från nervdatan följer en thomasprocess. (20) Genom att använda oss av ekvationerna (17), (18) i avsnitt 2.6 kan vi bara sätta in PCF istället för S(t). Vi får då Gi = ∑ r∈[0,min(a,b)] [ĝi(r)− gi(r : κ̂i, µ̂i)] 2 (21) och efter insättning av det teoretiska uttrycket för g(r) erhålles Gi = ∑ r∈[0,min(a,b)] [ ĝi(r)− 1− r2 4πκ̂iσ̂2 i exp ( − r2 4σ̂2 i )]2 . (22) För Gobs har vi Gobs = ∑ r∈[0,min(a,b)] [ĝ(r)− g(r : κ̂, µ̂)]. (23) Slutligen beräknar vi p-värdet genom ekvation (19) där vi förkastar H0 utifrån signifikansnivån α = 0.05. Resultaten presenteras nedan i tabell 5 och 6. Det är tydligt att vi i samtliga fall inte kan förkasta H0. Grupp 1 p-värde Patient 1 0.41 Patient 2 0.57 Patient 3 0.50 Patient 4 0.50 Patient 5 0.77 Patient 6 0.49 Patient 7 0.63 Patient 8 0.61 Tabell 5: p-värden för grupp 1. Grupp 2 p-värde Patient 1 0.87 Patient 2 0.64 Patient 3 0.45 Patient 4 0.69 Patient 5 0.50 Patient 6 0.63 Tabell 6: p-värden för grupp 2. 18 4.5 Sammanfattning av resultaten från analysen av nervdata Figurerna nedan jämför de poolade teststatistikorna mot varandra för de två grupperna av pati- enter. (a) Ripleys K-funktion. (b) PCF. Figur 19: Jämförelse av de två sammanfattande statistikorna, Ripleys K-funktion och PCF, för de två grupperna. Genom att observera dessa kurvor kan vi även sammanfatta det tidigare delar av detta avsnitt konstaterat. Båda grupperna är klart klustrade och det är inte ett klart felaktigt antagande att säga att de kommer från en thomasprocess. Vi kan även konstatera att grupp 2 verkar har betydligt färre basnerver, föräldrapunkter, än grupp 1 samt mer tätare och separerade kluster. Det är dock svårt att få fram några intervall som skulle kunna användas för att avgöra om ett nervmönster kommer från en frisk eller sjuk patient eftersom vi har blivit försedda med en väldigt liten datamängd. 5 Sammanfattning Som resultatet av simuleringsstudien visar i sektion 3.2 ger Ripleys K-funktion och PCF olika intervall för vilka värden på σ som är acceptabla. Vid simulering med paremetervärdena (κ, µ) = (22.9, 4) så bör σ ligga i IK,σ eller IP,σ för att MCM-algoritmerna thomas.estK och thomas.estpcf ska konvergera. För Ripleys K-funktion är intervallet Iσ,K = [0.02, 0.09] och för PCF är intervallet Iσ,P = [0.02, 0.03], vid simulering i enhetskvadraten. Att RipleysK-funktion har ett längre intervall för acceptabla värden för σ betyder att mer stabila och säkra uppskattningar kan göras med Ripleys K-funktion och är därför att föredra framför PCF. När vi varierar σ så varieras det förväntade avståndet mellan föräldrapunkt och berörda dot- terpunkter. När små värden på σ används, givet fixerade κ, µ, kommer varje kluster att vara tätare och hela punktmönstret kommer att se mer klumpat ut. Anledningen till att algoritmerna börjar konvergera dåligt och ge sämre skattningar när vi ökar σ över ett visst värde är att dotterpunk- terna till varje kluster börjar överlappa. I ett idealt fall där varje kluster är en cirkel där radien är avståndet från föräldrapunkten till den yttersta dotterpunkten så betyder detta att dessa cirklar överlappar när σ blir stort. Då fås en yta där två cirklar överlappar, en så kallad vesica piscis, se figur 20. När dotterpunkterna befinner sig i denna yta så har algoritmen svårt att veta till vilken förälderpunkt given dotterpunkt tillhör, vilket ger upphov till felskattningar. Dessutom kan algoritmerna heller inte identifiera vilka av punkterna som är dotter- respektive föräldrapunkter och då ger det fel i skattningar för κ. När man får fel skattning på κ skattas även µ fel eftersom µ̂ = λ/κ̂. Figur 14 i sektion 3, samt figurerna i appendix D.3 visar att simuleringarna ger en hel del avvikande värden. Vid skattning av parametrarna genom att använda ett medelvärde jämfört med ett medianvärde ger detta stora skillnader, vilket kan ses i appendix C. Detta är ett välkänt fenomen inom statistiken och generellt blir skattningarna mycket bättre när ett medianvärde används istället 19 Figur 20: Överlappande kluster som delar på dotterpunkter. De ifyllda punkterna är dotterpunkter och ringarna är föräldrapunkter. Notera att det inte är möjligt att veta vilken föräldrapunkt de blå dotterpunkterna tillhör. för ett medelvärde, just på grund av att det finns avvikande värden med extremt höga värden. Detta är något som är viktigt att ha i åtanke när man till exempel vill hitta de gemensamma skattningarna i grupp 1 respektive grupp 2, vid användning av medelvärdet eller medianvärdet kan resultatet skilja sig avsevärt. Som nämns i 4.4 är thomasprocessen ett bra sätt att representera nervdatan. Detta trots att vi bara haft data på antal dotterpunkter, utan någon information angående föräldrapunkterna. Resultatet i sektion 4.5 visar att grupp 2, som var patienter drabbade av diabetesneuropati, generellt hade tätare kluster av nervtrådar än de patienter från grupp 1 samt färre baspunkter vilket minskar antalet totala nervändpunkter. Detta stämmer väl överens med den tidigare forskning som gjorts på diabetesneuropati och som nämns i sektion 1. Utifrån de resultaten vi kommit fram till i denna studie är det dock väldigt svårt att kunna dra några slutsatser huruvida en enskild individ lider av diabetesneuropati eller inte, något som också diskuteras i sektion 4.5. Detta beror främst på att den givna datamängden är väldigt liten, och det är därför svårt att göra en noggrann dataanalys, eftersom det i nuläget är svårt att kvantifiera värden på κ, µ och σ för att kunna åtskilja individer med och utan sjukdomen diabetesneuropati. I simuleringsstudien valde vi värdena för intensiteten av föräldrapunkter κ och det genomsnitt- liga antalet dotterpunter per kluster µ så att de liknar värdena i nervdatan så att det är möjligt att använda resultaten av simuleringsstudien för tolkning av resultaten baserat på nervdatan. Det skulle vara intressant att undersöka hur skattningarna beter sig om andra värden hade valts för κ och µ. Om vi skulle fortsätta detta arbete skulle vi utföra en större simuleringsstudie med flera värden för de tre parametrarna, och eventuellt kunnat ta fram exakt hur mycket σ får skilja sig i förhållande till de andra parametrarna. Det skulle också ha varit intressant att använda de verkliga antalen föräldrapunkter och se om detta hade ändrat skattningarna av de två andra parametrarna i nervdatan. 20 Referenser [1] Ulf Adamsson. Diabetesneuropati. 2013. url: https://www.netdoktorpro.se/diabetes/ medicinska-oversikter/diabetesneuropati/. Accessed: 2020-02-08. [2] Claes Andersson, Peter Guttorp och Aila Särkkä. “Discovering early diabetic neuropathy from epidermal nerve fiber patterns”. I: (2016). doi: 10.1002/sim.7009. [3] Adrian Baddeley.mincontrast. url: https://www.rdocumentation.org/packages/spatstat/ versions/1.61-0/topics/mincontrast. Accessed: 2020-02-24. [4] Julian Besag och Peter J. Diggle. “Simple Monte Carlo Tests for Spatial Pattern”. I: Journal of the Royal Statistical Society 26.3 (1977), s. 327–333. doi: 10.1080/01621459.1975. 10480272. [5] Markus J. Buehler. Property calculation II. 2011. url: https://ocw.mit.edu/courses/ materials- science- and- engineering/3- 021j- introduction- to- modeling- and- simulation-spring-2012/part-i-lectures-readings/MIT3_021JS12_P1_L4.pdf. [6] Peter J. Diggle. Statistical Analysis of Spatial and Spatio-Temporal Point Patterns. 2014, s. 92. [7] Philip M. Dixon. “Ripley’s K-function”. I: Encyclopedia of Environmetrics 3.3 (2002), s. 1796– 1803. doi: https://doi.org/10.1002/9780470057339.var046. [8] Giuseppe Espa, Diego Giuliani och Giuseppe Arbia. “Weighting Ripley’s K-function to ac- count for the firm dimension in the analysis of spatial concentration”. I: (2010). [9] Tim Johnson. Introduction to Spatial Point Processes. 2010. url: https://warwick.ac. uk/fac/sci/statistics/staff/academic-research/nichols/research/spatbayes/ Johnson_SpatialPointProc.pdf. [10] Tim Johnson. Introduction to Spatial Point Processes. 2010. url: https://warwick.ac. uk/fac/sci/statistics/staff/academic-research/nichols/research/spatbayes/ Johnson_SpatialPointProc.pdf. Accessed: 2020-03-24. [11] Dale W. Jorgenson. “Multiple Regression Analysis of a Poisson Process”. I: Journal of the American Statistical Association 56.294 (1961), s. 235–245. doi: https://doi.org/10. 1080/01621459.1961.10482106. [12] Christian Källgren, Hampus Lane och Simon Eriksson. “Modellering av nervmönster med spatiala punkt- processer”. I: (2017). [13] Günter Last och Mathew Penrose. Lectures on the Poisson Process. 2018. [14] Jesper Moller och Rasmus Plenge Waagepetersen. Statistical Inference and Simulation for Spatial Point Processes. 2003. [15] MVEX01-20-13 Punktprocesser och nervtrådar. url: https : / / www . chalmers . se / sv / institutioner/math/utbildning/grundutbildning-chalmers/examensarbeten-och- kandidatprojekt/Kandidatprojekt_2020/Sidor/MVEX01-20-13.aspx. Accessed: 2020- 02-08. [16] Lance Waller, Aila Särkkä, Viktor Olsbo, Mari Myllymäki, Ioanna Panoutsopoulou, William Kennedy och Gwen Wendelschafer-Crabb. “Second-order spatial analysis of epidermal nerve fibers”. I: Statistics in medicine 30 (okt. 2011), s. 2827–41. doi: 10.1002/sim.4315. 21 https://www.netdoktorpro.se/diabetes/medicinska-oversikter/diabetesneuropati/ https://www.netdoktorpro.se/diabetes/medicinska-oversikter/diabetesneuropati/ https://doi.org/10.1002/sim.7009 https://www.rdocumentation.org/packages/spatstat/versions/1.61-0/topics/mincontrast https://www.rdocumentation.org/packages/spatstat/versions/1.61-0/topics/mincontrast https://doi.org/10.1080/01621459.1975.10480272 https://doi.org/10.1080/01621459.1975.10480272 https://ocw.mit.edu/courses/materials-science-and-engineering/3-021j-introduction-to-modeling-and-simulation-spring-2012/part-i-lectures-readings/MIT3_021JS12_P1_L4.pdf https://ocw.mit.edu/courses/materials-science-and-engineering/3-021j-introduction-to-modeling-and-simulation-spring-2012/part-i-lectures-readings/MIT3_021JS12_P1_L4.pdf https://ocw.mit.edu/courses/materials-science-and-engineering/3-021j-introduction-to-modeling-and-simulation-spring-2012/part-i-lectures-readings/MIT3_021JS12_P1_L4.pdf https://doi.org/https://doi.org/10.1002/9780470057339.var046 https://warwick.ac.uk/fac/sci/statistics/staff/academic-research/nichols/research/spatbayes/Johnson_SpatialPointProc.pdf https://warwick.ac.uk/fac/sci/statistics/staff/academic-research/nichols/research/spatbayes/Johnson_SpatialPointProc.pdf https://warwick.ac.uk/fac/sci/statistics/staff/academic-research/nichols/research/spatbayes/Johnson_SpatialPointProc.pdf https://warwick.ac.uk/fac/sci/statistics/staff/academic-research/nichols/research/spatbayes/Johnson_SpatialPointProc.pdf https://warwick.ac.uk/fac/sci/statistics/staff/academic-research/nichols/research/spatbayes/Johnson_SpatialPointProc.pdf https://warwick.ac.uk/fac/sci/statistics/staff/academic-research/nichols/research/spatbayes/Johnson_SpatialPointProc.pdf https://doi.org/https://doi.org/10.1080/01621459.1961.10482106 https://doi.org/https://doi.org/10.1080/01621459.1961.10482106 https://www.chalmers.se/sv/institutioner/math/utbildning/grundutbildning-chalmers/examensarbeten-och-kandidatprojekt/Kandidatprojekt_2020/Sidor/MVEX01-20-13.aspx https://www.chalmers.se/sv/institutioner/math/utbildning/grundutbildning-chalmers/examensarbeten-och-kandidatprojekt/Kandidatprojekt_2020/Sidor/MVEX01-20-13.aspx https://www.chalmers.se/sv/institutioner/math/utbildning/grundutbildning-chalmers/examensarbeten-och-kandidatprojekt/Kandidatprojekt_2020/Sidor/MVEX01-20-13.aspx https://doi.org/10.1002/sim.4315 A Härledning av statistikor för thomasprocessen A.1 Ripleys K-funktion Här presenteras en härledning för den teoretiska K-funktionen för thomasprocessen [12]. Ripleys K-funktion definieras enligt följande ekvation, K(r) = λ−1E[N0(r)], (24) där λ är antalet punkter per enhetsarea, intensiteten, och N0(r) är antalet punkter hos processen inom ett avstånd r till en godtycklig punkt av processen. Antag att vi har en punkt, x, som är placerad i origo och tillhör ett kluster som vi benämner c. Väntevärdet för området med centrum i x och radie r beror på antalet förväntade punkter från kluster c, E[N0,c(r)], och antalet förväntade punkter från andra kluster, E[N0,−c(r)]. Därmed är det möjligt att skriva väntevärdet som en summa av två väntevärden, E[N0(r)] = E[N0,c(r)] + E[N0,−c(r)]. (25) Den andra termen beror inte på något specifikt kluster utan beror bara på arean, πr2, och inten- siteten, λ = κµ. Den andra termen i (25) kan därmed skrivas E[N0,−c(r)] = πr2κµ. För den första termen i (25) använder vi lagen om sammanlagd förväntan E[N0,c(r)] = ∞∑ k=2 E[N0,c(r)|Kc = k]P (Kc = k). (26) Vi summerar väntevärdet av N0,c(r) givet att klustret är av storlek k multiplicerad med sannolik- heten för att klustret är av storlek k. Anledningen till att k börjar från 2 är att ett kluster med 0 punkter är inte rimligt, och när k = 1 så är väntevärdet lika med 0, då vi bara räknar ytterligare punkter. N0,c(r)|Kc = k är därmed antalet ytterligare punkter som tillhör kluster c inom en radie r givet att det finns totalt k punkter i klustret. Sannolikheten för att dessa k− 1 ytterligare punkter befinner sig inom den radien är givet av en fördelningsfunktion, F (r), för avståndet mellan två punkter som befinner sig i samma kluster. Eftersom att varje punkt är oberoende av varandra så är N0,c(r)|Kc = k binomialfördelad med k− 1 utfall med sannolikheten F (r) för varje utfall. Vän- tevärdet av en variabel med binomialfördelning är antalet utfall multiplicerat med sannolikheten för varje utfall, E[N0,c(r)|Kc = k] = (k − 1)F (r). (27) Den andra faktorn i (26), P (Kc = k), är samma som att multiplicera k med sannolikheten att ett kluster är av storlek k, P (Dc = k). Men en sannolikhetsfördelning kräver att ∑∞ k=1 P (Kc = k) = 1. Därför delar vi P (Kc = k) för k = 1, 2, ... med ∑∞ l=1 lP (Dc = l), vilket ger ∞∑ k=1 P (Kc = k) = ∑∞ k=1 kP (Dc = k)∑∞ l=1 lP (Dc = l) . Men då ∑∞ l=1 lP (Dc = l) = E[Dc] = µ är det möjligt att förenkla uttrycket på följande sätt P (Kc = k) = kP (Dc = k) µ (28) 22 Givet (27), (28) och faktumet att Dc ∼ P(µ) kan (26) utvecklas, E[N0,c(r)] = ∞∑ k=2 E[N0,c(r)|Kc = k]P (Kc = k) = ∞∑ k=2 (k − 1)F (r) kP (Dc = k) µ = F (r) µ ∞∑ k=2 (k − 1)kP (Dc = k) = F (r) µ E[Dc(Dc − 1)] = F (r) µ (µ2 + µ− µ) = µF (r). (29) För att utveckla F (r), avståndet mellan två oberoende och normalfördelade variabler, krävs ytter- ligare en härledning. Anta att X = (x1, x2) och Y = (y1, y2) är två punkter i R2, där x1, x2, y1 och y2 är oberoende och N (a, ω2), som ger F (r) = P (||X − Y || ≤ r) = P (||X − Y ||2 ≤ r2) = P ((x1 − y1)2 + (x2 − y2)2 ≤ r2). (30) Då x1 − y1 och x2 − y2 båda har väntevärdet 0 och variansen 2ω2, kan vi multiplicera dem med inversen av variansen, som då ger en ny fördelning, N (0, 1) innanför kvadraterna som beskrivs av (30). P ((x1 − y1)2 + (x2 − y2)2 ≤ r2) = P (( x1 − y1√ 2ω )2 + ( x2 − y2√ 2ω )2 ≤ r2 2ω2 ) . (31) Summan av kvadraterna i högerledet av ekvation (31) visar en χ2-fördelad variabel med två fri- hetsgrader, ( x1 − y1√ 2ω )2 + ( x2 − y2√ 2ω )2 = Q ∼ χ2 2. Fördelningsfunktionen för en χ2-fördelad variabel med två frihetsgrader är följande, P ( Q ≤ r2 2ω2 ) = 1− exp ( −r2 4ω2 ) . (32) Sista steget är att koppla samman E[N0,−c(r)] = πr2κµ, (29), (32), med ω = σ, K(r) = 1 κµ (E[N0,c(r)] + E[N0,−c(r)]) = πr2 + 1 κ F (r) = πr2 + 1 κ ( 1− exp ( − r2 4σ2 )) . (33) 23 A.2 Parkorrelationsfunktionen Här presenteras än härledning av den teoretiska PCF för thomasprocessen. PCF definieras enligt g(r) = K ′(r) 2πr , (34) där K(r) = πr2 + 1 κ (1− exp ( −r2 4σ2 ) ). (35) Detta leder till att K ′(r) = 2πr + 2r 4κσ2 exp ( −r2 4σ2 ) , (36) vilket ger oss den teoretiska funktionen för PCF i en thomasprocess enligt: g(r) = K ′(r) 2πr = 2πr + 2r 4κσ2 exp ( −r2 4σ2 ) 2πr = 1 + 1 4πκσ2 exp ( −r2 4σ2 ) . (37) 24 B Simulering B.1 Ändring av data och skalning av simuleringsfönster I genereringen av en realisation av en thomasprocess kräver algoritmen dimensionerna av det fönster i vilken punkterna genereras. Som standard tar funktionen rThomas in parametern c(0,1), c(0,1), vilket betyder att punkterna genereras inom enhetskvadraten [0, 1]× [0, 1]. För att på ett lämpligt sätt kunna simulera ny data som liknar den givna nervdatan vill vi ha dem i ungefär samma storlek. Den givna datan har flera brister som vi måste åtgärda för att kunna uppnå detta. Först och främst måste vi se till att alla koordinater ligger i den första kvadranten, dvs. endast positiva koordinater. Därför skalar vi de nervmönster med negativa koordinater med −1, vilket resulterar i att koordinaterna speglas i motsvarande axel. Detta kommer inte påverka resultatet av våra skattningar. Sedan vill vi se till att alla fönster börjar i origo och sträcker sig till samma x- respektive y-gräns. Vi väljer dessa gränser och skjuter koordinaterna i x- och y-led så att alla punktmönster får plats i samma fönster, nämligen [0, 433]× [0, 330]. Slutligen skalar vi fönstret och koordinaterna med 1 330 längs båda axlarna. Resultatet blir då att nervdatan nu ligger i fönstret [0, 1.31]× [0, 1]. Detta är så nära den enhetskvadraten vi kan komma, eftersom att skala om nervdatan till en kvadrat påverkar skattningen av σ på ett sätt som vi vill undvika. Vi kommer hädanefter låta λ̃ beteckna det totala antalet dotterpunkter i fönstret [0, 1.31]×[0, 1] och därmed, eftersom λ̃ = κ̃µ, kommer vi låta κ̃ beteckna antalet föräldrapunkter i fönstret [0, 1.31]× [0, 1]. För µ̃ = µ blir det ingen skillnad då det inte är relaterat till storleken av fönstret. Detta betyder att när vi ska mata in värden för κ och σ i rThomas tvingas vi dividera de parametrar vi uppskattat utifrån datan med 1.31. Dessa skalade parametrar kommer refereras till som κ och σ. För att sammanfatta detta stycke matematiskt låter vi den generiska rektangeln i vilken punk- terna genereras vara Ω och dess area |Ω|. Då gäller κ = κ̃ |Ω| , Ω = [0, A]× [0, B], |Ω| = AB, A,B ∈ R, (38) Nedan i Figur B.21 kan vi se bilder på nervdatan före och efter skalning. (a) Innan skalning. (b) Efter skalning. Figur B.21: Data från person 1, grupp 1, innan och efter skalning. Avslutningsvis måste noteras att när vi ska utföra de slutgiltiga analyserna av parametrarna så kommer vi vilja skala om de tillbaka till det ursprungliga dimensionen av datan. Exempelvis kommer intensiteten av föräldrapunkter i det stora fönstret [0, 433] × [0, 330] vara κ∗ = κ 433×330 , där κ är totala antalet föräldrapunkter i nervmönstret. 25 B.2 Framtagning av lämpliga parametrar för simuleringsstudien I simuleringsstudien användes det totala medelvärdet för båda grupper, det vill säga n̄ = n̄1. + n̄2. 2 ≈ 125.6. (39) Detta är också det förväntade antalet dotterpunkter λ̃ för thomasprocessen i det givna fönstret. Eftersom vi simulerar i fönstret [0, 1]×[0, 1] vill vi använda oss av λ istället. Detta kan åstadkommas genom att först bestämma κ̃. För att bättre förstå hur fönsterstorleken påverkar skattningarna, se appendix B.1. Eftersom λ̃ = κ̃µ gäller det att hitta lämpliga κ̃ och µ sådana att deras produkt kommer nära λ̃ = 125.6. Parametern κ̃ representerar antalet föräldrapunkter och µ representerar antalet dotterpunkter till varje föräldrapunkt. Varje föräldrapunkt är i själva verket en baspunkt utifrån vilken nervtrådar förgrenar sig från, alltså måste det i de flesta fall gälla att antal föräldrapunkter är färre än det totala antalet dotterpunkter. Ett rimligt värde på µ ligger mellan 2 och 5 [2]. Genom att välja ett µ ∈ [2, 5] så kan vi ta fram ett rimligt κ̃ utifrån λ̃ = κ̃µ. Vi har valt µ = 4 och får därmed κ̃ = 30 =⇒ κ = 30 1.31 ≈ 22.9, vilket slutligen ger λ = 22.9 · 4 ≈ 91.6. Notera att parametrarna λ = 91.6, κ = 22.9 och µ = 4 är de parametrar som nervdatan genomsnittligt hade givit i thomasprocessen om det istället varit ett [0, 1]× [0, 1] fönster och det är dessa vi matade in algoritmen rThomas som in-parametrar under simuleringsstudien. 26 C Tabeller och resultat Resultat av skattningar från simuleringsstudien i tabellform. Verkliga värden κ = 22.9 µ = 4 σ = 0.01 : 0.15 Skattningar κ̂ µ̂ σ̂ σ = 0.01 24.15359 24.417317 0.03083419 σ = 0.02 25.45370 9.843589 0.02875202 σ = 0.03 26.15853 11.495781 0.04218288 σ = 0.04 26.98548 7.130926 0.04560641 σ = 0.05 29.92780 12.640477 0.06679906 σ = 0.06 31.18656 17.959941 0.09168543 σ = 0.07 34.23076 29.757148 0.16791205 σ = 0.08 56.22653 77.464241 0.49188860 σ = 0.09 91.61842 131.678098 0.99224079 σ = 0.10 155.52405 187.217000 2.03331761 σ = 0.11 308.54949 206.823284 3.84776921 σ = 0.12 471.68262 254.867052 5.92067473 σ = 0.13 666.46080 266.121395 8.18679975 σ = 0.14 944.84811 267.339053 10.97474268 σ = 0.15 1298.63805 222.187319 15.07363424 Tabell C.7: Parameterskattningar med K-funktionen då medelvärde används. 27 Verkliga värden κ = 22.9 µ = 4 σ = 0.01 : 0.15 Skattningar κ̂ µ̂ σ̂ σ = 0.01 25.62470 5.681309 0.08201874 σ = 0.02 24.38073 6.661124 0.03834710 σ = 0.03 25.96076 5.521146 0.03024794 σ = 0.04 28.54149 3.887598 0.03760578 σ = 0.05 36.60449 4.273384 0.10420623 σ = 0.06 49.78715 4.388362 0.11217682 σ = 0.07 102.24685 11.059789 0.37786868 σ = 0.08 214.59063 20.681136 0.78503836 σ = 0.09 295.34809 35.965411 1.08958946 σ = 0.10 457.87318 71.349599 2.50063923 σ = 0.11 540.95338 109.523013 3.87994809 σ = 0.12 935.75361 105.090919 7.60277767 σ = 0.13 1054.10587 131.930662 8.17112979 σ = 0.14 1410.07767 144.104112 11.82033020 σ = 0.15 1558.85421 147.439841 13.94213757 Tabell C.8: Parameterskattningar med PCF då medelvärde används. Verkliga värden κ = 22.9 µ = 4 σ = 0.01 : 0.15 Skattningar κ̂ µ̂ σ̂ σ = 0.01 24.026030 3.760774 0.009809099 σ = 0.02 24.619311 3.638875 0.019179392 σ = 0.03 25.414828 3.587787 0.028593227 σ = 0.04 24.690064 3.567788 0.039147788 σ = 0.05 26.064837 3.513402 0.048603808 σ = 0.06 25.082254 3.611746 0.060528256 σ = 0.07 24.479339 3.621957 0.070765998 σ = 0.08 22.003507 4.047569 0.083422713 σ = 0.09 20.338054 4.404147 0.098905060 σ = 0.10 15.130373 5.892428 0.129802895 σ = 0.11 13.353137 6.674042 0.175497927 σ = 0.12 11.785156 7.465746 0.284439483 σ = 0.13 11.636988 7.545923 0.581603998 σ = 0.14 7.048926 12.079488 0.981315831 σ = 0.15 8.092667 10.399924 1.410125572 Tabell C.9: Parameterskattningar med K-funktionen då median används. 28 Verkliga värden κ = 22.9 µ = 4 σ = 0.01 : 0.15 Skattningar κ̂ µ̂ σ̂ σ = 0.01 22.91234 3.945210 0.009904673 σ = 0.02 23.23429 3.830248 0.018518627 σ = 0.03 24.59410 3.733664 0.027270609 σ = 0.04 26.42808 3.421231 0.035787127 σ = 0.05 28.58298 3.114991 0.042898650 σ = 0.06 29.63581 3.055611 0.050125561 σ = 0.07 31.29329 2.865039 0.057859330 σ = 0.08 35.63664 2.507547 0.063570632 σ = 0.09 34.40629 2.667989 0.073551943 σ = 0.10 38.09664 2.321308 0.077499501 σ = 0.11 42.82919 2.128642 0.082632651 σ = 0.12 40.33463 2.124513 0.092541965 σ = 0.13 46.68532 1.835994 0.102312497 σ = 0.14 59.41088 1.525845 0.103795704 σ = 0.15 83.45480 1.105411 0.115755801 Tabell C.10: Parameterskattningar med PCF då median används. 29 D Figurer D.1 Punktmönster från nervdatan I sektion D.1.1 och D.1.2 presenteras en visualisering av nervdatan från endpoints.rds för gruppen utan diabetesneuropati (grupp 1) respektive gruppen med diabetesneuropati (grupp 2). D.1.1 Grupp 1 (a) Punktmönster för Grupp 1 - Patient 1. (b) Punktmönster för Grupp 1 - Patient 2. (c) Punktmönster för Grupp 1 - Patient 3. (d) Punktmönster för Grupp 1 - Patient 4. 30 (e) Punktmönster för Grupp 1 - Patient 5. (f) Punktmönster för Grupp 1 - Patient 6. (g) Punktmönster för Grupp 1 - Patient 7. (h) Punktmönster för Grupp 1 - Patient 8. 31 D.1.2 Grupp 2 (i) Punktmönster för Grupp 2 - Patient 1. (j) Punktmönster för Grupp 2 - Patient 2. (k) Punktmönster för Grupp 2 - Patient 3. (l) Punktmönster för Grupp 2 - Patient 4. (m) Punktmönster för Grupp 2 - Patient 5. (n) Punktmönster för Grupp 2 - Patient 6. 32 D.2 Grafisk undersökning för anpassningen av thomasprocessen I sektion D.2.1 och D.2.2 presenteras resultatet över den skattade PCF för patienterna i grupp 1 repsektive grupp 2, tillsammans med envelopes från den simulerade datan och en röd linje som visar den genomsnittliga kurvan för thomasprocessen med motsvarande parametrar. D.2.1 Grupp 1 (a) Grupp 1 - Patient 1. (b) Grupp 1 - Patient 2. (c) Grupp 1 - Patient 3. (d) Grupp 1 - Patient 4. (e) Grupp 1 - Patient 5. (f) Grupp 1 - Patient 6. 33 (g) Grupp 1 - Patient 7. (h) Grupp 1 - Patient 8. 34 D.2.2 Grupp 2 (i) Grupp 2 - Patient 1. (j) Grupp 2 - Patient 2. (k) Grupp 2 - Patient 3. (l) Grupp 2 - Patient 4. (m) Grupp 2 - Patient 5. (n) Grupp 2 - Patient 6. 35 D.3 Violindiagram I sektion D.3.1 och D.3.2 presenteras de skattade parametrarna κ̂, µ̂ och σ̂ från Ripleys K-funktion respektive PCF för att visa spridningen på skattningarna. Notera att y-axlarna är logaritmerade. D.3.1 Ripleys K-funktion Figur D.22: Violindiagram över uppskattade värden för κ för olika värden på σ med Ripleys K- funktion. Figur D.23: Violindiagram över uppskattade värden för µ för olika värden på σ med Ripleys K- funktion. 36 Figur D.24: Violindiagram över uppskattade värden för σ för olika värden på σ med Ripleys K- funktion. D.3.2 Parkorrelationsfunktionen Figur D.25: Violindiagram över uppskattade värden för κ för olika värden på σ med PCF. 37 Figur D.26: Violindiagram över uppskattade värden för µ för olika värden på σ med PCF. Figur D.27: Violindiagram över uppskattade värden för σ för olika värden på σ med PCF. 38 E Kod E.1 Kod för histogram för skattningar då medelvärde används 1 library("latex2exp") 2 library("spatstat") 3 4 kappa = 22.9 5 mu = 4 6 7 Hist_mean = function(kappa , mu) { 8 9 timesteps = 100 10 kappaVector = vector("numeric", timesteps) 11 sigmaVector = vector("numeric", timesteps) 12 muVector = vector("numeric", timesteps) 13 numberOfPoints = vector("numeric", timesteps) 14 estimates = matrix(0, nrow = 15, ncol = 4) 15 sigmas = seq(from = 0.03, to = 0.07, by = 0.04) 16 17 for (ind in seq_along(sigmas)){ 18 19 for (i in 1: timesteps) { 20 21 thomasSim = rThomas(kappa = kappa , scale = sigmas[ind], mu = mu, win = owin(c (0,1),c(0,1))) 22 23 numberOfPoints[i] = thomasSim [["n"]] 24 25 result = thomas.estK(thomasSim) 26 27 kappaVector[i] = result$modelpar [[1]] 28 sigmaVector[i] = result$modelpar [[2]] 29 muVector[i] = numberOfPoints[i]/kappaVector[i] 30 31 } 32 33 hist(kappaVector , 34 main = TeX(’$\\ sigma$ = 0.07$’), 35 xlab = TeX(’$\\ kappa$’), 36 ylab = "Frekvens", 37 border = "blue", 38 col = "grey", 39 las = 1, 40 breaks = 50) 41 42 hist(muVector , 43 main = TeX(’$\\ sigma = 0.07$’), 44 xlab = TeX(’$\\mu$’), 45 ylab = "Frekvens", 46 border = "blue", 47 col = "grey", 48 las = 1, 49 breaks = 50) 50 51 } 52 } 53 54 Hist_mean (22.9 ,4) 39 E.2 Skillnaden i avvikande datapunkter, median kontra medelvärde 1 library("latex2exp") 2 library("spatstat") 3 4 kappa = 22.9 5 mu = 4 6 7 ripleysMean = function(kappa , mu) { 8 9 timesteps = 100 10 kappaVector = vector("numeric", timesteps) 11 sigmaVector = vector("numeric", timesteps) 12 muVector = vector("numeric", timesteps) 13 numberOfPoints = vector("numeric", timesteps) 14 estimates = matrix(0, nrow = 15, ncol = 4) 15 sigmas = seq(from = 0.03, to = 0.07, by = 0.04) 16 17 for (ind in seq_along(sigmas)){ 18 19 for (i in 1: timesteps) { 20 21 thomasSim = rThomas(kappa = kappa , scale = sigmas[ind], mu = mu, win = owin(c (0,1),c(0,1))) 22 23 numberOfPoints[i] = thomasSim [["n"]] 24 25 result = thomas.estK(thomasSim) 26 27 kappaVector[i] = result$modelpar [[1]] 28 sigmaVector[i] = result$modelpar [[2]] 29 muVector[i] = numberOfPoints[i]/kappaVector[i] 30 31 } 32 33 hist(kappaVector , 34 main = TeX(’$\\ sigma$ = 0.07$’), 35 xlab = TeX(’$\\ kappa$’), 36 ylab = "Frekvens", 37 border = "blue", 38 col = "grey", 39 las = 1, 40 breaks = 50) 41 42 hist(muVector , 43 main = TeX(’$\\ sigma = 0.07$’), 44 xlab = TeX(’$\\mu$’), 45 ylab = "Frekvens", 46 border = "blue", 47 col = "grey", 48 las = 1, 49 breaks = 50) 50 51 } 52 print(mean(kappaVector)) 53 print(median(kappaVector)) 54 print(mean(muVector)) 55 print(median(muVector)) 56 } 57 58 ripleysMean (22.9 ,4) 59 60 61 62 63 64 65 66 67 40 68 PCFMean = function(kappa , mu) { 69 70 timesteps = 100 71 kappaVector = vector("numeric", timesteps) 72 sigmaVector = vector("numeric", timesteps) 73 muVector = vector("numeric", timesteps) 74 numberOfPoints = vector("numeric", timesteps) 75 estimates = matrix(0, nrow = 15, ncol = 4) 76 sigmas <- seq(from = 0.01, to = 0.15, by = 0.01) 77 78 for (ind in seq_along(sigmas)){ 79 80 for (i in 1: timesteps) { 81 82 thomasSim = rThomas(kappa = kappa , scale = sigmas[ind], mu = mu, win = owin(c (0,1),c(0,1))) 83 84 numberOfPoints[i] = thomasSim [["n"]] 85 86 result = thomas.estpcf(thomasSim) 87 88 kappaVector[i] = result$modelpar [[1]] 89 sigmaVector[i] = result$modelpar [[2]] 90 muVector[i] = numberOfPoints[i]/kappaVector[i] 91 92 } 93 94 meanSigma = mean(sigmaVector) 95 meanKappa = mean(kappaVector) 96 meanMu = mean(muVector) 97 98 estimates[ind ,4] = meanSigma 99 estimates[ind ,3] = meanMu 100 estimates[ind ,2] = meanKappa 101 estimates[ind ,1] = sigmas[ind] 102 103 } 104 105 return(estimates) 106 107 } 108 109 PCFMean (22.9 ,4) 41 E.3 Kod för parameterskattningar då medelvärde används 1 ############################ Run these first 2 3 require("spatstat") 4 require(ggplot2) 5 require(reshape2) 6 7 ############################ Estimations for Ripleys with mean 8 9 rm(list = ls()) 10 kappa = 22.9 11 mu = 4 12 timesteps = 100 13 kappaVector = vector("numeric", timesteps) 14 sigmaVector = vector("numeric", timesteps) 15 muVector = vector("numeric", timesteps) 16 numberOfPoints = vector("numeric", timesteps) 17 estimates = matrix(0, nrow = 15, ncol = 4) 18 sigmas = seq(from = 0.01, to = 0.15, by = 0.01) 19 20 for (ind in seq_along(sigmas)){ 21 22 for (i in 1: timesteps) { 23 24 thomasSim = rThomas(kappa = kappa , scale = sigmas[ind], mu = mu, win = owin(c (0,1),c(0,1))) 25 26 numberOfPoints[i] = thomasSim [["n"]] 27 28 result = thomas.estK(thomasSim) 29 30 kappaVector[i] = result$modelpar [[1]] 31 sigmaVector[i] = result$modelpar [[2]] 32 muVector[i] = numberOfPoints[i]/kappaVector[i] 33 34 } 35 36 meanSigma = mean(sigmaVector) 37 meanKappa = mean(kappaVector) 38 meanMu = mean(muVector) 39 40 estimates[ind ,4] = meanSigma 41 estimates[ind ,3] = meanMu 42 estimates[ind ,2] = meanKappa 43 estimates[ind ,1] = sigmas[ind] 44 45 } 46 47 print(estimates) 48 df = as.data.frame(estimates) 49 col_names = paste(c("Varying_sigma", "Kappa", "Mu", "Sigma"), sep = "") 50 names(df) = col_names 51 df_mod = melt(df , id.vars = ’Varying_sigma ’, variable.name = ’Skattningar ’) 52 ggplot(df_mod , aes(Varying_sigma , value)) + geom_line(aes(colour = Skattningar)) + ggtitle("Skattningar f r varierande sigma (Ripleys/ M e d e l v rde )") + xlab("Sigma" ) + ylab(" V rde ") 53 54 ############################ Estimations for PCF with mean 55 56 rm(list = ls()) 57 kappa = 22.9 58 mu = 4 59 timesteps = 1000 60 kappaVector = vector("numeric", timesteps) 61 sigmaVector = vector("numeric", timesteps) 62 muVector = vector("numeric", timesteps) 63 numberOfPoints = vector("numeric", timesteps) 64 estimates = matrix(0, nrow = 15, ncol = 4) 65 sigmas = seq(from = 0.01, to = 0.15, by = 0.01) 42 66 67 for (ind in seq_along(sigmas)){ 68 69 for (i in 1: timesteps) { 70 71 thomasSim = rThomas(kappa = kappa , scale = sigmas[ind], mu = mu, win = owin(c (0,1),c(0,1))) 72 73 numberOfPoints[i] = thomasSim [["n"]] 74 75 result = thomas.estpcf(thomasSim) 76 77 kappaVector[i] = result$modelpar 78 sigmaVector[i] = result$modelpar [[2]] 79 muVector[i] = numberOfPoints[i]/kappaVector[i] 80 81 } 82 83 meanSigma = mean(sigmaVector) 84 meanKappa = mean(kappaVector) 85 meanMu = mean(muVector) 86 87 estimates[ind ,4] = meanSigma 88 estimates[ind ,3] = meanMu 89 estimates[ind ,2] = meanKappa 90 estimates[ind ,1] = sigmas[ind] 91 92 } 93 94 print(estimates) 95 df = as.data.frame(estimates) 96 col_names = paste(c("Varying_sigma", "Kappa", "Mu", "Sigma"), sep = "") 97 names(df) = col_names 98 df_mod = melt(df , id.vars = ’Varying_sigma ’, variable.name = ’Skattningar ’) 99 ggplot(df_mod , aes(Varying_sigma , value)) + geom_line(aes(colour = Skattningar)) + ggtitle("Skattningar f r varierande sigma (PCF/ M e d e l v rde )") + xlab("Sigma") + ylab(" V rde ") 43 E.4 Kod för visualisering av statistikorna 1 library("spatstat") 2 3 # Generate CSR point pattern 4 PoissonCSR <- rThomas(mu = 4, kappa = 25, scale = 0.1, win = owin(c(0,1), c(0,1))) 5 png(filename = "pointpattern_CSR.png", bg = "white") 6 par(mar=c(0,0,0,0)) 7 plot(PoissonCSR , pch = 20, main = "") 8 dev.off() 9 10 # Kest for CSR 11 CSR_Kest <- Kest(PoissonCSR , correction = "Ripley") 12 png(filename = "Kest_CSR.png", width = 800, bg = "white") 13 par(mar=c(