Matematiska modeller av läkemedelsprojekt Känslighetsanalys som verktyg i beslutsfattning Examensarbete för kandidatexamen i matematik vid Göteborgs universitet Kandidatarbete inom civilingenjörsutbildningen vid Chalmers Abraham Deniz Elijah Ferreira Erik Johansson Hanna Johansson Jacob Lindbäck Linus Wiskman Institutionen för matematiska vetenskaper Chalmers tekniska högskola Göteborgs universitet Göteborg 2017 Matematiska modeller av läkemedelsprojekt Känslighetsanalys som verktyg i beslutsfattning Examensarbete för kandidatexamen i matematik vid Göteborgs universitet Abraham Deniz Elijah Ferreira Kandidatarbete i matematik inom civilingenjörsprogrammet Bioteknik vid Chalmers Hanna Johansson Kandidatarbete i matematik inom civilingenjörsprogrammet Teknisk fysik vid Chalmers Linus Wiskman Kandidatarbete i matematik inom civilingenjörsprogrammet Teknisk matematik vid Chalmers Erik Johansson Jacob Lindbäck Handledare: David Bolin Matematiska vetenskaper Magnus Ytterstad Captario Examinator: Maria Roginskaya Marina Axelson-Fisk Institutionen för matematiska vetenskaper Chalmers tekniska högskola Göteborgs universitet Göteborg 2017 Populärvetenskaplig presentation Kostnaden för att ta fram ett nytt läkemedel har tredubblats under de senaste 20 åren, samtidigt som endast en liten andel av alla upptäckta läkemedel når marknaden. Detta beror huvudsakligen på att många nya läkemedel inte klarar av de höga krav som ställs. Det är därför viktigt att veta hur läkemedelsutveckling kan effektiviseras. Vi har med hjälp av matematik undersökt om hur man kan åstadkomma just detta. Läkemedelsutveckling Läkemedelsutvecklingen innefattar hela processen från upptäckten av ett nytt läkemedel till att det börjar säljas på marknaden. I processen utför läkemedelsföretaget ett flertal tester för att undersöka om läkemedlet är säkert och effektivt. Det är dessa tester som tar lång tid och är kostsamma. Efter varje test måste läkemedelsutvecklaren bestämma om det är värt att fortsätta med kommande tester, eller om risken att projektet inte blir lönsamt är för stor. Även om ett läkemedel klarar alla tester så finns det andra anledningar som kan påverka lönsamheten av projektet. Exempelvis är lanseringstiden betydelsefull. Om någon konkur- rent börjar sälja ett liknande läkemedel först så riskerar man att få väsentligt mycket lägre försäljningsvolym. Detta kan i värsta fall resultera i att man inte täcker upp kostnaderna för framtagningen av läkemedlet. Dessutom måste vinsterna för lyckade projekt täcka förlusterna av de misslyckade, vilket drar upp vinstkraven. Därför är det viktigt att identifiera vad som är viktigt för att få upp vinsten i läkemedelsprojekt. Det är här matematiken kommer in. Matematiska modeller Med matematikens hjälp kan man ställa upp matematiska modeller för utvecklingen av läke- medel. När man skapar en modell försöker man efterlikna verkligheten för att kunna förutspå vad som kan hända. Vidare kan en modell ses som en förenkling av någonting komplicerat, till någonting mer överskådligt. Rent matematiskt innebär detta att man försöker uttrycka det man är intresserad av med hjälp av en formel, som tar in värden som har att göra med det man försöker beskriva. Dessa värden kallar vi för variabler. En variabel skulle kunna vara hur många konkurrenter som utvecklar liknande läkemedel, eller hur mycket ett test kostar. Ett klassiskt sätt att beskriva denna formel är som en svart låda. Den svarta lådan tar in flera variabler, och ger ut exempelvis vinsten på ett projekt. Problemet är att man i allmänhet inte vet i förväg om t.ex. hur många konkurrenter man har, eller hur höga testkostnaderna kommer bli. Därför måste man göra flera kvalificerade gissningar, och beräkna med hjälp av formeln vinsten för varje gissning. Matematiken förser oss med ett systematiskt sätt att göra detta. Med hjälp av så kallade Monte Carlo-metoder kan en förväntad vinst fås fram på ett effektivt sätt genom datorberäkningar. En bild på en matematisk modell för ett läkemedelsprojekt. Den stora lådan är vår modell. Captario och vårt arbete Naturligtvis kan en matematisk modell för ett läkemedelsprojekt se ut på flera olika sätt. Företaget Captario har tagit fram ett nät-baserat verktyg där användaren själv kan skapa en modell och använda Monte Carlo-metoder för att uppskatta exempelvis förväntad vinst. Det är där vi kommer in i bilden. I vårt arbete har vi gjort en undersökning av vilka variabler som är viktigast för att driva upp vinsten på ett läkemedelsprojekt. De mest informativa metoder- na kan vara väldigt tidskrävande att genomföra, även om snabba datorer finns tillgängliga. Dock sparas beräkningstid om modellen har färre variabler. Därför började vi med att sortera bort variabler som inte påverkade vinst tillräckligt mycket med hjälp enklare metoder. Efter att vi hade färre variabler gjorde vi en undersökning av vad som spelade störst roll i vad man får för vinst med hjälp av de mer sofistikerade metoderna. Ett tydligt resultat vi fick var att det allra viktigaste, ur lönsamhetsperspektiv är att se till att läkemedlet lanseras så tidigt som möjligt. En spontan tanke är kanske att det är kostnaden som kommer vara avgörande för att garantera vinst. Skälet till att tiderna är såpass viktiga är för att försäljningsvolymen sjunker om konkurrenterna lanserar läkemedlet först. Om man då inte hade haft hjälp av ma- tematiken för att ta reda på detta, skulle man möjligtvis kommit fram till helt fel lösningar på sitt problem. Sammanfattning Läkemedelsutveckling styrs idag av att aktörer investerar i de läkemedelsprojekt som tros kunna ge stor avkastning. Hela processen från att ha hittat ett potentiellt läkemedel till eventuell försäljning på marknaden innefattar en rad olika studier och kontroller med höga kostnader. Dessutom svarar de allra flesta läkemedelsprojekt som genomgår denna process ej mot de höga kraven som ställs på nya läkemedel och läggs därmed ned. Detta resulterar sedermera i att läkemedelsföretaget inte får tillbaka det investerade kapitalet. Mot denna bakgrund har företaget Captario utvecklat ett nätbaserat verktyg, Captario SUM. I denna programvara kan användaren konstruera en matematisk modell över sitt läkemedelsprojekt och simulera utfallet för att erhålla en prognos över lönsamheten och andra intressanta aspekter av projektet. Vi har undersökt vilken betydelse variationer i variablerna till en sådan modell har på prognosen som erhålls. Dessutom har vi begränsat oss till att enbart fokusera på pro- gnosen över nettonuvärde, som är ett mått på lönsamhet. Till vår hjälp använder vi oss utav metoder inom fältet känslighetsanalys som har implementerats på en representativ modell av ett läkemedelsprojekt. Ett sekundärt mål under projektets gång har varit att försöka hitta lämpliga metoder som Captario i framtiden ska kunna integrera med sitt verktyg. I förlängningen vill Captario kunna erbjuda sina kunder en känslighetsanalys av sina modeller. Vi har efter denna undersökning hittat flera känslighetsanalysmetoder som lämpar sig väl för Captarios verktyg. En genomgående trend i samtliga metoder är att rekryterings- hastigheten för en av de större studierna var den variabel som hade störst betydelse för lönsamheten. Vidare kan denna rapport ses som ett första steg i analysen av läkeme- delsmodeller, och det diskuteras hur man skulle kunna gå vidare med en djupare analys. Abstract Drug development is today controlled by investing in projects or drugs that are believed to provide a big return on investment. The entire process of finding a potential drug until it is possibly sold on the market involves a variety of tests and controls resulting in huge expenses. The vast majority of drug projects undergoing this process do not meet the high demands set on new drugs and is therefore discontinued. This subsequently results in a loss for the pharmaceutical company in that project. It is against this background that Captario has developed a cloud-based tool, Captario SUM. In Captario SUM, an user is able to construct a mathematical model over their pharmaceutical project and simulate the outcome to obtain a forecast of profitability and other interesting aspects of the project. In this report we investigate the impact variations in the variables, of such a model, have on the forecast that is obtained. We have chosen to focus only on the forecast of net present value which is a measure of profitability. To our help, we use methods in the field emph sensitivity analysis which have been implemented on a prototypical model of a pharmaceutical project. A secondary goal during the project has been to try to find suitable methods that Captario in the future will be able to integrate with the tool and further offer their customers a sensitivity analysis on their models. After this investigation we found various sensitivity analysis methods that are well suitable for Captarios tool. All methods identified that the recruitment rate for one of the larger studies was the variable that had the greatest impact for profitability. Fur- thermore, this report can be seen as a first step in the analysis of drug development models and it is discussed how one could proceed with a deeper analysis. Förord Vi vill tacka Captario för att ha gett oss möjligheten att fördjupa oss i detta problem. Framförallt vill vi rikta ett stort tack till Magnus Ytterstad som har gett oss mycket insikter i arbetet, framförallt om hur det är att jobba i projekt. Vi vill även tacka vår handledare på Chalmers, David Bolin, som har gett oss mycket vägledning genom hela projektet och varit tillmötesgående när vi stött på problem. Bidragsrapport Under det här projektet har en dagbok och en individuell tidslogg förts. Båda dessa redovisar vad varje gruppmedlem bidragit med, men även vad vi har gjort varje vecka och de problem som har uppstått. Den större delen av projektet har bestått av individuellt arbete eller arbete i mindre grupper. Fördelning av arbetsuppgifter har skett löpande på veckovisa möten där vi även har diskuterat problem som har uppstått. Vi har i regel haft möten en gång i veckan med vår handledare på Captario för att diskutera projektets upplägg och mål. MATLAB-funktionen som ligger till grund för implementeringen har skapats av Jacob och Linus. De har även implementerat Sobols index, FAST och Morris metod. Morris metod har implementerats med hjälp av Elijah. Den större delen av litteraturstudier och implemente- ring av Monte Carlo-filtrering har också utförts av Elijah. Den populärvetenskapliga artikeln som följer med rapporten har skrivits av Erik och Abraham. De har även implementerat den grafiska känslighetsanalysen. Varje kapitel i rapporten tillskrivs en eller flera författare. Alla i gruppen har under ar- betet bidragit med korrekturläsning. Syfte: Jacob Bakgrund: Hanna / Elijah Teori: Elijah / Jacob / Hanna / Abraham Metod: Jacob / Linus / Hanna Resultat: Samtliga i gruppen Diskussion: Jacob / Hanna / Erik Appendix: Abraham / Jacob Innehåll 1 Bakgrund 1 1.1 Beslutsfattning inom läkemedelsutveckling . . . . . . . . . . . . . . . . . . . . 1 1.2 Captario och programvaran SUM . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.3 Känslighetsanalys . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Syfte 2 3 Teori 3 3.1 Notation för allmänna modeller . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3.2 Metoder inom känslighetsanalys . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3.2.1 One-at-a-time-experiment . . . . . . . . . . . . . . . . . . . . . . . . . 3 3.2.2 Morris metod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 3.2.3 Monte Carlo-filtrering . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3.2.4 Sobols index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3.2.5 Fourier Amplitude Sensitivity Test – FAST . . . . . . . . . . . . . . . 7 4 Metod 10 4.1 Evaluering av NPV för modeller i SUM . . . . . . . . . . . . . . . . . . . . . 10 4.1.1 Marknadsmodell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4.1.2 Val av responsfunktion . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 4.2 Val och implementering av känslighetsanalysmetoder . . . . . . . . . . . . . . 12 5 Resultat 14 5.1 Monte Carlo-filtrering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 5.2 Morris metod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 5.3 Sobols index och FAST . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 5.4 Grafisk känslighetsanalys . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 5.5 Relevans för Captario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 6 Diskussion 18 6.1 Hantering av modeller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 6.2 Applicering i SUM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 6.3 Framtida frågeställningar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 7 Referenser 21 Bilagor 23 A Morris metod på modell med Bernoullivariabler 23 B Morris Implementeringsalgoritm 24 C Modellbeskrivning av Accelerated 25 C.1 Variabelförteckning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 C.2 Fasernas kostnader och tider . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 C.3 Marknadsmodellen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 C.4 Andra projektvärdesstorheter . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 1 Bakgrund 1.1 Beslutsfattning inom läkemedelsutveckling Läkemedelsutveckling innefattar hela processen från upptäckten av en läkemedelskandidat till lansering på marknaden. I processen sker ett flertal aktiviteter, bland annat toxiska studier, kliniska studier och rekrytering av försökspersoner. Processen från den första kliniska studien till lansering på marknaden brukar delas upp i fyra faser. Detta illustreras i figur 1. Faserna kan karaktäriseras av tester på helt friska individer, tester på ett fåtal patienter som har sjukdomen mot vilken läkemedelskandidaten verkar, långvariga tester på större grupper som bär på sjukdomen och till sist registrering av läkemedelskandidaten [1]. Utvecklingsprocessen är i regel lång och kostsam. Det finns exempel på fall då det tagit mer än ett decennium för ett läkemedel att ta sig igenom alla faser med kostnader på flera miljarder kronor som följd [2][3]. Dessutom är det endast några få procent av alla läkemedelskandidaterna som i slutändan lanseras på marknaden. Detta beror på att en övervägande del av alla substanser misslyckas i någon aktivitet. Start Fas 1 Fas 2 Fas 3 Registrering Lansering Stopp Figur 1: Schematisk bild över de fyra faserna inom läkemedelsutveckling. Efter varje fas tas ett beslut om utvecklingen ska fortsätta eller om projektet ska läggas ner. Mot den bakgrunden måste läkemedelsföretagen ständigt fatta beslut om de skall fortsätta med kommande aktiviteter eller om risken att projektet inte blir lönsamt är för stor. Även om ett läkemedel skulle ta sig igenom alla faser finns det fortfarande andra faktorer som skulle kunna göra projektet olönsamt, en sådan faktor är tiden. Vilken tidpunkt som lanseringen sker i förhållande till eventuella konkurrenter är avgörande för lönsamheten på ett projekt. Lanserar konkurrenterna ett liknande läkemedel först riskerar man att förlora en stor del av marknadsandelarna vilket resulterar i lägre vinst (Magnus Ytterstad, Analytics Director, Captario). Under de senaste 20 åren har den totala kostnaden för att framställa och utveckla ett läke- medel tredubblats [3] och sannolikheten att nå marknaden är fortsatt låg med något negativ utveckling [2][3]. Med de fortsatta stora riskerna och den kostnadsmässigt negativa utveckling- en i branschen är det av stor vikt för läkemedelsföretag att kunna bygga trovärdiga modeller över läkemedelsprojekt och göra pålitliga riskanalyser av utvecklingsprocessen. Utan effek- tiv läkemedelsutveckling blir det ännu svårare och mer kostsamt för nya läkemedel att nå marknaden och sedermera patienter. Samtidigt tyder mycket på att behovet av läkemedel för kroniska sjukdomar kommer öka [3]. Dessa läkemedel är förhållandevis dyra, och därmed är frågan om effektiva läkemedelsprocesser inte bara viktiga för läkemedelsföretagen utan även samhället i stort. 1.2 Captario och programvaran SUM Captario är ett Göteborgsbaserat företag som grundades 2012. De hjälper läkemedelsut- vecklare att fatta strategiska beslut genom att modellera risker och uppskatta projektvär- desstorheter för ett läkemedelsprojekt. Captario har utvecklat ett nätbaserat verktyg, Cap- tario SUM, där användaren kan bygga en modell över utvecklingsprocessen av en läkeme- 1 delskandidat. I verktyget kan användaren specificera antaganden kring tid och kostnad för olika aktiviteter samt konstruera en marknadsmodell. Under dessa antaganden uppskattar SUM projektvärdesstorheter, exempelvis nettonuvärde och tid tills projektet är avslutat. Nettonuvärde, eller NPV från det engelska ordet net present value, är kanske den viktigaste projektvärdesstorheten och därför kommer vår analys av en modell från SUM i huvudsak kretsa kring denna. NPV definieras som nuvärdet av alla framtida intäkter och kostnader[4] vilket fungerar som ett mått på lönsamhet. Resterande projektvärdesstorheter som SUM be- räknar samt förklaringar till dessa finns i bilaga C. Sammanfattningsvis kan användaren få en prognos på hur lönsamt ett projekt kan bli, vilket sedermera kan användas som underlag vid beslutsfattning. Förutom en prognos över projektvärdesstorheter är det av intresse för användaren att få information om vad i modellen och ens antaganden som ger ett lönsamt projekt. Mer specifikt, eftersom det mesta i modellen varierar efter antaganden från använ- daren, är det intressant att ta reda på vilka osäkerheter från dessa antaganden som har stor påverkan på responsen från en modell. 1.3 Känslighetsanalys Matematiska modeller används, bland annat inom läkemedelsbranschen, i allt större utsträck- ning för att simulera utfall när traditionella experiment är för dyra eller omöjliga att genom- föra [5]. Modeller kan vara väldigt komplexa och bero på många variabler där fördelningarna oftast är okända och måste därför uppskattas. För att förstå hur dessa antaganden påverkar ens resultat, krävs att en så kallad känslighetsanalys utförs. Känslighetsanalys bygger på att analysera hur robusta modeller är under störning i indata. Metoderna inom detta fält delas in i två kategorier; lokal samt global känslighetsanalys, vilka förklaras ingående i kapitel 3. En definition av känslighetsanalys given av Andrea Saltelli, en av de ledande matematikerna inom området, lyder: ”The study of how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input.” [6] Den metod som Captario använder idag för att analysera vilka variabler som bidragit till osäkerhet i utdata är att göra ett så kallat tornadodiagram. Dessa har emellertid visat sig kräva många simuleringar för att ge tillförlitliga resultat. Dessutom kan denna metod inte påvisa interaktiva effekter mellan variabler. Vidare är det intressant för Captario att testa mer sofistikerade metoder inom känslighetsanalys, för att på så sätt kunna ge kunden bättre verktyg för att förstå sin konstruerade modell på ett bättre sätt. 2 Syfte Syftet med detta projekt är att hjälpa Captario att få en djupare förståelse över hur variatio- ner i indata påverkar de projektvärdesstorheter SUM returnerar. För att åstadkomma detta använder vi metodik inom känslighetsanalys för att identifiera de variabler som bidrar med mest känslighet i modellen. Detta kan i förlängningen användas som ett välgrundat underlag för Captarios kunder om vilka variabler som bidrar mest till variation i utvariablerna. Detta är dels intressant ur modelleringssynpunkt för att förstå hur osäkerheter i antagandena på- verkar osäkerheten i utvariabler, men också vilka variabler som behöver kontrolleras för att öka sannolikheten att projekt blir lönsamma. Vårt fokus ligger på att analysera känsligheten i projektvärdesstorheten NPV, och vår analys kommer begränsas till en specifik modell som speglar hur läkemedelsframtagning ser ut generellt. På så sätt testar vi vilka metoder som lämpar sig bäst för Captarios produkt. 2 3 Teori Känslighetsanalys är ett brett område och innefattar många olika metoder. De metoder som används i det här projektet, och presenteras i avsnitt 3.2, är bland de vanligaste inom områ- det. I avsnitt 3.1 beskrivs hur vi definierar en modell som känslighetsanalyser kan utföras på. I styckena nedan förklaras skillnaden mellan två typer av känslighetsanalyser; lokal och global. Lokal känslighetsanalys är den enklaste typen av känslighetsanalys. Det bygger på att en variabel Xi varieras med en förbestämd störning ∆ kring ett nominellt värde. Hur denna störning påverkar responsen undersöks genom att olika partiella derivator uppskattas i det nominella värdet. Metoderna ger information om vilka variabler som har en linjär effekt på responsen kring det nominella värdet [7]. Om modellen består av ett större antal variabler så är fördelarna med en lokal känslighetsanalys dess tidseffektivitet och kan därmed behöva relativt få beräkningar. Nackdelarna är däremot att endast förändringar kring ett nominellt värde utforskas och att eventuella interaktionseffekter mellan variabler förbises. För att kunna upptäcka interaktionseffekter kan man använda global känslighetsanalys. Ge- mensamt för dessa metoder är att de undersöker störningar i hela indatarummet och ger på så sätt information kring hur känsligheten i responsen kan förklaras i form av osäkerheten i modellens ingående variabler [8]. I förlängningen ger detta mer detaljerad information över vilka variabler som bidrar mest till känsligheten i responsen. Begränsningarna med global känslighetsanalys ligger i att metoderna kan vara väldigt beräkningsmässigt tidskrävande, framförallt om antalet variabler är stort [8]. 3.1 Notation för allmänna modeller Metoder för att utföra en känslighetsanalys beskrivs i avsnitt 3.2 och där används en allmän modell som exempel. I det här avsnittet beskrivs den notation som har valts för att beskriva en sådan modell. När en modell analyseras antas det att en deterministisk funktion existerar och är tillgänglig, vilken evaluerar responsen av intresse Y för en given vektor X av variabler. Alltså gäller det att Y = h(X) (3.1) därX = (X1, X2, . . . , Xp) ∈ Rp är vektorn innehållande alla variabler. Denna vektor kommer vara stokastisk med någon viss fördelning. För modeller i SUM kan det utan inskränkning antas att alla variabler är oberoende. Av bekvämlighetsskäl begränsas analysen till det fall då utfallsrummet av projektets variabler är hyperkuben Kp = [0, 1]p och slumpvariablerna är likformigt fördelade över denna mängd. Detta är ingen inskränkning, ty Y kan omskrivas till Y = (h ◦ g)(Z) (3.2) där Z ∼ unif([0, 1]p) och är oberoende samt gi(Z) = F−1Xi (Zi) 1. I fortsättningen kommer vi låta f = h ◦ g. 3.2 Metoder inom känslighetsanalys De metoder som presenteras i det här avsnittet är uteslutande globala känslighetsanalysme- toder. Varför just dessa metoder valts ut motiveras i avsnitt 4.2. 3.2.1 One-at-a-time-experiment One-at-a-time-expermient, eller OAT-experiment, är en kategori av metoder inom känslig- hetsanalys. I dessa metoder undersöker man vilken effekt man får på responsen Y när man varierar någon variabel Xi kring ett värde medan man håller de resterande variablerna fixe- rade [9]. Detta görs i tur och ordning för alla variabler där praxis i OAT-experiment är att 1I det diskreta fallet ska F−1 Xi tolkas som den generaliserade inversen. 3 variera Xi i ett intervall kring dess nominella- eller standardvärde [10]. Den låga beräknings- kostnaden är den stora fördelen med OAT-experiment men de tar dessvärre inte hänsyn till eventuella interaktionseffekter [9]. 3.2.2 Morris metod Morris metod [11] är en typ av OAT-experiment, vilken kan betraktas som en generalise- ring av de lokala känslighetsmetoderna med partiella derivator till en metod som besvarar vilka variabler som bidrar till den globala känsligheten. Metoden är förhållandevis beräk- ningsmässigt effektiv, samtidigt som den identifierar de variabler som bidrar till den globala känsligheten. Dessutom ger metoden en indikation av i vilka variabler responsfunktionen är icke-linjär och/eller vilka variabler som interagerar med andra. I och med att Morris metod tillhör gruppen OAT-experiment, kan inte interaktiva effekter mellan variabler hittas med hjälp av denna metod. Däremot säger metoden vilka variabler som inte interagerar, eller interagerar svagt med andra. Metoden bygger på att storheter associerade till fördelningen av den elementära effekten för den i:te variabeln beräknas. Den elementära effekten defineras – för något utfall Z enligt (3.3) Ei(Z) = f(Z + ∆ei)− f(Z) ∆ (3.3) där ∆ är en förbestämd fixerad störning och ei är den i:te kanoniska basvektorn. Vida- re är (3.3) en numerisk uppskattning av ∂f ∂Zi . Fördelningen av Ei ger information om den globala känsligheten av den i:te variabeln. Om m stycken simuleringar av Z genomförts (Z1, Z2 . . .Zm) kan följande storheter, (3.4) och (3.5), beräknas µ∗i = 1 m m∑ j=1 |Ei(Zj)| (3.4) σi = √√√√ 1 m m∑ j=1 (Ei(Zj)− Ēi)2 (3.5) där Ēi är medelvärdet av de elementära effekterna för variabel i. Morris [11] föreslog att µi = Ēi ska betraktas istället för µ∗i . Dock har µ∗i blivit allt mer vanligt förekommande i vetenskapliga sammanhang under senare år. Bland annat introduceras denna av Compolon- go [12]. Vi har mot den bakgrunden valt att studera µ∗i istället för µi. Om µ∗i är relativt stor innebär det att en störning i den i:te variabeln fortplantar sig, ty µ∗i uppskattar den genomsnittliga absoluta avvikelsen under en störning av storlek ∆. σi, standardavvikelsen av den i:te elementära effekten, ger en indikation om den i:te variabeln interagerar med andra, eller om responsfunktionen f är icke-linjär i just denna variabel. Vidare bör denna metod användas för att identifiera variabler som är okänsliga och inte interagerar med andra. Ofta illustreras dessa index i en graf med µ∗ på ena axeln och σ på andra. På så sätt kan använ- daren på ett kvalitativt sätt identifiera vilka variabler som är mer känsliga än andra, genom att analysera vilka punkter (µ∗i , σi) som är längst bort ifrån origo. En implementeringsmetod som föreslogs av Morris bygger på slumpmässiga banor i hyperku- ben [0, 1]p. En närmare beskrivning av denna finns beskriven i appendix B. Denna metod har emellertid under senare år förfinats av bland andra Compolongo [12], och bör vidare endast betraktas som en skiss över hur Morris metod implementeras på ett effektivt sätt. Morris implementeringsmetod [11] bygger på att hyperkuben diskretiseras likformigt med n steg i varje dimension. I detta diskretiseringsnät evalueras f i varje nod, vilka sedermera används i metoden. På så sätt kommer responsen inte beräknas fler gånger än nödvändigt, vilket är beräkningsmässigt fördelaktigt om f är tidskrävande att evaluera. 4 3.2.3 Monte Carlo-filtrering Monte Carlo-filtrering (MCF) är en metod för att kartlägga hur variabler är fördelade i oli- ka områden av responsen Y . Metoden går ut på att man delar upp simuleringsdatan i två disjunkta delmängder: acceptabla, G, och icke-acceptabla, B, efter förutbestämda villkor satt på utfallet av responsen. Exempelvis kan vissa utfall klassas som bra om responsen överstiger ett visst tröskelvärde ξ, vilket sedermera utgör villkoret vid indelning. MCF används fram- förallt inom diagnostisering av modeller för att ta reda på vilka variabler som bidrar till ett visst utfall för responsen [9]. För att implementera metoden genererasN stycken realiseringar av variablerX1, X2 . . .XN , vilka sedan används för att beräkna N värden av responsen Yi = h(Xi), i = 1, 2, ..., N . Där- efter delas observationerna upp i respektive delmängd utifrån det förutbestämda värdet på ξ. G ={(Xi, Yi); i = 1, ..., N, Yi > ξ} B ={(Xi, Yi); i = 1, ..., N, Yi ≤ ξ} (3.6) För varje variabel i modellen kan nu två stycken empiriska fördelningar fås, F̂i(Xi|G) och F̂i(Xi|B), där Xi|G ∈ G och Xi|B ∈ B. Ifall de två empiriska fördelningarna avviker mycket från varandra kan slutsatsen dras att variabeln Xi har stor påverkan på Y [8]. För att mäta denna avvikelse kan ett Kolmogorov-Smirnovtest (KS-test) tillämpas som testar ifall två stickprov är dragna från samma underliggande fördelning. Teststatistikan är det största avståndet mellan två fördelningar och definieras enligt di = sup Xi | F̂i(Xi|G)− F̂i(Xi|B) |. (3.7) Hypotesen att Xi|G har samma fördelning som Xi|B förkastas med signifikansnivå α om di > √ −1 2 ln ( α 2 ) · √ n+m nm (3.8) där n och m är lika med antalet Xi i mängden G respektive B. MCF är ett bra sätt att kvantifiera vilka variabler i modellen som har störst påverkan på Y . Metoden ihop har många egenskaper gemensamt med global känslighetsanalys däribland att hela indatarummet tas i beaktning samt att alla variabler varieras samtidigt. Dock så tas ingen hänsyn till interaktionseffekter överhuvudtaget [13]. Vilka variabler som kommer ge utslag i ett KS-test kommer att vara högst beroende på vad man sätter tröskelvärdet till. Det bör även poängteras att om en variabel Xi har någon inverkan på Y , kommer ett KS-test att påfinna denna om stickprovsstorleken är tillräckligt stor. Därav bör man studera graferna av de empiriska fördelningarna för att undersöka om det finns en praktisk signifikant avvikelse. 3.2.4 Sobols index Den mest centrala känslighetsanalysmetoden som vi valt att betrakta i denna rapport bygger på variansdekomposition. Detta innebär att variansen av responsvariablen additativt delas upp i funktioner som är beroende av färre variabler. Dessa termer utnyttjas sedan för att uppskatta det bidrag variabler och interaktioner mellan variabler ger till den totala variansen i responsen. Låt f vara responsfunktionen till modellen definierad på samma vis som i avsnitt 3.2. Vi antar dessutom att f är kvadratiskt integrerbar över hyperkuben Kp. Under detta antagan- det följer det att f kan uttryckas som en summa av kvadratiskt integrerbara funktioner av färre variabler [15]: f(X) = f0 + p∑ i=1 fi(Xi) + ∑ i ∑ j>i fij(Xi, Xj) + · · ·+ f12...p(X) (3.9) 5 Denna expansion innehåller totalt 2p−1 stycken summor och är unik under villkoret i ekvation (3.9) [15]: ∫ 1 0 fi1...is(Xi1 , ..., Xis) dXik = 0, 1 ≤ k ≤ s, {i1, ..., is} ⊆ {1, ..., p} (3.10) Det följer från ekvation (3.10) att termerna i ekvation (3.9) är ortogonala och kan uttryckas som integraler av f(X) genom [15]∫ f(X) dX = f0∫ f(X) ∏ k 6=i dXk = f0 + fi(Xi)∫ f(X) ∏ k 6=i,j dXk = f0 + fi(Xi) + fj(Xj) + fij(Xi, Xj) (3.11) där varje integral ska tolkas som en bestämd multipelintegral över en hyperkub med lämplig dimension. En följd av ekvation (3.11) är att f0 är konstant lika med väntevärdet av f(X). Genom att kvadrera båda leden i ekvation (3.9) och integrera över Kp så följer det att∫ f2(X) dX − f20 = n∑ s=1 n∑ i1<...i Vij(Y ) + · · ·+ V12...p(Y ), (3.13) Vij(Y ) = V (fij(Xi, Xj)) = VXiXj (EX∼ij (Y |XiXj))− VXi (EX∼i (Y |Xi))− VXj (EX∼j (Y |Xj)) = ∫ f2ij dXi dXj . (3.14) Vi(Y ) kallas första ordningens varians, vilken är den förväntade minskningen av variansen om vi håller Xi fixerat. Om vi dividerar båda leden i ekvationen ovan med V (Y ) så får vi att 1 = p∑ i=1 Vi(Y ) V (Y ) + ∑ i ∑ j>i Vij(Y ) V (Y ) + · · ·+ V12...p(Y ) V (Y ) . (3.15) Sobols känslighetsindex definieras nu på följande vis Si = V (E[Y |Xi]) V (Y ) Sij = V (E[Y |Xi, Xj ])− Vi − Vj V (Y ) ... STi = 1− V (E[Y |X∼i]) V (Y ) , (3.16) där X∼i = (X1, X2, . . . , Xi−1, Xi+1, . . . , Xp). Si är Sobols index för första ordningens effekt och ger information om hur stor del av variansen i responsfunktionen som kan tillskrivas den 6 i:te variabeln isolerad från de andra variablerna. Sij är index för första ordningens interak- tionseffekt mellan den i:te och den j:te variabeln, och i en modell med p stycken variabler finns ∑p i=2 ( p i ) möjliga interaktionseffekter. Vidare kommer det för ett stort p finnas ett stort antal interaktionstermer, och av beräkningsmässiga skäl uppskattas vanligen endast index av låg ordning. Till följd av ekvation (3.15) kommer summan av alla känslighetsindex vara 1. STi är index för den totala effekten, vilken är den förväntade variansen om alla variabler utom Xi hålls fixerade och kommer att innehålla information om första ordningens effekt samt alla möjliga interaktionseffekter. Detta index ger information om hur influerande den i:te variabeln är i modellen genom att samla alla interaktionseffeter i en storhet. Genom att slå samman alla interaktionseffekter på detta sätt till en term gör att man inte förlorar någon information då man empiriskt har sett att interaktioner mellan tre variabler eller fler oftast är obefintlig (Kerstin Wiklander, Uni.lektor Chalmers). Vidare uppskattas i regel endast första ordningens samt totala variansen för att hålla nere beräkningstiden, ty dessa två index ger god information om känsligheten i modellen. Variansen som används vid beräkning av första ordningens index kan omskrivas [15]: Vi = V (E[Y |Xi]) = ∫ f(X)f(Xi,X ′ ∼i) dX dX ′∼i − f20 (3.17) där Xi är det i:te elementet i X, X∼i är en vektor innehållandes alla variabler utom den i:te och (Xi,X∼i) utgör en komplett parameteruppsättning. Sobol visade också att den totala variansen kan omskrivas: V∼i = V (E[Y |X∼i]) = V (Y )− 1 2 ∫ ( f(X)− f(X ′i,X∼i) )2 dX dX ′i. (3.18) Integralerna i (3.17) och (3.18) uppskattas lämpligen med hjälp av Monte Carlo-simuleringar. Om x(k) ∈ Rp, y(k) ∈ Rp−1 och z(k) ∈ R (k = 1, 2 . . . N) är simulerade vektorer av oberoende likformigt fördelade slumpvariabler, skattas variansena med följande summor: f̂0 = 1 N N∑ k=1 f(x(k)) V̂ = 1 N N∑ k=1 f2(x(k))− f̂0 2 V̂i = 1 N ( N∑ k=1 f(x(k))f(x (k) i ,y(k)) ) − f̂0 2 V̂∼i = V̂ − 1 2N N∑ k=1 ( f(x(k))− f(z(k),x (k) ∼i ) )2 . (3.19) Uppskattade sobolindex Ŝi och ŜTi erhålls om varianserna i (3.17) byts mot de skattade i (3.19). Till följd av att felet i Monte Carlo-simuleringar avtar som N−1/2 fordras ett stort N för att få bra resultat. Mer effektiva skattningar har emellertid utvecklats under de senare åren, vilka presenteras av Saltelli i Computer Physics Communication [14]. 3.2.5 Fourier Amplitude Sensitivity Test – FAST De integraler som behöver beräknas för Sobols index uppskattas vanligen med Monte Carlo- metoder, genom simuleringar av variabler som antas ligga i hyperkuben Kp. Ett alternativt tillvägagångssätt, som först presenterades av Cukier [16] är att generera slumptal längs en kurva som fyller hyperkuben likformigt. Alltså, antas att alla variabler kan uttryckas xi = ri(sinωis) för någon funktion ri och vinkelhastighet ωi, och s tillåts variera i intervallet [−π, π]. Då kan f0 approximeras som 7 f0 = E(Y ) = ∫ f(X) dX ≈ 1 2π ∫ π −π f(r(s)) ds (3.20) där r(s) = (r1(sinω1s), r2(sinω2s) . . . rp(sinωps)). Denna approximation gäller endast om ri och ωi är väl specificerade, därbland att r utrymmer hyperkuben likformigt, och att vinkelfrekvenserna är linjärt oberoende i heltalsmening, det vill säga att ∑p i=1 λiωi = 0 saknar heltalslösningar i λi [17]. Vidare kan variansen av responsen approximeras på följande sätt: V (Y ) = ∫ f2(X) dX − (E(Y ))2 ≈ 1 2π ∫ π −π f2(r(s)) ds− ( 1 2π ∫ π −π f(r(s)) ds )2 . (3.21) Den första termen i högerledet kan enligt Parseval’s formel uttryckas genom 1 2π ∫ π −π f2(r(s)) ds = ∞∑ k=−∞ |ck|2 (3.22) där ck är den k:te Fourierseriekoefficienten till f ◦ r: ck = 1 2π ∫ π −π f(r(s))e−jks ds, j2 = −1. (3.23) I och med att c−k = c̄k =⇒ |c−k| = |ck|, och att sista termen i högerledet av (3.21) är |c0|2, kan (3.21) omskrivas: V (Y ) ≈ 2 ∞∑ k=1 |ck|2 = D. (3.24) Dessutom, om ωi, för i = 1, 2 . . . p, är heltal, kan V (E(Y |Xi)) uppskattas genom Vi(Y ) ≈ 2 ∞∑ k=1 |ckωi |2 = Di. (3.25) I praktiken måste dessa serier trunkeras, men på grund av konvergens hos fourierserier, kan detta med fördel göras så länge tillräckligt många termer används. Heuristik visar att M ≥ 4 termer brukar vanligtvis vara tillräckligt [17]. Dessutom kan ck effektivt beräknas med hjälp av den diskreta fouriertransformen eller FFT, vilket innebär att det finns beräkningsmässiga fördelar med att tillämpa FAST. Dessutom har det visats att SFAST i = Di/D är ekvivalent med första ordningens Sobol index. Saltelli föreslog att åtminstone Ns = 2Mωmax + 1 simu- leringar måste genomföras för att erhålla ett tillförlitligt resultat. För att uppskatta den totala effekten för den i:te variabeln tilldelas den, i analogi med uppskattningen av första ordningens index, en frekvens ωi. Till skillnad från tidigare tilldelas alla andra variabler en och samma vinkelfrekvens ω∼i. Vidare kan man med hjälp av detta beräkna den varians ”alla utom i” ger upphov till genom V∼i ≈ 2 ∞∑ k=1 |ckω∼i |2 = D∼i (3.26) vilken används för att uppskatta totala effekten: SFAST Ti = 1− D∼i D . (3.27) Saltelli [17] föreslog att kurvor på formen ri(·) = 1 2 + 1 π arcsin(·) kan användas som utfyllande banor i hyperkuben. Dessa har visat sig ge förhållandevis rättvisa resultat men nackdelen med denna kurva är att den alltid kommer börja i samma punkt. För att undvika potentiell sys- tematisk avvikelse, kan istället kurvan xi = ri(sin(ωis+ φi)) användas, där φi ∼ unif(0, 2π). Detta bör emellertid upprepas med Nr realiseringar av (φ1, φ2, . . . , φr), för att få ner den 8 systematiska avvikelsen φi ger upphov till. Då fås indexen med hjälp av aritmetiska medel- värden SFAST i = ∑Nr j=1D (j) i∑Nr j=1D (j) (3.28) där D(j) i och D(j) är de uppskattade varianserna associerade till den j:te realiseringen av (φ1, φ2, . . . , φr). Vidare kommer totala antalet simuleringar med denna kurva bli (2Mωmax + 1)Nr. 9 4 Metod En användare i SUM bygger sin egen modell och gör antaganden kring sitt läkemedelsprojekt och därför finns det, i teorin, oändligt många unika modeller. För att utföra en känslighets- analys på läkemedelsprojekt begränsades undersökningen till en modell som Captario byggt som kallas Accelerated. Den ansågs vara representativ för de flesta läkemedelsprojekt och tar vara på den sekvensiella strukturen av ett projekt som finns beskriven i avsnitt 1.1. Model- lerna i SUM uppskattar flera projektvärdesstorheter men vi valde att bara fokusera på den mest centrala av dem; NPV. Andra projektvärdesstorheter kan analyseras helt analogt. De flesta metoder som har valts för det här projektet, och som finns beskrivna i avsnitt 3.2, kan inte direkt använda sig av simuleringsdata från SUM. Metoderna behöver en funktion för Accelerated som givet en uppsättning variabler kan evaluera NPV. Evalueringsfunktionen som SUM använder var inte tillgänglig under projektet så vi beslutade att skriva en funktion i MATLAB som imiterar beräkningen. Funktionen återfinns i [18] tillsammans med all källkod för projektet. Modeller i SUM kan vara komplexa och imitationen kräver grundlig förståelse för hur de är uppbyggda. Accelerated är en icke-linjär modell med 40 invariabler varav sju stycken är Bernoullifördelade. Fördelningarna av de resterande variablerna är specificerade i bilaga C, vilken innehåller en mer detaljerad beskrivning av accelerated. Vid beräkning av NPV nuvärdesjusteras dessutom alla kostnader över tid. Idéen bakom modeller i SUM finns beskriven i avsnitt 4.1. Med en funktion på plats kunde ett antal olika metoder för känslighetsanalys appliceras. Valet av dessa metoder och hur de implementerades motiveras i avsnitt 4.2 och resultaten av implementeringarna finns i avsnitt 5. De förenklingar och vägval som gjorts under projektet är sammanställda i figur 2. Captario SUM Accelerated NPV(X) Känslighetsanalys- metoder Resultat Bygg representativ modell Imitation MATLAB Implementering Figur 2: Schematisk figur över problemlössningsstrategin i projektet. I programvaran SUM byggs den representativa modellen Accelerated, vars beräkning av NPV imiteras i MATLAB. Känslighets- analysmetoderna implementeras sedan genom att anropa MATLAB-funktionen, och resultat erhålls. 4.1 Evaluering av NPV för modeller i SUM För att kunna konstruera en funktion från modellen Accelerated behövde vi först förstå i allmänhet hur modeller i SUM fungerar och dessutom hur vi skulle hantera de Bernoullivari- abler som finns. I vårt fall, när en projektvärdesstorhet i ett läkemedelsprojekt är av intresse att modellera, kommer variablerna i modellen vara händelser associerade till bland annat kliniska studier samt händelser i marknaden som påverkar försäljningsvolymen. Exempelvis kan elementen i X vara kostnader, tider för olika kliniska studier och lanseringstider för kon- 10 kurrenter. För NPV, kan responsfunktionen för en allmän projektmodell explicit uttryckas enligt (4.1) Y = h1(X) = −C1 − I1C2 − I1I2C3 − I1I2I3Creg. + I1I2I3I4(R− Cm) (4.1) där Ci och R = R är den nuvärdesjusterade kostnaden för den i:te fasen respektive inkomsten givet variablerna i X. Dessutom är Ii ∼ Bernoulli(pi) och antas vara oberoende, vilken antar värdet 1 om läkemedelskandidaten tar sig igenom den i:te fasen och 0 annars. Vidare om Ii = 0 för något i kommer avkastningen bli negativ. Ett naturligt tillvägagångssätt att uppskatta eNPV = E(Y ) är genom att simulera att stort antal utfall X1,X2 . . .XN och beräkna medelvärdet êNPV1 = 1 N N∑ i=1 h1(Xi). (4.2) Ett problem med skattningen (4.2) är att en stor del simuleringsdata inte används till följd av Bernoullivariablerna i (4.1). Exempelvis, om I1 = 0 kommer alla variabler som inte finns i X(1) att simuleras, men inte användas. En alternativ punktskattning som kringgår detta utnyttjar antagandet att Bernoullivariablerna är oberoende eNPV = E[h1(X)] = E [ − C1 − I1C2 − I1I2C3 − I1I2I3Creg. + I1I2I3I4(R− Cm) ] = = −E(C1)−p1E(C2)−p1p2E(C3)−p1p2p3E(C4)+p1p2p3p4E(R−Cm) = E[h2(X)] (4.3) där h2 defineras som h2(X) = −C1 − p1C2 − p1p2C3 − p1p2p3Creg. + p1p2p3p4(R− Cm). (4.4) Vidare kan eNPV uppskattas med följande punktskattning, vilken använder all simulerings- data: êNPV2 = 1 N N∑ i=1 h2(Xi). (4.5) 4.1.1 Marknadsmodell I många marknadsmodeller kan det vara intressant att modellera konkurrenternas anspråk på de tillgängliga marknadsandelarna. Vår analys begränsas till att en konkurrents beteende enbart bestäms av huruvida de lyckas lansera produkten överhuvudtaget, samt konkurren- tens lanseringstid. Det förstnämnda modelleras med en oberoende Bernoullivariabel, och det senare med en oberoende kontinuerlig slumpvariabel med specificerad fördelning. Mot samma bakgrund som tidigare, är det behändigt att få bort Bernoullivariablerna ur modellen, ty en simulering av en lanseringstid kommer inte till användning om konkurrenten i fråga inte lanserar. I allmänhet kommer ordningen aktörerna lanserar att ha betydelse för hur många marknadsandelar som vi erhåller och därför är det ofta svårt att formulera den årliga avkastningen RA(t) som en additiv modell, liksom (4.1). Vidare måste andra verktyg tillämpas. Vi låter lanseringstiden för projektet betecknas TRE , vektorn med alla Bernoullivariabler för alla konkurrenter I = (Ic1 , Ic2 , . . . , Icn) och K alla möjliga utfall av I. Då gäller att E[RA(t)|TRe] = ∑ κ∈K E((RA(t)|TRe)|I = κ)P (I = κ) = = ∑ κ∈K ( E((RA(t)|TRe)|I = κ) n∏ i=1 pκi i (1− pi)1−κi ) (4.6) 11 där κi är den i:te komponenten i κ och pi är parametern för bernoullivariabeln som svarar mot den i:te konkurrenten. Om intäkterna diskonteras årsvis med räntan r gäller att R|TRE = ∫ T TRE RA(t)r−([t]+1) dt (4.6) =⇒ (4.7) E(R|TRE) = ∫ T TRE (∑ κ∈K pκi i (1− pi)1−κiE((RA(t)|TRE)|I = κ)r−([t]+1) ) dt (4.8) där T är den tid då projektet avslutas och [t] är heltalsdelen av t. Används (4.8) vid upp- skattningen av R i (4.4) kommer det inte att finnas några Bernoullivariabler kvar i modellen över huvudtaget, vilket utnyttjas nedan. 4.1.2 Val av responsfunktion Det kanske mest naturliga valet av modell är Y = h1(X), ty simuleringarna av X kommer efterlikna ett verkligt scenario, där någon aktivitet antingen misslyckas, eller alla studier påvisar positiva effekter. Om h1 (4.1) studeras i en känslighetsanalys, så kommer i allmänhet Bernoullivariablerna bidra mest till osäkerhet i utdata. Detta på grund av att dessa variabler ger upphov till diskontinuitet i funktionen h1 ◦ g, vilket är definierat i analogi med ekvation (3.2). Av den anledningen har vi valt att genomföra en känslighetsanalys på responsfunktionen h2 definerad i ekvation (4.4), i vilken samtliga Bernoullivariabler är borttagna. 4.2 Val och implementering av känslighetsanalysmetoder I det här projektet implementerade vi metoderna Morris metod, Monte Carlo-filtrering och Sobols index samt genomförde en grafisk känslighetsanalys. Den mest centrala av dessa me- toder är Sobols index eftersom den ger oss både information från hela indata-rummet samt kvantifierar interaktionseffekter. Vår modell består av många variabler och indatan är alltid stokastisk och på grund av komplexiteten i modellen kan vi anta att det finns många inter- aktionseffekter. Nackdelen med att använda Sobol är att den är beräkningsmässigt krävande och för att behandla detta valde vi att avskärma insignifikanta variabler. Att avskärma variabler innebär att fixera okänsliga variabler, oftast kring dess nominella värde. De variabler som avskärmades fixerades till dess väntevärde. Den vanligaste metoden för detta är Morris metod som dessutom kvalificerar som en typ av global känslighetsana- lysmetod. Fördelarna med metoden i vårt fall är att vi både får ett mått på vilken effekt variablerna utgör på NPV, men även om de kan ingå i någon form av interaktionseffekter. Vid implementeringen utfördes transformen (3.2) till hyperkuben, innan en algoritm från en artikel i tidningen Acta Astronautica [19] tillämpades. Monte Carlo-filtrering skulle potentiellt kunna användas som en avskärmningsmetod men metoden förbiser helt interaktionseffekter. Det finns dock fördelar med MCF; den undersöker hela indata-rummet och den är enklare att genomföra eftersom realisationsdatan kan använ- das direkt från simuleringar i SUM. I det här projektet användes MCF som ett komplement till Morris metod för att jämföra resultaten. Vid implementering av MCF gjordes 10 0000 realiseringar för att få tillförlitliga resultat och tröskelvärdet sattes till eNPV från vår modell. Även om det går att använda realiseringar direkt från SUM användes realiseringar från vår MATLAB-funktion för att få helt konsekventa resultat. För att avgöra om en variabel var signifikant eller inte användes KS-testet. Eftersom resultaten i MCF är så beroende av vilket värde tröskelvärdet har, gjordes ett test där KS-testet genomfördes för olika tröskelvärden som varierade likformigt från 2/3 eNPV till 4/3 eNPV med 100 värden. Efter att flera variabler avskärmats med hjälp av Morris metod kunde Sobols index och So- bols totala index räknas ut under en väsentligt kortare beräkningstid. Vid implementeringen användes ett MATLAB-program som använder skattningarna presenterade i [14]. Det beräk- ningsmässiga felet uppskattdes med hjälp av normal approximation med konfidensgrad 0,95. 12 Ett alternativt sätt att räkna ut Sobols index är att använda metoden FAST vilket innebär mindre beräkningskostnader men potentiellt mindre pålitliga resultat. I vår implementering- en av FAST skrev vi ett eget MATLAB-program som använde MATLAB:s funktion FFT. Vi valde att genomföra båda metoderna för att räkna ut Sobols index av första ordningen och sedan jämföra om FAST kan vara ett lika tillförlitligt alternativ. Som komplement till Sobols index valde vi att genomföra en grafisk analys, som tillsam- mans med skalärerna från Sobols index kunde ge en ännu tydligare bild av hur en variabel påverkar responsen. Vi utförde detta genom att beräkna hur eNPV ändrades när vi fixerade en variabel Xi med ett antal punkter från dess teoretiska fördelning. Sedan illustrerade vi detta grafiskt med E[Y |Xi] − E[Y ] och Xi. Att utföra detta på alla variabler i modellen är beräkningsmässigt tidskrävande. Men eftersom vi redan hade avskärmat bort ett antal variabler med Morris metod kunde vi utföra detta på de variabler som var kvar. 13 5 Resultat I avsnitten nedan presenteras de figurer och tabeller som är relevanta för vårt projekt enligt den typ av modell som finns beskriven i avsnitt 4.1. I 4.1 motiveras även varför känslighetsana- lysen har genomförts på modellen Accelerated utan Bernoullivariabler. Några av metoderna utfördes på modellen med Bernoullivariabler bara för att visa på att de dominerar känslig- heten. Dessa resultat presenteras i bilaga A. Samtliga metoder under känslighetsanalysen identifierade POC_RecruitmentRate som den dominerande variabeln. Dess värde beskriver rekryteringsfarten per klinik per månad för de tre kliniker som antas vara med och rekrytera patienter till POC-studien i fas 2. I bilaga C finns en mer detaljerad beskrivning av hur POC_RecruitmentRate samverkar i Accelerated. 5.1 Monte Carlo-filtrering MCF utfördes på simuleringsdata innehållande 10 000 realiseringar av variablerna och be- räkningar av NPV där tröskelvärdet ξ är satt till eNPV. Resultaten från MCF visas i figur 3. x 0.5 0.6 0.7 0.8 0.9 F (x ) 0 0.5 1 cadTime x 0 2 4 6 8 F (x ) 0 0.5 1 POC RecruitmentRate x 0.5 1 1.5 2 F (x ) 0 0.5 1 Ph3StudyTime x 0.9 1 1.1 1.2 F (x ) 0 0.5 1 RegTime x ×10 -3 2.4 2.6 2.8 3 3.2 F (x ) 0 0.5 1 prevalence x 0 2 4 6 8 F (x ) 0 0.5 1 RampTimeInYears NPV<ξ NPV>ξ Figur 3: De empiriska fördelningarna F (x) för de variabler som fick ett signifikant utfall med KS- testet där signifikansnivån α sattes till 5 % och tröskelvärdet ξ sattes till eNPV. Den blå kurvan är den empiriska fördelningen för variabeln när responsen var lägre än eNPV och den orangea när utfallet på responsen var större än eNPV. Av modellens 33 variabler så är det fem stycken som ger ett signifikant resultat med KS- testet. Utav dessa är POC_RecruitmentRate den variabel som visar på störst skillnad mellan de två olika empiriska fördelningarna. MCF utfördes även på simuleringsdata när ξ varierades i ett symmetriskt intervall kring eNPV. För varje unikt ξ inom intervallet genomfördes ett KS-test på alla variabler där resultatet av denna analys kan ses i tabell 1, i vilken andelen signifikanta test för en variabel anges. 14 Tabell 1: Tabellen illustrerar andelen av KS-testet visade på signifikant avvikelse i de empiriska fördelningarna, med signifikansnivå α satt till 5%. Testet utfördes när ξ varierade från 2/3 av eNPV till 4/3 av eNPV med 100 värden. Enbart de variabler som någon gång visade på ett signifikant utfall med KS-testet återfinns i tabellen. Variabelnamn ξ ≤ eNPV ξ > eNPV batchForPH3Time 0.12 0.60 POC_RecruitmentRate 1 1 Ph3StudyTime 1 1 RegTime 0.42 0.02 POC_CostPerCenter 0.14 0 Development_costs3 0.06 0.38 prevalence 1 1 RampTimeInYears 1 1 Comp1LaunchTime 0.12 0 Comp3LaunchTime 0.02 0 När tröskelvärdet varierar så börjar även en del andra variabler att påvisa skillnader i sina empiriska fördelningar. Men det återfinns fyra variabler som förblir signifikanta för varje test som görs, POC_RecruitmentRate , Ph3StudyTime, prevalence, RampTimeInYears. Något som även går att observera är att batchForPH3Time påvisar signifikant avvikelse i mer än hälften av testen när ξ var större än eNPV. 5.2 Morris metod I figur 4 illustreras resultaten från implementeringen av Morris metod. De variabler som återfinns tillräckligt nära origo kommer inte betraktas i beräkningen av Sobols index för att minska den beräkningsmässiga kostnaden. Valet av kritiskt område i figur 4 är en ellips kring origo. Observera att Morris metod och MCF enas på vissa variablers inverkan på NPV, och att Morris uteslutande identifierar tider. 0 5 10 15 20 25 30 35 0 20 40 60 12 3 4 5 10 11 6 26 27 µ∗ σ 1 cadtime 2 Time 3 3mToxTime 4 POC_SetupTime 5 POC_RecruitmentRate 6 POC_AnalysisTime 10 Ph3StudyTime 11 RegTime 26 prevalence 27 RampTimeInYears Figur 4: Diagram med absoluta medelvärde µ∗ och standardavvikelse σ för de elementära effek- terna associerade till alla variabler. Endast de variabler med relativt stora index är markerade med variabelnamn. 5.3 Sobols index och FAST Beräknade värden av Sobols index för första ordningens effekter visas i tabell 2. Resultatet är framräknat från 1 000 simuleringar. Värden för första ordningens effekter framräknade med FAST visas även i tabellen. De variabler som visas är de som Morris metod identifierade som mest inflytelserika. 15 Tabell 2: Beräknade värden på Sobols index för första ordningens effekter, Si med tillhörande felmarginaler, ∆Si, samt resultat från implementation av FAST, SFAST i . Variabelnamn Si ±∆Si SFAST i POC_RecruitmentRate 0, 9312± 0.0591 0, 9136 RampTimeInYears 0, 0158± 0.0093 0, 0163 prevalence 0, 0061± 0.0057 0, 0065 Ph3StudyTime 0, 0046± 0.0057 0, 0054 Time 0, 0016± 0.0018 0, 0005 3mToxTime −0, 0001± 0.0016 0, 0005 RegTime −0, 0012± 0.0028 0, 0003 POC_AnalysisTime 0, 0009± 0.002 0, 0001 POC_SetupTime 0, 0001± 0.0027 0, 0004 cadtime −0, 0009± 0.0027 0, 0010 POC_RecruitmentRate är förklarar cirka 93 procent av variansen i NPV för modellen utan Bernoullivariabler, och prevalence bidrar näst mest med känslighet. Övriga variabler har relativt liten inverkan. POC_RecruitmentRate var dominant även vid implementation av FAST, men den inbördes rangordningen skiljer sig något för övriga variabler. I tabell 3 visas motsvarande beräkningar för Sobols totala index, STi. Notera de stora beräkningsmässiga felmarignalerna. Tabell 3: Beräknade värden på Sobols totala index, STi, med tillhörande felmarginal, ∆STi, för de variabler som Morris metod identifierade som mest inflytelserika. Variabelnamn STi ±∆STi POC_RecruitmentRate 0, 9546± 0.0246 RampTimeInYears 0, 0246± 0.008 Ph3StudyTime 0, 0095± 0.0004 prevalence 0, 0091± 0.0033 3mToxTime 0, 0007± 0.001 Time 0, 0009± 0.001 RegTime 0, 0022± 0.0002 POC_AnalysisTime 0, 0011± 0.0001 POC_SetupTime 0, 0021± 0.002 cadtime 0, 0021± 0.002 5.4 Grafisk känslighetsanalys Som komplement till variablerna i Sobols index utfördes en grafisk analys. I figur 5 återfinns resultatet från denna. Notera att värden på POC_RecruitmentRate bidrar till en markant skillnad jämfört med ändringar för de andra variablerna. Det kan även ses i graferna att skill- naden i väntevärde av NPV är nästan konstant för variablerna ph3studyCosts och RegTime medan andra variabler i graferna har större inverkan på väntevärdet. 16 0.8 1 1.2 1.4 1.6 1.8 −10 −5 0 5 10 Ph3StudyTime 0.95 1 1.05 1.1 1.15 −4 −2 0 2 4 Regtime 70 80 90 100 −4 −2 0 2 4 Ph3StudyCosts 2.5 2.6 2.7 2.8 2.9 3 ·10−3 −10 −5 0 5 10 prevalence 2 3 4 5 6 7 −20 −10 0 10 20 RampTimeInYears 0.65 0.7 0.75 0.8 0.85 0.9 −2 −1 0 1 2 Time 0.4 0.5 0.6 0.7 0.8 0.9 −1 −0.5 0 0.5 1 3mToxTime 0 0.1 0.2 0.3 0.4 −4 −2 0 2 4 POC_SetupTime 0 1 2 3 4 5 −100 −50 0 50 100 POC_RecruitmentRate 0.08 0.1 0.12 0.14 0.16 −1 −0.5 0 0.5 1 POC_AnalysisTime Figur 5: I graferna ovan visas skillnaden E[NPV|Xi]−E[NPV] för ett stickprov av storlek 100 från variabelns fördelning. För varje värde av Xi har 10.000 beräkningar av NPV gjorts för att uppskatta E[NPV|Xi] 5.5 Relevans för Captario En del av frågeställningen i det här projektet är huruvida Captario kan använda de käns- lighetsanalysmetoder vi har använt. Rent praktiskt innebär detta att för att Captario skall kunna använda sig av metoderna bör de var effektiva. Den mest centrala metoden, Sobols index, har visat sig vara beräkningsmässigt tidskrävande och tog flera timmar att beräkna. Indexen beräknades betydligt snabbare, under en minut, när FAST implementerades. Den grafiska känslighetsanalysen tog även den flera timmar att beräkna. Morris metod och MCF beräknades snabbt. 17 6 Diskussion I samtliga presenterade resultat är POC_RecruitmentRate den mest influerande variabeln. Denna variabel är approximativt inverst proportionell till tiden för en aktivitet i den andra fasen och vidare kommer NPV vara högst olinjär i detta argument. Dessutom är ett antagan- de i modellen att samtliga tre deltagande kliniker i denna modell har samma rekryteringstakt. Vidare vore ett lämpligt sätt att få ned känsligheten i modellen att modellera rekryterings- takten på ett annat sätt, exempelvis genom att låta varje klinik ha en egen rekryteringstakt. En annan tydlig genomgående tendens är att variabler associerade till tider är mer påver- kande än de andra. Detta har att göra med att en sänkt projekttid har flera effekter på lönsamheten. Dels kommer en tidig lansering innebära att produkten säljs på marknaden under en längre tid, och vidare blir den totala försäljningsvolymen högre. Dessutom kommer en kortare projekttid göra det mer sannolikt att företaget lanserar läkemedlet före konkur- renterna, vilket resulterar i att en större marknadsandel erhålls. En tredje effekt är att en tidig lansering gör att de första intäkterna inkasseras närmare projektstarten, och om alla intäkter och kostnader nuvärdesjusteras kommer detta göra att eNPV ökar. Därmed kommer eNPV vara icke-linjär, eller svagt icke-linjär i dessa variabler, till skillnad från exempelvis kostnadsvariabler, som förväntas öka eNPV linjärt. Detta framgår tydligt i den grafiska käns- lighetsanalysen, men även i graferna som erhölls från Morris metod, där σ << 1 för samtliga kostnadsvariabler. Den variabel som avgör hur stor andel av populationen som insjuknar och behöver vår lä- kemedelskandidat,Prevalence, kommer trots dess låga varians ha ett väsentligt bidrag till känsligheten i modellen. Detta beror på att prevalence skalas med en faktor av storleksord- ning 108 på samtliga ställen där den används. Det bör poängteras att känsligheter bör hanteras på olika sätt av en modellskapare bero- ende på vad den känsliga variabeln i fråga representerar. Variabler såsom prevalence och POC_RecruitmentRate kan vara omöjliga att kontrollera. Vidare bör den som konstruerar modellen se till att de antagna fördelningarna av dessa variabler är så exakta som möjligt. Annars måste man överväga ett annat sätt att modellera vad variabeln i fråga representerar. Exempelvis, i fallet med POC_RecruitmentRate vore det fördelaktigt att dela upp denna hastighet, och låta varje klinik ha en egen rekryteringstakt. Däremot finns variabler som går att kontrollera rent praktiskt. De resultat som presenterats tyder på att det är viktigare att lansera tidigt än att få ned kostnaden för vardera studie. Alltså är det avgörande för att öka sannolikhet för vinst på ett läkemedelsprojekt, att planera projektet på ett sådant sätt att den förväntade lanseringstiden sker så tidigt som möjligt. 6.1 Hantering av modeller Genom att skapa vår egen modell, med Accelerated i SUM som utgångspunkt, har vi kun- nat välja hur vi bäst kan behandla de Bernoullivariabler som finns. Istället för att inkludera Bernoullivariablerna i modellen valde vi att integrera dess sannnolikheter i funktionen. Re- sultatet blev att vi tydligare kunde se vilka andra variabler som bidrog till känsligheten i modellen. Detta är möjligt tack vare antagandet att variablerna i Accelerated är oberoende. Detta kommer i regel att gälla, eftersom det enda som utnyttjas i avsnitt 4.1, ekvation (4.2) är att Bernoullivariablerna är oberoende av varandra. Dessutom att Bernoullivariablerna som avgör ifall projektet går vidare till nästa fas är oberoende av kostnaden i nästa fas. Alltså gäller att, även om modellen hade utvecklats så att Bernoullivariablerna är beroende av ak- tiviteter i fasen innan, hade vår uppdelning fortfarande varit giltig. Vi begränsade känslighetsanalysen till endast NPV som utvariabel. SUM returnerar fler pro- jektvärdesstorheter som inte behandlas alls i detta kandidatarbete. Identifieringen av de mest känsliga variablerna hade möjligtvis blivit annorlunda om vi hade valt att undersöka en an- nan projektvärdesstorhet. Trots detta kan samtliga metoder vi studerat användas, oberoende av projektvärdesstorhet, eftersom vår hantering av Bernoullivariabler är applicerbar på al- 18 la projektvärdesstorheter. Exempelvis om sannolikheten att göra förlust undersöktes skulle följande kunna utnyttjas: P (NPV < 0) = P (NPV < 0|I1I2I3I4 = 1)p+ P (NPV < 0|I1I2I3I4 = 0)︸ ︷︷ ︸ =1 (1− p) = = P (NPV < 0|I1I2I3I4 = 1)p+ (1− p) (6.1) där p = P (I1I2I3I4 = 1). Vidare, om (6.1) utnyttjas exkluderas alla Bernoullivariablerna ur modellen. Det bör även poängteras att en känslighetsanalys kan utföras på flera utvariabler samtidigt. Detta hade möjligen gett en bättre beskrivning av känsligheten i de ingående variablerna. Nackdelen med detta är naturligtvis att många fler beräkningar krävs. 6.2 Applicering i SUM Captario kan använda våra resultat för att förstå sina modeller bättre och till och med välja att bygga dem annorlunda. En känslighetsanalys i deras fall behöver inte bara vara för att göra prognoser utan även för att förbättra modeller. Ett exempel som nämndes tidigare är att POC_RecruitmentRate är den mest influerande variabeln. Detta beror troligtvis på att den används i alla tre kliniker i modellen, och om variablens fördelning skattas individuellt för varje klinik så kommer känsligheten i den att minska. Från våra resultat har vi kunnat hjälpa Captario att komma till insikter kring en modell i SUM, men en viktig fråga är om Captario, tillsammans med deras kunder, kan använda de känslighetsanalysmetoder vi har använt oss av i framtida modeller. Både Sobols index och Morris metod kräver evaluering av modellen i specificerade simulerade variabler. I exempel- vis Morris metod krävs detta för att beräkna varje elementär effekt. Idag har inte SUM den funktionen, och vidare måste Captario utveckla nya applikationer till SUM för att genom- föra nya simuleringar. Den metoden som vi ansåg var mest informativ – Sobol index – är för beräkingsmässigt tidskrävande för Captario om programmen som vi tillämpade används. Här finns utrymme för effektivisering av koden, och vid integrering av denna metod i SUM, skulle potentiellt de simuleringar kunna anpassas på ett sådant sätt att simuleringar kan återanvändas i flera metoder. Ett annat alternativ vore att använda FAST, vilket drar ned beräkningstiden väsentligt, framförallt för första ordningens Sobol index. Förutom att MCF kan användas med simuleringsdata direkt från SUM är metoden snabb och dessutom anpassningsbar. Genom att göra ändringar i SUM och justera tröskelvärdet som delar upp simuleringsdatan kan en användare i SUM använda MCF för att kolla på olika fall. Exempelvis kan en användare välja att bara kolla på projekt som når lansering och göra en känslighetsanalys för att ta reda på vilka ingående variabler som driver upp ett väldigt lönsamt projekt. 6.3 Framtida frågeställningar Många metoder vi har använt har varit beräkningsmässigt tidskrävande, framförallt beräk- ningen av Sobols index och den grafiska känslighetsanalysen. Detta beror huvudsakligen på att funktionen vi har använt i implementeringen inte är beräkningsmässigt effektiv. För att göra metoderna lättare att använda för Captario hade vår kod kunnat effektiviseras så att beräkningstiderna hade minskats. Vårt huvudsakliga fokus har varit på variationer i ingående variabler. Eftersom många av de ingående variablerna tar sina värden från skattade fördelningar kan det vara intressant att ut- föra en känslighetsanalys på parametrarna från dessa fördelningar. På så sätt kan man bygga 19 mer robusta modeller genom att förstå hur responsen påverkas av de antaganden man gjort. En sådan analys hade varit mer beräkningsmässigt krävande då eNPV behöver uppskattas för varje parameteruppsättning. Av de metoder som behandlats i detta kandidatarbete vore Morris metod den mest lämpade för en sådan typ av analys, eftersom den kräver förhållan- devis få funktionsanrop. Dessutom kan man anta att variabler som inte ger något utslag på en känslighetsanalys inte har någon känslighet i sin parameteruppsättning. Man kan alltså avskärma icke-signifikanta variabler inför en sådan analys. I Accelerated antas att alla ingående variabler är oberoende, vilket dessutom varit ett ge- nomgående antagande i samtliga metoder vi betraktat. I mer verklighetstrogna modeller bör detta inte vara fallet. Exempelvis om en aktivitet tar lång tid bör detta innebära en högre kostnad. För att få mer verklighetstrogna modeller och kunna utföra en känslighetsanalys kan modellerna byggas så att exempelvis tid och kostnad för en viss aktivitet slås samman till variabler som är oberoende i sig. På detta vis skulle känslighetsanalyser kunna utföras på mer verklighetstrogna modeller i SUM. 20 7 Referenser [1] Ciociola AA, Cohen LB, Kulkarni P. How Drugs are Developed and Approved by the FDA: Current Process and Future Directions, The American Journal of Gastroente- rology. 2014 Maj;109(5):620-3. doi:10.1038/ajg.2013.407 [2] DiMasi JA, Grabowski HG, Hansen RW. Innovation in the pharmaceutical industry: New estimates of R&D costs. Journal of Health Economics. 2016 Maj;47:20-33. doi:10.1016/j.jhealeco.2016.01.012 [3] Paul SM, Mytelka DS, Dunwiddle CT, Persinger CC, Munos BH, Lindborg SR, et al. How to improve R&D productivity: the pharmaceutical industry’s grand challenge. Nat Rev Drug Discov. 2010 mars;9(3):203-14. doi:10.1038/nrd3078 [4] Investopedia. Net Present Value - NPV [Internet]. New York: In- vestopedia, LLC; 2017 [citerad 30 mars 2017]. Hämtad från: http://www.investopedia.com/terms/n/npv.asp [5] Iooss B, Saltelli A. Introduction to Sensitivity Analysis. I: Ghanem R, Higdon D, Owhadi H, redaktörer. Handbook of Uncertainty Quantification. 1:a uppl. Zürich. Springer International Publishing; 2015. 1-17. [6] Saltelli A, Tarantola S, Campolongo F, Ratto M. Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models. 1:a uppl. Chichester, England: John Wiley Sons Ltd; 2002. 45. [7] Iooss B, Lema P. A Review on Global Sensitivity Analysis Methods. I: Dellino G, Me- loni C, redaktörer. Uncertainty Management in Simulation Optimization of Complex Systems. 59 uppl. New York. Springer; 2015. 101-22. [8] Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, et al.Global Sensitivity Analysis. 1:a uppl. Chichester, England: John Wiley & Sons Ltd; 2008. [9] Saltelli A, Chan K, Scott EM. Sensitivity analysis. 4:e uppl. Chichester, England: John Wiley Sons Ltd; 2006. [10] Daniel C. One-at-a-Time Plans. Journal of the American Statistical Association. 1973 jun;68(342):353-60. [11] Morris MD. Factorial Sampling Plans for Preliminary Computational Experiments. Technometrics. 1991 maj;33(2):161-74. [12] Campolongo F, Cariboni J, Saltelli A.An effective screening design for sensitivity analysis of large models. Environmental Modelling and Software. 2007;22(10):1509–18. Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, et al. Global Sensitivity Analysis. 1:a uppl. Chichester, England: John Wiley & Sons Ltd; 2008. [13] Saltelli A, Tarantola S, Campolongo F, Rattor M. Sensitivity analysis in practice. 1:a uppl. Chichester, England: John Wiley & Sons Ltd; 2004. [14] Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M, Tarantola S. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index.Computer Physics Communications. 2010 feb;181(2):259-70. [15] Sobol IM. Sensitivity Estimates for Nonlinear Mathematical Models. Matemathics and Computers in Simulation 2001;55(1):271-280. [16] Cukiert RI, Fortuin CM, Shuler KE, Petschek AG, Schaibly JH. Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. I Theory. The Journal of Chemical Physics. 1973;59(8):3873-8. 21 [17] Saltelli A, Tarantola S, Chan K. P.-S. A Quantitative Model-Independent Method for Global Sensitivity Analysis of Model Output. Techometrics. 1999;41(1):39-56. [18] Deniz A, Ferreira E, Johansson E, Johansson H, Lindbäck J, Wiskman L. sa-accelerated GitHub repository, 2017. https://github.com/wiskman/sa-accelerated [19] Sohier H, Piet-Lahanier H, JL Farges. Analysis and optimization of an air-launch-to- orbit separation Acta Astronautica. 2015 mar/apr ;108:18-29. http://doi.org/10.1016/j.actaastro.2014.11.043 [20] Saltelli A. Sensitivity Analysis for importance Assessment. Risk Analysis. 2002;22(3):579-90. 22 https://github.com/wiskman/sa-accelerated http://doi.org/10.1016/j.actaastro.2014.11.043 Bilagor A Morris metod på modell med Bernoullivariabler 0 10 20 30 40 50 60 70 µ ∗ 0 20 40 60 80 100 120 140 160 180 σ 1 5 10 26 27 28 30 32 Figur A.1: Diagram med absoluta medelvärde µ∗ och standardavikelse σ för de elementära effek- terna associerade till alla variabler. Endast de variabler med höga index är markerade med variabel- namn. 23 B Morris Implementeringsalgoritm Morris implementeringsmetod [11] bygger på att hyperkuben diskretiseras likformigt med n steg i varje dimension. Alltså betraktas den diskreta delmängden Ω = { 1 n+1 , 2 n+1 , . . . n n+1} p istället för hela parameterrummet. ∆ förbestäms till en heltalsmultipel av 1 n+1 . I nätet Ω evalueras f i varje nod. För att skatta en elementär effekt för varje variabler föreslog Morris följande tillvägagångsätt: i) Specificera en designmatris B av storlek l×p där l = p+1 och B innehåller enbart ettor och nollor. Denna ska även uppfylla att i varje kolonn i = 1, 2, . . . p så måste de i första elementen vara samma - vilka ofta sätts till 0. Exempelvis kan designmatrisen väljas till en nedåt triangulär matris med bij = 0 om i ≤ j och bij = 1 annars. ∆B kan i princip användas för att skatta en elementär effekt för (∆B)i+1 och (∆B)i skiljer sig med ∆ i en komponent. Men i och med att ∆B inte är slumpmässig, måste denna uppdateras. ii) Generera en slumpmässig p×p matris D∗ med 1 och −1 på diagonalen och de resterande elementen satta till 0. Diagonalelementen generas oberoende under villkoret att P (d∗ii = 1) = P (d∗ii = −1) = 1 2 . Låt sedan N = 1 2 ((2B − 1l×p)D∗ + 1l×p). Den första faktorn (2B − 1l×p) kommer ha samma värden som B i alla dess nollskilda element. De andra elementen är blir −1. Om sedan D∗ verkar på denna uppdaterade matris byts tecken på dess kolonner enligt diagonalelementen i D∗. Utförs sedan de sista operationena för att erhålla N kommer de element i (2B−1l×p)D∗ som är 1 bevaras, medan de element som är negativa bli 0. Vidare kommer varje kolonn i N vara lika med motsvarande kolonn i B eller i 1l×p − B enligt simuleringen av D∗. Två rader på varandra i N kommer endast skiljas i ett element med ±1. Detta är på grund av att två rader på varandra i B skiljer sig åt i exakt en kolonn. Efter uppdateringen kommer dessa element fortfarande skiljas åt, potentiellt i omvänd ordning. I de andra kolonnerna, där båda raderna är lika med 1 eller 0 samtidigt, kommer båda dessa element bevaras, eller ändras samtidigt efter uppdateringen till N . iii) Skapa en startvektor x∗ ∈ [0, 1]p där alla element är slumpade från { 1 n+1 , 2 n+1 , . . . n n+1 − ∆}. iv) En slumpmässig inverterbar permutationsmatris P ∗ av storlek p × p genereras i vilken varje kolonn innehåller exakt en etta och resten nollor. Med hjälp av denna beräknas en så kallad slumpmässig orientering, vilken är definerad: B∗ = (X∗+ ∆N)P ∗, där X∗ är en l×p-matris vars rader är x∗. Till följd av N :s egenskap, kommer skillnaden mellan två rader i (X∗+ ∆N) endast skilja med ±∆. Om P ∗ verkar på denna matris kommer dess kolonner att byta plats slumpmässigt. Vidare kommer även B∗ att ha denna egenskap. En slumpmässig orientering kan användas för att beräkna en elementäreffekt för varje varia- bel, i vilken raderna parvis utgör den data som fordras i (3.3). Detta genom Ei(x (i)) = ±f(x(i+1))− f(x(i)) ∆ (B.1) där x(i) är den i:te raden i B∗ och tecknet motsvarar det tecken som det enda nollskilda element i x(i+1) − x(i) har. Om man är intresserad av m stycken elementära effekter, för att uppskatta storheterna (3.4) och (3.5), genereras m stycken olika slumpmässiga orienteringar oberoende från varandra, vilka används för att beräkna m stycken elementära effekter för p variabler. 24 C Modellbeskrivning av Accelerated Nedan följer en mer detaljerad beskrivning av Captarios exempelmodel Accelerated som vår känslighetsanalys är utförd på. Accelerated-modellen består av 4 faser, varav varje fas består av olika aktiviteter som pågår parallellt. Projektet antas starta i mars 2016 (år 2016.3) och avslutas i början av 2047. Nedan är en bild som illustrerar starten till slutet av ett läkemedelsprojekt med accelerated-modellen: där varje rektangel representerar en aktivitet, och varje romb är ett vägskäl. En ifylld romb innebär att flera aktivitet startas alternativt avslutas och de tommma romberna avgör om läkemedlet går igenom fasen till vänster. Här nedan kommer vi beskriva de olika faserna med tillhörande variabler, såsom kostnader och tider. Vi låter LogN(a, b) beteckna en LogNormal- fördelning med a som 5%-kvantil och b som 95%-kvantil. C.1 Variabelförteckning Första fasen består av 4 olika aktiviteter: SMAD Study w biomarkers, Batch for ph2 studies, 3M TOX Study och Other ph1 activities. De variabler som tillhör vardera aktivitet är. Vi låter dessutom Tij vara tiden det tar att genomföra aktivitet j i den t:te fasen. Cij och Ĉij betecknar den diskonterade, respektive den icke-diskonterade kostnaden för just denna aktivitet. Aktivitet Variabelnamn Fördelning Förklaring SMAD Study w cadtime LogN(2/3, 5/6) tiden biomarkers Development_costs1 LogN(2, 11/5) kostnaden Batch for ph2 Time LogN(3/5, 4/5) tiden studies Development_costs2 LogN(3/10, 2/5) kostnaden 3M TOX Study 3mToxTime LogN(1/2, 4/5) tiden 3mToxCost LogN(1/4, 3/10) kostnaden Other ph1 activities cadtime – tiden (se SMAD) Other_Ph1_Costs LogN(3/2, 9/5) kostnader Dessutom gäller att sannolikheten att projektet går vidare till nästa fas 50%. Det vill säga I1 ∼ Bernoulli(0.5). Andra fasen består av 4 olika aktiviteter: 3M POC Study 400 patients, Batch for ph3, 6M TOX Studies och Other ph2 Costs. Associerade variabler till dessa är: 25 Aktivitet Variabelnamn Fördelning Förklaring POC_SetupTime LogN(1/10, 1/5) POC_AnalysisTime LogN(1/10, 3/20) 3M POC Study POC_RecruitmentRate LogN(1, 3) 400 patients POC_SetupCost LogN(1/10, 3/10) POC_AnalysisCost LogN(1/10, 1/5) POC_CostPerPatient LogN(3/100, 1/25) POC_CostPerCenter LogN(1, 6/5) Batch for ph3 batchForPH3Time LogN(0.6, 1) tiden batchForPH3Cost LogN(1, 3) kostnaden 6M TOX Studies 6MToxTime LogN(1, 1.1) tiden 6MToxCost LogN(0.3, 0.4) kostnaden Other ph2 Costs 3MPOCTime LogN(0.75, 1) tiden Development_costs_3 LogN(2, 3) kostnaden Tiden för 3M POC Study 400 patients fås genom T2,1 = POC_SetupTime + 400 3× 12× POC_RecruitmentRate + POC_AnalysisTime + 0.25 där 400 är antalet patienter i studien, 3 är antalet kliniker och 12 konverterar enheten till år, eftersom POC_RecruitmentRate är rekryteringstakt i patient/(klinik×månad). Den sista termen är en fix behandlingstid. Kostnaden för denna aktivitet fås av: Ĉ2,1 = POC_SetupCost + POC_AnalysisCost + POC_FSI-LSO_Cost, POC_FSI-LSO_Cost = 400× POC_CostPerPatient + 3× POC_CostPerCenter. Projektet går vidare till nästa fas med sannolikheten 0.5. Alltså gäller att I2 ∼ Bernoulli(0.5). Tredje fasen består av 2 aktiviteter: 6M Ph3 Study 2500 patients och Other Ph3 Costs. Här är en lista på variablerna som tillhör samtliga aktiviteter Aktivitet Variabelnamn Fördelning Förklaring 6M Ph3 Study Ph3StudyTime LogN(1, 1.5) tiden 2500 Patients ph3studyCosts LogN(75, 95) kostnaden Other Ph3 Costs Ph3StudyTime – tiden (se 6M Ph3) ph3OtherCosts LogN(10, 20) kostnaden Sannolikheten att gå igenom fas 3 är 70%, vidare gäller att I3 ∼ Bernoulli(0.7). Registreringsfasen är bara en aktivitet. Vi går alltid igenom registreringsfasen. Här är en lista på variablerna som tillhör registreringsfasen: 26 Variabelnamn Fördelning Förklaring RegTime LogN(1, 1.1) tiden RegCosts LogN(4, 6) kostnaden Sannolikheten att gå igenom registreringsfasen är 90%, och därför blir den sista bernoulliva- riabeln I4 ∼ Bernoulli(0.9). C.2 Fasernas kostnader och tider Tiden att genomföra den i:te fasen är Ti := max{Ti,1, Ti,2, . . . Ti,n}. Den diskonterade kostna- den för en fas fås genom att summera alla diskonterade kostnader för alla aktiviteter i fasen, det vill säga Ci = ∑ j Ci,j . För att diskontera kostnaden för en aktivitet antas att kostnaden är likformigt fördelad över hela tidspannet aktiviteten äger rum. För årsvis diskontering blir den nuvärdesjusterade kostnaden Ci,j = ∫ Ti,j+T (s) i T (s) i Ĉi,j Ti,j r[t]+1 dt (C.1) där fasens starttid T (s) i = ∑ k 16.9 där 16.9 tiden till patentet av produkten går ut. Under antagandet att PeakMarketShare är konstant, modelleras alltså att företagets marknadsandelar ökar linjärt upp till en den max- imala marknadsandelen PeakMarketShare. Efter att patentet gått ut, förloras 90% av alla marknadsandelar. I allmänhet kommer PeakMarketShare inte vara konstant, utan kommer vara beroende på vilka konkurrenter som lanserat. Denna variabel fås av PeakMarketShare NumberOfEntries och OrderOfEntry 0.5 NumberOfEntries = 1, OrderOfEntry = 1 0.4 NumberOfEntries = 2, OrderOfEntry = 1 0.2 NumberOfEntries = 2, OrderOfEntry = 2 0.3 NumberOfEntries = 3, OrderOfEntry = 1 0.2 NumberOfEntries = 3, OrderOfEntry = 2 0.15 NumberOfEntries = 3, OrderOfEntry = 3 där OrderOfEntry och NumberOfEntries fås av: OrderOfEntry = 1 + ∑ i Comp(i)DoesLaunch× IAi(T4) NumberOfEntries = 1 + ∑ i Comp(i)DoesLaunch× IAi(t) IAi(t) = { 1, Comp(i)LaunchTime− 2016.3 ≤ t 0, Comp(i)LaunchTime− 2016.3 > t . Ett sätt utvecklingen i erhållna marknadsandelar kan bli illustreras i C.1 t M a r k e t S h a r e s Figur C.1: Den tjockare linjen är ett scenario av vilka marknadsandelar som erhålls över tid, i vilken en konkurrent lanserar innan (RampTimeInYears + T4) och en efter. De streckade graferna är den utveckling som hade skett om inga konkurrenter lanserat efter företaget lanserat läkemedlet. 28 C.4 Andra projektvärdesstorheter Tabell C.1: Lista över andra projektvärdesstorheter Captario SUM beräknar Förklaring NPV Summan av all framtida nuvärdesjusterad kassaflöde. eNPV Väntevärdet av NPV. Nominell kostnad Summan av alla kostnader under projektets gång. Nominell avkastning Summan av alla intäkter under projektets gån. ROINom Nominell avkastning på investering: Nominell avkastning−Nominell kostnad Nominell kostnad . ROIPV Nuvärdesjusterad avkastning på investering. pf Risken att göra förlust: P (NPV < 0). TRE Lanseringstid. TA Tid tills kliniska studier avslutar. Speciellt är TA = TRE om läkemedlet lanserar. 29 Bakgrund Beslutsfattning inom läkemedelsutveckling Captario och programvaran SUM Känslighetsanalys Syfte Teori Notation för allmänna modeller Metoder inom känslighetsanalys One-at-a-time-experiment Morris metod Monte Carlo-filtrering Sobols index Fourier Amplitude Sensitivity Test – FAST Metod Evaluering av NPV för modeller i SUM Marknadsmodell Val av responsfunktion Val och implementering av känslighetsanalysmetoder Resultat Monte Carlo-filtrering Morris metod Sobols index och FAST Grafisk känslighetsanalys Relevans för Captario Diskussion Hantering av modeller Applicering i SUM Framtida frågeställningar Referenser Bilagor Morris metod på modell med Bernoullivariabler Morris Implementeringsalgoritm Modellbeskrivning av Accelerated Variabelförteckning Fasernas kostnader och tider Marknadsmodellen Andra projektvärdesstorheter