Nothing
#' Derivative of Tukey's bi-square loss function.
#'
#' This function evaluates the first derivative of Tukey's bi-square loss function.
#'
#' This function evaluates the first derivative of Tukey's bi-square loss function.
#'
#' @param r A vector of real numbers
#' @param k A positive tuning constant.
#'
#' @return A vector of the same length as \code{r}.
#'
#' @author Matias Salibian-Barrera, \email{matias@stat.ubc.ca}, Alejandra Martinez
#'
#' @examples
#' x <- seq(-2, 2, length=10)
#' psi.tukey(r=x, k = 1.5)
#'
#' @export
#' @import stats graphics
#' @useDynLib rmargint, .registration = TRUE
psi.tukey <- function(r, k=4.685){
u <- abs(r/k)
w <- r*((1-u)*(1+u))^2
w[u>1] <- 0
return(w)
}
# 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)
}
#Tukey's loss function
rho.tukey <- function(r, k=4.685){
if(abs(r)<=k){
return( 1-(1-(r/k)^2)^3 )
}else{
return(1)
}
}
#' Derivative of Huber's loss function.
#'
#' This function evaluates the first derivative of Huber's loss function.
#'
#' This function evaluates the first derivative of Huber's loss function.
#'
#' @param r A vector of real numbers.
#' @param k A positive tuning constant.
#'
#' @return A vector of the same length as \code{r}.
#'
#' @author Matias Salibian-Barrera, \email{matias@stat.ubc.ca}, Alejandra Martinez
#'
#' @examples
#' x <- seq(-2, 2, length=10)
#' psi.huber(r=x, k = 1.5)
#'
#' @export
psi.huber <- function(r, k=1.345)
pmin(k, pmax(-k, r))
#Huber's weight function "Psi(r)/r"
psi.huber.w <- function(r, k=1.345)
pmin(1, k/abs(r))
#Huber's loss function
rho.huber <- function(r, k=1.345){
if(abs(r)<=k){
return(r^2)
}else{
return(2*k*abs(r)-k^2)
}
}
#' Euclidean norm of a vector
#'
#' This function calculates the Euclidean norm of a vector.
#'
#' @param x A real vector.
#'
#' @return The Euclidean norm of the input vector.
#'
#' @author Matias Salibian-Barrera, \email{matias@stat.ubc.ca}, Alejandra Martinez
#'
#' @examples
#' x <- seq(-2, 2, length=10)
#' my.norm.2(x)
#'
#' @export
my.norm.2 <- function(x) sqrt(sum(x^2))
#' Epanechnikov kernel
#'
#' This function evaluates an Epanechnikov kernel
#'
#' This function evaluates an Epanechnikov kernel.
#'
#' @param x a vector of real numbers
#'
#' @return A vector of the same length as \code{x} where each entry is
#' \code{0.75 * (1 - x^2)} if \code{x < 1} and 0 otherwise.
#'
#' @author Matias Salibian-Barrera, \email{matias@stat.ubc.ca}, Alejandra Martinez
#'
#' @examples
#' x <- seq(-2, 2, length=10)
#' k.epan(x)
#'
#' @export
k.epan<-function(x) {
a <- 0.75*(1-x^2)
tmp <- a*(abs(x)<=1)
return(tmp)
}
#Order 2 kernel = Epanechnikov kernel
kernel2<-function(t){
nucleo<- (3/4)*(1-t^2)*(abs(t)<=1)
nucleo
}
#- Higher order kernels -#
#' Order 4 kernel
#'
#' This function evaluates a kernel of order 4.
#'
#' This function evaluates a kernel of order 4. 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.
#'
#' @param x A vector of real numbers.
#'
#' @return A vector of the same length as \code{x} where each entry is \code{( 15/32 ) * ( 1 - x^2 ) * ( 3 - 7 * x^2 )} if \code{abs(x) < 1} and 0 otherwise.
#'
#' @author Alejandra Martinez, \email{ale_m_martinez@hotmail.com}, Matias Salibian-Barrera
#'
#' @examples
#' x <- seq(-2,2,length=10)
#' kernel4(x)
#'
#' @export
kernel4<-function(x) {
a <- (15/32)*(1-x^2)*(3-7*x^2)
tmp <- a*(abs(x)<= 1)
return(tmp)
}
#' Order 6 kernel
#'
#' This function evaluates a kernel of order 6.
#'
#' This function evaluates a kernel of order 6. 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.
#'
#'
#' @param x A vector of real numbers.
#'
#' @return A vector of the same length as \code{x} where each entry is \code{( 105/256 ) * ( 1 - x^2 ) * ( 5 - 30 * x^2 + 33 * x^4 )} if \code{abs(x) < 1} and 0 otherwise.
#'
#' @author Alejandra Martinez, \email{ale_m_martinez@hotmail.com}, Matias Salibian-Barrera
#'
#' @examples
#' x <- seq(-2,2,length=10)
#' kernel6(x)
#'
#' @export
kernel6<-function(x) {
a <- (105/256)*(1-x^2)*(5-30*x^2+33*x^4)
tmp <- a*(abs(x)<=1)
return(tmp)
}
#' Order 8 kernel
#'
#' This function evaluates a kernel of order 8.
#'
#' This function evaluates a kernel of order 8. 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.
#'
#' @param x A vector of real numbers.
#'
#' @return A vector of the same length as \code{x} where each entry is \code{( 315/4096 ) * ( 1 - x^2 ) * ( 35 - 385 * x^2 + 1001 * x^4 - 715 * x^6 )} and 0 otherwise.
#'
#' @author Alejandra Martinez, \email{ale_m_martinez@hotmail.com}, Matias Salibian-Barrera
#'
#' @examples
#' x <- seq(-2,2,length=10)
#' kernel8(x)
#' @export
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)
}
#' Order 10 kernel
#'
#' This function evaluates a kernel of order 10. A kernel of order 10.
#'
#' This function evaluates a kernel of order 10. 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.
#'
#' @param x A vector of real numbers.
#'
#' @return A vector of the same length as \code{x} where each entry is \code{0.75 * ( 1 - x^2 ) * ( 315/128 - 105/32 * x^2 + 63/64 * x^4 - 3/32 * x^6 - 1/384 * x^8 )} and 0 otherwise.
#'
#' @author Alejandra Martinez, \email{ale_m_martinez@hotmail.com}, Matias Salibian-Barrera
#'
#' @examples
#' x <- seq(-2,2,length=10)
#' kernel10(x)
#'
#' @export
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)
}
#' Classic marginal integration procedures for additive models
#'
#' This function computes the standard marginal integration procedures for additive models.
#'
#' This function computes three types of classical marginal integration procedures for additive models, that is, considering a squared loss function.
#'
#' @param formula an object of class \code{formula} (or one that can be coerced to that class): a symbolic description of the model to be fitted.
#' @param data an optional data frame, list or environment (or object coercible by \link{as.data.frame} to a data frame) containing the variables in the model. If not found in \code{data}, the variables are taken from \code{environment(formula)}, typically the environment from which the function was called.
#' @param subset an optional vector specifying a subset of observations to be used in the fitting process.
#' @param point a matrix of points where predictions will be computed and returned.
#' @param windows a vector or a squared matrix of bandwidths for the smoothing estimation procedure.
#' @param epsilon convergence criterion.
#' @param prob a vector of probabilities of observing each response (n). Defaults to \code{NULL}.
#' @param type three different type of estimators can be selected: type \code{'0'} (local constant on all the covariates), type \code{'1'} (local linear smoother on all the covariates), type \code{'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 \code{'alpha'}. Defaults to \code{NULL} for the case when using estimators of type \code{'0'} or \code{'1'}.
#' @param orderkernel order of the kernel used in the nuisance directions when using the estimator of type \code{'alpha'}. Defaults to \code{2}.
#' @param qderivate if TRUE, it calculates \code{g^(q+1)/(q+1)!} for each component only for the type \code{'alpha'} method. Defaults to \code{FALSE}.
#' @param Qmeasure a matrix of points where the integration procedure ocurrs. Defaults to \code{NULL} for calcuting the integrals over the sample.
#'
#' @return A list with the following components:
#' \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 \code{'alpha'} method.}
#' \item{g.derivative }{Matrix of estimated derivatives of the additive components (only when qderivate is \code{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 \code{TRUE}).}
#' \item{Xp}{Matrix of explanatory variables.}
#' \item{yp}{Vector of responses.}
#' \item{formula}{Model formula}
#'
#' @references
#' Chen R., Hardle W., Linton O.B. and Severance-Lossin E. (1996). Nonparametric estimation of additive separable regression models. Physica-Verlag HD, Switzerland.
#' Linton O. and Nielsen J. (1995). A kernel method of estimating structured nonparametric regression based on marginal integration. Biometrika, 82(1), 93-101.
#' Severance-Lossin E. and Sperlich S. (1999). Estimation of derivatives for additive separable models. Statistics, 33(3), 241-265.
#' Tjostheim D. and Auestad B. (1994). Nonparametric identification of nonlinear time series: Selecting significant lags. Journal of the American Statistical Association, 89(428), 1410-1430.
#'
#' @author Alejandra Martinez, \email{ale_m_martinez@hotmail.com}, Matias Salibian-Barrera
#'
#' @examples
#' function.g1 <- function(x1) 24*(x1-1/2)^2-2
#' function.g2 <- function(x2) 2*pi*sin(pi*x2)-4
#' n <- 150
#' x1 <- runif(n)
#' x2 <- runif(n)
#' X <- cbind(x1, x2)
#' eps <- rnorm(n,0,sd=0.15)
#' regresion <- function.g1(x1) + function.g2(x2)
#' y <- regresion + eps
#' bandw <- matrix(0.25,2,2)
#' set.seed(8090)
#' nQ <- 80
#' Qmeasure <- matrix(runif(nQ*2), nQ, 2)
#' fit.cl <- margint.cl(y ~ X, windows=bandw, type='alpha', degree=1, Qmeasure=Qmeasure)
#'
#' @export
margint.cl <- function(formula, data, subset, point=NULL, windows, epsilon=1e-6, prob=NULL,
type='0', degree=NULL, qderivate=FALSE, orderkernel=2,
Qmeasure=NULL) {
#AUX <- get_all_vars(formula)
#yp <- AUX[,1]
#Xp <- AUX[,-1]
cl <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset"), names(mf), 0)
mf <- mf[c(1, m)]
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms") # allow model.frame to update it
yp <- model.response(mf, "numeric")
Xp <- model.matrix(mt, mf, NULL)
# above line typically is "model.matrix(mt, mf, contrasts)" and contrasts equals NULL
# remove the default intercept, it is estimated separately
if( all( Xp[, 1] == 1 ) ) Xp <- Xp[, -1]
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){stop("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,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]
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))) {
if(q==1){
prediccion <- mpunto <- as.matrix(punto)
}else{
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))){
if(q==1){
prediccion <- mpunto <- as.matrix(punto)
}else{
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, formula=formula)
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, formula=formula)
class(object) <- c("margint.cl", "margint", "list")
return(object)
}
} else {
object <- list(mu=alpha,g.matrix=g.matriz, prediction=prediccion, Xp=Xp, yp=yp, formula=formula)
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, formula=formula)
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, formula=formula)
class(object) <- c("margint.cl", "margint", "list")
return(object)
}
} else {
object <- list(mu=alpha,g.matrix=g.matriz, Xp=Xp, yp=yp, formula=formula)
class(object) <- c("margint.cl", "margint", "list")
return(object)
}
}
}
#' Robust marginal integration procedures for additive models
#'
#' This function computes robust marginal integration procedures for additive models.
#'
#' This function computes three types of robust marginal integration procedures for additive models.
#'
#' @param formula an object of class \code{formula} (or one that can be coerced to that class): a symbolic description of the model to be fitted.
#' @param data an optional data frame, list or environment (or object coercible by \link{as.data.frame} to a data frame) containing the variables in the model. If not found in \code{data}, the variables are taken from \code{environment(formula)}, typically the environment from which the function was called.
#' @param subset an optional vector specifying a subset of observations to be used in the fitting process.
#' @param point a matrix of points where predictions will be computed and returned.
#' @param windows a vector or a squared matrix of bandwidths for the smoothing estimation procedure.
#' @param prob a vector of probabilities of observing each response (n). Defaults to \code{NULL}.
#' @param sigma.hat estimate of the residual standard error. If \code{NULL} we use the mad of the residuals obtained with local medians.
#' @param win.sigma a vector of bandwidths for estimating sigma.hat. If \code{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 \code{'0'} (local constant on all the covariates), type \code{'1'} (local linear smoother on all the covariates), type \code{'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 \code{'alpha'}. Defaults to \code{NULL} for the case when using estimators of type \code{'0'} or \code{'1'}.
#' @param typePhi one of either \code{'Tukey'} or \code{'Huber'}.
#' @param k.h tuning constant for a Huber-type loss function. Defaults to \code{1.345}.
#' @param k.t tuning constant for a Tukey-type loss function. Defaults to \code{4.685}.
#' @param max.it maximum number of iterations for the algorithm.
#' @param qderivate if TRUE, it calculates \code{g^(q+1)/(q+1)!} for each component only for the type \code{'alpha'} method. Defaults to \code{FALSE}.
#' @param orderkernel order of the kernel used in the nuisance directions when using the estimator of type \code{'alpha'}. Defaults to \code{2}.
#' @param Qmeasure a matrix of points where the integration procedure ocurrs. Defaults to \code{NULL} for calcuting the integrals over the sample.
#'
#' @return A list with the following components:
#' \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 \code{'alpha'} method.}
#' \item{g.derivative }{Matrix of estimated derivatives of the additive components (only when qderivate is \code{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 \code{TRUE}).}
#' \item{Xp}{Matrix of explanatory variables.}
#' \item{yp}{Vector of responses.}
#' \item{formula}{Model formula}
#'
#' @author Alejandra Martinez, \email{ale_m_martinez@hotmail.com}, Matias Salibian-Barrera
#'
#' @references Boente G. and Martinez A. (2017). Marginal integration M-estimators for additive models. TEST, 26(2), 231-260. https://doi.org/10.1007/s11749-016-0508-0
#'
#' @examples
#' function.g1 <- function(x1) 24*(x1-1/2)^2-2
#' function.g2 <- function(x2) 2*pi*sin(pi*x2)-4
#' set.seed(140)
#' n <- 150
#' x1 <- runif(n)
#' x2 <- runif(n)
#' X <- cbind(x1, x2)
#' eps <- rnorm(n,0,sd=0.15)
#' regresion <- function.g1(x1) + function.g2(x2)
#' y <- regresion + eps
#' bandw <- matrix(0.25,2,2)
#' set.seed(8090)
#' nQ <- 80
#' Qmeasure <- matrix(runif(nQ*2), nQ, 2)
#' fit.rob <- margint.rob(y ~ X, windows=bandw, type='alpha', degree=1, Qmeasure=Qmeasure)
#'
#' @export
margint.rob <- function(formula, data, subset, 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){
#AUX <- get_all_vars(formula)
#yp <- AUX[,1]
#Xp <- AUX[,-1]
cl <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset"), names(mf), 0)
mf <- mf[c(1, m)]
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms") # allow model.frame to update it
yp <- model.response(mf, "numeric")
Xp <- model.matrix(mt, mf, NULL)
# above line typically is "model.matrix(mt, mf, contrasts)" and contrasts equals NULL
# remove the default intercept, it is estimated separately
if( all( Xp[, 1] == 1 ) ) Xp <- Xp[, -1]
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){stop("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)
}
sigma.hat <- mad(yp - ab)
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))){
if(q==1){
prediccion <- mpunto <- as.matrix(punto)
}else{
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))){
if(q==1){
prediccion <- mpunto <- as.matrix(punto)
}else{
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, formula=formula, typePhi=typePhi)
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, formula=formula, typePhi=typePhi)
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, formula=formula, typePhi=typePhi)
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, formula=formula, typePhi=typePhi)
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, formula=formula, typePhi=typePhi)
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, formula=formula, typePhi=typePhi)
class(object) <- c("margint.rob", "margint", "list")
return(object)
}
}
}
pos.estimate <- function(y,ini=NULL,sigma.hat, epsilon=1e-6, iter.max=10,typePhi){
yp <- y[ tmp<-!is.na(y) ]
if(is.null(ini)){
ini <- median(yp)
}
corte <- 10
iter <- 0
n <- length(yp)
prob <- rep(1,n)
k.h <- 1.345
k.t <- 4.685
if(typePhi=='Huber'){
beta <- .C("huber_pos", as.integer(n), as.double(yp), as.double(ini), as.double(epsilon),
as.double(sigma.hat), as.double(prob), as.double(k.h), as.integer(iter.max), salida=as.double(0) )$salida
}
if(typePhi=='Tukey'){
beta <- .C("tukey_pos", as.integer(n), as.double(yp), as.double(ini), as.double(epsilon),
as.double(sigma.hat), as.double(prob), as.double(k.t), as.integer(iter.max), salida=as.double(0) )$salida
}
return(beta)
}
#- S3 Methods -#
#' Residuals for objects of class \code{margint}
#'
#' 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
#' \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. Currently ignored.
#'
#' @return A vector of residuals.
#'
#' @author Alejandra Mercedes Martinez \email{ale_m_martinez@hotmail.com}
#'
#' @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. Currently ignored.
#'
#' @return A vector of fitted values.
#'
#' @author Alejandra Mercedes Martinez \email{ale_m_martinez@hotmail.com}
#'
#' @export
predict.margint <- function(object, ...){
return( 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{margint.cl} or \code{margint.rob}.
#' @param ... additional other arguments. Currently ignored.
#'
#' @return A vector of fitted values.
#'
#' @author Alejandra Mercedes Martinez \email{ale_m_martinez@hotmail.com}
#'
#' @method fitted.values margint
#' @export
fitted.values.margint <- function(object,...){
UseMethod("fitted")
}
#' @export
fitted.margint <- function(object,...){
return(predict(object))
}
#' Diagnostic plots for objects of class \code{margint}
#'
#' Plot method for class \code{margint}.
#'
#' @param x an object of class \code{margint}, a result of a call to \code{\link{margint.cl}} or \code{\link{margint.rob}}.
#' @param derivative if TRUE, it plots the q-th derivatives. Defaults to FALSE.
#' @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{ale_m_martinez@hotmail.com}
#'
#' @examples
#' function.g1 <- function(x1) 24*(x1-1/2)^2-2
#' function.g2 <- function(x2) 2*pi*sin(pi*x2)-4
#' set.seed(140)
#' n <- 150
#' x1 <- runif(n)
#' x2 <- runif(n)
#' X <- cbind(x1, x2)
#' eps <- rnorm(n,0,sd=0.15)
#' regresion <- function.g1(x1) + function.g2(x2)
#' y <- regresion + eps
#' bandw <- matrix(0.25,2,2)
#' set.seed(8090)
#' nQ <- 80
#' Qmeasure <- matrix(runif(nQ*2), nQ, 2)
#' fit.rob <- margint.rob(y ~ X, windows=bandw, type='alpha', degree=1, Qmeasure=Qmeasure)
#' plot(fit.rob, which=1)
#'
#' @export
plot.margint <- function(x, derivative=FALSE, which=1:np, ask=FALSE,...){
object <- x
Xp <- object$Xp
np <- dim(Xp)[2]
opar <- par(ask=ask)
on.exit(par(opar))
these <- rep(FALSE, np)
these[ which ] <- TRUE
if(!derivative){
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[,i], res, pch=20,col='gray45',main="",xlab=x_name,ylab=y_name, ylim=lim_cl,cex.lab=0.8)
lines(Xp[ord,i],object$g.matrix[ord,i],lwd=3)
}
}
}else{
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('g'[.(i)]^('q')))
plot(Xp[ord,i], object$g.derivate[ord,i], type='l',lwd=3, main="",xlab=x_name,ylab=y_name,cex.lab=0.8)
}
}
}
}
#Multiple R-squared
R2 <- function(object,...){
yp <- object$yp
y <- yp[ tmp<-!is.na(yp) ]
n <- length(y)
res <- residuals(object)
S02 <- sum((y-mean(y))^2)
S2 <- sum(res^2)
R2 <- (S02-S2)/S02
return(R2)
}
#Robust multiple R-squared
R2.rob <- function(object,...){
yp <- object$yp
y <- yp[ tmp<-!is.na(yp) ]
n <- length(y)
S02 <- 0
S2 <- 0
res <- residuals(object)
sigma.hat <- object$sigma.hat
typePhi <- object$typePhi
pos <- pos.estimate(y,ini=NULL,sigma.hat, epsilon=1e-6, iter.max=50,typePhi=typePhi)
if(typePhi=='Tukey'){
for(i in 1:n){
S02 <- S02 + rho.tukey((y[i]-pos)/sigma.hat)
S2 <- S2 + rho.tukey(res[i]/sigma.hat)
}
}
if(typePhi=='Huber'){
for(i in 1:n){
S02 <- S02 + rho.huber((y[i]-pos)/sigma.hat)
S2 <- S2 + rho.huber(res[i]/sigma.hat)
}
}
R2.rob <- (S02-S2)/S02
return(R2.rob)
}
#' 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{ale_m_martinez@hotmail.com}
#'
#' @export
#' @aliases summary.margint summary.margint.cl summary.margint.rob
summary.margint <- function(object,...){
NextMethod()
}
#' @export
summary.margint.cl <- function(object,...){
message("Estimate of the intercept: ", round(object$mu,5))
message("Multiple R-squared: ", round(R2(object),5))
res <- residuals(object)
message("Residuals:")
summary(res)
}
#' @export
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))
message("Robust multiple R-squared: ", round(R2.rob(object),5))
res <- residuals(object)
message("Residuals:")
summary(res)
}
#' Deviance for objects of class \code{margint}
#'
#' This function returns the deviance of the fitted additive model using one of the three
#' classical or robust marginal integration estimators, as 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. Currently ignored.
#'
#' @return A real number.
#'
#' @author Alejandra Mercedes Martinez \email{ale_m_martinez@hotmail.com}
#'
#' @export
deviance.margint <- function(object, ...){
NextMethod()
}
deviance.margint.cl <- function(object, ...){
return( sum( (residuals(object))^2) )
}
deviance.margint.rob <- function(object, ...){
yp <- object$yp
y <- yp[ tmp<-!is.na(yp) ]
n <- length(y)
S2 <- 0
res <- residuals(object)
sigma.hat <- object$sigma.hat
typePhi <- object$typePhi
if(typePhi=='Tukey'){
for(i in 1:n){
S2 <- S2 + rho.tukey(res[i]/sigma.hat)
}
}
if(typePhi=='Huber'){
for(i in 1:n){
S2 <- S2 + rho.huber(res[i]/sigma.hat)
}
}
return( S2 )
}
#' Additive model formula
#'
#' Description of the additive model formula extracted from an object of class \code{margint}.
#'
#' @param x 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. Currently ignored.
#'
#' @return A model formula.
#'
#' @author Alejandra Mercedes Martinez \email{ale_m_martinez@hotmail.com}
#'
#' @export
formula.margint <- function(x, ...){
return(x$formula )
}
#' Print a Marginal Integration procedure
#'
#' The default print method for a \code{margint} object.
#'
#' @param x 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. Currently ignored.
#'
#' @return A real number.
#'
#' @author Alejandra Mercedes Martinez \email{ale_m_martinez@hotmail.com}
#'
#' @export
print.margint <- function(x, ...){
cat("Formula:\n")
print(x$formula)
#cat("\n")
}
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.