# R/RMI-fn.R In alemermartinez/RMI-GitHub: Robust marginal integration

#' @useDynLib RMI, .registration = TRUE
#' @import stats graphics

#' @title Derivative of Tukey's bi-square loss function
#' @description Derivative of Tukey's bi-square loss function.
#' @param r A real number.
#' @param k Tuning constant.
#' @author Alejandra Martinez, Matias Salibian-Barrera
#' @return A real number.
#' @export

psi.tukey <- function(r, k=4.685){
u <- abs(r/k)
w <- r*((1-u)*(1+u))^2
w[u>1] <- 0
return(w)
}

#' @title Tukey bi-square weight function
#' @description Tukey bi-square weight function.
#' @param r A real number.
#' @param k Tuning constant.
#' @details he weight function used in the re-weighted least squares iterations.
#' @return A real number.
#' @author Alejandra Martinez, Matias Salibian-Barrera
#' @export

#Tukey's weight function "Psi(r)/r"
psi.w <- function(r, k= 4.685){
u <- abs(r/k)
w <- ((1 + u) * (1 - u))^2
w[u > 1] <- 0
return(w)
}

#' @title Derivative of Huber's loss function
#' @description Derivative of Huber's loss function.
#' @param r A real number.
#' @param k Tuning constant.
#' @details The weight function used in the re-weighted least squares iterations.
#' @return A real number.
#' @author Alejandra Martinez, Matias Salibian-Barrera
#' @export

# Huber's Psi
psi.huber <- function(r, k=1.345)
pmin(k, pmax(-k, r))

#' @title Huber weight function
#' @description Huber weight function.
#' @param r A real number.
#' @param k Tuning constant.
#' @details The weight function used in the re-weighted least squares iterations.
#' @return A real number.
#' @author Alejandra Martinez, Matias Salibian-Barrera
#' @export

#Huber's weight function "Psi(r)/r"
psi.huber.w <- function(r, k=1.345)
pmin(1, k/abs(r))

#' @title Euclidean norm of a vector
#' @description Euclidean norm of a vector.
#' @param x A real vector.
#' @return The Euclidean norm of the input vector.
#' @author Alejandra Martinez, Matias Salibian-Barrera
#' @export

#Euclidean norm
my.norm.2 <- function(x) sqrt(sum(x^2))

#' @title Epanechnikov kernel
#' @description Epanechnikov kernel.
#' @usage k.epan(x)
#' kernel2(x)
#' @param x A real number.
#' @author Alejandra Martinez, Matias Salibian-Barrera
#' @return 0 if abs(x) > 1 and 0.75 * (1 - x^2) otherwise.
#' @export

#Epanechnikov kernel
k.epan<-function(x) {
a <- 0.75*(1-x^2)
tmp <- a*(abs(x)<=1)
return(tmp)
}

#' @export

#Order 2 kernel = Epanechnikov kernel
kernel2<-function(t){
nucleo<- (3/4)*(1-t^2)*(abs(t)<=1)
nucleo
}

#- Higher order kernels -#

#' @title Order 4 kernel
#' @description A kernel of order 4.
#' @param x A real number.
#' @author Alejandra Martinez, Matias Salibian-Barrera
#' @return 0 if abs(x) > 1 and ( 15/32 ) * ( 1 - x^2 ) * ( 3 - 7 * x^2 ) otherwise.
#' @details A kernel L is a kernel of order 4 if it integrates 1, the integrals of u^j L(u) are 0 for 1 <= j < 4 (j integer) and the integral of u^4 L(u) is different from 0.
#' @export

#Order 4
kernel4<-function(x) {
a <- (15/32)*(1-x^2)*(3-7*x^2)
tmp <- a*(abs(x)<= 1)
return(tmp)
}

#' @title Order 6 kernel
#' @description A kernel of order 6.
#' @param x A real number.
#' @author Alejandra Martinez, Matias Salibian-Barrera
#' @return 0 if abs(x) > 1 and ( 105/256 ) * ( 1 - x^2 ) * ( 5 - 30 * x^2 + 33 * x^4 ) otherwise.
#' @details A kernel L is a kernel of order 6 if it integrates 1, the integrals of u^j L(u) are 0 for 1 <= j < 6 (j integer) and the integral of u^6 L(u) is different from 0.
#' @export

#Order 6
kernel6<-function(x) {
a <- (105/256)*(1-x^2)*(5-30*x^2+33*x^4)
tmp <- a*(abs(x)<=1)
return(tmp)
}

#' @title Order 8 kernel
#' @description A kernel of order 8.
#' @param x A real number.
#' @author Alejandra Martinez, Matias Salibian-Barrera
#' @return 0 if abs(x) > 1 and ( 315/4096 ) * ( 1 - x^2 ) * ( 35 - 385 * x^2 + 1001 * x^4 - 715 * x^6 ) otherwise.
#' @details A kernel L is a kernel of order 8 if it integrates 1, the integrals of u^j L(u) are 0 for 1 <= j < 8 (j integer) and the integral of u^8 L(u) is different from 0.
#' @export

#Order 8
kernel8<-function(x) {
a <- (315/4096)*(1-x^2)*(35-385*x^2+1001*x^4-715*x^6)
tmp <- a*(abs(x)<=1)
return(tmp)
}

#' @title Order 10 kernel
#' @description A kernel of order 10.
#' @param x A real number.
#' @author Alejandra Martinez, Matias Salibian-Barrera
#' @return 0 if abs(x) > 1 and 0.75 * ( 1 - x^2 ) * ( 315/128 - 105/32 * x^2 + 63/64 * x^4 - 3/32 * x^6 - 1/384 * x^8 ) otherwise.
#' @details A kernel L is a kernel of order 10 if it integrates 1, the integrals of u^j L(u) are 0 for 1 <= j < 10 (j integer) and the integral of u^10 L(u) is different from 0.
#' @export

#Order 10
kernel10<-function(x) {
a <- 0.75*(1-x^2)*(315/128-105/32*x^2+63/64*x^4-3/32*x^6-1/384*x^8)
tmp <- a*(abs(x)<=1)
return(tmp)
}

