I'm trying to understand this code from R and I am running the code step by step. When I simply call the function, the code works perfectly but if I do it step by step, there is an error in the line with red color. I don't understand the usage of dots before functions in R. The question is about this kind of lines x_sort <- .sort.missing(x, x_nonmiss). Thank you.
> GSE
function (x, tol = 1e-04, maxiter = 150, method = c("bisquare",
"rocke"), init = c("emve", "qc", "huber", "imputed", "emve_c"),
mu0, S0, ...)
{
xcall <- match.call()
method <- match.arg(method)
init <- match.arg(init)
if (is.data.frame(x) | is.matrix(x))
x <- data.matrix(x)
else stop("Data matrix must be of class matrix or data.frame")
p <- ncol(x)
if (p > 200 | p < 2)
stop("Column dimension of 'x' must be in between 2 and 200.")
x_nonmiss <- is.na(x) * -1 + 1
pp <- rowSums(x_nonmiss)
pp_col <- colSums(x_nonmiss)
if (all(pp == 0))
stop("All observations have missing values!")
if (any(pp_col == 0))
stop("Data matrix cannot contain column(s) with completely missing data!")
ok <- which(pp > 0)
x_orig <- x
x <- x[ok, ]
x_nonmiss <- x_nonmiss[ok, ]
x_sort <- .sort.missing(x, x_nonmiss)
EM.mu0 <- colMeans(x_sort$x, na.rm = T)
EM.S0 <- diag(apply(x_sort$x, 2, var, na.rm = T))
x_sort <- c(x_sort, .CovEM.setparam(p, EM.mu0, EM.S0))
n <- nrow(x_sort$x)
p <- ncol(x_sort$x)
if (n <= p + 1)
stop(if (n <= p)
"n <= p -- you can't be serious!"
else "n == p+1 is too small sample size")
if (n < 2 * p)
warning("n < 2 * p, i.e., possibly too small sample size")
if (xor(missing(mu0), missing(S0)))
warning("Both 'mu0' and 'S0' must be provided. Default 'init' is used...")
if (missing(mu0) || missing(S0)) {
init.res <- switch(init, emve = {
emve(x, sampling = "uniform", ...)
}, qc = {
HuberPairwise(x, psi = "sign", computePmd = FALSE)
}, huber = {
HuberPairwise(x, psi = "huber", computePmd = FALSE,
...)
}, imputed = {
ImpS(x, ...)
}, emve_c = {
emve(x, sampling = "cluster", ...)
})
S0 <- init.res@S
mu0 <- init.res@mu
}
S0.chol <- tryCatch(chol(S0), error = function(e) NA)
if (!is.matrix(S0.chol))
stop("Estimated initial covariance matrix 'S0' is not positive definite.")
bdp <- 0.5
print.step <- 0
tol.scale <- 1e-09
miter.scale <- 300
res <- with(x_sort, .GSE.init(x, x_nonmiss, bdp, pu, n, p,
miss.group.unique, miss.group.counts, mu0, S0, tol, maxiter,
tol.scale, miter.scale, print.step = print.step, method))
pmd <- pmd.adj <- rep(NA, nrow(x_orig))
pmd[ok] <- res$pmd[x_sort$id.ro]
pmd.adj[ok] <- res$pmd.adj[x_sort$id.ro]
wgts <- rep(NA, nrow(x_orig))
wgts[ok] <- res$weights[x_sort$id.ro]
wgtsp <- rep(NA, nrow(x_orig))
wgtsp[ok] <- res$weightsprm[x_sort$id.ro]
ximp <- matrix(NA, nrow(x_orig), ncol(x_orig))
ximp[ok, ] <- res$ximp[x_sort$id.ro, ]
res <- new("GSE", call = xcall, S = res$S, mu = res$mu, sc = res$stilde0,
mu0 = mu0, S0 = S0, iter = res$iter, eps = res$ep, estimator = "Generalized S-Estimator",
x = x_orig, ximp = ximp, weights = wgts, weightsp = wgtsp,
pmd = pmd, pmd.adj = pmd.adj, p = p, pu = pp)
res
}
This is another example with the same situation.
> partial.mahalanobis
function (x, mu, Sigma)
{
xcall <- match.call()
if (is.data.frame(x) | is.matrix(x))
x <- data.matrix(x)
else stop("Data matrix must be of class matrix or data.frame")
n <- nrow(x)
p <- ncol(x)
x_nonmiss <- is.na(x) * -1 + 1
pp <- rowSums(x_nonmiss)
pp_col <- colSums(x_nonmiss)
if (all(pp == 0))
stop("All observations have missing values!")
if (any(pp_col == 0))
stop("Data matrix cannot contain column(s) with completely missing data!")
ok <- which(pp > 0)
not.ok <- which(pp == 0)
x.orig <- x
x <- x[ok, ]
x_nonmiss <- x_nonmiss[ok, ]
x_sort <- .sort.missing(x, x_nonmiss)
if (any(x_nonmiss == 0)) {
x_cent <- sweep(x_sort$x, 2, mu, "-")
pmd.tmp <- .partial.mahalanobis.Rcpp(x_cent, Sigma, x_sort$miss.group.unique,
x_sort$miss.group.counts)
pmd.tmp <- pmd.tmp[x_sort$id.ro]
}
else {
pmd.tmp <- mahalanobis(x_sort$x, mu, Sigma)
pmd.tmp <- pmd.tmp[x_sort$id.ro]
}
pmd <- rep(NA, nrow(x.orig))
pmd[ok] <- pmd.tmp
pu <- rowSums(!is.na(x.orig))
pmd.adj <- qchisq(pchisq(pmd, df = pu, log.p = T, lower.tail = F),
df = p, log.p = T, lower.tail = F)
pmd.adj[which(pu == p)] <- pmd[which(pu == p)]
new("CovRobMiss", mu = mu, S = Sigma, call = xcall, estimator = "unknown",
x = x.orig, pmd = pmd, pmd.adj = pmd.adj, p = p, pu = pu)
}
<environment: namespace:GSE>
