Simulatsioonid

Simulatsioonid on läbi viidud Tartu Ülikooli teadusarvutuste keskuse arvutusklastris, kasutades järgmisi R-funktsioone.

library(logistf)
library(dplyr)
# N - valimimaht
# p - diagnoosi osakaal, p0 - kontrollide CNV osakaal,
# OR - tegelik šansside suhe
# p01 - valepositiivsuse määr, p10 - valenegatiivsuse määr
simulatsioon1 <- function(N, p, p0, OR, p01, p10){
    n1 = N * p # juhtude arv
    n0 = N - n1 # kontrollide arv
    p1 = p0 * OR / (p0 * (OR - 1) + 1) # juhtude CNV osakaal
    p11 = 1 - p10 #õigepositiivsuse määr
    p1_penn = (p1*p11+(1-p1)*p01) # juhtude CNV osakaal PennCNV tunnuse kaudu
    p0_penn = (p0*p11+(1-p0)*p01) # kontrollide CNV osakaal PennCNV tunnuse kaudu
    T0 = seq(0,n0) # CNV kandjate võimalikud arvud kontrollide hulgas
    T1 = seq(0,n1) # CNV kandjate võimalikud arvud juhtude hulgas
    pT0 = dbinom(T0, n0, p0_penn) # tõenäosused P(T_0=i)
    names(pT0) = T0
    pT1 = dbinom(T1, n1, p1_penn) # tõenäosused P(T_1=j)
    names(pT1) = T1
    m = outer(pT0, pT1) # paaride (i,j) tõenäosused P(T_0=i,T_1=j)
    # jätame alles need paarid  (i,j), mille tõenäosus on piirist suurem
    piir = 1e-10
    v = which(m > piir, arr.ind = T)
    paarid = data.frame(i = as.numeric(rownames(m)[v[, 1]]),
                        j = as.numeric(colnames(m)[v[, 2]]))
    tulemused = NULL
    for (k in 1:nrow(paarid)){
        i = paarid[k,1] # kontrollid
        j = paarid[k,2] # juhud
        x = c(c(rep(1, i), rep(0, n0-i), c(rep(1, j), rep(0, n1-j)))) # CNV tunnus
        y = c(rep(0, n0), rep(1, n1)) #diagnoosi tunnus
        m1 = glm(y~x, family = binomial(link = "logit"))
        m2 = logistf(y~x, pl=F, firth=F)
        tulemused = rbind(tulemused,data.frame(glm.pvalue=last(coef(summary(m1))),
                                                glm.OR = last(exp(coef(m1))),
                                                firth.pvalue = last(m2$prob),
                                                firth.OR = last(exp(coef(m2))),
                                                p_ij = m[as.character(i),as.character(j)],
                                                p=p, p0=p0, OR=OR, p01=p01, p10=p10))
    }
    return(tulemused)
}

library(logistf)
library(dplyr)
# valepositiivsete jaotus TÜ EGV andmete põhjal
S0 = function(n=1){
    u = runif(n,0,1)
    s1 = c(runif(sum((u<=0.02*25.5)),0,0.02),runif(sum(u>0.02*25.5),0.02,1))}
# õigepositiivsete jaotus TÜ EGV andmete põhjal
S1 = function(n=1){
    u = runif(n,0,1)
    s0 = c(runif(sum(u<=0.98*0.8),0,0.98),runif(sum(u>0.98*0.8),0.98,1))}
# N - valimimaht
# p - diagnoosi osakaal, p0 - kontrollide CNV osakaal,
# OR - tegelik šansside suhe
# p01 - valepositiivsuse määr, p10 - valenegatiivsuse määr
simulatsioon2 <- function(N, p, p0, OR, p01, p10, arv){
    n1 = N * p # juhtude arv
    n0 = N - n1 # kontrollide arv
    p1 = p0 * OR / (p0 * (OR - 1) + 1) # juhtude CNV osakaal
    p11 = 1 - p10 # õigepositiivsuse määr
    tulemused = NULL
    for (i in 1:arv) {
# kvaliteediskoori simuleerimine
QS = function(n,cnv_freq){
    u = runif(n,0,1)
    qs = c(S1(sum(u <= cnv_freq*p11)),
        S0(sum((u > cnv_freq*p11) & (u <= (cnv_freq*p11 + (1-cnv_freq)*p01)))),
        rep(0, sum(u>cnv_freq*p11 + (1-cnv_freq)*p01)))}
x = c(QS(n0, p0), QS(n1, p1)) # CNV tunnus
y = c(rep(0, n0), rep(1, n1)) # diagnoosi tunnus
m1 = glm(y~x, family = binomial(link = "logit"))
m2 = logistf(y~x)
tulemused = rbind(tulemused, data.frame(glm.pvalue=last(coef(summary(m1))),
                                        glm.OR = last(exp(coef(m1))),
                                        firth.pvalue = last(m2$prob),
                                        firth.OR = last(exp(coef(m2))),
                                        p=p, p0=p0, OR=OR, p01=p01, p10=p10))
}
  return(tulemused)
}