#' @title Classic marginal integration procedures for additive models
#' @description Standard marginal integration procedures for additive models.
#' @param Xp Matrix of explanatory variables (n by p).
#' @param yp  Vector of responses (missing values are allowed).
#' @param point Matrix of points where predictions will be computed and returned.
#' @param windows Vector or a squared matrix of bandwidths for the smoothing estimation procedure.
#' @param epsilon Convergence criterion.
#' @param prob Probabilities of observing each response (n). Defaults to NULL''.
#' @param type Three different type of estimators can be selected: type '0' (local constant on all the covariates), type '1' (local linear smoother on all the covariates), type 'alpha' (local polynomial smoother only on the direction of interest).
#' @param degree Degree of the local polynomial smoother in the direction of interest when using the estimator of type 'alpha'. Defaults to NULL'' for the case when using estimators of type '0' or '1'.
#' @param orderkernel Order of the kernel used in the nuisance directions when using the estimator of type 'alpha'. Defaults to 2''.
#' @param qderivate If TRUE, it calculates g^(q+1)/(q+1)! for each component only for the type 'alpha' method. Defaults to FALSE''.
#' @param Qmeasure A matrix of points where the integration procedure ocurrs. Defaults to NULL'' for calcuting the integrals over the sample.
#' @details Three types of classical marginal integration procedures for additive models, that is, considering a squared loss function.
#' @return
#' \item{mu}{Estimate for the intercept.}
#' \item{g.matrix}{Matrix of estimated additive components (n by p).}
#' \item{prediction }{Matrix of estimated additive components for the points listed in the argument point.}
#' \item{mul}{A vector of size p showing in each component the estimated intercept that considers only that direction of interest when using the type 'alpha' method.}
#' \item{g.derivative }{Matrix of estimated derivatives of the additive components (only when qderivate is TRUE'') (n by p).}
#' \item{prediction.derivate }{Matrix of estimated derivatives of the additive components for the points listed in the argument point (only when qderivate is TRUE'').}
#' \item{Xp}{Matrix of explanatory variables.}
#' \item{yp}{Vector of responses.}
#'
#' @author Alejandra Martinez, Matias Salibian-Barrera
#' @export

## Classic Marginal Integration
margint.cl <- function(Xp, yp, point=NULL, windows, epsilon=1e-6, prob=NULL,
type='0', degree=NULL, qderivate=FALSE, orderkernel=2,
Qmeasure=NULL) {
# Xp = covariance matrix (n x q).
# yp = response vector (NA's are allowed).
# point = vector of length q or a matrix with q columns where predictions
#  are computed. If missing, predictions are computed for each row of Xp.
# windows = vector of length q or a q times q matrix with kernel windows.
#  The matrix is used for the 'alpha' procedure.
# epsilon = convergence criterion.
# prob = probabilities of observing each response (n).
# type= '0' (local constant), '1' (local linear), 'alpha' (local
#  polynomial using 'degree').
# degree = degree of the local polynomial smoother in the direction of
# interest when using the estimator of type 'alpha'.
# orderkernel = order of the kernel used in the nuisance directions when
#  using the estimator of type 'alpha'.
# qderivate = if TRUE, it calculates g^(q+1)/(q+1)! for each component
#  only for the type 'alpha' method.
# Qmeasure = if NULL, the integration procedure is over the sample.

if(!is.null(dim(Xp))){
if(type=='alpha'){
if(is.null(degree)){
stop("Degree of local polynomial missing")
}else{
if( length(windows)==1 ){
stop("Windows should be a vector o a matrix")
}
}
}else{
if( (is.matrix(windows)) | (length(windows)==1)  ){
stop("Windows should be a vector")
}
}
}else{
if(!is.null(dim(windows))){
stop("Windows should be a number")
}
}

n <- length(yp)
Xp <- as.matrix(Xp)
q <- dim(Xp)[2]
corte <- 10*epsilon
punto <- point

# Remove observations with missing responses
yy <- yp
XX <- Xp
yp <- yp[ tmp<-!is.na(yp) ]
Xp <- Xp[tmp, , drop=FALSE]
n.miss <- length(yp)
if(is.null(prob)){prob <- rep(1,n.miss)
} else {
prob <- prob[tmp]
}
if(dim(t(as.matrix(windows)))[2]!=q){return("Error Dimension of Bandwidths")}

#Initializations
puntoj <- rep(0,q)
pesosi <- rep(1,n.miss)

# If a Qmeasure is provided
if(!is.null(Qmeasure)){

grilla <- rbind(Qmeasure, XX)
nQ <- dim(grilla)[1]

alphal.aux <- matrix(0,nQ,q)
g.matriz.bis <- matrix(0,nQ,q)
g.matriz <- matrix(0,n,q)
aux.b <- rep(0,nQ)
alpha.aux <- rep(0, nQ)
nq <- dim(Qmeasure)[1]
aa <- rep(0,n)

if(qderivate){
g.derivate.bis <- matrix(0,nQ,q)
g.derivate <- matrix(0,n,q)
aa.derivate <- rep(0,nq)
}else{
g.derivate <- NULL
}

for(k in 1:nQ){
for(j in 1:q){
for(i in 1:nq){
puntoj <- Qmeasure[i,]
puntoj[j]<- grilla[k,j]
if(type=='0'){
aa[i] <- .C("kernel_cl_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(windows), as.double(epsilon), as.double(prob),salida=as.double(0) )$salida if(i==k){alpha.aux[k] <- aa[i]} } if(type=='1'){ aa[i] <- .C("kernel_cl_lin_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(windows), as.double(epsilon), as.double(prob),salida=as.double(rep(0, q+1)))$salida[1]
if(i==k){alpha.aux[k] <- aa[i]}
}
if(type=='alpha'){
if(is.null(dim(windows))){
windows <- t(matrix(windows,q,q))
}
AUX <- .C("kernel_cl_alpha_multi", puntoj, as.double(Xp), as.integer(j), as.integer(degree), as.integer(q), as.integer(n.miss), as.integer(orderkernel),
as.double(yp), as.double(windows[j,]), as.double(epsilon), as.double(prob), salida=as.double(rep(0, degree+1)) )$salida aa[i] <- AUX[1] if(qderivate){ aa.derivate[i] <- AUX[degree+1] } if(i==k){alphal.aux[k,j] <- aa[i]} } } g.matriz.bis[k,j]<-mean(aa,na.rm=TRUE) if(qderivate){ g.derivate.bis[k,j]<-mean(aa.derivate,na.rm=TRUE) } } } alphal <- NULL if(type=='alpha'){ alpha.aux <- colMeans(alphal.aux[1:nq,],na.rm=TRUE) alphal <- alpha.aux } alpha <- mean(alpha.aux,na.rm=TRUE) g.matriz <- g.matriz.bis[(nq+1):nQ,] - alpha if(qderivate){ g.derivate <- g.derivate.bis[(nq+1):nQ,] } }else{ #That is, if a Qmeasure is not provided alphal.aux <- matrix(0,n,q) g.matriz <- matrix(0,n,q) if(qderivate){ g.derivate <- matrix(0,n,q) A.derivate <- rep(0,n) }else{ g.derivate <- NULL } aux.b <- rep(0,n) aa <- rep(0,n) alpha.aux <- rep(0, n) for(k in 1:n){ for(j in 1:q){ for(i in 1:n){ puntoj <- XX[i,] puntoj[j]<-XX[k,j] if(type=='0'){ aa[i] <- .C("kernel_cl_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(windows), as.double(epsilon), as.double(prob),salida=as.double(0) )$salida
if(i==k){alpha.aux[k] <- aa[i]}
}
if(type=='1'){
aa[i] <- .C("kernel_cl_lin_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(windows), as.double(epsilon), as.double(prob),salida=as.double(rep(0, q+1)))$salida[1] if(i==k){alpha.aux[k] <- aa[i]} } if(type=='alpha'){ if(is.null(dim(windows))){ windows <- t(matrix(windows,q,q)) } AUX <- .C("kernel_cl_alpha_multi", puntoj, as.double(Xp), as.integer(j), as.integer(degree), as.integer(q), as.integer(n.miss), as.integer(orderkernel), as.double(yp), as.double(windows[j,]), as.double(epsilon), as.double(prob), salida=as.double(rep(0, degree+1)) )$salida
aa[i] <- AUX[1]
if(qderivate){
A.derivate[i] <- AUX[degree+1]
}
if(i==k){alphal.aux[k,j] <- aa[i]}
}
}
g.matriz[k,j]<-mean(aa,na.rm=TRUE)
if(qderivate){
g.derivate[k,j]<-mean(A.derivate,na.rm=TRUE)
}
}
}

