Randomizačné procedúry
Lenka Filova
Generovanie randomizačnych postupností
Na generovanie a porovnávanie randomizačných procedúr máme v R k dispozícii balík randomizeR.
Na začiatok si vygenerujme postupnosť dĺžky \(N=10\) pomocou kompletnej randomizácie, ktorá je ekvivalentná hádzaniu mincou. Takáto randomizácia vedie k \(2^N\) postupnostiam, z ktorých sú všetky rovnako pravdepodobné. crPar nižie je tzv. konštruktor postupnosti a param obsahuje informácie o randomizačnej postupnosti.
library(randomizeR)## Warning: package 'randomizeR' was built under R version 3.6.3## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 3.6.3## Loading required package: plotrix## Loading required package: survivalN <- 10
(params <- crPar(N))##
## Object of class "crPar"
##
## design = CR
## N = 10
## groups = A BKeď chceme vygenerovať len jednu postupnosť, postupujeme nasledovne (seed sa ukladá kvôli reprodukovateľnosti).
(R <- genSeq(params))##
## Object of class "rCrSeq"
##
## design = CR
## seed = 1442293831
## N = 10
## groups = A B
##
## The sequence M:
##
## 1 B B A A B A B B A BgetRandList(R)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] "B" "B" "A" "A" "B" "A" "B" "B" "A" "B"plotSeq(R, plotAllSeq = T)
Všetky postupnosti kompletnej randomizácie s N=10:
(allSeqs <- getAllSeq(params))##
## Object of class "crSeq"
##
## design = CR
## N = 10
## groups = A B
##
## The first 3 of 1024 sequences of M:
##
## 1 A A A A A A A A A A
## 2 B A A A A A A A A A
## 3 A B A A A A A A A A
## ...Ak by bolo postupností priveľa alebo generovanie by bolo príliš náročné, môžeme simulačne vygenerovať nejaký menší počet postupností pomocou genSeq
N <- 50
params <- crPar(N)
(randomSeqs <- genSeq(params, r = 10000))##
## Object of class "rCrSeq"
##
## design = CR
## seed = 124808508
## N = 50
## groups = A B
##
## The first 3 of 10000 sequences of M:
##
## 1 A A A B A B A B B A ...
## 2 B B B A B B B A B A ...
## 3 B B B A B A B B A B ...
## ...Výpočet pravdepodobnosti nastatia danej postupnosti:
p <- getProb(allSeqs)
head(data.frame(Sequences = cbind(getRandList(allSeqs)), Probability = round(p, 6)))## Sequences.1 Sequences.2 Sequences.3 Sequences.4 Sequences.5 Sequences.6
## 1 A A A A A A
## 2 B A A A A A
## 3 A B A A A A
## 4 B B A A A A
## 5 A A B A A A
## 6 B A B A A A
## Sequences.7 Sequences.8 Sequences.9 Sequences.10 Probability
## 1 A A A A 0.000977
## 2 A A A A 0.000977
## 3 A A A A 0.000977
## 4 A A A A 0.000977
## 5 A A A A 0.000977
## 6 A A A A 0.000977Posudzovanie kvality randomizačných metód
Predpokladajme, že odozvy pri oboch ošetreniach sú normálne rozdelené s nulovými strednými hodnotami \(\mu_A=\mu+B=0\) a jednotkovými disperziami \(\sigma_A^2=\sigma_B^2=1\).
muA <- muB <- 0
sigmaA <- sigmaB <- 1
normalEndpoint <- normEndp(mu = c(muA, muB), sigma = c(sigmaA, sigmaB))Vygenerujeme objekt reprezentujúci pravdepodobnosť zamietnutia za prítomnosti chronologickej odchýlky vzniknutej vplyvom lineárneho časového trendu s \(\theta=1\):
(cb <- chronBias(type = "linT", theta = 1, method = "exact"))##
## Object of class "chronBias"
##
## TYPE = linT
## THETA = 1
## METHOD = exact
## ALPHA = 0.05Chceme vyhodnotiť Big Stick Design s N=12 a toleranciou na nevyváženosť mti=2 vzhľadom na chronologickú a vyberovú odchýlku a stratu sily:
N <- 12
mti <- 2
bsdSeq <- getAllSeq(bsdPar(N, mti))
d <- 1.796 #rozdiel vyberovych priemerov
sb <- selBias("CS", eta = d/4, method = "exact")
cb <- chronBias("linT", theta = 1/N, method = "exact")
pw <- setPower(d, method = "exact")Funkcia assess potom vypočíta pravdepodobnosť zamietnutia pre každé kritérium a každú postupnosť. V našom prípade máme 972 možných postupností.
prvý stĺpec: postupnosť
druhý stĺpec: pravdepodobnosť nastatia
ďalšie stĺpce: pravdepodobnosť zamietnutia za prítomnosti daného typu odchýlky
stĺpec
power(exact)popisuje silu každej postupnosti
Súhrnné charakteristiky rozdelení alokačných postupností získame pomocou summary:
(A <-assess(bsdSeq, sb, cb, pw, endp = normalEndpoint))##
## Assessment of a randomization procedure
##
## design = BSD(2)
## N = 12
## K = 2
## groups = A B
##
##
## The first 3 rows of 972 rows of D:
##
## Sequence Probability P(rej)(CS) P(rej)(linT) power(exact)
## 1 BBABABABA ... 0.004 0.045 0.05 0.789
## 2 BABBABABA ... 0.002 0.042 0.05 0.789
## 3 ABBBABABA ... 0.002 0.039 0.05 0.789
## ...summary(A)## P(rej)(CS) P(rej)(linT) power(exact)
## mean 0.056 0.05 0.795
## sd 0.013 0.00 0.006
## max 0.109 0.05 0.800
## min 0.034 0.05 0.789
## x05 0.037 0.05 0.789
## x25 0.048 0.05 0.789
## x50 0.054 0.05 0.789
## x75 0.062 0.05 0.800
## x95 0.079 0.05 0.800Porovnanie metód randomizácie
Chceme porovnať Big Stick Design, maximálnu procedúru s mti = 2 a blokovo-permutovanú randomizáciu s blokom veľkosti 4 vzhľadom na náchylnosť k selekčnej odchýlke.
mpSeq <- getAllSeq(mpPar(N, mti))
bc <- rep(4, N/4)
pbrSeq <- getAllSeq(pbrPar(bc))
(C <- compare(sb, bsdSeq, mpSeq, pbrSeq, endp = normalEndpoint))##
## Comparison for P(rej)(CS)
##
## BSD.2. MP.2. PBR.4.
## mean 0.056 0.072 0.082
## sd 0.013 0.015 0.015
## max 0.109 0.109 0.109
## min 0.034 0.040 0.050
## x05 0.037 0.050 0.061
## x25 0.048 0.061 0.072
## x50 0.054 0.072 0.079
## x75 0.062 0.079 0.099
## x95 0.079 0.100 0.103plot(C)
plot(C, y = "boxplot")## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique