# Quantile Regression for Nonlinear Mixed-Effects Models

### Description

Performs a quantile regression for a NLMEM using the Stochastic-Approximation of the EM Algorithm (SAEM) for an unique or a set of quantiles.

### Usage

 1 2 QRNLMM(y,x,groups,initial,exprNL,covar=NA,p=0.5,precision=0.0001,MaxIter=500, M=20,cp=0.25,beta=NA,sigma=NA,Psi=NA,show.convergence=TRUE,CI=95) 

### Arguments

 y the response vector of dimension N where N is the total of observations. x vector of longitudinal (repeated measures) covariate of dimension N. For example: Time, location, etc. groups factor of dimension N specifying the partitions of the data over which the random effects vary. initial an numeric vector, or list of initial estimates for the fixed effects. It must be provide adequately (see details section) in order to ensure a proper convergence. exprNL expression containing the proposed nonlinear function. It can be of class character or expression. It must have a defined structure defined in the details section in order to be correctly read by the derivate R function deriv. covar a matrix of dimension N \times r where r represents the number of covariates. p unique quantile or a set of quantiles related to the quantile regression. precision the convergence maximum error. MaxIter the maximum number of iterations of the SAEM algorithm. Default = 500. M Number of Monte Carlo simulations used by the SAEM Algorithm. Default = 20. For more accuracy we suggest to use M=20. cp cut point (0 ≤ cp ≤ 1) which determines the percentage of initial iterations with no memory. beta fixed effects vector of initial parameters, if desired. sigma dispersion initial parameter for the error term, if desired. Psi Variance-covariance random effects matrix of initial parameters, if desired. show.convergence if TRUE, it will show a graphical summary for the convergence of the estimates of all parameters for each quantile in order to assess the convergence. CI Confidence to be used for the Confidence Interval when a grid of quantiles is provided. Default=95.

### Details

This algorithm performs the SAEM algorithm proposed by Delyon et al. (1999), a stochastic version of the usual EM Algorithm deriving exact maximum likelihood estimates of the fixed-effects and variance components. Covariates are allowed, the longitudinal (repeated measures) coded x and a set of covariates covar.

About initial values: Estimation for fixed effects parameters envolves a Newton-Raphson step. In adition, NL models are highly sensitive to initial values. So, we suggest to set of intial values quite good, this based in the parameter interpretation of the proposed NL function.

About the nonlinear expression: For the NL expression exprNL just the variables x, covar, fixed and random can be defined. Both x and covar represent the covariates defined above. The fixed effects must be declared as fixed[1], fixed[2],..., fixed[d] representing the first, second and dth fixed effect. Exactly the same for the random effects and covariates where the term fixed should be replace for random and covar respectively.

For instance, if we use the exponential nonlinear function with two parameters, each parameter represented by a fixed and a random effect, this will be defined by

y_{ij} = (β_1 + b_1)\exp^{-(β_2 + b_2)x_{ij}}

and the exprNL should be a character or and expression defined by

exprNL = "(fixed[1]+random[1])*exp(-(fixed[2]+random[2])*x)"

or

exprNL = expression((fixed[1]+random[1])*exp(-(fixed[2]+random[2])*x)).

If we are interested in adding two covariates in order to explain on of the parameters, the covariates covar[1] and covar[2] must be included in the model. For example, for the nonlinear function

y_{ij} = (β_1 + β_3*covar1_{ij} + b_1)\exp^{-(β_2 + β_4* covar2_{ij} + b_2)x_{ij}}

the exprNL should be

exprNL = "(fixed[1]+fixed[3]*covar[1]+random[1])*exp(-(fixed[2]+fixed[4]*covar[2]+random[2])*x)"

or

exprNL = expression((fixed[1]+fixed[3]*covar[1]+random[1])*exp(-(fixed[2]+ fixed[4]*covar[2]+random[2])*x)).

Note that the mathematical function exp was used. For derivating the deriv R function recognizes in the exprNL expression the arithmetic operators +, -, *, / and ^, and the single-variable functions exp, log, sin, cos, tan, sinh, cosh, sqrt, pnorm, dnorm, asin, acos, atan, gamma, lgamma, digamma and trigamma, as well as psigamma for one or two arguments (but derivative only with respect to the first).

General details: When a grid of quantiles is provided, a graphical summary with point estimates and Confidence Intervals for model parameters is shown and also a graphical summary for the convergence of these estimates (for each quantile), if show.convergence=TRUE.

If the convergence graphical summary shows that convergence has not be attained, it's suggested to increase the total number of iterations MaxIter.

About the cut point parameter cp, a number between 0 and 1 (0 ≤ cp ≤ 1) will assure an initial convergence in distribution to a solution neighborhood for the first cp*MaxIter iterations and an almost sure convergence for the rest of the iterations. If you do not know how SAEM algorithm works, these parameters SHOULD NOT be changed.

This program uses progress bars that will close when the algorithm ends. They must not be closed before, if not, the algorithm will stop.

### Value

