I am using a modified quantile regression function found in a scientific paper. (source:

Here is what I would like to plot more or less from the models made with @Parfait solution, the intercept and the explanatory variable evolution according to Tau (the quantiles) on the x-axis.

CodePudding user response：

Consider encapsulating the author's entire processing in a user-defined method that receives data (swiss from authors) as input parameter along with other variables including formula (Fertility ~ . from authors) and response column ("Fertility" from authors).

Then, pass data subsets with group_by. Also, the authors for loop can be refactored to a vectorized loop such as sapply or vapply since the return is a numeric vector.

Generalized Functions

minimize.logcosh <- function(par, X, y, tau) {
    diff <- y-(X %*% par)
    check <- (tau-0.5)*diff (0.5/0.7)*logcosh(0.7*diff) 0.4
    return(sum(check))
}

smrq <- function(X, y, tau){
    p <- ncol(X)
    op.result <- optim(
        rep(0, p),
        fn = minimize.logcosh,
        method = 'BFGS',
        X = X,
        y = y,
        tau = tau
    )
    beta <- op.result$par
    return(beta)
}

run_smrq <- function(data, fml, response) {
    x <- model.matrix(fml, data)[,-1]
    y <- data[[response]]
    X <- cbind(x, rep(1,nrow(x)))

    n <- 99
    betas <- sapply(1:n, function(i) smrq(X, y, tau=i/(n 1)))
    # betas <- vapply(1:n, function(i) smrq(X, y, tau=i/(n 1)), numeric(1))
    return(betas)        
}

Callers

Test authors' example

swiss <- datasets::swiss
smrq_models <- run_smrq(data=swiss, fml=Fertility~., response="Fertility")

dplyr (using group_map)

smrq_models <- data %>%
  group_by(latitude) %>%
  group_map(~ run_smrq(data=., fml=julian_day~year, response="julian_day")

base (using by, the object-oriented wrapper to tapply)

smrq_models <- by(
   data, 
   data$latitude, 
   function(sub) run_smrq(data=sub, fml=julian_day~year, response="julian_day")
)