R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

SIR-malli markkinoinnissa

Johdanto

Epidemiologinen SIR-malli pyrkii mallintamaan infektion leviämistä populaatiossa. Mallissa populaatio on jaettu kolmeen luokkaan, alttiisiin (susceptible eli S), infektoituneisiin (infected eli I) ja tervehtyneisiin (recovered eli R). Yksilöt voivat siirtyä näiden luokkien välillä tietyillä nopeuksilla (S -> I = beta eli infection rate ja I -> R = recovery rate), mutta recovered-luokan oletetaan olevan siinä mielessä pysyvä, ettei siihen kerran juotunut yksilö voi enää myöhemmin muuttua uudelleen alttiiksi tai tartuttavaksi.

Mallia voidaan laajentaa monin tavoin, esimerkiksi lisäämällä siihen syntyvyys ja kuolleisuus, jolloin voidaan mallintaa myös pidempien ajanjaksojen kuluessa toistuvia epidemioita. Lisäksi malliin voidaan lisätä vaikkapa rokotusten (asiakkaan omistuksessa on jo toinen, kilpaileva tuote) vaikutus epidemioiden kulkuun.

Markkinointimielessä yllä mainitut suureet voisi ehkä koodata seuraavasti:

S = potentiaalisten asiakkaiden lukumäärä
I = nykyisten asiakkaiden lukumäärä
R = entisten asiakkaiden lukumäärä

infection rate = uusien asiakkaiden kertymisnopeus
recovery rate = tuotteesta luopuvien asiakkaiden kertymisnopeus

R-koodi

Seuraavassa käytetään R:n deSolve-pakettia tarvittavien yhtälöiden ratkaisemiseen. Alkuperäinen, hieman muuttelemani koodi on kopioitu muualta.

Itse simulaatiokoodi koostuu kahdesta funktiosta: Funktio sir() tekee simuloinnin, ja plot.sir() puolestaan hoitaa tuloksen visualisoinnin:

sir <- function(beta=1.5, recovery=1/12, death=0.00005, birth=0.00005, steps=365, step.by=1, susceptible=0.048, infected=0.002, recovered=0.95){
    require(deSolve)
    pars<-c(beta = beta, recovery = recovery, death = death, birth = birth)
    print(pars)
    times <- seq(from = 0, to = steps, by = step.by)                 # we want to run the model for 3000 time steps
    yinit <- c(Susc = susceptible, Infected = infected, Recovered = recovered) # this parameter sets the initial conditions, susceptible = 0.98*0.02+(1-0.974)=0.046
    SIR_model <- function(times, yinit, pars){
    with(as.list(c(yinit,pars)), {
            dSusc      <- birth - beta*Infected*Susc                     - death*Susc
            dInfected  <-         beta*Infected*Susc - recovery*Infected - death*Infected
            dRecovered <-                              recovery*Infected - death*Recovered
            return(list(c(dSusc, dInfected, dRecovered)))})
    }
    results <- ode(func = SIR_model, times = times, y = yinit, parms = pars)
    results <- as.data.frame(results)
    return(results)
}
 
plot.sir <- function(results, N=1000) {
   matplot(results[, 1], results[, 2:4]*N, type="l", lty=1, lwd=2, las=1, xlab="Aika (pv.)", ylab="Henkilöiden lukumäärä")
   abline(v=seq(from=0, to=nrow(results), by=28), col="grey95")
   abline(v=seq(from=0, to=nrow(results), by=365), col="grey65")
   abline(h=seq(from=0, to=max(N*results[,2:4]), by=N/10), col="grey95")
   matplot(results[, 1], results[, 2:4]*N, type="l", lty=1, lwd=2, add=TRUE)
   legend("topright", col=1:3, legend=c("S", "I", "R"), lwd=2, bty="n")
   box()
}

Parametrien arviointi

Funktioiden oletusarvot vastaavat karkeasti tuhkarokkoepidemian simulointiin tarvittavia parametrejä. Parametrit voidaan kuitenkin arvioida myös muille ”epidemioille”, kuten uuden tuotteen leviämiselle populaatioon. Jos lisääntymisluku R0 on se luku, joka kuvaa jokaisen ”infektoituneen” tartuttamien toisten henkilöiden keskimääräistä lukumäärää, voidaan mallin beta arvioida kaavasta:

beta = R0 * recovery rate

”Infektoituneiden” parantumisvauhti (recovery rate) puolestaan voidaan laskea keskimääräisen ”infektion” kestosta:

recovery rate = 1 / duration of infection

Ajatellaan vaikkapa uuden tietoteknisen laitteen tuloa Suomen markkinoille. Jos kiinnitetään markkinoinnille alttiin porukan muutosnopeudet death = 0 ja birth = 0 eli uusia asiakkaita ei saavu markkinoille, eikä vanhoja toisaalta myöskään poistu niiltä, ja tarkastellaan 5 vuoden aikajännettä eli steps=365*5 ja tarkastellaan muutoksia päivittäin eli step=1 on jo perusoletukset saatu arvioitua. Jos arvioidaan, että vaikkapa 10 % väestöstä on kiinnostunut uudesta tuotteesta saadaan, ja aluksi tuotteen on ottanut käyttöön 0.1% väestöstä saadaan eri luokkien osuuksia kuvaaville parametreille seuraavat arvot: susceptible = 100/1000, infected = 1/1000 ja recovered = 899/1000.

Jos edelleen ajatellaan, että laitteesta ollaan valmiita luopumaan esimerkiksi kahden vuoden syklillä (infektio siis kestää 2 vuotta), saadaan recovery rate = 1/(365*2). Siitä voidaan sitten arvioida erilaisia betan arvoa riippuen siitä, kuinka monta muuta kunkin tuotteen hankkineen uskotaan saavan myös hankkimaan tuotteen. Jos ajatellaan, että lisääntymisluku olisi vaikkapa 200, voitaisiin malli simuloida seuraavalla komennolla:

plot.sir(results<-sir(beta=200*1/(365*2), recovery=1/(365*2), steps=365*5, death=0, birth=0, step.by=1, susceptible=100/1000, infected=1/1000, recovered=899/1000), N=1000000)

Tuloksena on seuraava kuva, jos oletetaan, että populaation koko markkinointialueella on 1 000 000:

Näillä arvoilla tuote saavuttaisi suurimman mahdollisen sen ostaneiden asiakkaiden määrän jo alle vuodessa. Toisaalta, miltei kaikki olisivat jo luopuneet siitä alle viidessä vuodessa.

Mallin parametrit voisi varmasti arvioida käytännössä joidenkin tuotteiden myyntiluvuista, mutten onnistunut tähän hätän löytämään soveltuvia myyntilukuja, joilla tätä olisi voinut yrittää. Ehdotuksia otetaan vastaan….