Statistiline jõud

Seekordse kirjatüki eesmärk on rääkida statistilisest jõust. Annan väga põgusa ülevaate, mida see mõiste endast kujutab ning näitan ühte võimalikku lähenemist sele arvutamiseks. Artikkel ei pretendeeri kõikehõlmavuse ning ainulaadsuse tiitlile. Suurem osa materjali põhineb raamatul Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan 2nd Edition. Kui ei tea suurt midagi statistilise testimise kohta, loe seda artiklit.

Jõud

Statisitiline jõud näitab klassikaliste statistiliste testide puhul, kui suure tõenäosusega test lükkab ümber nullhüpoteesi, kui nullhüpotees pole tõene. Võib sõnastada ka nii, et tõenäosus, et test aktsepteerib alternatiivse hüpoteesi, kui see on tõsi (kuigi reaalsuse ei ole need kaks ekvivalentsed, võime siin sellise käsitlusega leppida). Kordame üle, mida null- ja alternatiivne hüpotees tähendasid:

nullhüpotees: nagu ka nimi reedab, väidab see hüpotees, et nähtuste vahel ei ole seoseid
alternatiivne hüpotees: vastpidiselt nullhüpoteesile väidab see, et nähtuste vahel on seos (võib anda ka konreetsema väite, nt mis suunaline see seos võiks olla)

Kuna olen viimasel ajal rohkem tegelenud Bayes’i andmeanalüüsi meetoditega, ei lähtu ma klassikalisest lähenemisest. Laiemalt tähendab statistiline jõud, et planeeritud uuring saavutab oma eesmärgi, kui selle planeerimise aluseks olevad andmed kirjeldavad maailma piisavalt hästi. Uuringute eesmärgiks on info saamine maailma kohta. Uuringu eesmärgiks võivad olla: - mingi parameetri väärtuse (nullhüpoteesi) ümberlükkamine - mingi parameetri väärtuse (nullhüpoteesi) kinnitamine - parameetri väärtuse võimalikult täpne mõõtmine

Statistilise jõu hindamiseks kasutatakse Bayes’i andmeanalüüsi meetodite puhul järgnevat lähenemist:

hüpoteetilisest parameetrite jaotusest genereeri juhuslikud parameetrite väärtused. Siin võib aluseks võtta eelnevate uuringute andmed, mis võivad anda infot parameetrite võimalike väärtuste kohta.
eelnevas punktis genereeritud parameetrite väärtustest genereeri valim andmeid. Oluline on, et andmete kogumine jälgiks protsessi, kuidas pärast uuringus tegelikult andmeid kogutakse (kas fikseeritud on valimi suurus, kogumise aeg jne)
simuleeritud andmete põhjal arvuta nö tuletatud hinnang, kasutades Bayes’i reeglit. Analüüs peab olema läbi viidud nii, et see veenaks ka skeptilisi vaatlejaid.
eelmises punktis läbi viidud analüüsi põhjal hinda, kas eesmärk saavutati või mitte (eesmärgiks võib olla siis ükskõik milline kolmest eelmainitud eesmärgist).
korda kõiki eelnevaid samme mitu korda, et hinnata uuringu statistilist jõudu. Statistiline jõud näitab eesmärkide saavutamist pikemas perspektiivis, mistõttu ühe korra sammude läbiviimine võib juhuslikkuse tõttu anda vale tulemuse.

Näide

Esiteks laen sisse funktsioonid, mis on näite jaoks vajalikud. Nende töötamise loogikast arusaamiseks pead pühendama aega, kuna nende selgitamine siin, võtaks palju ruumi ja aega. Kui asja vastu huvi, loe eelmainitud raamatut.

