library("DFIT") library("xtable") library("ggplot2") ## Dichotomous item parameters taken from Wright 2011 data("dichotomousItemParameters", package = "DFIT") ## Polytomous item parameters taken from Raju 2009 data("polytomousItemParameters", package = "DFIT") ## dataSubset: Rasch model raschParameters <- lapply(dichotomousItemParameters, function(x) x[, 2, drop = FALSE]) raschParameters <- as.list(unique(as.data.frame(raschParameters))) raschParameters <- lapply(raschParameters, function(x) matrix(x, ncol = 1)) ## Three parameter logistic model items3Pl <- c(2, 20, 22, 8, 10, 28, 46, 32) threePlParameters <- lapply(dichotomousItemParameters, function(x) x[items3Pl, ]) ## itemParameters: lapply(threePlParameters, head, 5) lapply(polytomousItemParameters, head, 5) ## difExample: One parameter logistic model normal metric ## NCDIF, CDIF and DFT ncdif1pl <- Ncdif(itemParameters = raschParameters, irtModel = "1pl", focalAbilities = NULL, focalDistribution = "norm", subdivisions = 5000, logistic = FALSE) cdif1pl <- Cdif(itemParameters = raschParameters, irtModel = "1pl", focalAbilities = NULL, focalDistribution = "norm", subdivisions = 5000, logistic = FALSE) (dtf1plWithCdif <- Dtf(cdif = cdif1pl)) (dtf1plWithoutCdif <- Dtf(cdif = NULL, itemParameters = raschParameters, irtModel = "1pl", focalAbilities = NULL, focalDistribution = "norm", subdivisions = 5000, logistic = FALSE)) ## RajuExample: Raju's Area Measures (SA: Signed Area, UA: Unsigned Area) sam1pl <- SignedArea(itemParameters = raschParameters, irtModel = "1pl", subdivisions = 5000, logistic = FALSE) uam1pl <- UnsignedArea(itemParameters = raschParameters, irtModel = "1pl", subdivisions = 5000, logistic = FALSE) ## Mantel-Haenszel mh1pl <- IrtMh(itemParameters = raschParameters, irtModel = "1pl", focalDistribution = "norm", referenceDistribution = "norm", focalDistrExtra = list(mean = 0, sd = 1), referenceDistrExtra = list(mean = 0.5, sd = 1), groupRatio = 1, logistic = FALSE) delta1pl <- DeltaMhIrt(mh1pl) ## table1pl items <- seq(nrow(raschParameters[["focal"]])) tab1pl <- cbind(Item = items, `Reference b` = raschParameters[["reference"]][items, ], `Focal b` = raschParameters[["focal"]][items, ], MH = mh1pl, `ETS $\\Delta$` = delta1pl, SA = sam1pl[items, ], UA = uam1pl[items, ], NCDIF = ncdif1pl, CDIF = cdif1pl) xtab1pl <- xtable(tab1pl) digits(xtab1pl) <- c(0, 0, 1, 1, 3, 2, 2, 2, 5, 5) align(xtab1pl) <- rep("r", length(digits(xtab1pl))) caption(xtab1pl) <- "DIF statistics for example items under the 1pl model." label(xtab1pl) <- "tab:1plItems" ## tab1plPrint print(xtab1pl, sanitize.colnames.function = function(x) { x }, include.rownames = FALSE) ## PlotDIF: No DIF item it5PlotD <- PlotNcdif(iiItem = 5, itemParameters = raschParameters, irtModel = "1pl", plotDensity = TRUE, logistic = FALSE, focalDensityText = "Focal group density (theoretical)", main = "Item 5. NO DIF. 1PL model") ## DIF item close to focal distribution it7PlotS <- PlotNcdif(iiItem = 7, itemParameters = raschParameters, irtModel = "1pl", plotDensity = FALSE, logistic = FALSE, main = "Item 7. Uniform DIF. 1PL model") ## otherPlotsRasch: Plot examples for item with and without uniform DIF under 1pl ## model DIF item close to focal distribution it7PlotD <- PlotNcdif(iiItem = 7, itemParameters = raschParameters, irtModel = "1pl", plotDensity = TRUE, logistic = FALSE, main = "Item 7. Uniform DIF. 1PL model") ## DIF item far from focal distribution it3PlotS <- PlotNcdif(iiItem = 3, itemParameters = raschParameters, irtModel = "1pl", plotDensity = FALSE, logistic = FALSE, main = "Item 3. Uniform DIF. 1PL model") ## plotNCDIFRaschNoDif it5PlotD ## plotNCDIFRaschDifClose it7PlotS it3PlotS ## difExampleOtherModels: Three parameter logistic model normal metric ## NCDIF, CDIF and DFT ncdif3pl <- Ncdif(itemParameters = threePlParameters, irtModel = "3pl", focalAbilities = NULL, focalDistribution = "norm", subdivisions = 5000, logistic = FALSE) cdif3pl <- Cdif(itemParameters = threePlParameters, irtModel = "3pl", focalAbilities = NULL, focalDistribution = "norm", subdivisions = 5000, logistic = FALSE) (dtf3plWithCdif <- Dtf(cdif = cdif3pl)) ## Raju's Area Measures (SA: Signed Area, UA: Unsigned Area) sam3pl <- SignedArea(itemParameters = threePlParameters, irtModel = "3pl", subdivisions = 5000, logistic = TRUE) uam3pl <- UnsignedArea(itemParameters = threePlParameters, irtModel = "3pl", subdivisions = 5000, logistic = FALSE) ## Mantel-Haenszel mh3Pl <- IrtMh(itemParameters = threePlParameters, irtModel = "3pl", focalDistribution = "norm", referenceDistribution = "norm", focalDistrExtra = list(mean = 0, sd = 1), referenceDistrExtra = list(mean = 0, sd = 1), groupRatio = 1, logistic = FALSE) delta3Pl <- DeltaMhIrt(mh3Pl) ## Non Uniform - != guess DIF item close to focal distribution it46PlotS <- PlotNcdif(iiItem = 7, itemParameters = threePlParameters, irtModel = "3pl", plotDensity = FALSE, logistic = FALSE, main = "Item 46. Non uniform DIF with different guessing") ## GRM parameter logistic model normal metric ## NCDIF, CDIF and DFT ncdifGrm <- Ncdif(itemParameters = polytomousItemParameters, irtModel = "grm", focalAbilities = NULL, focalDistribution = "norm", subdivisions = 5000, logistic = FALSE) cdifGrm <- Cdif(itemParameters = polytomousItemParameters, irtModel = "grm", focalAbilities = NULL, focalDistribution = "norm", subdivisions = 5000, logistic = FALSE) (dtfGrmWithCdif <- Dtf(cdif = cdifGrm)) ## Raju's Area Measures (SA: Signed Area, UA: Unsigned Area) samGrm <- SignedArea(itemParameters = polytomousItemParameters, irtModel = "grm", subdivisions = 5000, logistic = TRUE) uamGrm <- UnsignedArea(itemParameters = polytomousItemParameters, irtModel = "grm", subdivisions = 5000, logistic = FALSE) ## Mixed DIF item it4PlotS <- PlotNcdif(iiItem = 4, itemParameters = polytomousItemParameters, irtModel = "grm", plotDensity = FALSE, logistic = FALSE, main = "Item 4. Mixed DIF. GRM", ylab = "Expected score") ## table3pl items <- items3Pl reference <- apply(threePlParameters[["reference"]], 1, function(x) paste("(", x[1], ", ", x[2], ", ", x[3], ")", sep = "")) focal <- apply(threePlParameters[["focal"]], 1, function(x) paste("(", x[1], ", ", x[2], ", ", x[3], ")", sep = "")) tab3Pl <- cbind(Item = items, Reference = reference, Focal = focal, MH = sprintf("%.2f", mh3Pl), `ETS $\\Delta$` = sprintf("%.3f", delta3Pl), SA = sprintf("%.3f", sam3pl), UA = sprintf("%.3f", uam3pl), NCDIF = sprintf("%.4f", ncdif3pl), CDIF = sprintf("%.3f", cdif3pl)) xtab3Pl <- xtable(tab3Pl) align(xtab3Pl) <- rep("r", length(digits(xtab3Pl))) caption(xtab3Pl) <- "DIF statistics for example items under the 3Pl model." label(xtab3Pl) <- "tab:threePlItems" ## tab3PlPrint print(xtab3Pl, sanitize.colnames.function = function(x) { x }, include.rownames = FALSE) ## plotNCDIF3pl it46PlotS ## tablegrm items <- seq(nrow(polytomousItemParameters[["reference"]])) reference <- apply(polytomousItemParameters[["reference"]], 1, function(x) paste("(", round(x[1], 2), ", ", x[2], ", ", x[3], ", ", x[4], ", ", x[5], ")", sep = "")) focal <- apply(polytomousItemParameters[["focal"]], 1, function(x) paste("(", round(x[1], 2), ", ", x[2], ", ", x[3], ", ", x[4], ", ", x[5], ")", sep = "")) tabGrm <- cbind(Item = items, Reference = reference, Focal = focal, SA = sprintf("%.1f", samGrm), UA = sprintf("%.2f", uamGrm), NCDIF = sprintf("%.3f", ncdifGrm), CDIF = sprintf("%.3f", cdifGrm)) xtabGrm <- xtable(tabGrm) align(xtabGrm) <- rep("r", length(digits(xtabGrm))) caption(xtabGrm) <- "DIF statistics for example items under the GRM model." label(xtabGrm) <- "tab:GRM" ## tabGrmPrint print(xtabGrm, sanitize.colnames.function = function(x) { x }, include.rownames = FALSE, scalebox = "0.8") ## plotNCDIFgrm it4PlotS ## cutoffIprWithNoAse: Cutoff points directly from item parameters and calculating ## asymptotic variance-covariance set.seed(89334828) cutoffRaschNcdif1 <- CutoffIpr(quantiles = 0.95, itemParameters = raschParameters, itemCovariances = "asymptotic", nullGroup = "focal", irtModel = "1pl", focalSampleSize = 500, referenceSampleSize = 1500, referenceDistrExtra = list(mean = 0.5, sd = 1), logistic = TRUE, statistic = "ncdif", nReplicates = 1000) ## ASERasch: Asymptotic variance-covariance matrices Rasch nullParameters <- list() nullParameters[["focal"]] <- raschParameters[["focal"]] nullParameters[["reference"]] <- raschParameters[["focal"]] raschAse <- list() raschAse[["focal"]] <- AseIrt(itemParameters = nullParameters[["focal"]], distribution = "norm", sampleSize = 500, irtModel = "1pl", distributionParameters = list(mean = 0, sd = 1), logistic = TRUE) raschAse[["reference"]] <- AseIrt(itemParameters = nullParameters[["reference"]], distribution = "norm", sampleSize = 1500, irtModel = "1pl", distributionParameters = list(mean = 0.5, sd = 1), logistic = TRUE) ## cutoffIprDirect: Cutoff points directly from item parameters and ## variance-covariance set.seed(29834328) cutoffRaschNcdif2 <- CutoffIpr(quantiles = 0.95, itemParameters = nullParameters, itemCovariances = raschAse, irtModel = "1pl", logistic = TRUE, statistic = "ncdif", nReplicates = 1000) ## iprFromIPRstatistics: Cutoff points from DIF statistics on generated item ## replications set.seed(29833326) raschIpr <- Ipr(itemParameters = nullParameters, itemCovariances = raschAse, nReplicates = 1000) raschNcdifIpr <- IprNcdif(itemParameterList = raschIpr, irtModel = "1pl", logistic = TRUE) raschUamIpr <- IprUam(itemParameterList = raschIpr, irtModel = "1pl", logistic = TRUE) raschSamIpr <- IprSam(itemParameterList = raschIpr, irtModel = "1pl", logistic = TRUE) raschMhIpr <- IprMh(itemParameterList = raschIpr, irtModel = "1pl", logistic = TRUE) ## cutoffIprFromIPRstatistics cutoffRaschNcdif3 <- CutoffIpr(quantiles = 0.95, iprStatistics = raschNcdifIpr, itemParameterList = raschIpr, itemParameters = nullParameters, itemCovariances = raschAse, irtModel = "1pl", statistic = "ncdif") cutoffRaschUam <- CutoffIpr(quantiles = 0.95, iprStatistics = raschUamIpr, itemParameterList = raschIpr, itemParameters = nullParameters, itemCovariances = raschAse, irtModel = "1pl", statistic = "uam") cutoffRaschSam <- CutoffIpr(quantiles = c(0.025, 0.95), iprStatistics = raschSamIpr, itemParameterList = raschIpr, itemParameters = nullParameters, itemCovariances = raschAse, irtModel = "1pl", statistic = "sam") cutoffRaschMh <- CutoffIpr(quantiles = c(0.025, 0.975), iprStatistics = raschMhIpr, itemParameterList = raschIpr, itemParameters = nullParameters, itemCovariances = raschAse, irtModel = "1pl", statistic = "mh") ## tableNcdif items <- seq(length(ncdif1pl)) ncdifRasch <- Ncdif(itemParameters = raschParameters, irtModel = "1pl", focalAbilities = NULL, focalDistribution = "norm", subdivisions = 5000, logistic = TRUE) samRasch <- SignedArea(itemParameters = raschParameters, irtModel = "1pl", subdivisions = 5000, logistic = TRUE) uamRasch <- UnsignedArea(itemParameters = raschParameters, irtModel = "1pl", subdivisions = 5000, logistic = TRUE) mhRasch <- IrtMh(itemParameters = raschParameters, irtModel = "1pl", focalDistribution = "norm", referenceDistribution = "norm", focalDistrExtra = list(mean = 0, sd = 1), referenceDistrExtra = list(mean = 0.5, sd = 1), groupRatio = 1, logistic = TRUE) tabNcdif <- cbind(Item = items, `Reference b` = raschParameters[["reference"]][items, ], `Focal b` = raschParameters[["focal"]][items, ], `True NCDIF` = sprintf("%.5f", ncdifRasch), `Cutoff 1` = sprintf("%.5f", cutoffRaschNcdif1$quantiles), `Cutoff 2` = sprintf("%.5f", cutoffRaschNcdif2$quantiles), `Cutoff 3` = sprintf("%.5f", cutoffRaschNcdif3$quantiles)) xtabNcdif <- xtable(tabNcdif) align(xtabNcdif) <- rep("r", length(digits(xtabNcdif))) caption(xtabNcdif) <- "NCDIF cutoff points under the Rasch model through the IPR approach." label(xtabNcdif) <- "tab:ncdif" ## tabNcdifPrint print(xtabNcdif, sanitize.colnames.function = function(x) { x }, include.rownames = FALSE) #, scalebox = '0.8') ## tableIpr items <- seq(length(ncdif1pl)) tabIpr <- cbind(Item = items, `MH lower` = sprintf("%.3f", cutoffRaschMh$quantiles[1, ]), `True MH ` = sprintf("%.3f", mhRasch), `MH upper` = sprintf("%.3f", cutoffRaschMh$quantiles[2, ]), `SA lower` = sprintf("%.4f", cutoffRaschSam$quantiles[1, ]), `True SA` = samRasch[items, ], `SA upper` = sprintf("%.4f", cutoffRaschSam$quantiles[2, ]), `True UA` = uamRasch[items, ], UA = sprintf("%.4f", cutoffRaschUam$quantiles)) xtabIpr <- xtable(tabIpr) align(xtabIpr) <- rep("r", length(digits(xtabIpr))) caption(xtabIpr) <- "Cut-off points for different DIF statistics under the Rasch model through the IPR approach." label(xtabIpr) <- "tab:ipr" ## tabIprPrint print(xtabIpr, sanitize.colnames.function = function(x) { x }, include.rownames = FALSE, scalebox = "0.9") ## powerCalculationSetUp nReplicates <- 3000 nFocal <- 800 nReference <- 2500 kRatio <- nReference/nFocal focalParam <- list(mean = 0, sd = 1) referenceParam <- list(mean = 0, sd = 1) ## UniformPowerSetup Uniform power example itemParameters <- list(focal = cbind(rep(1, 51), seq(-0.5, 0.5, length = 51)), reference = cbind(rep(1, 51), rep(0, 51))) nullFocal <- which(itemParameters[["focal"]][, 2] == itemParameters[["reference"]][, 2]) itemParametersNull <- lapply(itemParameters, function(x) x[nullFocal, , drop = FALSE]) names(itemParametersNull) <- c("focal", "reference") ## trueNCDIFUnif: True NCDIF twoPlUniNcdifTrue <- Ncdif(itemParameters, irtModel = "2pl", focalDistribution = "norm", focalDistrExtra = focalParam, logistic = FALSE) ## powerUniformCutoffProposedAse: Calculate cutoff point proposed algorithm twoPlUniAse <- list() twoPlUniAse[["focal"]] <- AseIrt(itemParameters = itemParameters[["focal"]], distribution = "norm", distributionParameters = focalParam, sampleSize = nFocal, irtModel = "2pl", logistic = FALSE) twoPlUniAse[["reference"]] <- AseIrt(itemParameters = itemParameters[["reference"]], distribution = "norm", distributionParameters = referenceParam, sampleSize = nReference, irtModel = "2pl", logistic = FALSE) ## powerUniformCutoffProposedIpr set.seed(29834328) twoPlUniIpr <- Ipr(itemParameters = itemParameters, itemCovariances = twoPlUniAse, nReplicates = nReplicates) twoPlUniNcdif <- IprNcdif(itemParameterList = twoPlUniIpr, irtModel = "2pl", logistic = FALSE, subdivisions = 1000) ## powerUniformCutoffProposed cutoffPointEachSZUni <- CutoffIpr(quantiles = 0.95, iprStatistics = twoPlUniNcdif[nullFocal, , drop = FALSE]) ## powerUniformcutoffCurrent: Calculate cutoff point current algorithm set.seed(29834328) twoPlUniAseCurrent <- twoPlUniAse twoPlUniAseCurrent[["focal"]] <- twoPlUniAseCurrent[["focal"]][nullFocal] twoPlUniAseCurrent[["reference"]] <- lapply(twoPlUniAseCurrent[["reference"]][nullFocal], "*", kRatio) twoPlUniIprCurrent <- Ipr(itemParameters = itemParametersNull, itemCovariances = twoPlUniAseCurrent, nReplicates = nReplicates) twoPlUniNcdifCurrent <- IprNcdif(itemParameterList = twoPlUniIprCurrent, irtModel = "2pl", logistic = FALSE, subdivisions = 1000) cutoffPointUni <- CutoffIpr(quantiles = 0.95, iprStatistics = matrix(twoPlUniNcdifCurrent, nrow = length(nullFocal))) ## powerUniformCalculation: Calculate power powerUni <- data.frame(Type = "Cutoff with only focal group", b = itemParameters[["focal"]][, 2], NCDIF = (twoPlUniNcdifTrue * sign(itemParameters[["focal"]][, 2] - itemParametersNull[["focal"]][, 2])), Power = rowMeans(twoPlUniNcdif > cutoffPointUni$quantiles)) powerUniEach <- data.frame(Type = "Cutoff with both groups", b = itemParameters[["focal"]][, 2], NCDIF = (twoPlUniNcdifTrue * sign(itemParameters[["focal"]][, 2] - itemParametersNull[["focal"]][, 2])), Power = rowMeans(twoPlUniNcdif > cutoffPointEachSZUni$quantiles)) power80P <- which(abs(powerUniEach[, "Power"] - 0.8) == min(abs(powerUniEach[, "Power"] - 0.8))) powerUniJ <- rbind(powerUni, powerUniEach) ## plotUniPowerPrep powerPlotUniDif <- ggplot(powerUniJ, aes(x = b, y = Power)) powerPlotUniDif <- powerPlotUniDif + xlim(-0.5, 0.5) powerPlotUniDif <- powerPlotUniDif + geom_line(aes(colour = Type, linetype = Type), size = 1) powerPlotUniNcdif <- ggplot(powerUniJ, aes(x = NCDIF, y = Power)) powerPlotUniNcdif <- powerPlotUniNcdif + xlim(-0.012, 0.012) powerPlotUniNcdif <- powerPlotUniNcdif + geom_line(aes(colour = Type, linetype = Type), size = 1) ## plotpowerUnidif powerPlotUniDif <- powerPlotUniDif + geom_hline(yintercept = 0.05, linetype = "longdash", colour = "grey") powerPlotUniDif <- powerPlotUniDif + geom_hline(yintercept = powerUniEach[power80P, "Power"], linetype = "longdash", colour = "grey") powerPlotUniDif <- powerPlotUniDif + geom_hline(yintercept = powerUni[power80P, "Power"], linetype = "longdash", colour = "grey") powerPlotUniDif <- powerPlotUniDif + geom_vline(xintercept = powerUni[power80P, "b"], linetype = "longdash", colour = "grey") powerPlotUniDif <- powerPlotUniDif + geom_vline(xintercept = -powerUni[power80P, "b"], linetype = "longdash", colour = "grey") powerPlotUniDif ## plotpowerUnincdif powerPlotUniNcdif <- powerPlotUniNcdif + geom_hline(yintercept = 0.05, linetype = "longdash", colour = "grey") powerPlotUniNcdif <- powerPlotUniNcdif + geom_hline(yintercept = powerUniEach[power80P, "Power"], linetype = "longdash", colour = "grey") powerPlotUniNcdif <- powerPlotUniNcdif + geom_hline(yintercept = powerUni[power80P, "Power"], linetype = "longdash", colour = "grey") powerPlotUniNcdif <- powerPlotUniNcdif + geom_vline(xintercept = powerUniEach[power80P, "NCDIF"], linetype = "longdash", colour = "grey") powerPlotUniNcdif <- powerPlotUniNcdif + geom_vline(xintercept = -powerUniEach[power80P, "NCDIF"], linetype = "longdash", colour = "grey") powerPlotUniNcdif ## powerNonUniCalculationSetUp NonUniform power example itemParameters <- list(focal = cbind(unique(c(seq(0.5, 1, length = 26), 1, seq(1, 1/0.5, length = 26))), rep(0, 51)), reference = cbind(rep(1, 51), rep(0, 51))) nullFocal <- which(itemParameters[["focal"]][, 1] == itemParameters[["reference"]][, 1]) itemParametersNull <- lapply(itemParameters, function(x) matrix(x[nullFocal, ], nrow = length(nullFocal))) names(itemParametersNull) <- c("focal", "reference") ## True NCDIF twoPlNonUniNcdifTrue <- Ncdif(itemParameters, irtModel = "2pl", focalDistribution = "norm", focalDistrExtra = focalParam, logistic = FALSE) ## powerNonUniformCutoffProposed: Calculate cutoff point proposed algorithm twoPlNonUniAse <- list() twoPlNonUniAse[["focal"]] <- AseIrt(itemParameters = itemParameters[["focal"]], distribution = "norm", distributionParameters = focalParam, sampleSize = nFocal, irtModel = "2pl", logistic = FALSE) twoPlNonUniAse[["reference"]] <- AseIrt(itemParameters = itemParameters[["reference"]], distribution = "norm", distributionParameters = referenceParam, sampleSize = nReference, irtModel = "2pl", logistic = FALSE) set.seed(29834328) twoPlNonUniIpr <- Ipr(itemParameters = itemParameters, itemCovariances = twoPlNonUniAse, nReplicates = nReplicates) twoPlNonUniNcdif <- IprNcdif(itemParameterList = twoPlNonUniIpr, irtModel = "2pl", logistic = FALSE, subdivisions = 1000) cutoffPointEachSZNonUni <- CutoffIpr(quantiles = 0.95, iprStatistics = matrix(twoPlNonUniNcdif[nullFocal, ], nrow = length(nullFocal))) ## powerNonUniformcutoffCurrent: Calculate cutoff point current algorithm set.seed(29834328) twoPlNonUniAseCurrent <- twoPlNonUniAse twoPlNonUniAseCurrent[["focal"]] <- twoPlNonUniAseCurrent[["focal"]][nullFocal] twoPlNonUniAseCurrent[["reference"]] <- lapply(twoPlNonUniAseCurrent[["reference"]][nullFocal], "*", kRatio) twoPlNonUniIprCurrent <- Ipr(itemParameters = itemParametersNull, itemCovariances = twoPlNonUniAseCurrent, nReplicates = nReplicates) twoPlNonUniNcdifCurrent <- IprNcdif(itemParameterList = twoPlNonUniIprCurrent, irtModel = "2pl", logistic = FALSE, subdivisions = 1000) cutoffPointNonUni <- CutoffIpr(quantiles = 0.95, iprStatistics = matrix(twoPlNonUniNcdifCurrent, nrow = length(nullFocal))) ## powerNonUniformCalculation: Calculate power powerNonUni <- data.frame(Type = "Cutoff with only focal group", a = itemParameters[["focal"]][, 1], NCDIF = (twoPlNonUniNcdifTrue * sign(itemParameters[["focal"]][, 1] - itemParametersNull[["focal"]][, 1])), Power = rowMeans(twoPlNonUniNcdif > cutoffPointNonUni$quantiles)) powerNonUniEach <- data.frame(Type = "Cutoff with both groups", a = itemParameters[["focal"]][, 1], NCDIF = (twoPlNonUniNcdifTrue * sign(itemParameters[["focal"]][, 1] - itemParametersNull[["focal"]][, 1])), Power = rowMeans(twoPlNonUniNcdif > cutoffPointEachSZNonUni$quantiles)) power80N <- which(abs(powerNonUniEach[, "Power"] - 0.8) == min(abs(powerNonUniEach[, "Power"] - 0.8))) powerNonUniJ <- rbind(powerNonUni, powerNonUniEach) ## plotNonUniPowerPrep powerPlotNonUniDif <- ggplot(powerNonUniJ, aes(x = a, y = Power)) powerPlotNonUniDif <- powerPlotNonUniDif + xlim(0.5, 2) powerPlotNonUniDif <- powerPlotNonUniDif + geom_hline(yintercept = 0.05, linetype = "longdash", colour = "grey") powerPlotNonUniDif <- powerPlotNonUniDif + geom_hline(yintercept = powerNonUniEach[power80N, "Power"], linetype = "longdash", colour = "grey") powerPlotNonUniDif <- powerPlotNonUniDif + geom_hline(yintercept = powerNonUni[power80N, "Power"], linetype = "longdash", colour = "grey") powerPlotNonUniDif <- powerPlotNonUniDif + geom_vline(xintercept = powerNonUni[power80N, "a"], linetype = "longdash", colour = "grey") powerPlotNonUniDif <- powerPlotNonUniDif + geom_vline(xintercept = -powerNonUni[power80N, "a"], linetype = "longdash", colour = "grey") powerPlotNonUniDif <- powerPlotNonUniDif + geom_line(aes(colour = Type, linetype = Type), size = 1) powerPlotNonUniNcdif <- ggplot(powerNonUniJ, aes(x = NCDIF, y = Power)) powerPlotNonUniNcdif <- powerPlotNonUniNcdif + xlim(-0.015, 0.015) powerPlotNonUniNcdif <- powerPlotNonUniNcdif + geom_hline(yintercept = 0.05, linetype = "longdash", colour = "grey") powerPlotNonUniNcdif <- powerPlotNonUniNcdif + geom_hline(yintercept = powerNonUniEach[power80N, "Power"], linetype = "longdash", colour = "grey") powerPlotNonUniNcdif <- powerPlotNonUniNcdif + geom_hline(yintercept = powerNonUni[power80N, "Power"], linetype = "longdash", colour = "grey") powerPlotNonUniNcdif <- powerPlotNonUniNcdif + geom_vline(xintercept = powerNonUniEach[power80N, "NCDIF"], linetype = "longdash", colour = "grey") powerPlotNonUniNcdif <- powerPlotNonUniNcdif + geom_vline(xintercept = -powerNonUniEach[power80N, "NCDIF"], linetype = "longdash", colour = "grey") powerPlotNonUniNcdif <- powerPlotNonUniNcdif + geom_line(aes(colour = Type, linetype = Type), size = 1) ## valuesNCDIF80P power80Pa <- (abs(powerUniEach[, "Power"] - 0.8)) # == powerUniEach[c(15:19, 33:37), ] power80Na <- (abs(powerNonUniEach[, "Power"] - 0.8)) # == powerNonUniEach[c(11:15, 33:37), ] lowUni <- powerUniEach[17, "NCDIF"] lowNonUni <- powerNonUniEach[13, "NCDIF"] ## plotpowerNonUnidif powerPlotNonUniDif ## plotpowerNonUnincdif powerPlotNonUniNcdif ## powerFixed: Global variables nFocal <- 500 nReference <- 1500 distriParam <- list(mean = 0, sd = 1) nReplicates <- 3000 ## asePower nullParameters <- list() nullParameters[["focal"]] <- raschParameters[["focal"]] nullParameters[["reference"]] <- raschParameters[["focal"]] ## Variance-covariance matrices under null hypothesis nullAse <- list() nullAse[["focal"]] <- AseIrt(itemParameters = nullParameters[["focal"]], distribution = "norm", sampleSize = nFocal, distributionParameters = distriParam, irtModel = "1pl", logistic = FALSE) nullAse[["reference"]] <- AseIrt(itemParameters = nullParameters[["reference"]], distribution = "norm", sampleSize = nReference, distributionParameters = distriParam, irtModel = "1pl", logistic = FALSE) ## cutoffItems: Cut-off points set.seed(29834528) cutoffPoints <- CutoffIpr(quantiles = 0.95, itemParameters = nullParameters, itemCovariances = nullAse, irtModel = "1pl", logistic = FALSE, statistic = "ncdif", nReplicates = nReplicates) ## alternativeAse: Variance-covariance matrices under alternative hypothesis altAse <- list() altAse[["focal"]] <- AseIrt(itemParameters = raschParameters[["focal"]], distribution = "norm", sampleSize = nFocal, distributionParameters = distriParam, irtModel = "1pl", logistic = FALSE) altAse[["reference"]] <- AseIrt(itemParameters = raschParameters[["reference"]], distribution = "norm", sampleSize = nReference, distributionParameters = distriParam, irtModel = "1pl", logistic = FALSE) ## IPR sample under alternative hypothesis altIPR <- Ipr(itemParameters = raschParameters, itemCovariances = altAse, nReplicates = nReplicates) ## NCDIF sample under alternative altNcdif <- IprNcdif(itemParameterList = altIPR, irtModel = "1pl", logistic = FALSE) ## powerAPriori: Power calculation powerRasch <- rowMeans(altNcdif > cutoffPoints$quantiles) ## tablePower items <- seq(length(ncdif1pl)) tabPower <- cbind(Item = items, `Reference b` = raschParameters[["reference"]][items, ], `Focal b` = raschParameters[["focal"]][items, ], `True NCDIF` = sprintf("%.5f", ncdif1pl), Cutoff = sprintf("%.5f", cutoffPoints$quantiles), Power = sprintf("%.2f", powerRasch)) xtabPower <- xtable(tabPower) align(xtabPower) <- rep("r", length(digits(xtabPower))) caption(xtabPower) <- "NCDIF Power calculated for particular item difficulties under the 1PL model using the IPR approach." label(xtabPower) <- "tab:power" ## tabPowerPrint print(xtabPower, sanitize.colnames.function = function(x) x, include.rownames = FALSE) #, scalebox = '0.8')