library(mcr) library(parallel) ## Algoritma. # 1. istenilen da??l?m ve ?zelliklerde K adet veri seti ?ret. # 2. Her k. veri seti i?in # 2.1 'mcreg' fonksiyonu yard?m? ile ilgilenilen y?nteme ait g?ven aral?klar ile elde et. # 2.2 Bootstrap y?ntemi ile g?ven aral?klar ile elde edilmesi i?in B adet bootstrap ?rneklemi # olu?tur ve ilgilenilen parametre i?in bootstrap g?ven aral?klar ile elde et. # 3. 2.1 ve 2.2 ad?mlar?ndann elde edilen g?ven aral?klar?n? coverage ve/veya tip 2 hata d?zeylerini # K tekrar ?zerinden hesapla. cl <- makeCluster(6, "PSOCK") set.seed(2064) #generateData fonksiyonu b0 sabit hatalar?, b1 oransal hatal??, #sd.resid aral?klar? st.sapmas??, x referans y?ntem, y test y?ntemi #mean referans y?ntemin ortalamas?, sd ise referans y?ntemin st.sapmas?? #olmak ?zere n g?zlemli x ve y verisini t?retir. generateData <- function(n = 100, b0 = 0, b1 = 1, sd.xerr = 5, sd.yerr = sd.x, xdist = c("symmetric", "skewed"), errdist = c("symmetric", "skewed"), xdist.param = list(mean = 0, sd = 1), ...){ # n = 5000 # b0 = 0 # b1 = 1 # sd.xerr = 5 # sd.yerr = sd.xerr # xdist = "skewed" #c("symmetric", "skewed") # errdist = "symmetric" #c("symmetric", "skewed") # xdist.param = list(mean = 135.5, sd = 30) # errdist.param = list(mean = 0, sd = 1) xdist <- match.arg(xdist) errdist <- match.arg(errdist) # Generate true values of reference test (X) --- Target Values # X'in dağılımı skewed olduğu durumda veriler ki-kare dağılımından üretiliyor. # Ki-kare dağılımından üretilen veriler daha sonra istenilen ortalamaya eşit olması için # shift ediliyor. if (xdist == "symmetric"){ x <- rnorm(n, xdist.param$mean, xdist.param$sd) } else { v <- xdist.param$sd^2 df <- v / 2 x <- rchisq(n, df = df) + (xdist.param$mean - df) } # Generate true values of new test Y using b0 and b1 parameters --- Target Values y <- b0 + b1 * x # Add Random Error to X and Y (Analytical Variances) if (xdist == "symmetric"){ x <- x + rnorm(n = n, mean = 0, sd = sd.xerr) y <- y + rnorm(n = n, mean = 0, sd = sd.yerr) } else { vx <- sd.xerr^2 dfx <- vx / 2 vy <- sd.yerr^2 dfy <- vy / 2 x <- x + (rchisq(n = n, df = dfx) - dfx) y <- y + (rchisq(n = n, df = dfy) - dfy) } return(list(x = x, y = y)) } ################### Örnek veri üretimi. ############### # x <- rnorm(n, xdist.param$mean, xdist.param$sd) # y <- b0 + b1 * x # # x <- x + rnorm(n = n, mean = 0, sd = sd.xerr) # y <- y + rnorm(n = n, mean = 0, sd = sd.yerr) # # # Var(y) = b1^2 * Var(x) # Var(y) = b1^2 * Var(x) + Var(e.y) # # Var(y) = b1^2 * Var(x + e.x) + Var(e.y) # = b1^2 * {Var(x) + Var(e.x)} + Var(e.y) set.seed(1) df <- generateData(n = 10000, sd.xerr = 0.5, sd.yerr = 0.5, xdist = "symmetric", errdist = "symmetric", xdist.param = list(mean = 135.5, sd = 3.8)) max(df$x) / min(df$x) max(df$y) / min(df$y) a <- 135.5 + c(-1, 1) * 3 * 30 a[2] / a[1] plot(df) dat <- df tmp <- mcreg(dat$x, dat$y, error.ratio = 1, method.reg = "Deming", method.ci = "analytical", method.bootstrap.ci = "quantile", mref.name = "Reference.Test", mtest.name = "New.Test", na.rm = TRUE, rng.seed = 3627, nsamples = 499) tmp@para max(df$x) / min(df$x) max(df$y) / min(df$y) ######################################################### ## Simulation parameters: n <- seq(40, 200, 40) SDx_SDy_Ratio <- c(0.5, 1, 2) # Error ratio b0 <- c(0, 5, 10) b1 <- c(1, 1.1, 1.2) # range <- c(1, 2, 3) # 1: low, 2: medium, 3: large outlier <- c(0, 1, 2) # number of outliers distX_distErr <- c("nn", "ns", "sn", "ss") simCombs <- expand.grid(N = n, SDx_SDy_Ratio = SDx_SDy_Ratio, b0 = b0, b1 = b1, range = range, outlier = outlier, distX_distErr = distX_distErr) simCombs <- simCombs[1:3, ] #Benzetim say?s?? N <- 250 #generateData fonksiyonu argumanlar?? ile N say?da benzetim verisi t?retir. allData <- lapply(1:N, function(xx){ generateData(n = 100, b0 = 0, b1 = 1, sd.xerr = 5, sd.yerr = 5, xdist = "symmetric", errdist = "symmetric", xdist.param = list(mean = 135, sd = 25)) }) ### Simulasyon kombinasyonları # Sample Size: {} #b0 ve b1 için sonuçlar. clusterSetRNGStream(cl, 3627) allRes <- parLapply(cl, X = allData, fun = function(X, b0.null = 0, b1.null = 1, ...){ library(mcr) dat <- X tmp <- mcreg(dat$x, dat$y, error.ratio = 1, method.reg = "PaBa", method.ci = "analytical", method.bootstrap.ci = "quantile", mref.name = "Reference.Test", mtest.name = "New.Test", na.rm = TRUE, rng.seed = 3627, nsamples = 499) coeff.est <- getCoefficients(tmp) ci.b0 <- coeff.est["Intercept", c("EST", "LCI", "UCI")] if (ci.b0["LCI"] <= b0.null & ci.b0["UCI"] >= b0.null){ ci.b0 <- c(ci.b0, 0) } else { ci.b0 <- c(ci.b0, 1) } ci.b1 <- coeff.est["Slope", c("EST", "LCI", "UCI")] if (ci.b1["LCI"] <= b1.null & ci.b1["UCI"] >= b1.null){ ci.b1 <- c(ci.b1, 0) } else { ci.b1 <- c(ci.b1, 1) } names(ci.b1) <- names(ci.b0) <- c("EST", "LCI", "UCI", "Decision") return(list(b0.res = ci.b0, b1.res = ci.b1)) }) b1Res <- b0Res <- NULL for (i in 1:length(allRes)){ tmp0 <- allRes[[i]]$b0.res tmp1 <- allRes[[i]]$b1.res b0Res <- rbind(b0Res, tmp0) b1Res <- rbind(b1Res, tmp1) } prop.table(table(b0Res[ ,"Decision"])) prop.table(table(b1Res[ ,"Decision"]))