#funktsioonid eelmainitud raamatust, et eist aru saada, pead natukene aega pühendama ja soovituslikult ka raamatu läbi lugema
HDIofICDF = function( ICDFname , credMass=0.95 , tol=1e-8 , ... ) {
  # Arguments:
  #   ICDFname is R's name for the inverse cumulative density function
  #     of the distribution.
  #   credMass is the desired mass of the HDI region.
  #   tol is passed to R's optimize function.
  # Return value:
  #   Highest density iterval (HDI) limits in a vector.
  # Example of use: For determining HDI of a beta(30,12) distribution, type
  #   HDIofICDF( qbeta , shape1 = 30 , shape2 = 12 )
  #   Notice that the parameters of the ICDFname must be explicitly named;
  #   e.g., HDIofICDF( qbeta , 30 , 12 ) does not work.
  # Adapted and corrected from Greg Snow's TeachingDemos package.
  incredMass =  1.0 - credMass
  intervalWidth = function( lowTailPr , ICDFname , credMass , ... ) {
    ICDFname( credMass + lowTailPr , ... ) - ICDFname( lowTailPr , ... )
  }
  optInfo = optimize( intervalWidth , c( 0 , incredMass ) , ICDFname=ICDFname ,
                      credMass=credMass , tol=tol , ... )
  HDIlowTailPr = optInfo$minimum
  return( c( ICDFname( HDIlowTailPr , ... ) ,
             ICDFname( credMass + HDIlowTailPr , ... ) ) )
}

minNforHDIpower = function( genPriorMode , genPriorN ,
                            HDImaxwid=NULL , nullVal=NULL , 
                            ROPE=c(max(0,nullVal-.02),min(1,nullVal+.02)) ,
                            desiredPower=0.8 , audPriorMode=0.5 , audPriorN=2 ,
                            HDImass=0.95 , initSampSize=20 , verbose=TRUE ) {
  # Check for argument consistency:
  if ( !xor( is.null(HDImaxwid) , is.null(nullVal) ) ) {
    stop("One and only one of HDImaxwid and nullVal must be specified.")
  }
    # Convert prior mode and N to a,b parameters of beta distribution:
  genPriorA = genPriorMode * (genPriorN-2) + 1
  genPriorB = ( 1.0 - genPriorMode ) * (genPriorN-2) + 1
  audPriorA = audPriorMode * (audPriorN-2) + 1
  audPriorB = ( 1.0 - audPriorMode ) * (audPriorN-2) + 1
  # Initialize loop for incrementing sampleSize:
  sampleSize = initSampSize
  notPowerfulEnough = TRUE
  # Increment sampleSize until desired power is achieved:
  while( notPowerfulEnough ) {
    zvec = 0:sampleSize # vector of all possible z values for N flips.
    # Compute probability of each z value for data-generating prior:
    pzvec = exp( lchoose( sampleSize , zvec )
                 + lbeta( zvec + genPriorA , sampleSize-zvec + genPriorB )
                 - lbeta( genPriorA , genPriorB ) )
    # For each z value, compute posterior HDI: 
    # hdiMat will hold HDI limits for each z:
    hdiMat = matrix( 0 , nrow=length(zvec) , ncol=2 )
    for ( zIdx in 1:length(zvec) ) {
      z = zvec[zIdx]
      hdiMat[zIdx,] = HDIofICDF( qbeta , 
                                 shape1 = z + audPriorA ,
                                 shape2 = sampleSize - z + audPriorB ,
                                 credMass = HDImass )
    }
    # Compute HDI widths:
    hdiWid = hdiMat[,2] - hdiMat[,1]
    # Sum the probabilities of outcomes with satisfactory HDI widths:
    if ( !is.null( HDImaxwid ) ) {
      powerHDI = sum( pzvec[ hdiWid < HDImaxwid ] )
    }
    # Sum the probabilities of outcomes with HDI excluding ROPE:
    if ( !is.null( nullVal ) ) {
      powerHDI = sum( pzvec[ hdiMat[,1] > ROPE[2] | hdiMat[,2] < ROPE[1] ] )
    }
    if ( verbose ) {
      cat( " For sample size = ", sampleSize , ", power = " , powerHDI ,
           "\n" , sep="" ) ; flush.console()
    }
    if ( powerHDI > desiredPower ) {  # If desired power is attained,
      notPowerfulEnough = FALSE       # set flag to stop,
    } else {                          # otherwise
      sampleSize = sampleSize + 1     # increment the sample size.
    }
  } # End while( notPowerfulEnough ).
  # Return the sample size that achieved the desired power:
  return( sampleSize )
} # End of function.

