# R/estIndepModel.R In combi: Compositional omics model based visual integration

#### Documented in estIndepModel

#' Estimate the independence model belonging to one view
#'
#' @param data a list of data matrices with the same number of samples n in the rows.
#' Also phyloseq objects are acceptable
#' @param distribution a character string describing which distributional assumption should be used.
#' @param compositional A logical indicating if the dataset should be treated as compositional
#' @param maxIt an integer, the maximum number of iterations
#' @param tol A small scalar, the convergence tolerance
#' @param meanVarFit mean variance model
#' @param newtonRaphson a boolean, should newton-raphson be used
#' @param dispFreq An integer, frequency of dispersion estimation
#' @param ... passed on to the estOff() function
#'
#' @return A list with elements
#' \item{rowOff}{The row offsets}
#' \item{colOff}{The column offsets}
#' \item{converged}{A logical flag, indicating whether the fit converged}
#' \item{iter}{An integer, the number of iterations}
#' @importFrom nleqslv nleqslv
estIndepModel = function(data, distribution, compositional, maxIt, tol, link,
converged = FALSE
iter = 1L
#Disregard NA's
NArows = apply(data, 1, function(x) all(is.na(x)))
#Create starting values
rowOff = link(rowSums(data, na.rm = TRUE))
colOff = link(colSums(data, na.rm = TRUE)/sum(data, na.rm = TRUE))
while(iter < maxIt && !converged){
#Estimate the mean-variance trend if needed
if(((iter-1) %% dispFreq == 0) && distribution == "quasi"){
#Mean variance trend
meanVarTrend = suppressWarnings(estMeanVarTrend(data = data,
colOff, "+")),
meanVarFit = meanVarFit,
}
rowOffOld = rowOff; colOffOld = colOff

rowOff[!NArows] = estOff(colOff = colOff, data = data[!NArows,], distribution = distribution,
col = FALSE, rowOff = rowOff[!NArows], meanVarTrend = meanVarTrend,
newtonRaphson = newtonRaphson, libSizes = invLink(rowOff),...)
if(distribution == "quasi"){
}
colOff = estOff(colOff = colOff, data = data[!NArows,], distribution = distribution,
rowOff = rowOff[!NArows], meanVarTrend = meanVarTrend, col = TRUE,
newtonRaphson = newtonRaphson, libSizes = invLink(rowOff),...)
if(distribution == "quasi"){
}
#Enforce compositionality here, avoids trouble with the jacobian,
#amongst others
if(compositional) {
if(distribution == "quasi"){
colOff = colOff - log(sum(exp(colOff)))
} else if(distribution == "gaussian"){
colOff = colOff/sum(colOff)
}
}
iter = iter + 1L
converged = sqrt(mean((rowOff[!NArows] - rowOffOld[!NArows])^2)) < tol &&
sqrt(mean((colOff - colOffOld)^2)) < tol
}
if(!converged){
warning(immediate. = TRUE,
"Fit of independence model did not converge. Please investigate cause.\n")
}
rowOff[NArows] = 0
return(list(rowOff = rowOff, colOff = colOff, converged = converged, iter = iter))
}


## Try the combi package in your browser

Any scripts or data that you put into this service are public.

combi documentation built on Nov. 8, 2020, 5:34 p.m.