Nothing
### calibrateType1.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: apr 5 2018 (10:23)
## Version:
## Last-Updated: Jan 11 2022 (16:51)
## By: Brice Ozenne
## Update #: 852
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
## * Documentation - calibrateType1
##' @title Simulation Study Assessing Bias and Type 1 Error
##' @description Perform a simulation study over one or several sample size
##' to assess the bias of the estimate
##' and the type 1 error of the Wald test and robust Wald test
##' @name calibrateType1
##'
##' @param object a \code{lvm} object defining the model to be fitted.
##' @param warmup [list of lvm] a list of \code{lvm} objects that will be sequentially fitted with for
##' starting values the parameter of the previous model in the list (if any). The parameters of the final
##' model of the list are used to initialize the fit of the model of interest (i.e. object).
##' @param param [character vector] names of the coefficient whose value will be tested.
##' @param null [numeric vector] vector of null hypotheses, one for each model coefficient.
##' By default a vector of 0.
##' @param n [integer vector, >0] sample size(s) considered in the simulation study.
##' @param n.rep [integer, >0] number of simulations per sample size.
##' @param cluster [integer vector] the grouping variable relative to which the observations are iid.
##' Will be passed to \code{lava::estimate}.
##' @param correction [logical] should the type 1 error after correction be computed?
##' @param generative.object [lvm] object defining the statistical model generating the data.
##' @param generative.coef [name numeric vector] values for the parameters of the generative model.
##' Can also be \code{NULL}: in such a case the coefficients are set to default values decided by lava (usually 0 or 1).
##' @param true.coef [name numeric vector] expected values for the parameters of the fitted model.
##' @param n.true [integer, >0] sample size at which the estimated coefficients will be a reliable approximation of the true coefficients.
##' @param round.true [integer, >0] the number of decimal places to be used for the true value of the coefficients. No rounding is done if \code{NULL}.
##' @param bootstrap [logical] should bootstrap resampling be performed?
##' @param n.bootstrap [integer, >0] the number of bootstrap sample to be used for each bootstrap.
##' @param checkType1 [logical] returns an error if the coefficients associated to the null hypotheses do not equal 0.
##' @param checkType2 [logical] returns an error if the coefficients associated to the null hypotheses equal 0.
##' @param dir.save [character] path to the directory were the results should be exported.
##' Can also be \code{NULL}: in such a case the results are not exported.
##' @param F.test [logical] should a multivariate Wald test be perform testing simultaneously all the null hypotheses?
##' @param label.file [character] element to include in the file name.
##' @param seed [integer, >0] value that will be set before adjustment for multiple comparisons to ensure reproducible results.
##' Can also be \code{NULL}: in such a case no seed is set.
##' @param cpus [integer >0] the number of processors to use.
##' If greater than 1, the simulations are performed in parallel.
##' @param trace [integer] should the execution of the function be trace. Can be 0, 1 or 2.
##' @param ... [internal] Only used by the generic method.
##'
##' @return An object of class \code{calibrateType1}.
##' @seealso \code{link{autoplot.calibrateType1}} for a graphical display of the bias or of the type 1 error.
##'
##' @author Brice Ozenne
##'
##' @examples
##' \dontrun{
##' #### simulate data ####
##' m.Sim <- lvm(c(Y1[mu1:sigma]~1*eta,
##' Y2[mu2:sigma]~1*eta,
##' Y3[mu3:sigma]~1*eta,
##' eta~beta1*Group+beta2*Gender))
##' latent(m.Sim) <- ~eta
##' categorical(m.Sim, labels = c("M","F")) <- ~Gender
##'
##' d <- lava::sim(m.Sim, 1e2)
##'
##' #### calibrate type 1 error on the estimated model ####
##' m <- lvm(Y1~eta,
##' Y2~eta,
##' Y3~eta,
##' eta~Group+Gender)
##' e <- lava::estimate(m, data = d)
##' res <- calibrateType1(e, param = "eta~Group", n.rep = 100)
##' res <- calibrateType1(e, param = c("eta~Group","Y1~eta"), F.test = TRUE, n.rep = 100)
##' res <- calibrateType1(e, param = "eta~Group", n.rep = 100, cpus = 4)
##' summary(res)
##' }
##'
##' @export
`calibrateType1` <-
function(object, param, n.rep, ...) UseMethod("calibrateType1")
## * calibrateType1.lvm
##' @rdname calibrateType1
##' @export
calibrateType1.lvm <- function(object, param, n.rep, n, correction = TRUE, warmup = NULL, null = NULL,
F.test = FALSE, cluster = NULL,
generative.object = NULL, generative.coef = NULL,
true.coef = NULL, n.true = 1e6, round.true = 2,
bootstrap = FALSE, n.bootstrap = 1e3,
checkType1 = FALSE, checkType2 = FALSE,
dir.save = NULL, label.file = NULL,
seed = NULL, cpus = 1, trace = 2, ...){
### ** test
if(cpus>1){
test.package <- try(requireNamespace("foreach"), silent = TRUE)
if(inherits(test.package,"try-error")){
stop("There is no package \'foreach\' \n",
"This package is necessary when argument \'cpus\' is greater than 1 \n")
}
if(cpus > parallel::detectCores()){
stop("Argument \'cpus\' is greater than the number of available CPU cores \n",
"available CPU cores: ",parallel::detectCores(),"\n")
}
}
### ** prepare
n.n <- length(n)
## *** generative model
if(is.null(generative.object)){
generative.object <- object
}
if(any(lava::manifest(object) %in% lava::manifest(generative.object) == FALSE)){
missingVar <- lava::manifest(object)[lava::manifest(object) %in% lava::manifest(generative.object) == FALSE]
stop("The object contains manifest variables that are not in the generative model \n",
"missing manifest variables: \"",paste0(missingVar, collapse = "\" \""),"\" \n")
}
df.type_generative <- coefType(generative.object, as.lava = FALSE)
name.param_generative <- df.type_generative[!is.na(df.type_generative$lava),"param"]
n.param_generative <- length(name.param_generative)
## *** coef generative
if(!is.null(generative.coef) && any(names(generative.coef) %in% name.param_generative == FALSE)){
extraParam <- names(generative.coef)[names(generative.coef) %in% name.param_generative == FALSE]
stop("Invalid argument \'generative.coef\': some of the coefficient names do not match those of the generative model \n",
"extra coefficients: \"",paste0(extraParam, collapse = "\" \""),"\"\n")
}
## *** coef of the fitted model
if(is.null(true.coef)){
if(trace>1){
cat(" Estimate true coefficients using a sample size of n=",n.true," ", sep="")
}
e.true <- lava::estimate(object,
data = lava::sim(generative.object, n = n.true, p = generative.coef, latent = FALSE))
coef.true <- coef(e.true)
if(!is.null(round.true)){
coef.true <- round(coef.true, digits = round.true)
}
if(trace>1){
cat("- done \n")
}
}else{
if(trace>1){
cat(" Check true coefficients ")
}
n.true <- n[n.n]
e.true <- suppressWarnings(lava::estimate(object, cluster = cluster,
data = lava::sim(generative.object, n = n.true, p = generative.coef, latent = FALSE)))
name.test <- names(coef(e.true))
if(e.true$opt$convergence>0 && trace>1){
cat("(incorrect convergence of the model) ")
}
if(!identical(sort(name.test),sort(names(true.coef)))){
extraNames <- setdiff(names(true.coef),name.test)
missingNames <- setdiff(name.test,names(true.coef))
stop("Names of the coefficients in argument \'true.coef\' does not matches those of the estimated coefficients \n",
"missing names: \"",paste0(missingNames, collapse = "\" \""),"\" \n",
"extra names: \"",paste0(extraNames, collapse = "\" \""),"\" \n")
}
coef.true <- true.coef
if(trace>1){
cat("- done \n")
}
}
name.coef <- names(coef.true)
n.coef <- length(name.coef)
## *** type of the coef of the fitted model
df.type <- coefType(e.true, as.lava = FALSE)
df.type <- df.type[df.type$name %in% name.coef,]
type.coef <- stats::setNames(df.type$detail, df.type$name)
## *** null hypothesis
n.param <- length(param)
if(any(param %in% name.coef == FALSE)){
incorrect.name <- param[param %in% name.coef == FALSE]
possible.name <- setdiff(name.coef, param)
ls.name <- lapply(incorrect.name, function(iN){
dist.tempo <- utils::adist(x = iN, y = possible.name)
return(possible.name[which.min(dist.tempo)])
})
ex.name <- unique(unlist(ls.name))
stop("Invalid argument \'param\': some of the coefficient names does not match those of the estimate model \n",
"incorrect names: \"",paste(incorrect.name, collapse = "\" \""),"\" \n",
"example of valid names: \"",paste(ex.name, collapse = "\" \""),"\"\n")
}
if(checkType1 && any(coef.true[param]!=0)){
txtCoef <- paste(param[coef.true[param]!=0], collapse = "\" \"")
stop("Control type 1 error: coefficients \"",txtCoef,"\" are not 0 while their belong to the param hypothesis\n")
}
if(checkType2 && any(coef.true[param]==0)){
txtCoef <- paste(param[coef.true[param]==0], collapse = "\" \"")
stop("Control type 2 error: coefficients \"",txtCoef,"\" are 0 while their belong to the param hypothesis\n")
}
res.C <- .createContrast(param,
name.param = name.coef,
add.rowname = TRUE,
rowname.rhs = FALSE)
contrast <- res.C$contrast
if(is.null(null)){
rhs <- res.C$null
}else{
if(length(null)!=n.param){
stop("Argument \'null\' have the same length as argument \'param\' \n")
}
if(!is.null(names(null)) && any(names(null)!=param)){
stop("If named, argument \'null\' must have for name argument \'param\' \n")
}
rhs <- null
}
## *** warmup
if(!is.null(warmup) && inherits(warmup,"lvm")){
stop("Argument \'warmup\' must be a list of lvm objects \n")
}
### ** display
if(trace>1){
cat(" Settings: \n")
cat(" > simulation for n=",paste(n,collapse = " "),"\n",sep="")
cat(" > model: \n")
print(object)
cat(" > expected coefficients: \n")
print(coef.true)
cat(" > null hypotheses: \n")
print(rhs)
if(bootstrap){
cat(" > bootstrap: ",bootstrap,"\n")
}
if(!is.null(seed)){
cat(" > seed: ",seed,"\n")
}
}
### ** loop
store.coef <- paste0(rownames(contrast)," == ",rhs)
if(F.test){
store.coef <- c(store.coef, "global")
}
n.store <- length(store.coef)
if(trace>1){
cat("\n")
cat(" Perform simulation: \n")
}
if(cpus>1){
## *** define cluster
if(trace>0){
cl <- suppressMessages(parallel::makeCluster(cpus, outfile = ""))
pb <- utils::txtProgressBar(max = n.rep, style = 3)
}else{
cl <- parallel::makeCluster(cpus)
}
## *** link to foreach
doParallel::registerDoParallel(cl)
## *** export package
parallel::clusterCall(cl, fun = function(x){
suppressPackageStartupMessages(requireNamespace("lava", quietly = TRUE))
suppressPackageStartupMessages(requireNamespace("lavaSearch2", quietly = TRUE))
})
## *** seed
cpus.name <- unlist(parallel::clusterApply(cl = cl, 1:cpus, function(x){
myName <- paste(Sys.info()[["nodename"]], Sys.getpid(), sep="-")
return(myName)
}))
if(length(seed)==0){
seed <- rep(as.numeric(NA), cpus)
}else if(length(seed)==1){
set.seed(seed)
seed <- sample(1:1e5,size=cpus,replace=FALSE)
}else{
if(length(seed)!=cpus){
stop("Length of argument \'seed\' does not match argument \'cpus\' \n")
}
}
names(seed) <- cpus.name
## *** parallel computation
toExport <- c(".warperType1", "cpus.name")
iRep <- NULL # [:for CRAN check] foreach
resSim <- foreach::`%dopar%`(
foreach::foreach(iRep = 1:n.rep,
.export = toExport,
.packages = c("lavaSearch2")),{ # iRep <- 1
if(trace>0){utils::setTxtProgressBar(pb, iRep)}
myName <- paste(Sys.info()[['nodename']], Sys.getpid(), sep='-')
iSeed <- as.double(seed[myName])
ls.pvalue <- vector(mode = "list", length = n.n)
ls.estimate <- vector(mode = "list", length = n.n)
iIndex <- 1
for(iN in n){ # iN <- n[1]
iRes <- .warperType1(iRep,
n = iN,
correction = correction,
generative.object = generative.object,
generative.coef = generative.coef,
object = object,
warmup = warmup,
cluster = cluster,
coef.true = coef.true,
type.coef = type.coef,
name.coef = name.coef,
store.coef = store.coef,
n.coef = n.coef,
n.store = n.store,
F.test = F.test,
param = param,
contrast = contrast,
rhs = rhs,
bootstrap = bootstrap,
n.bootstrap = n.bootstrap,
seed = iSeed,
...)
ls.pvalue[[iIndex]] <- iRes$pvalue
ls.estimate[[iIndex]] <- iRes$estimate
iIndex <- iIndex + 1
}
return(list(pvalue = do.call("rbind",ls.pvalue),
estimate = do.call("rbind",ls.estimate)))
})
parallel::stopCluster(cl)
if(trace>0){close(pb)}
## *** post process
dt.pvalue <- do.call("rbind",lapply(resSim,"[[","pvalue"))
dt.pvalue <- dt.pvalue[order(dt.pvalue$n,dt.pvalue$rep),,drop=FALSE]
dt.estimate <- do.call("rbind",lapply(resSim,"[[","estimate"))
dt.estimate <- dt.estimate[order(dt.estimate$n,dt.estimate$rep),,drop=FALSE]
}else{
## *** filename
if(is.null(label.file)){label.file <- seed}
filename_tempo.pvalue <- paste0("type1error-S",label.file,"(tempo).rds")
filename_tempo.estimate <- paste0("estimate-S",label.file,"(tempo).rds")
filename.pvalue <- gsub("\\(tempo\\)","",filename_tempo.pvalue)
filename.estimate <- gsub("\\(tempo\\)","",filename_tempo.estimate)
if(!is.null(dir.save)){
if (dir.exists(dir.save) == FALSE) {
stop("Argument \'dir.save\' does not lead to an existing directory \n")
}
}
if(!is.null(dir.save)){
cat(" > export results in ",dir.save,"\n")
}
if(trace>0){
test.package <- try(requireNamespace("pbapply"), silent = TRUE)
if(inherits(test.package,"try-error")){
stop("There is no package \'pbapply\' \n",
"This package is necessary when argument \'trace\' is TRUE \n")
}
FCTapply <- pbapply::pblapply
}else{
FCTapply <- lapply
}
if(!is.null(seed)){
set.seed(seed)
}else{
seed <- as.numeric(NA)
}
## *** sequential simulation
dt.pvalue <- NULL
dt.estimate <- NULL
for(iN in n){ ## iN <- n[1]
if(trace>0){
cat(" > sample size=",iN,"\n", sep = "")
}
resSim <- do.call(FCTapply, args = list(X = 1:n.rep, FUN = function(iRep){
.warperType1(iRep,
n = iN,
correction = correction,
generative.object = generative.object,
generative.coef = generative.coef,
object = object, warmup = warmup, cluster = cluster,
coef.true = coef.true, type.coef = type.coef, name.coef = name.coef,
store.coef = store.coef, n.coef = n.coef, n.store = n.store,
F.test = F.test, param = param, contrast = contrast, rhs = rhs,
bootstrap = bootstrap,
n.bootstrap = n.bootstrap,
seed = seed,
...)
}))
dt.pvalue <- rbind(dt.pvalue,
do.call("rbind",lapply(resSim,"[[","pvalue"))
)
dt.estimate <- rbind(dt.estimate,
do.call("rbind",lapply(resSim,"[[","estimate")))
## export (tempo)
if(!is.null(dir.save)){
saveRDS(dt.pvalue, file = file.path(dir.save,filename_tempo.pvalue))
saveRDS(dt.estimate, file = file.path(dir.save,filename_tempo.estimate))
}
}
}
## ** export
if(!is.null(dir.save)){
saveRDS(dt.pvalue, file = file.path(dir.save,filename.pvalue))
saveRDS(dt.estimate, file = file.path(dir.save,filename.estimate))
}
out <- list(estimate = dt.estimate,
p.value = dt.pvalue,
e.true = e.true,
param = param)
class(out) <- append("calibrateType1",class(out))
return(out)
}
## * calibrateType1.lvmfit
##' @rdname calibrateType1
##' @export
calibrateType1.lvmfit <- function(object, param, n.rep, correction = TRUE, F.test = FALSE,
bootstrap = FALSE, n.bootstrap = 1e3,
seed = NULL, trace = 2, cpus = 1, ...){
## ** Prepare
## *** model
object.model <- object$model
## *** coef
coef.true <- coef(object)
name.coef <- names(coef.true)
if(any(param %in% name.coef == FALSE)){
txt <- param[param %in% name.coef == FALSE]
txt2 <- setdiff(name.coef, param)
stop("Argument \'param\' does not match the names of the model coefficients \n",
"Incorrect param: \"",paste(txt, collapse = "\" \""),"\" \n",
"Possible param: \"",paste(txt2, collapse = "\" \""),"\" \n")
}
coef.true[param] <- 0
## *** data
n <- object$data$n
## ** Run
out <- calibrateType1(object.model,
param = param,
n.rep = n.rep,
n = n,
correction = correction,
F.test = F.test,
generative.object = object.model,
generative.coef = coef.true,
true.coef = coef.true,
bootstrap = bootstrap,
n.bootstrap = n.bootstrap,
checkType1 = FALSE,
checkType2 = FALSE,
dir.save = NULL,
label.file = NULL,
seed = seed,
cpus = cpus,
trace = trace,
...)
## ** Export
return(out)
}
## * .warperType1
.warperType1 <- function(iRep, n, correction, generative.object, generative.coef,
object, warmup, cluster,
coef.true, type.coef, name.coef, store.coef, n.coef, n.store,
F.test, param, contrast, rhs,
bootstrap, n.bootstrap,
seed,
...){
ls.iE <- list(estimate.truth = coef.true) ## temporary
ls.iP <- list() ## temporary
out <- list()
dots <- list(...)
if("control" %in% names(dots) == FALSE){
dots$control <- list()
}
## ** simulation
dt.sim <- lava::sim(generative.object, n = n, p = generative.coef, latent = FALSE)
## ** initialisation
if(!is.null(warmup)){
n.warmup <- length(warmup)
for(iW in 1:n.warmup){ ## iW <- 1
iE.lvm <- suppressWarnings(do.call(lava::estimate, args = c(list(warmup[[iW]], data = dt.sim, cluster = cluster), dots)))
if(("convergence" %in% names(iE.lvm$opt)) && (iE.lvm$opt$convergence==1)){return(list(pvalue=NULL,estimate=NULL))} ## exclude lvm that has not converged
dots$control$start <- coef(iE.lvm)
}
}
## ** model fit
e.lvm <- suppressWarnings(do.call(lava::estimate, args = c(list(object, data = dt.sim, cluster = cluster), dots)))
## check correct convergence of the latent variable model
if(("convergence" %in% names(e.lvm$opt)) && (e.lvm$opt$convergence==1)){return(list(pvalue=NULL,estimate=NULL))} ## exclude lvm that has not converged
if(any(eigen(getVarCov2(e.lvm, ssc = FALSE, df = FALSE))$values<=0)){return(list(pvalue=NULL,estimate=NULL))} ## exclude lvm where the residual covariance matrix is not semipositive definite
ratio_sd_beta <- sqrt(diag(vcov(e.lvm)))/(abs(coef(e.lvm))+1)
if(max(stats::na.omit(ratio_sd_beta))>1e3){return(list(pvalue=NULL,estimate=NULL))} ## exclude if standard error much larger than coefficient
if(inherits(suppressWarnings(try(summary(e.lvm)$coef, silent = TRUE)), "try-error")){return(list(pvalue=NULL,estimate=NULL))} ## exclude lvm where we cannot compute the summary
## ** corrections
e.lvm <- estimate2(e.lvm, df = "none", ssc = "none")
if(correction){
e.lvm.Satt <- try(estimate2(e.lvm, df = "satterthwaite", ssc = "none"), silent = TRUE)
e.lvm.KR <- try(suppressWarnings(estimate2(e.lvm, df = "satterthwaite", ssc = "residuals")), silent = TRUE)
}else{
e.lvm.Satt <- 1
class(e.lvm.Satt) <- "try-error"
e.lvm.KR <- 1
class(e.lvm.KR) <- "try-error"
}
## ** extract p.values
eS.ML <- summary2(e.lvm, robust = FALSE)$coef
F.ML <- summary(compare2(e.lvm, robust = FALSE,
linfct = contrast, rhs = rhs, F.test = F.test, as.lava = FALSE))$table2
eS.robustML <- summary2(e.lvm, robust = TRUE)$coef
F.robustML <- summary(compare2(e.lvm, robust = TRUE,
linfct = contrast, rhs = rhs, F.test = F.test, as.lava = FALSE))$table2
if(!inherits(e.lvm.Satt,"try-error")){
eS.Satt <- summary2(e.lvm.Satt, robust = FALSE)$coef
F.Satt <- summary(compare2(e.lvm.Satt, robust = FALSE,
linfct = contrast, rhs = rhs, F.test = F.test,
as.lava = FALSE))$table2
eS.robustSatt <- summary2(e.lvm.Satt, robust = TRUE)$coef
F.robustSatt <- summary(compare2(e.lvm.Satt, robust = TRUE,
linfct = contrast, rhs = rhs, F.test = F.test,
as.lava = FALSE))$table2
}
if(!inherits(e.lvm.KR,"try-error")){
eS.SSC <- summary2(e.lvm.KR, robust = FALSE)$coef
F.SSC <- summary(compare2(e.lvm.KR, robust = FALSE,
linfct = contrast, rhs = rhs, F.test = F.test,
as.lava = FALSE))$table2
eS.robustSSC <- summary2(e.lvm.KR, robust = TRUE)$coef
F.robustSSC <- summary(compare2(e.lvm.KR, robust = TRUE,
linfct = contrast, rhs = rhs, F.test = F.test,
as.lava = FALSE))$table2
eS.KR <- summary2(e.lvm.KR, robust = FALSE)$coef
F.KR <- summary(compare2(e.lvm.KR, robust = FALSE,
linfct = contrast, rhs = rhs, F.test = F.test,
as.lava = FALSE))$table2
eS.robustKR <- summary2(e.lvm.KR, robust = TRUE)$coef
F.robustKR <- summary(compare2(e.lvm.KR, robust = TRUE,
linfct = contrast, rhs = rhs, F.test = F.test,
as.lava = FALSE))$table2
}
## ** store
## *** estimates
ls.iE$estimate.ML <- eS.ML[name.coef,"Estimate"]
if(!inherits(e.lvm.KR,"try-error")){
ls.iE$estimate.MLcorrected <- eS.KR[name.coef,"Estimate"]
}else{
ls.iE$estimate.MLcorrected <- rep(as.numeric(NA), n.coef)
}
## *** standard errors
if(is.null(cluster)){
ls.iE$se.ML <- eS.ML[name.coef,"Std. Error"]
if(!inherits(e.lvm.KR,"try-error")){
ls.iE$se.MLcorrected <- eS.KR[name.coef,"Std. Error"]
}else{
ls.iE$se.MLcorrected <- rep(as.numeric(NA), n.coef)
}
ls.iE$se.robustML <- eS.robustML[name.coef,"robust SE"]
if(!inherits(e.lvm.KR,"try-error")){
ls.iE$se.robustMLcorrected <- eS.robustKR[name.coef,"robust SE"]
}else{
ls.iE$se.robustMLcorrected <- rep(as.numeric(NA), n.coef)
}
}
## *** degree of freedom
if(!inherits(e.lvm.Satt,"try-error")){
if(is.null(cluster)){
ls.iE$df.ML <- eS.Satt[name.coef,"df"]
}
ls.iE$df.robustML <- eS.robustSatt[name.coef,"df"]
}
if(!inherits(e.lvm.KR,"try-error")){
if(is.null(cluster)){
ls.iE$df.MLcorrected <- eS.KR[name.coef,"df"]
}
ls.iE$df.robustMLcorrected <- eS.robustKR[name.coef,"df"]
}
## *** p-value
if(is.null(cluster)){
ls.iP$p.Ztest <- F.ML[store.coef,"p.value"]
if(!inherits(e.lvm.Satt,"try-error")){
ls.iP$p.Satt <- F.Satt[store.coef,"p.value"]
}
if(!inherits(e.lvm.KR,"try-error")){
ls.iP$p.SSC <- F.SSC[store.coef,"p.value"]
ls.iP$p.KR <- F.KR[store.coef,"p.value"]
}
}
ls.iP$p.robustZtest <- F.robustML[store.coef,"p.value"]
if(!inherits(e.lvm.Satt,"try-error")){
ls.iP$p.robustSatt <- F.robustSatt[store.coef,"p.value"]
}
if(!inherits(e.lvm.KR,"try-error")){
ls.iP$p.robustSSC <- F.robustSSC[store.coef,"p.value"]
ls.iP$p.robustKR <- F.robustKR[store.coef,"p.value"]
}
## *** niter.correct / warning
if(!inherits(e.lvm.KR,"try-error")){
test.warning <- inherits(attr(e.lvm.KR$sCorrect,"warning"),"try-error")
niter.correct <- as.double(e.lvm.KR$sCorrect$ssc$iter)
}else{
test.warning <- NA
niter.correct <- as.double(NA)
}
## ** bootstrap
if(bootstrap>0){
e.boot <- eval(parse(text = "butils::bootReg(e.lvm, type = \"coef\", n.boot = n.bootstrap"))
index.coef.boot <- match(param, name.coef)
boot.perc <- summary(e.boot, p.value = TRUE, type = "perc", print = FALSE, index = index.coef.boot)
boot.stud <- summary(e.boot, p.value = TRUE, type = "stud", print = FALSE, index = index.coef.boot)
boot.bca <- summary(e.boot, p.value = TRUE, type = "bca", print = FALSE, index = index.coef.boot)
ls.iP$p.bootPerc <- boot.perc[store.coef,"p.value"]
ls.iP$p.bootStud <- boot.stud[store.coef,"p.value"]
ls.iP$p.bootBca <- boot.bca[store.coef,"p.value"]
}else{
n.bootstrap <- as.numeric(NA)
}
## ** collect results
## estimates
## use rep to avoid warning
## return(list(n = n,
## rep = iRep,
## seed = seed,
## ninter = niter.correct,
## warning = test.warning,
## name = names(coef.true),
## type = type.coef[name.coef]))
df1 <- data.frame(n = n,
rep = iRep, n.coef,
seed = seed, n.coef,
niter = niter.correct,
warning = test.warning,
name = names(coef.true),
type = unname(type.coef[name.coef]),
stringsAsFactors = FALSE)
rownames(df1) <- NULL
df2 <- do.call(cbind,ls.iE)
rownames(df2) <- NULL
out$estimate <- cbind(df1, df2)
## p.value
df1 <- data.frame(n = n,
rep = iRep,
seed = seed,
nboot = n.bootstrap,
niter = niter.correct,
warning = test.warning,
link = unname(store.coef),
stringsAsFactors = FALSE)
rownames(df1) <- NULL
df2 <- do.call(cbind,ls.iP)
rownames(df2) <- NULL
out$pvalue <- cbind(df1, df2)
## ** export
return(out)
}
######################################################################
### calibrateType1.R ends here
## * .mergeCalibrateType1
.mergeCalibrateType1 <- function(list){
if(any(sapply(list,function(iL){inherits(iL,"calibrateType1")})==FALSE)){
stop("Argument \'list\' must only contain \"calibrateType1\" objects \n")
}
newObject <- list[[1]]
newObject$estimate <- do.call(rbind,lapply(list,"[[","estimate"))
newObject$p.value <- do.call(rbind,lapply(list,"[[","p.value"))
return(newObject)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.