The function returns a list with two objects

 conv A two elements list with the matrices teta and se containing the point estimates and standard error estimate for all parameters along all iterations.

The second element of the list is res, a list of 12 elements detailed as

 iter number of iterations. criteria attained criteria value. nlmodel the proposed nonlinear function. beta fixed effects estimates. sigma scale parameter estimate for the error term. Psi Random effects variance-covariance estimate matrix. SE Standard Error estimates. table Table containing the inference for the fixed effects parameters. loglik Log-likelihood value. AIC Akaike information criterion. BIC Bayesian information criterion. HQ Hannan-Quinn information criterion. time processing time.

### Note

If a grid of quantiles was provided, the result is a list of the same dimension where each element corresponds to each quantile as detailed above.

### Author(s)

Christian E. Galarza <cgalarza88@gmail.com> and Victor H. Lachos <hlachos@ime.unicamp.br>

### References

Galarza, C.M., Castro, LM. Louzada, F. and V. H. Lachos (2016) Quantile Regression for Nonlinear Mixed Effects Models: A Likelihood Based Perspective. Technical Report 8, Universidade Estadual de Campinas. http://www.ime.unicamp.br/sites/default/files/rp08-2016.pdf

Delyon, B., Lavielle, M. & Moulines, E. (1999). Convergence of a stochastic approximation version of the EM algorithm. Annals of Statistics, pages 94-128.

Soybean, HIV, lqr, QRLMM, group.plots
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 ## Not run: #Using the Soybean data data(Soybean) attach(Soybean) ################################# #A full model (no covariate) y = weight #response x = Time #time #Expression for the three parameter logistic curve exprNL = expression((fixed[1]+random[1])/ (1 + exp(((fixed[2]+random[2])- x)/(fixed[3]+random[3])))) #Initial values for fixed effects initial = c(max(y),0.6*max(y),0.73*max(y)) #A median regression (by default) median_reg = QRNLMM(y,x,Plot,initial,exprNL) #Assing the fit fxd = median_reg$res$beta nlmodel = median_reg$res$nlmodel seqc = seq(min(x),max(x),length.out = 500) group.plot(x = Time,y = weight,groups = Plot,type="l", main="Soybean profiles",xlab="time (days)", ylab="mean leaf weight (gr)",col="gray") lines(seqc,nlmodel(x = seqc,fixed = fxd,random = rep(0,3)), lwd=2,col="blue") ######################################### #A model for compairing the two genotypes y = weight #response x = Time #time covar = c(Variety)-1 #factor genotype (0=Forrest, 1=Plan Introduction) #Expression for the three parameter logistic curve with a covariate exprNL = expression((fixed[1]+(fixed[4]*covar[1])+random[1])/ (1 + exp(((fixed[2]+random[2])- x)/(fixed[3]+random[3])))) #Initial values for fixed effects initial = c(max(y),0.6*max(y),0.73*max(y),3) # A quantile regression for the three quartiles box_reg = QRNLMM(y,x,Plot,initial,exprNL,covar,p=c(0.25,0.50,0.75)) #Assing the fit for the median (second quartile) fxd = box_reg[[2]]$res$beta nlmodel = box_reg[[2]]$res$nlmodel seqc = seq(min(x),max(x),length.out = 500) group.plot(x = Time[Variety=="P"],y = weight[Variety=="P"], groups = Plot[Variety=="P"],type="l",col="light blue", main="Soybean profiles by genotype",xlab="time (days)", ylab="mean leaf weight (gr)") group.lines(x = Time[Variety=="F"],y = weight[Variety=="F"], groups = Plot[Variety=="F"],col="gray") lines(seqc,nlmodel(x = seqc,fixed = fxd,random = rep(0,3),covar=1), lwd=2,col="blue") lines(seqc,nlmodel(x = seqc,fixed = fxd,random = rep(0,3),covar=0), lwd=2,col="black") ######################################### #A simple output example --------------------------------------------------- Quantile Regression for Nonlinear Mixed Model --------------------------------------------------- Quantile = 0.5 Subjects = 48 ; Observations = 412 - Nonlinear function function(x,fixed,random,covar=NA){ resp = (fixed[1] + random[1])/(1 + exp(((fixed[2] + random[2]) - x)/(fixed[3] + random[3]))) return(resp)} ----------- Estimates ----------- - Fixed effects Estimate Std. Error z value Pr(>|z|) beta 1 18.80029 0.53098 35.40704 0 beta 2 54.47930 0.29571 184.23015 0 beta 3 8.25797 0.09198 89.78489 0 sigma = 0.31569 Random effects Variance-Covariance Matrix matrix b1 b2 b3 b1 24.36687 12.27297 3.24721 b2 12.27297 15.15890 3.09129 b3 3.24721 3.09129 0.67193 ------------------------ Model selection criteria ------------------------ Loglik AIC BIC HQ Value -622.899 1265.798 1306.008 1281.703 ------- Details ------- Convergence reached? = FALSE Iterations = 300 / 300 Criteria = 0.00058 MC sample = 20 Cut point = 0.25 Processing time = 22.83885 mins ## End(Not run)