alphal <- NULL
if(type=='alpha'){
alpha.aux <- colMeans(alphal.aux,na.rm=TRUE)
alphal <- alpha.aux
}

alpha <- mean(alpha.aux,na.rm=TRUE)
g.matriz <- g.matriz - alpha
}

#Predictions:

prediccion <- NULL

#If a Qmeasure is provided
if(!is.null(Qmeasure)){
aa <- rep(0,nq)
aa.deri <- rep(0,nq)
if(!is.null(punto)) {
if(is.null(dim(punto))) {
prediccion <- mpunto <- t(as.matrix(punto))
if(qderivate){
prediccion.deri <- prediccion
}
} else {
prediccion <- mpunto <- punto
if(qderivate){
prediccion.deri <- prediccion
}
}
np <- dim(mpunto)[1]
for(k in 1:np){
for(j in 1:q){
for(i in 1:nq){
puntoj <- Qmeasure[i,]
puntoj[j] <- mpunto[k,j]

if(type=='0'){
aa[i] <- .C("kernel_cl_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(windows), as.double(epsilon), as.double(prob),salida=as.double(0) )$salida } if(type=='1'){ aa[i] <- .C("kernel_cl_lin_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(windows), as.double(epsilon), as.double(prob),salida=as.double(rep(0, q+1)))$salida[1]
}

if(type=='alpha'){
AUX <- .C("kernel_cl_alpha_multi", puntoj, as.double(Xp), as.integer(j), as.integer(degree), as.integer(q), as.integer(n.miss), as.integer(orderkernel),
as.double(yp), as.double(windows[j,]), as.double(epsilon), as.double(prob), salida=as.double(rep(0, degree+1)) )$salida aa[i] <- AUX[1] if(qderivate){ aa.deri[i] <- AUX[degree+1] } } } prediccion[k,j] <- mean(aa,na.rm=TRUE)-alpha if(qderivate){ prediccion.deri[k,j] <- mean(aa.deri,na.rm=TRUE) } } } } #end if }else{ #If a Qmeasure is not provided aa <- rep(0,n) aa.deri <- rep(0,n) if(!is.null(punto)){ if(is.null(dim(punto))){ prediccion <- mpunto <- t(as.matrix(punto)) if(qderivate){ prediccion.deri <- prediccion } } else { prediccion <- mpunto <- punto if(qderivate){ prediccion.deri <- prediccion } } np <- dim(mpunto)[1] for(k in 1:np){ for(j in 1:q){ for(i in 1:n){ puntoj <- XX[i,] puntoj[j] <- mpunto[k,j] if(type=='0'){ aa[i] <- .C("kernel_cl_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(windows), as.double(epsilon), as.double(prob),salida=as.double(0) )$salida
}
if(type=='1'){
aa[i] <- .C("kernel_cl_lin_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(windows), as.double(epsilon), as.double(prob),salida=as.double(rep(0, q+1)))$salida[1] } if(type=='alpha'){ AUX <- .C("kernel_cl_alpha_multi", puntoj, as.double(Xp), as.integer(j), as.integer(degree), as.integer(q), as.integer(n.miss), as.integer(orderkernel), as.double(yp), as.double(windows[j,]), as.double(epsilon), as.double(prob), salida=as.double(rep(0, degree+1)) )$salida
aa[i] <- AUX[1]
if(qderivate){
aa.deri[i] <- AUX[degree+1]
}
}
}
prediccion[k,j] <- mean(aa,na.rm=TRUE)-alpha
if(qderivate){
prediccion.deri[k,j] <- mean(aa.deri,na.rm=TRUE)
}
}
}
} #End of if
} #End of else

