#' @include internals.R #' @include model.procmod.R #' #' @title ProcMod #' @description blabla #' @author Christelle Gonindard-Melodelima #' NULL pm.fit = function(covmat,y,xs, singular.ok = singular.ok, tol) { xy.cov = covmat[xs,y,drop=FALSE] xx.cov = covmat[xs,xs,drop=FALSE] #xx.qr = qr(xx.cov,tol=tol) #if (!singular.ok && xx.qr$rank < ncol(xx.cov)) # stop("singular fit encountered") coef = as.numeric(solve(xx.cov) %*% xy.cov) names(coef)=rownames(xx.cov) #coef = qr.coef(xx.qr,xy.cov) #effects= qr.qty(xx.qr,xy.cov) return(list(coefficients=coef #rank=xx.qr$rank + 1 #qr=xx.qr, ) ) } .ctrace <- function(MAT) sum(MAT^2) #' Tests that the object is as `pm` instance. #' #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export is.pm = function(obj) { inherits(obj, "pm") } #' Performs a procruste model on a set of coordinate matrices #' #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export pm = function (formula,data, subset, weights, na.action, method = "qr", singular.ok = TRUE, tol = 1e-07, model = TRUE, x = TRUE, y = TRUE, A = TRUE) { ret.x <- x ret.y <- y cl <- match.call() mf <- match.call(expand.dots = FALSE) m <- match(c("formula", "data", "subset", "weights", "na.action"), table = names(mf), nomatch = 0L) mf <- mf[c(1L, m)] mf$drop.unused.levels <- TRUE mf[[1L]] <- quote(ProcMod::model.procmod.default) mf <- eval(mf, parent.frame()) if (method == "model.frame") return(mf) else if (method != "qr") warning(gettextf("method = '%s' is not supported. Using 'qr'", method), domain = NA) mt <- attr(mf, "terms") Y <- model.response(mf, "numeric") w <- as.vector(model.weights(mf)) if (!is.null(w) && !is.numeric(w)) stop("'weights' must be a numeric vector") if (is.empty.model(mt)) { x <- NULL z <- list(coefficients = NULL, anova = NULL, SCT = SCT, residuals = Y, fitted.values = 0 * Y, weights = w, rank = 0L, df.residual = if (!is.null(w)) sum(w != 0) else nrow(Y)) if (ret.x) z$x <- NULL if (!A) z$A <- NULL } else { vars <- vars.procmod(mt, mf) nvars=ncol(vars) irep = attr(vars,"response") data.cov = mvar(vars) variances= diag(data.cov) std.dev = sqrt(variances) vars.norm = as.procmod.frame(mapply(function(x) scale(x,scale = FALSE), vars, SIMPLIFY = FALSE)) if (is.null(w)) { subset.w=rep(TRUE,nvars) } else { sw = sqrt(w) vars.norm = as.procmod.frame(lapply(vars.norm, function(v) v * sw)) subset.w = sw > 0 vars.norm = vars.norm[subset.w,] } data.cov = mvar(vars.norm) z = pm.fit(data.cov, y = irep, xs = seq_len(nvars)!=irep, singular.ok = singular.ok, tol = tol) z$cov = data.cov y.sd = std.dev[irep] xs.sd= std.dev[seq_len(nvars)!=irep] z$varpart = (z$coefficients * xs.sd / y.sd) ^ 2 Y = vars.norm[[irep]] Xs= vars.norm[seq_len(nvars)!=irep,drop=FALSE] Ymean <- colMeans(Y) Xmeans<- lapply(Xs,colMeans) # Shall we compute rotation on weigthed matrix or not ... ? XYs <- lapply(Xs,function(x) crossprod(Y, x)) sol_yxs <- lapply(XYs,svd) A_xys <- lapply(sol_yxs, function(x) x$v %*% t(x$u)) Xrots <- mapply(function(x,a) x %*% a, Xs,A_xys, SIMPLIFY = FALSE) # Intercept estimation scales = z$coefficients b = Ymean - Reduce(function(a,b) a+b, mapply(function(d,coef,rot) coef * (d %*% rot), Xmeans, scales, A_xys, SIMPLIFY = FALSE)) yhat = mapply(function(coef,xrot) coef * xrot, scales, Xrots, SIMPLIFY = FALSE) init=matrix(rep(0,nrow(yhat[[1]])*ncol(yhat[[1]])), nrow=nrow(yhat[[1]])) yhat = Reduce(function(a,b) a+b, yhat, init) yhat = sweep(yhat,2,b,"+") z$residuals = Y - yhat z$fitted.values = yhat z$weights = w z$df.residual = (if (!is.null(w)) sum(w != 0) else nrow(Y)) - nvars + 1 if (ret.x) z$x <- vars.norm[-irep] ## if (!A) z$A <- A_xys } class(z) <- "pm" z$na.action <- attr(mf, "na.action") # z$contrasts <- attr(x, "contrasts") # z$xlevels <- .getXlevels(mt, mf) z$call <- cl z$terms <- mt z$SCT = .ctrace(scale(Y,scale=FALSE)) if (model) z$model <- mf if (ret.y) z$y <- Y return(z) } #' Default display function for `pm` objects. #' #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export print.pm = function(m,...) { cat(sprintf("\nCall:\n%s\n\nCoefficients:\n", deparse(m$call))) print(m$coefficients) invisible(m) } summary.pm = function(object, correlation = FALSE, symbolic.cor = FALSE, ...) { results = list( call = object$call, terms= object$terms ) return(results) } #' Compute the Akaike Information criterion (AIC) for a procruste model. #' #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export extractAIC.pm = function(fit, scale, k = 2,...) { n = nrow(fit$residuals) RSS = sum(fit$residuals^2) p = length(fit$coefficients) return(c(n,(n * log(RSS/n)+ p * k))) } #' Compute the Akaike Information criterion (AIC) for a procruste model. #' #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export AIC.pm = function(object, ..., k = length(object$coefficients)) { RSS = sum(object$residuals^2) n = nrow(object$fitted.values) return(n * log(RSS/n)+ 2 * k) } #' Extract Model Residuals. #' #' @description residuals is a generic function which extracts model #' residuals from objects returned by modeling functions. #' The abbreviated form resid is an alias for residuals. #' It is intended to encourage users to access object #' components through an accessor function rather than by #' directly referencing an object slot. #' #' @param obj R object of class `pm` #' @param type the type of residuals which should be returned. #' Can be abbreviated. #' #' @return Numeric matrix of n rows, where n is the number #' of observations. #' #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export residuals.pm = function (object, type = c("working", "response", "deviance", "pearson", "partial"), ...) { type <- match.arg(type) r <- object$residuals res <- switch(type, working = , response = r, deviance = , pearson = if (is.null(object$weights)) r else sweep(r,MARGIN = 1, STATS = sqrt(object$weights), FUN = "*"), partial = r) res <- naresid(object$na.action, res) if (type == "partial") { stop("'partial' residuals not yet implemented for procrustean models ") res <- res + predict(object, type = "terms") } res } #' Compute Weighted Residuals. #' #' @description Computed weighted residuals from a procrustean model fit. #' #' @param obj R object of class `pm` #' @param drop0 logical. If `TRUE``, drop all cases with `weights` == 0. #' #' @return Numeric matrix of n' rows, where n' is the number #' of of non-0 weights (drop0 = TRUE) or #' the number of observations, otherwise. #' #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export weighted.residuals = function (obj, drop0 = TRUE) { w <- weights(obj) r <- residuals(obj, type = "deviance") if (drop0 && !is.null(w)) { if (is.matrix(r)) r[w != 0, , drop = FALSE] else r[w != 0] } else r } #' Model Deviance for the procruste model. #' #' @description Returns the deviance of a fitted model object. #' #' @return The value of the deviance extracted from the object object. #' #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export deviance.pm = function (object, ...) sum(weighted.residuals(object)^2, na.rm = TRUE) #' Compute the Bayesian Information criterion (BIC) for a procruste model. #' #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export BIC.pm = function(object, ..., k = length(object$coefficients)) { RSS = sum(object$residuals^2) n = nrow(object$fitted.values) return(n * log(RSS/n)+ log(n) * k) }