mprocuste.R 8.86 KB
Newer Older
1 2 3
#' @include internals.R
#' @include model.procmod.R
#'
4 5 6 7
#' @title ProcMod
#' @description blabla
#' @author Christelle Gonindard-Melodelima
#'
Eric Coissac committed
8
NULL
9

10

11
pm.fit = function(covmat,y,xs,
12
                  singular.ok = singular.ok,
13
                  tol) {
14

15
  xy.cov = covmat[xs,y,drop=FALSE]
16 17
  xx.cov = covmat[xs,xs,drop=FALSE]

18
  #xx.qr  = qr(xx.cov,tol=tol)
19

20 21
  #if (!singular.ok && xx.qr$rank < ncol(xx.cov))
  #  stop("singular fit encountered")
22

23 24 25 26
  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)
27 28


29 30 31
  return(list(coefficients=coef
              #rank=xx.qr$rank + 1
              #qr=xx.qr,
32 33
             )
        )
34 35 36 37
}

.ctrace <- function(MAT) sum(MAT^2)

38 39 40 41 42 43 44 45 46
#' Tests that the object is as `pm` instance.
#' 
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
is.pm = function(obj) {
  inherits(obj, "pm")
}

47 48
#' Performs a procruste model on a set of coordinate matrices
#'
49 50 51
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
52 53
pm = function (formula,data, subset, weights, na.action, method = "qr",
               singular.ok = TRUE, tol = 1e-07,
54
               model = TRUE, x = TRUE, y = TRUE, A = TRUE) {
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103

  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
104
  }
105
  else {
106
    
107 108 109 110 111 112 113 114
    vars <- vars.procmod(mt, mf)
    nvars=ncol(vars)
    irep = attr(vars,"response")

    data.cov = mvar(vars)
    variances= diag(data.cov)
    std.dev = sqrt(variances)

115 116 117
    vars.norm = as.procmod.frame(mapply(function(x) scale(x,scale = FALSE),
                                        vars,
                                        SIMPLIFY = FALSE))
118 119
    
    
120 121 122 123 124 125

    if (is.null(w)) {
      subset.w=rep(TRUE,nvars)
      }
    else {
      sw = sqrt(w)
126 127 128

      vars.norm = as.procmod.frame(lapply(vars.norm, function(v) v * sw))
      subset.w  = sw > 0
129

130
      vars.norm = vars.norm[subset.w,]
131 132
    }

133
    data.cov = mvar(vars.norm)
134

135 136 137 138
    z = pm.fit(data.cov,
               y = irep, xs = seq_len(nvars)!=irep,
               singular.ok = singular.ok,
               tol = tol)
139

140 141
    z$cov = data.cov

142 143
    y.sd = std.dev[irep]
    xs.sd= std.dev[seq_len(nvars)!=irep]
144

145
    z$varpart = (z$coefficients * xs.sd / y.sd) ^ 2
146

147

148
    Y = vars.norm[[irep]]
149
    Xs= vars.norm[seq_len(nvars)!=irep,drop=FALSE]
150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167

    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

168 169
    scales = z$coefficients

170 171 172 173 174 175 176 177 178 179
    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)
180

181 182
    init=matrix(rep(0,nrow(yhat[[1]])*ncol(yhat[[1]])),
                nrow=nrow(yhat[[1]]))
183

184 185 186
    yhat = Reduce(function(a,b) a+b,
                  yhat,
                  init)
187

188
    yhat = sweep(yhat,2,b,"+")
189

190 191 192
    z$residuals = Y - yhat
    z$fitted.values = yhat
    z$weights = w
193
    z$df.residual = (if (!is.null(w))
194
                      sum(w != 0)
195
                     else nrow(Y)) - nvars + 1
196

197 198
    if (ret.x)
      z$x <- vars.norm[-irep]
199

200
##    if (!A)
201 202
      z$A <- A_xys
  }
203 204


205 206 207 208 209 210
  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
211

212
  z$SCT = .ctrace(scale(Y,scale=FALSE))
213

214 215
  if (model)
    z$model <- mf
216

217 218
  if (ret.y)
    z$y <- Y
219

220
  return(z)
221 222 223

}

224 225
#' Default display function for `pm` objects.
#'
226 227 228 229 230 231 232 233 234 235
#' @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)
}

236

237

238 239 240 241 242 243 244 245 246 247 248 249 250 251
summary.pm = function(object,
                      correlation = FALSE,
                      symbolic.cor = FALSE,
                      ...) {

  results = list(
    call = object$call,
    terms= object$terms

  )

  return(results)
}

252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276
#' 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)
}

Eric Coissac committed
277

278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312
#' 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
Eric Coissac committed
313
}
314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366


#' 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)
}