Nothing
#' Fitting Linear Regression Model with Missing Values
#'
#' This function is used inside \code{miss.lm} to fit linear regression model with missing values, by EM algorithm.
#' @param x design matrix with missingness \eqn{N \times p}{N * p}.
#' @param y response vector \eqn{N \times 1}{N * 1}.
#' @param control a list of parameters for controlling the fitting process. For \code{miss.lm.fit} this is passed to \code{\link{miss.lm.control}}.
#' @return a list with following components:
#' \item{coefficients}{Estimated \eqn{\beta}{\beta}.}
#' \item{ll}{Observed log-likelihood.}
#' \item{s.resid}{Estimated standard error for residuals.}
#' \item{s.err}{Standard error for estimated parameters.}
#' \item{mu.X}{Estimated \eqn{\mu}{\mu}.}
#' \item{Sig.X}{Estimated \eqn{\Sigma}{\Sigma}.}
#' @import norm MASS
#' @importFrom methods is
#' @examples
#' ## For examples see example(miss.lm)
miss.lm.fit <- function (x, y,
control = list())
{
control <- do.call("miss.lm.control", control)
maxruns = control$maxruns
tol_em = control$tol_em
print_iter = control$print_iter
xnames <- dimnames(x)[[2L]]
x = as.matrix(x[,-1])
p <- ncol(x)
#delete rows completely missing
if(any(apply(is.na(x),1,sum)==p)){
i_allNA=which(apply(is.na(x),1,sum)==p)
x = x[-i_allNA,]
y = y[-i_allNA]
}
if(any((is.na(y))==TRUE)){
i_YNA=which(is.na(y)==TRUE)
x = x[-i_YNA,]
y = y[-i_YNA]
}
# Check parameter consistency
if (!is(x, "matrix")){
x <- as.matrix(x)
}
if (sum(sapply(x, is.numeric)) < ncol(x) ||
sum(sapply(y, is.numeric)) < length(y)) {
stop("Error: parameters 'x' and 'y' should be numeric.")
}
# Estimate extended mean and covariance using EM
n <- nrow(x)
p <- ncol(x) + 1
s <- prelim.norm(cbind(y, x))
thetahat <- em.norm(s, showits=print_iter, maxits=maxruns, criterion=tol_em)
pars <- getparam.norm(s, thetahat)
ll <- loglik.norm(s, thetahat)
ptm <- Sys.time()
# Calculate regression estimates
b.est <- c(pars$mu[1] -
pars$sigma[1,2:p] %*% solve(pars$sigma[2:p,2:p]) %*% pars$mu[2:p],
pars$sigma[1,2:p] %*% solve(pars$sigma[2:p,2:p]))
esigma.est <- as.vector(sqrt(
pars$sigma[1,1] - b.est[2:p] %*% pars$sigma[2:p,2:p] %*% b.est[2:p]))
Gram.est <- rbind(rep(0, p), cbind(rep(0, p - 1), pars$sigma[2:p,2:p])) +
c(1, pars$mu[2:p]) %*% t(c(1, pars$mu[2:p]))
sdb.est <- sqrt(diag(esigma.est^2 * solve(Gram.est * n)))
time_run=Sys.time() - ptm
names(b.est) <- xnames
names(sdb.est) <- xnames
list(coefficients = b.est, s.err = sdb.est, s.resid = esigma.est, ll= ll, Sig.X = pars$sigma, mu.X = pars$mu)
}
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.