Oletame, et meil üks münt, mille “ausust” tahame hinnata. Meil on juba info, 10 eelneva mündiviske kohta (mis ei ole väga kindel), et visete keskmine tulemus on 0.8 (st et münt ei ole “aus”). Seda infot kasutame parameetrite genereerimiseks. Bayes’i arvutuste aluseks võtame lisaks info, mis rahuldab skeptilisi vaatlejaid. Oletame, et selleks on info, et eelduslikult on münt “aus” (keskmine tulemus on 0.5), kuid me ei ole selles väga kindel (visete arv, mis näitab ka konsentratsiooni, on 2). Meie eesmärgiks on parameetri väärtuse võimalikult täpne mõõtmine (see jääks vahemikku +/-0.2). Selleks, et saavutada statistiline jõud 80% (st 80% juhtudest saavutaksime oma eesmärgi), peab valim (visete arv) olema:

sampSize = minNforHDIpower( genPriorMode=0.80, genPriorN=10,
HDImaxwid=0.20, nullVal=NULL, ROPE=NULL,
desiredPower=0.8,
audPriorMode=0.5,#sektpikutele sobiv mood (prior)
audPriorN=2, #skeptikutele sobib kontsentratsioon, näitab kui tugevalt usume, et mood on selline
HDImass=0.95, initSampSize=80, verbose=TRUE )

##  For sample size = 80, power = 0.6925385
##  For sample size = 81, power = 0.7112756
##  For sample size = 82, power = 0.7290568
##  For sample size = 83, power = 0.7459331
##  For sample size = 84, power = 0.7619551
##  For sample size = 85, power = 0.7771725
##  For sample size = 86, power = 0.7916341
##  For sample size = 87, power = 0.8053871

Nagu näha peas valim olema 87 viset, et saavutaksime oma eesmärgi.

Meie eesmärk võib olla ka teistsugune. Näiteks tahame välistada vahemiku 0.4-0.46, kui mündi visete keskmise väärtuse. Sel juhul on valimi suuruseks vaja:

sampSize = minNforHDIpower( genPriorMode=0.80, genPriorN=10,
HDImaxwid=NULL, nullVal=0.43, ROPE=c(0.4,0.46),
desiredPower=0.8,
audPriorMode=0.5,#sektpikutele sobiv mood (prior)
audPriorN=2, #skeptikutele sobib kontsentratsioon, näitab kui tugevalt usume, et mood on selline
HDImass=0.95, initSampSize=30, verbose=TRUE )

##  For sample size = 30, power = 0.7400778
##  For sample size = 31, power = 0.7762192
##  For sample size = 32, power = 0.7570467
##  For sample size = 33, power = 0.789481
##  For sample size = 34, power = 0.7727936
##  For sample size = 35, power = 0.8019163

Sel juhul oleks meil vaja 35 viset.

See oli väga põgus sissejuhatus statistilisse jõusse. Miks seda arvutada? See aitab katset/uuringut paremini planeerida: saad teada, kui suurt valimit on vaja, et soovitud tulemust saavutada. Kui kogud liiga vähe andmeid, siis võid kogu uuringu põhja lasta. Teisipidi aitab see ka üle kulutamist vältida, sa ei kogu igaks juhuks liiga palju andmeid (mis võib olla kulukas).

Lõppu ka väike kriitika. Kui nähtus, mida uuritakse, on väga kiiresti muutuv või tema kohta pole andmeid, siis on statistilise jõu arvutamine suhteliselt keeruline. Siis võiks alustada lihtsalt andmete kogumisest ning vaadata, mis andmetest välja joonistub. Pärast seda saab juba põhjalikumat uuringut planeerida (kui nähtus on püsiv).

Kasutatud kirjandus

Kruschke, J., 2014, Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan 2nd Edition

Statistiline jõud

Risto Hinno

juuni 2016

Jõud

Näide

Kasutatud kirjandus