if(!is.null(point)){
if(type=='alpha'){
if(!qderivate){
object <- list(mu=alpha,g.matrix=g.matriz, prediction=prediccion, mul=alphal, Xp=Xp, yp=yp)
class(object) <- c("margint.cl", "margint", "list")
return(object)
} else {
object <- list(mu=alpha,g.matrix=g.matriz, prediction=prediccion, mul=alphal,g.derivate=g.derivate, prediction.derivate=prediccion.deri, Xp=Xp, yp=yp)
class(object) <- c("margint.cl", "margint", "list")
return(object)
}
} else {
object <- list(mu=alpha,g.matrix=g.matriz, prediction=prediccion, Xp=Xp, yp=yp)
class(object) <- c("margint.cl", "margint", "list")
return(object)
}
} else {
if(type=='alpha'){
if(!qderivate){
object <- list(mu=alpha,g.matrix=g.matriz, mul=alphal, Xp=Xp, yp=yp)
class(object) <- c("margint.cl", "margint", "list")
return(object)
} else {
object <- list(mu=alpha,g.matrix=g.matriz, mul=alphal,g.derivate=g.derivate, Xp=Xp, yp=yp)
class(object) <- c("margint.cl", "margint", "list")
return(object)
}
} else {
object <- list(mu=alpha,g.matrix=g.matriz, Xp=Xp, yp=yp)
class(object) <- c("margint.cl", "margint", "list")
return(object)
}
}
}

#' @title Robust marginal integration procedures for additive models
#' @description Robust marginal integration procedures for additive models.
#' @param Xp Matrix of explanatory variables (n by p).
#' @param yp  Vector of responses (missing values are allowed).
#' @param point Matrix of points where predictions will be computed and returned.
#' @param windows Vector or a squared matrix of bandwidths for the smoothing estimation procedure.
#' @param prob Probabilities of observing each response (n). Defaults to NULL''.
#' @param sigma.hat Estimate of the residual standard error. If NULL we use the mad of the residuals obtained with local medians.
#' @param win.sigma Vector of bandwidths for estimating sigma.hat. If NULL it uses the argument windows if it is a vector or its diagonal if it is a matrix.
#' @param epsilon Convergence criterion.
#' @param type Three different type of estimators can be selected: type '0' (local constant on all the covariates), type '1' (local linear smoother on all the covariates), type 'alpha' (local polynomial smoother only on the direction of interest).
#' @param degree Degree of the local polynomial smoother in the direction of interest when using the estimator of type 'alpha'. Defaults to NULL'' for the case when using estimators of type '0' or '1'.
#' @param typePhi One of either Tukey'' or Huber''.
#' @param k.h Tuning constant for a Huber-type loss function. Defaults to 1.345''.
#' @param k.t Tuning constant for a Tukey-type loss function. Defaults to 4.685''.
#' @param max.it Maximum number of iterations for the algorithm.
#' @param qderivate If TRUE, it calculates g^(q+1)/(q+1)! for each component only for the type 'alpha' method. Defaults to FALSE''.
#' @param orderkernel Order of the kernel used in the nuisance directions when using the estimator of type 'alpha'. Defaults to 2''.
#' @param Qmeasure A matrix of points where the integration procedure ocurrs. Defaults to NULL'' for calcuting the integrals over the sample.
#' @details Three types of robust marginal integration procedures for additive models.
#' @return
#' \item{mu }{Estimate for the intercept.}
#' \item{g.matrix }{Matrix of estimated additive components (n by p).}
#' \item{sigma.hat }{Estimate of the residual standard error.}
#' \item{prediction }{Matrix of estimated additive components for the points listed in the argument point.}
#' \item{mul }{A vector of size p showing in each component the estimated intercept that considers only that direction of interest when using the type 'alpha' method.}
#' \item{g.derivative }{Matrix of estimated derivatives of the additive components (only when qderivate is TRUE'') (n by p).}
#' \item{prediction.derivate }{Matrix of estimated derivatives of the additive components for the points listed in the argument point (only when qderivate is TRUE'').}
#' \item{Xp}{Matrix of explanatory variables.}
#' \item{yp}{Vector of responses.}
#'
#' @author Alejandra Martinez, Matias Salibian-Barrera
#' @export

