Nothing
imputeCA <-function (X, ncp = 2, threshold = 1e-08, maxiter = 1000, row.sup=NULL, col.sup=NULL, quanti.sup=NULL, quali.sup=NULL)
{
shrinkCA <- function(X, ncp = 2, row.sup=NULL, col.sup=NULL) {
X[row.sup,] <- X[row.sup,]*1e-08
X[,col.sup] <- X[,col.sup]*1e-08
P <- as.matrix(X/sum(X))
Rc <- apply(P, 2, sum)
Rr <- apply(P, 1, sum)
S <- t(t((P - Rr %*% t(Rc))/sqrt(Rr))/sqrt(Rc))
svdRes <- svd(S)
n <- nrow(X)-length(row.sup)
p <- ncol(X)-length(col.sup)
sigma2 <- sum(svdRes$d[-c(1:ncp)]^2)/((n - 1) * (p - 1) - (n - 1) * ncp - (p - 1) * ncp + ncp^2)
lambda.shrinked <- (svdRes$d[1:ncp]^2 - n * (p/min(p, (n - 1))) * sigma2)/svdRes$d[1:ncp]
if (ncp == 1) recon <- (svdRes$u[, 1] * lambda.shrinked) %*% t(svdRes$v[, 1])
else recon <- svdRes$u[, 1:ncp] %*% (t(svdRes$v[, 1:ncp]) * lambda.shrinked)
recon <- sum(X) * (t(t(recon * sqrt(Rr)) * sqrt(Rc)) + Rr %*% t(Rc))
recon[row.sup,] <- recon[row.sup,]*1e+08
recon[,col.sup] <- recon[,col.sup]*1e+08
rownames(recon) <- rownames(X)
colnames(recon) <- colnames(X)
res <- list(recon = recon)
return(res)
}
### main
if (!is.null(quanti.sup) | !is.null(quali.sup)){
Xtot <- X
if (!is.null(col.sup)) IndiceColSup <- colnames(Xtot)[col.sup]
X <- X[,-c(quali.sup,quanti.sup)]
if (!is.null(col.sup)) col.sup <- which(colnames(X)%in%IndiceColSup)
}
X <- as.matrix(X)
if (sum(is.na(X)) == 0) stop("No value is missing")
missing <- which(is.na(X))
Xhat <- X
Xhat[missing] <- sample(X[-missing], length(missing), replace = TRUE) + 1
nb.iter <- 1
old <- Inf
objective <- 0
recon <- Xhat
while (nb.iter > 0) {
Xhat[missing] <- recon[missing]
Xhat[missing][which((Xhat[missing]<=0))]=0
RXhat <- rowSums(Xhat)
CXhat <- colSums(Xhat)
if ((sum(RXhat > 1e-06) == nrow(Xhat)) & (sum(CXhat > 1e-06) == ncol(Xhat))) recon <- shrinkCA(Xhat, ncp = ncp, row.sup=row.sup, col.sup=col.sup)$recon
diff <- Xhat - recon
diff[missing] <- 0
objective <- sum((diff^2))
criterion <- abs(1 - objective/old)
old <- objective
nb.iter <- nb.iter + 1
if (!is.nan(criterion)) {
if ((criterion < threshold) && (nb.iter > 5))
nb.iter <- 0
if ((objective < threshold) && (nb.iter > 5))
nb.iter <- 0
}
if (nb.iter > maxiter) {
nb.iter <- 0
warning(paste("Stopped after ", maxiter, " iterations"))
}
}
completeObs <- X
completeObs[missing] <- Xhat[missing]
if (!is.null(quali.sup) | !is.null(quanti.sup)) {
Xtot[,-c(quali.sup,quanti.sup)] <- completeObs
completeObs <- Xtot
}
return(completeObs)
}
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.