## Robust Marginal Integration
margint.rob <- function(Xp, yp, point=NULL, windows, prob=NULL, sigma.hat=NULL,
win.sigma=NULL, epsilon=1e-06, type='0', degree=NULL, typePhi='Huber',
k.h=1.345, k.t = 4.685, max.it=20, qderivate=FALSE, orderkernel=2,
Qmeasure=NULL){

# Xp = covariance matrix (n x q).
# yp = response vector (NA's are allowed).
# point = vector of length q or a matrix with q columns where predictions
#  are computed. If missing, predictions are computed for each row of Xp.
# windows = vector of length q x q times q matrix with kernel windows.
#  The matrix is used for the 'alpha' procedure.
# epsilon = convergence criterion.
# prob = probabilities of observing each response (n).
# type= '0' (local constant), '1' (local linear), 'alpha' (local
#  polynomial using 'degree').
# degree = degree of the local polynomial smoother in the direction of
# interest when using the estimator of type 'alpha'.
# orderkernel = order of the kernel used in the nuisance directions when
#  using the estimator of type 'alpha'.
# qderivate = if TRUE, it calculates g^(q+1)/(q+1)! for each component (only for the
#  type 'alpha' method.
# Qmeasure = if NULL, the integration procedure is over the sample.
# sigma.hat = estimate of the residual standard error. If missing we use the
#   mad of the residuals obtained with local medians.
# initial = initial intercept estimate.
# k.h = tuning constant for the Huber function.
# k.t = tuning constant for the Tukey function.
# typePhi = 'Huber' or 'Tukey'
# max.it = max number of iterations for the robust estimation step.

if(!is.null(dim(Xp))){
if(type=='alpha'){
if(is.null(degree)){
stop("Degree of local polynomial missing")
}else{
if( length(windows)==1 ){
stop("Windows should be a vector or a matrix")
}
}
}else{
if( (is.matrix(windows)) | (length(windows)==1) ){
stop("Windows should be a vector")
}
}
}else{
if(!is.null(dim(windows))){
stop("Windows should be a number")
}
}

n <- length(yp)
Xp <- as.matrix(Xp)
q <- dim(Xp)[2]
corte <- 10*epsilon
punto <- point

if(is.null(win.sigma)){
if(is.null(dim(windows))){
win.sigma <- windows
}else{
win.sigma <- diag(windows)
}
}

# Remove observations with missing responses
yy<-yp
XX <- Xp
yp <- yp[ tmp<-!is.na(yp) ]
Xp <- Xp[tmp, ,drop=FALSE]
n.miss <- length(yp)
if(is.null(prob)){prob <- rep(1,n.miss)
}else{
prob <- prob[tmp]
}
if(dim(t(as.matrix(windows)))[2]!=q){return("Error Dimension of Bandwidths")}

#Initializations
puntoj <- rep(0,q)
aa <- rep(0,n)
pesosi <- rep(1,n.miss)

# Estimate residual standard error
if(is.null(sigma.hat)){
ab <- rep(0,n.miss)
for(i in 1:n.miss){
xtildebis <- scale(Xp, center=Xp[i,], scale=win.sigma)
a <- matrix(as.numeric(abs(xtildebis) < 1), n.miss, q)
a <- apply(a, 1, prod)
a[ a == 0 ] <- NA
ab[i] <- median( a*yp, na.rm = TRUE)
}
if( sigma.hat < 1e-10 ){sigma.hat <- 1e-10} # sigma.hat <- sd(yp-ab,na.rm=TRUE)
}

#If a Qmeasure is provided
if(!is.null(Qmeasure)){

grilla <- rbind(Qmeasure, XX)
nQ <- dim(grilla)[1]

alphal.aux <- matrix(0,nQ,q)
g.matriz.bis <- matrix(0,nQ,q)
g.matriz <- matrix(0,n,q)
aux.b <- rep(0,nQ)
alpha.aux <- rep(0, nQ)
nq <- dim(Qmeasure)[1]
aa <- rep(0,nq)
if(qderivate){
g.derivate.bis <- matrix(0,nQ,q)
g.derivate <- matrix(0,n,q)
aa.derivate <- rep(0,nq)
}else{
g.derivate <- NULL
}
for(k in 1:nQ){
for(j in 1:q){
for(i in 1:nq){
puntoj <- Qmeasure[i,]
puntoj[j]<-grilla[k,j]

#Inicializo el mu
if(type=='0'){
mu.ini <- .C("ini_mu_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(windows),
as.double(prob), salida=as.double(0) )$salida if(typePhi=='Huber'){ aa[i] <- .C("kernel_huber_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(mu.ini), as.double(windows), as.double(epsilon), as.double(sigma.hat), as.double(prob), as.double(k.h), as.integer(max.it),salida=as.double(0) )$salida
}
if(typePhi=='Tukey'){
aa[i] <- .C("kernel_tukey_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(mu.ini), as.double(windows),
as.double(epsilon), as.double(sigma.hat),
as.double(prob), as.double(k.t), as.integer(max.it),salida=as.double(0) )$salida } if(i==k){alpha.aux[k] <- aa[i]} } if(type=='1'){ mu.ini <- .C("ini_mu_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(windows), as.double(prob), salida=as.double(0) )$salida
beta.ini <- rep(0,q+1)
beta.ini[1] <- mu.ini
if(typePhi=='Huber'){
aa[i] <- .C("kernel_huber_lin_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(beta.ini), as.double(windows),
as.double(epsilon), as.double(sigma.hat),
as.double(prob), as.double(k.h), as.integer(max.it),salida=as.double(rep(0, q+1)) )$salida[1] } if(typePhi=='Tukey'){ aa[i] <- .C("kernel_tukey_lin_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(beta.ini), as.double(windows), as.double(epsilon), as.double(sigma.hat), as.double(prob), as.double(k.t), as.integer(max.it),salida=as.double(rep(0, q+1)) )$salida[1]
}
if(i==k){alpha.aux[k] <- aa[i]}
}
if(type=='alpha'){
mu.ini <- .C("ini_mu_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(windows),
as.double(prob), salida=as.double(0) )$salida beta.ini <- rep(0,degree+1) beta.ini[1] <- mu.ini if(is.null(dim(windows))){ windows <- t(matrix(windows,q,q)) } if(typePhi=='Huber'){ AUX <- .C("kernel_huber_alpha_multi", puntoj, as.double(Xp), as.integer(j), as.integer(degree), as.integer(q), as.integer(n.miss), as.integer(orderkernel), as.double(yp), as.double(beta.ini), as.double(windows[j,]), as.double(epsilon), as.double(sigma.hat), as.double(prob), as.double(k.h), as.integer(max.it),salida=as.double(rep(0, degree+1)) )$salida
aa[i] <- AUX[1]
if(qderivate){
aa.derivate[i] <- AUX[degree+1]
}
}
if(typePhi=='Tukey'){
AUX <- .C("kernel_tukey_alpha_multi", puntoj, as.double(Xp), as.integer(j), as.integer(degree), as.integer(q), as.integer(n.miss), as.integer(orderkernel),
as.double(yp), as.double(beta.ini), as.double(windows[j,]),
as.double(epsilon), as.double(sigma.hat),
as.double(prob), as.double(k.t), as.integer(max.it),salida=as.double(rep(0, degree+1)) )$salida aa[i] <- AUX[1] if(qderivate){ aa.derivate[i] <- AUX[degree+1] } } if(i==k){ alphal.aux[i,j] <- aa[i] } } } g.matriz.bis[k,j]<-mean(aa,na.rm=TRUE) if(qderivate){ g.derivate.bis[k,j]<-mean(aa.derivate,na.rm=TRUE) } } } alphal <- NULL if(type=='alpha'){ alpha.aux <- colMeans(alphal.aux[1:nq,],na.rm=TRUE) alphal <- alpha.aux } alpha <- mean(alpha.aux,na.rm=TRUE) g.matriz <- g.matriz.bis[(nq+1):nQ,] - alpha if(qderivate){ g.derivate <- g.derivate.bis[(nq+1):nQ,] } }else{ #If no Qmeasure is provided alpha.aux <- rep(0, n) g.matriz <- matrix(0,n,q) if(qderivate){ g.derivate <- matrix(0,n,q) aa.derivate <- rep(0,n) }else{ g.derivate <- NULL } alphal.aux <- matrix(0,n,q) for(k in 1:n){ for(j in 1:q){ for(i in 1:n){ puntoj <- XX[i,] puntoj[j]<-XX[k,j] #Inicializo el mu if(type=='0'){ mu.ini <- .C("ini_mu_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(windows), as.double(prob), salida=as.double(0) )$salida
if(typePhi=='Huber'){
aa[i] <- .C("kernel_huber_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(mu.ini), as.double(windows),
as.double(epsilon), as.double(sigma.hat),
as.double(prob), as.double(k.h), as.integer(max.it),salida=as.double(0) )$salida } if(typePhi=='Tukey'){ aa[i] <- .C("kernel_tukey_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(mu.ini), as.double(windows), as.double(epsilon), as.double(sigma.hat), as.double(prob), as.double(k.t), as.integer(max.it),salida=as.double(0) )$salida
}
if(i==k){alpha.aux[k] <- aa[i]}
}
if(type=='1'){
mu.ini <- .C("ini_mu_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(windows),
as.double(prob), salida=as.double(0) )$salida beta.ini <- rep(0,q+1) beta.ini[1] <- mu.ini if(typePhi=='Huber'){ aa[i] <- .C("kernel_huber_lin_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(beta.ini), as.double(windows), as.double(epsilon), as.double(sigma.hat), as.double(prob), as.double(k.h), as.integer(max.it),salida=as.double(rep(0, q+1)) )$salida[1]

}
if(typePhi=='Tukey'){
aa[i] <- .C("kernel_tukey_lin_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(beta.ini), as.double(windows),
as.double(epsilon), as.double(sigma.hat),
as.double(prob), as.double(k.t), as.integer(max.it),salida=as.double(rep(0, q+1)) )$salida[1] } if(i==k){alpha.aux[k] <- aa[i]} } if(type=='alpha'){ mu.ini <- .C("ini_mu_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(windows), as.double(prob), salida=as.double(0) )$salida
beta.ini <- rep(0,degree+1)
beta.ini[1] <- mu.ini
if(is.null(dim(windows))){
windows <- t(matrix(windows,q,q))
}
if(typePhi=='Huber'){
AUX <- .C("kernel_huber_alpha_multi", puntoj, as.double(Xp), as.integer(j), as.integer(degree), as.integer(q), as.integer(n.miss), as.integer(orderkernel),
as.double(yp), as.double(beta.ini), as.double(windows[j,]),
as.double(epsilon), as.double(sigma.hat),
as.double(prob), as.double(k.h), as.integer(max.it),salida=as.double(rep(0, degree+1)) )$salida aa[i] <- AUX[1] if(qderivate){ aa.derivate[i] <- AUX[degree+1] } } if(typePhi=='Tukey'){ AUX <- .C("kernel_tukey_alpha_multi", puntoj, as.double(Xp), as.integer(j), as.integer(degree), as.integer(q), as.integer(n.miss), as.integer(orderkernel), as.double(yp), as.double(beta.ini), as.double(windows[j,]), as.double(epsilon), as.double(sigma.hat), as.double(prob), as.double(k.t), as.integer(max.it),salida=as.double(rep(0, degree+1)) )$salida
aa[i] <- AUX[1]
if(qderivate){
aa.derivate[i] <- AUX[degree+1]
}
}
if(i==k){
alphal.aux[i,j] <- aa[i]
}
}
}
g.matriz[k,j]<-mean(aa,na.rm=TRUE)
if(qderivate){
g.derivate[k,j]<-mean(aa.derivate,na.rm=TRUE)
}
}
}
alphal <- NULL
if(type=='alpha'){
alpha.aux <- colMeans(alphal.aux,na.rm=TRUE)
alphal <- alpha.aux
}
alpha <- mean(alpha.aux,na.rm=TRUE)
g.matriz <- g.matriz - alpha
}#End of else

#Predictions:

prediccion <- NULL

#If a Qmeasure is provided
if(!is.null(Qmeasure)){
aa <- rep(0,nq)
aa.deri <- rep(0,nq)
if(!is.null(punto)) {
if(is.null(dim(punto))){
prediccion <- mpunto <- t(as.matrix(punto))
if(qderivate){
prediccion.deri <- prediccion
}
} else {
prediccion <- mpunto <- punto
if(qderivate){
prediccion.deri <- prediccion
}
}
np <- dim(mpunto)[1]
for(k in 1:np){
for(j in 1:q){
for(i in 1:nq){
puntoj <- grilla[i,]
puntoj[j] <- mpunto[k,j]
if(type=='0'){
#Inicializo el mu
mu.ini <- .C("ini_mu_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(windows),
as.double(prob), salida=as.double(0) )$salida if(typePhi=='Huber'){ aa[i] <- .C("kernel_huber_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(mu.ini), as.double(windows), as.double(epsilon), as.double(sigma.hat), as.double(prob), as.double(k.h), as.integer(max.it),salida=as.double(0) )$salida
}
if(typePhi=='Tukey'){
aa[i] <- .C("kernel_tukey_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(mu.ini), as.double(windows),
as.double(epsilon), as.double(sigma.hat),
as.double(prob), as.double(k.t), as.integer(max.it),salida=as.double(0) )$salida } } if(type=='1'){ mu.ini <- .C("ini_mu_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(windows), as.double(prob), salida=as.double(0) )$salida
beta.ini <- rep(0,q+1)
beta.ini[1] <- mu.ini
if(typePhi=='Huber'){
aa[i] <- .C("kernel_huber_lin_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(beta.ini), as.double(windows),
as.double(epsilon), as.double(sigma.hat),
as.double(prob), as.double(k.h), as.integer(max.it),salida=as.double(rep(0, q+1)) )$salida[1] } if(typePhi=='Tukey'){ aa[i] <- .C("kernel_tukey_lin_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(beta.ini), as.double(windows), as.double(epsilon), as.double(sigma.hat), as.double(prob), as.double(k.t), as.integer(max.it),salida=as.double(rep(0, q+1)) )$salida[1]
}
}
if(type=='alpha'){
mu.ini <- .C("ini_mu_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(windows),
as.double(prob), salida=as.double(0) )$salida beta.ini <- rep(0,degree+1) beta.ini[1] <- mu.ini if(typePhi=='Huber'){ AUX <- .C("kernel_huber_alpha_multi", puntoj, as.double(Xp), as.integer(j), as.integer(degree), as.integer(q), as.integer(n.miss), as.integer(orderkernel), as.double(yp), as.double(beta.ini), as.double(windows[j,]), as.double(epsilon), as.double(sigma.hat), as.double(prob), as.double(k.h), as.integer(max.it),salida=as.double(rep(0, degree+1)) )$salida
aa[i] <- AUX[1]
if(qderivate){
aa.deri[i] <- AUX[degree+1]
}
}
if(typePhi=='Tukey'){
AUX <- .C("kernel_tukey_alpha_multi", puntoj, as.double(Xp), as.integer(j), as.integer(degree), as.integer(q), as.integer(n.miss), as.integer(orderkernel),
as.double(yp), as.double(beta.ini), as.double(windows[j,]),
as.double(epsilon), as.double(sigma.hat),
as.double(prob), as.double(k.t), as.integer(max.it),salida=as.double(rep(0, degree+1)) )$salida aa[i] <- AUX[1] if(qderivate){ aa.deri[i] <- AUX[degree+1] } } } } prediccion[k,j] <- mean(aa,na.rm=TRUE)-alpha if(qderivate){ prediccion.deri[k,j] <- mean(aa.deri,na.rm=TRUE) } } } } }else{ #If no Qmeasure is provided aa <- rep(0,n) aa.deri <- rep(0,n) if(!is.null(punto)) { if(is.null(dim(punto))){ prediccion <- mpunto <- t(as.matrix(punto)) if(qderivate){ prediccion.deri <- prediccion } } else { prediccion <- mpunto <- punto if(qderivate){ prediccion.deri <- prediccion } } np <- dim(mpunto)[1] for(k in 1:np){ for(j in 1:q){ for(i in 1:n){ puntoj <- XX[i,] puntoj[j] <- mpunto[k,j] if(type=='0'){ #Inicializo el mu mu.ini <- .C("ini_mu_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(windows), as.double(prob), salida=as.double(0) )$salida

if(typePhi=='Huber'){
aa[i] <- .C("kernel_huber_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(mu.ini), as.double(windows),
as.double(epsilon), as.double(sigma.hat),
as.double(prob), as.double(k.h), as.integer(max.it),salida=as.double(0) )$salida } if(typePhi=='Tukey'){ aa[i] <- .C("kernel_tukey_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(mu.ini), as.double(windows), as.double(epsilon), as.double(sigma.hat), as.double(prob), as.double(k.t), as.integer(max.it),salida=as.double(0) )$salida
}
}
if(type=='1'){
mu.ini <- .C("ini_mu_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(windows),
as.double(prob), salida=as.double(0) )$salida beta.ini <- rep(0,q+1) beta.ini[1] <- mu.ini if(typePhi=='Huber'){ aa[i] <- .C("kernel_huber_lin_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(beta.ini), as.double(windows), as.double(epsilon), as.double(sigma.hat), as.double(prob), as.double(k.h), as.integer(max.it),salida=as.double(rep(0, q+1)) )$salida[1]
}
if(typePhi=='Tukey'){
aa[i] <- .C("kernel_tukey_lin_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss),
as.double(yp), as.double(beta.ini), as.double(windows),
as.double(epsilon), as.double(sigma.hat),
as.double(prob), as.double(k.t), as.integer(max.it),salida=as.double(rep(0, q+1)) )$salida[1] } } if(type=='alpha'){ mu.ini <- .C("ini_mu_pos_multi", puntoj, as.double(Xp), as.integer(q), as.integer(n.miss), as.double(yp), as.double(windows), as.double(prob), salida=as.double(0) )$salida
beta.ini <- rep(0,degree+1)
beta.ini[1] <- mu.ini
if(typePhi=='Huber'){
AUX <- .C("kernel_huber_alpha_multi", puntoj, as.double(Xp), as.integer(j), as.integer(degree), as.integer(q), as.integer(n.miss), as.integer(orderkernel),
as.double(yp), as.double(beta.ini), as.double(windows[j,]),
as.double(epsilon), as.double(sigma.hat),
as.double(prob), as.double(k.h), as.integer(max.it),salida=as.double(rep(0, degree+1)) )$salida aa[i] <- AUX[1] if(qderivate){ aa.deri[i] <- AUX[degree+1] } } if(typePhi=='Tukey'){ AUX <- .C("kernel_tukey_alpha_multi", puntoj, as.double(Xp), as.integer(j), as.integer(degree), as.integer(q), as.integer(n.miss), as.integer(orderkernel), as.double(yp), as.double(beta.ini), as.double(windows[j,]), as.double(epsilon), as.double(sigma.hat), as.double(prob), as.double(k.t), as.integer(max.it),salida=as.double(rep(0, degree+1)) )$salida
aa[i] <- AUX[1]
if(qderivate){
aa.deri[i] <- AUX[degree+1]
}
}
}
}
prediccion[k,j] <- mean(aa,na.rm=TRUE)-alpha
if(qderivate){
prediccion.deri[k,j] <- mean(aa.deri,na.rm=TRUE)
}
}
}
}
}

if(!is.null(point)){
if(type=='alpha'){
if(!qderivate){
object <- list(mu=alpha, g.matrix=g.matriz, sigma.hat=sigma.hat, prediction=prediccion, mul=alphal, Xp=Xp, yp=yp)
class(object) <- c("margint.rob", "margint", "list")
return(object)
} else {
object <- list(mu=alpha, g.matrix=g.matriz, sigma.hat=sigma.hat, prediction=prediccion, mul=alphal, g.derivate=g.derivate, prediction.derivate=prediccion.deri, Xp=Xp, yp=yp)
class(object) <- c("margint.rob", "margint", "list")
return(object)
}
} else {
object <- list(mu=alpha, g.matrix=g.matriz, sigma.hat=sigma.hat, prediction=prediccion, Xp=Xp, yp=yp)
class(object) <- c("margint.rob", "margint", "list")
return(object)
}
} else {
if(type=='alpha'){
if(!qderivate){
object <- list(mu=alpha, g.matrix=g.matriz, sigma.hat=sigma.hat, mul=alphal, Xp=Xp, yp=yp)
class(object) <- c("margint.rob", "margint", "list")
return(object)
} else {
object <- list(mu=alpha, g.matrix=g.matriz, sigma.hat=sigma.hat, mul=alphal,g.derivate=g.derivate, Xp=Xp, yp=yp)
class(object) <- c("margint.rob", "margint", "list")
return(object)
}
} else {
object <- list(mu=alpha, g.matrix=g.matriz, sigma.hat=sigma.hat, Xp=Xp, yp=yp)
class(object) <- c("margint.rob", "margint", "list")
return(object)
}
}
}

#S3 Methods

#' Residuals of a marginal integration fit in an additive model
#'
#' This function returns the residuals of the fitted additive model using one of the three
#' classical or robust marginal integration estimators, as computed with \code{\link{margint.cl}} or
#'
#' @param object an object of class \code{margint}, a result of a call to \code{\link{margint.cl}} or \code{\link{margint.rob}}.
#' @param ... additional other arguments.
#'
#' @return A vector of residuals.
#'
#' @author Alejandra Mercedes Martinez \email{alemartinez@unlu.edu.ar}
#'
#' @export

residuals.margint <- function(object, ...){
return( object$yp - rowSums(object$g.matrix) -object$mu ) } #' Fitted values for objects of class \code{margint} #' #' This function returns the fitted values given the covariates of the original sample under an additive model using a classical or robust marginal integration procedure estimator computed with \code{\link{margint.cl}} or \code{\link{margint.rob}}. #' #' @param object an object of class \code{margint}, a result of a call to \code{\link{margint.cl}} or \code{\link{margint.rob}}. #' @param ... additional other arguments. #' #' @return A vector of fitted values. #' #' @author Alejandra Mercedes Martinez \email{alemartinez@unlu.edu.ar} #' #' @export predict.margint <- function(object, ...){ return( rowSums(object$g.matrix) + object$mu ) } #' Diagnostic plots for objects of class \code{margint} #' #' Plot method for class \code{margint}. #' #' @param object an object of class \code{margint}, a result of a call to \code{\link{margint.cl}} or \code{\link{margint.rob}}. #' @param which vector of indices of explanatory variables for which partial residuals plots will be generated. Defaults to all available explanatory variables. #' @param ask logical value. If \code{TRUE}, the graphical device will prompt before going to the next page/screen of output. #' @param ... additional other arguments. #' #' @author Alejandra Mercedes Martinez \email{alemartinez@unlu.edu.ar} #' #' @export plot.margint <- function(object, which=1:np, ask=FALSE,...){ Xp <- object$Xp
np <- dim(Xp)[2]
on.exit(par(opar))
these <- rep(FALSE, np)
these[ which ] <- TRUE
for(i in 1:np) {
if(these[i]) {
ord <- order(Xp[,i])
if (is.null(colnames(Xp)) ){
x_name <- bquote(paste('x')[.(i)])
#paste("x",i,sep="")
} else {
x_name <- colnames(Xp)[i]
}
y_name <- bquote(paste(hat('g')[.(i)]))
res <- object$yp - rowSums(object$g.matrix[,-i, drop=FALSE])-object$mu lim_cl <- c(min(res), max(res)) plot(Xp[ord,i],object$g.matrix[ord,i],type="l",lwd=3,main="",xlab=x_name,ylab=y_name, ylim=lim_cl)
points(Xp[,i], res, pch=20,col='gray45')
}
}
}

#' Summary for additive models fits using a marginal integration procedure
#'
#' Summary method for class \code{margint}.
#'
#' This function returns the estimation of the intercept and also the five-number summary and the mean of the residuals for both classical and robust estimators. For the robust estimator it also returns the estimate of the residual standard error.
#'
#' @param object an object of class \code{margint}, a result of a call to \code{\link{margint.cl}} or \code{\link{margint.rob}}.
#' @param ... additional other arguments.
#'
#' @author Alejandra Mercedes Martinez \email{alemartinez@unlu.edu.ar}
#'
#' @export
#' @aliases summary.margint summary.margint.cl summary.margint.rob

summary.margint <- function(object,...){
NextMethod()
}

summary.margint.cl <- function(object,...){
message("Estimate of the intercept: ", round(object$mu,5)) res <- residuals(object) message("Residuals:") summary(res) } summary.margint.rob <- function(object,...){ message("Estimate of the intercept: ", round(object$mu,5))
message("Estimate of the residual standard error: ", round(object\$sigma,5))
res <- residuals(object)
message("Residuals:")
summary(res)
}

alemermartinez/RMI-GitHub documentation built on May 9, 2019, 2:21 a.m.