Nothing
R/ARpMMEC.est.R
In ARpLMEC: Censored Mixed-Effects Models with Different Correlation Structures
Defines functions
ARpMMEC.est
Documented in
ARpMMEC.est
#' @title Censored Mixed-Effects Models with Autoregressive Correlation Structure and DEC for Normal and t-Student Errors
#' @import numDeriv
#' @import TruncatedNormal
#' @import LaplacesDemon
#' @import tcltk
#' @import MASS
#' @import stats
#' @import relliptical
#' @import expm
#' @description This functino fits left, right or intervalar censored mixed-effects linear model, with autoregressive errors of order \code{p}, using the EM algorithm. It returns estimates, standard errors and prediction of future observations.
#' @param y Vector \code{1 x n} of censored responses, where \code{n} is the sum of the number of observations of each individual
#' @param x Design matrix of the fixed effects of order \code{n x s}, corresponding to vector of fixed effects.
#' @param z Design matrix of the random effects of order\code{n x b}, corresponding to vector of random effects.
#' @param cc Vector of censoring indicators of length \code{n}, where \code{n} is the total of observations. For each observation: \code{0} if non-censored, \code{1} if censored.
#' @param nj Vector \code{1 x m} with the number of observations for each subject, where \code{m} is the total number of individuals.
#' @param tt Vector \code{1 x n} with the time the measurements were made, where \code{n} is the total number of measurements for all individuals. Default it's considered regular times.
#' @param struc \code{UNC},\code{ARp},\code{DEC},\code{SYM} or \code{DEC(AR)} for uncorrelated ,autoregressive, DEC(phi1,phi2), DEC(phi1,phi2=1), DEC(DEC(phi1,phi2=1)) structure, respectively
#' @param order Order of the autoregressive process. Must be a positive integer value.
#' @param initial List with the initial values in the next orden: betas,sigma2,alphas,phi and nu. If it is not indicated it will be provided automatically. Default is \code{NULL}
#' @param nu.fixed Logical. Should estimate the parameter "nu" for the t-student distribution?. If is False indicates the value in the list of initial values. Default is \code{FALSE}
#' @param typeModel \code{Normal} for Normal distribution and \code{Student} for t-Student distribution. Default is \code{Normal}
#' @param cens.type \code{left} for left censoring, \code{right} for right censoring and \code{interval} for intervalar censoring. Default is \code{left}
#' @param LI Vector censoring lower limit indicator of length \code{n}. For each observation: \code{0} if non-censored, \code{-inf} if censored. It is only indicated for when \code{cens.type} is \code{both}. Default is \code{NULL}
#' @param LS Vector censoring upper limit indicator of length \code{n}. For each observation: \code{0} if non-censored, \code{inf} if censored.It is only indicated for when \code{cens.type} is \code{both}. Default is \code{NULL}
#' @param MaxIter The maximum number of iterations of the EM algorithm. Default is \code{200}
#' @param error The convergence maximum error. Default is \code{0.0001}
#' @param Prev Indicator of the prediction process. Available at the moment only for the \code{typeModel=normal} case. Default is \code{FALSE}
#' @param isubj Vector indicator of subject included in the prediction process. Default is \code{NULL}
#' @param step Number of steps for prediction. Default is \code{NULL}
#' @param xpre Design matrix of the fixed effects to be predicted. Default is \code{NULL}.
#' @param zpre Design matrix of the random effects to be predicted. Default is \code{NULL}.
#' @return returns list of class \dQuote{ARpMMEC}:
#' \item{FixEffect}{Data frame with: estimate, standar errors and confidence intervals of the fixed effects.}
#' \item{Sigma2}{Data frame with: estimate, standar errors and confidence intervals of the variance of the white noise process.}
#' \item{Phi}{Data frame with: estimate, standar errors and confidence intervals of the autoregressive parameters.}
#' \item{RandEffect}{Data frame with: estimate, standar errors and confidence intervals of the random effects.}
#' \item{nu}{the parameter "nu" for the t-student distribution}
#' \item{Est}{Vector of parameters estimate (fixed Effects, sigma2, phi, random effects).}
#' \item{SE}{Vector of the standard errors of (fixed Effects, sigma2, phi, random effects).}
#' \item{Residual}{Vector of the marginal residuals.}
#' \item{loglik}{Log-likelihood value.}
#' \item{AIC}{Akaike information criterion.}
#' \item{BIC}{Bayesian information criterion.}
#' \item{AICc}{Corrected Akaike information criterion.}
#' \item{iter}{Number of iterations until convergence.}
#' \item{Yfit}{Vector "y" fitted}
#' \item{MI}{Information matrix}
#' \item{Prev}{Predicted values (if xpre and zpre is not \code{NULL}).}
#' \item{time}{Processing time.}
#' \item{others}{The first and second moments of the random effect and vector Y}
#' @references Olivari, R. C., Garay, A. M., Lachos, V. H., & Matos, L. A. (2021). Mixed-effects
#' models for censored data with autoregressive errors. Journal of Biopharmaceutical Statistics, 31(3), 273-294.
#' \doi{10.1080/10543406.2020.1852246}
#' @examples
#' \dontrun{
#'p.cens = 0.1
#'m = 10
#'D = matrix(c(0.049,0.001,0.001,0.002),2,2)
#'sigma2 = 0.30
#'phi = 0.6
#'beta = c(1,2,1)
#'nj=rep(4,10)
#'tt=rep(1:4,length(nj))
#'x<-matrix(runif(sum(nj)*length(beta),-1,1),sum(nj),length(beta))
#'z<-matrix(runif(sum(nj)*dim(D)[1],-1,1),sum(nj),dim(D)[1])
#'data=ARpMMEC.sim(m,x,z,tt,nj,beta,sigma2,D,phi,struc="ARp",typeModel="Normal",p.cens=p.cens)
#'
#'teste1=ARpMMEC.est(data$y_cc,x,z,tt,data$cc,nj,struc="ARp",order=1,typeModel="Normal",MaxIter = 2)
#'teste2=ARpMMEC.est(data$y_cc,x,z,tt,data$cc,nj,struc="ARp",order=1,typeModel="Student",MaxIter = 2)
#'
#'xx=matrix(runif(6*length(beta),-1,1),6,length(beta))
#'zz=matrix(runif(6*dim(D)[1],-1,1),6,dim(D)[1])
#'isubj=c(1,4,5)
#'teste3=ARpMMEC.est(data$y_cc,x,z,tt,data$cc,nj,struc="ARp",order=1,typeModel="Normal",
#' MaxIter = 2,Prev=TRUE,step=2,isubj=isubj,xpre=xx,zpre=zz)
#'teste3$Prev
#'
#' }
#'
#'
#' @export
#'
#'
#'
ARpMMEC.est=function(y,x,z,tt,cc,nj,struc="UNC",order=1, initial=NULL,nu.fixed=TRUE,
typeModel="Normal",cens.type="left", LI=NULL,LS=NULL, MaxIter=200,
error=0.0001, Prev=FALSE,step=NULL,isubj=NULL,xpre=NULL,zpre=NULL)
{
m<-length(y); N<-sum(nj); p<-dim(x)[2]; q1<-dim(z)[2]; m1<-m*p; m2<-m*q1
if(is.matrix(y)) y <- y[as.vector(!is.na(as.vector(t(y))))]
if(is.matrix(cc)) cc <- cc[as.vector(!is.na(as.vector(t(cc))))]
if (!is.matrix(x)) x=as.matrix(x)
if (!is.matrix(z)) x=as.matrix(z)
if( is.matrix(nj)) nj <- nj[as.vector(!is.na(as.vector(t(nj))))]
if(!is.numeric(y)) stop("y must be a numeric vector. Check documentation!")
if(sum(is.na(y))>0) stop("Vector y does not support NA values.")
if(!is.vector(y)) stop("y must be a vector.Check documentation!")
if(length(y)!=nrow(as.matrix(x))) stop("x does not have the same number of lines than y.")
if(length(y)!=length(cc)) stop("cc does not have the same length than y.")
if(length(y)!=nrow(as.matrix(z))) stop("x does not have the same number of lines than y.")
if(length(y)!=sum(nj)) stop("not compatible sizes between the response y and the repetited measures nj")
if(length(y)==0) stop("The parameter y must be provided.")
if(length(y)!=length(tt)) stop("not compatible sizes between the response y and the vector time tt")
if(!is.numeric(x)) stop("x must be a numeric matrix. Check documentation!")
if(sum(is.na(x))>0) stop("There are some NA values in x.")
if(!is.matrix(x)) stop("x must be a matrix. Check documentation!")
if(det(t(x)%*%x)==0) stop("the columns of x must be linearly independent.")
if(length(x)==0) stop("The parameter x must be provided.")
if(!is.numeric(z)) stop("z must be a numeric matrix. Check documentation!")
if(!is.matrix(z)) stop("z must be a matrix. Check documentation!")
if(sum(is.na(z))>0) stop("There are some NA values in z.")
if(length(z)==0) stop("The parameter z must be provided.")
if(!is.numeric(cc)) stop("cc must be a numeric vector. Check documentation!")
if(!is.vector(cc)) stop("cc must be a vector.Check documentation!")
if(sum(is.na(cc))>0) stop("There are some NA values in cc.")
if(sum(cc%in%c(0,1))<length(cc)) stop("The elements of the vector cc must be 0 or 1.")
if(length(cc)==0) stop("The parameter cc must be provided.")
if(!is.numeric(nj)) stop("nj must be a numeric vector. Check documentation!")
if(!is.vector(nj)) stop("nj must be a vector. Check documentation!")
if(sum(is.na(nj))>0) stop("There are some NA values in nj")
if(length(nj)==0) stop("The parameter nj must be provided.")
if(struc!="DEC"&struc!="DEC(AR)"&struc!="SYM"&struc!="ARp"&struc!="UNC") stop("Struc must be UNC, DEC, DEC(AR), SYM or ARp. Check documentation!")
if(struc=="ARp"){
if(!is.numeric(order) ) stop("Orde must be a number. Check documentation!")
if(length(order)!=1) stop("Order must be a value.")
if(is.numeric(order))
{ if(order!=round(order)|order<=0) stop("Order must be a positive integer value.")}}
if(!is.null(initial))
{ if(!is.null(initial$betas))
{if(!is.numeric(initial$betas)) stop("betas must be a numeric vector. Check documentation!")
if(!is.vector(initial$betas)) stop("betas must be a vector. Check documentation!")
if(length(initial$betas)!=ncol(x)) stop("not compatible sizes between the matrix x and parameter betas.")}
if(!is.null(initial$sigma2))
{if(!is.numeric(initial$sigma2)) stop("sigma2 must be a scalar.")
if(length(initial$sigma2)>1) stop("sigma2 must be a scalar.")}
if(!is.null(initial$alphas))
{if(!is.matrix(initial$alphas)) stop("alphas must be a matrix.")
if(initial$alphas[upper.tri(initial$alphas)]!=initial$alphas[lower.tri(initial$alphas)])stop("alphas must be a simetric matrix.")
if(dim(initial$alphas)[2]!=ncol(z)) stop("not compatible sizes between the matrix z and parameter alphas.")}
if(struc=="ARp"){
if(!is.null(initial$phi))
{if(!is.numeric(initial$phi)) stop("phi must be a numeric vector. Check documentation!")
if(length(initial$phi)!=order) stop("not compatible sizes between the value Arp and parameter phi. Check documentation!")}
}
}
if(typeModel!='Normal'& typeModel!='Student') stop('typeModel must be Normal or Student. Check documentation!')
if(cens.type!="left" & cens.type!="right" & cens.type!="interval")stop('cens.type must be left, right or interval. Check documentation!')
if(cens.type=="interval"&is.null(LI)) stop("The parameter LI must be provided.. Check documentation!")
if(cens.type=="interval"&is.null(LS)) stop("The parameter LS must be provided.. Check documentation!")
if(!is.null(LI)&!is.numeric(LI)) stop("LI must be a numeric vector. Check documentation!")
if(!is.null(LS)&!is.numeric(LS)) stop("LS must be a numeric vector. Check documentation!")
if(length(LS)!=length(LI)) stop("not compatible sizes between the vectors LI and LS. Check documentation!")
if(cens.type=="interval") {
if(length(y)!=length(LI)) stop("not compatible sizes between the vectors y and LI. Check documentation!")
if(length(y)!=length(LS)) stop("not compatible sizes between the vectors y and LS. Check documentation!")
}
if (!is.numeric(MaxIter)) stop("MaxIter must be a positive number. Check documentation!")
if (length(MaxIter) > 1) stop("MaxIter parameter must be a scalar")
if (MaxIter <0) stop("MaxIter parameter must be positive number")
if (!is.numeric(error)) stop("error must be a positive number. Check documentation!")
if (length(error) > 1) stop("error parameter must be a scalar")
if (error <0) stop("error parameter must be positive number")
if (Prev) {
if(is.null(step)|is.null(xpre)|is.null(zpre)|is.null(isubj)) stop("step, isubj, xpre, zpre needs to be provided. Check documentation!")
if (!is.numeric(isubj)) stop("isubj must be a numeric vector. Check documentation!")
if (!is.numeric(step)) stop("step must be a positive number. Check documentation!")
if (step <0) stop("step parameter must be positive number")
if (length(step) > 1) stop("step parameter must be a scalar")
if (ncol(xpre)!=ncol(as.matrix(x))) stop("xpre must have the same number of columns than x")
if (sum(is.na(xpre))>0) stop("There are some NA values in xpre")
if (!is.numeric(xpre)) stop("xpred must be a numeric matrix")
if (ncol(zpre)!=ncol(as.matrix(z))) stop("zpre must have the same number of columns than z")
if (sum(is.na(zpre))>0) stop("There are some NA values in zpre")
if (!is.numeric(zpre)) stop("zpred must be a numeric matrix")
if(nrow(xpre)!=length(isubj)*step) stop("not compatible sizes between xpre and isubj. Check documentation!")
if(nrow(zpre)!=length(isubj)*step) stop("not compatible sizes between zpre and isubj. Check documentation!")
}
if(typeModel=="Normal"){
if(struc=="ARp"){
out<-EMCensArpN(cc=cc,y=y,x=x,z=z,tt=tt,nj=nj, Arp=order, initial=initial, cens.type=cens.type, LI=LI,LS=LS,MaxIter=MaxIter,ee=error,
Prev=Prev,step=step,isubj=isubj ,xpre=xpre,zpre=zpre)}
if(struc=="UNC"){
out<-EMCensArpN(cc=cc,y=y,x=x,z=z,tt=tt,nj=nj, Arp=struc, initial=initial, cens.type=cens.type, LI=LI,LS=LS,MaxIter=MaxIter,ee=error,
Prev=Prev,step=step,isubj=isubj ,xpre=xpre,zpre=zpre)}
if(struc=="DEC"|struc=="DEC(AR)"|struc=="SYM"){
out<-EMCensDECN(cc=cc,y=y,x=x,z=z,tt=tt,nj=nj, struc=struc, initial=initial, cens.type=cens.type, LI=LI,LS=LS,MaxIter=MaxIter,ee=error,
Prev=Prev,step=step,isubj=isubj ,xpre=xpre,zpre=zpre)}
}
if(typeModel=="Student"){
if(struc=="ARp"){
out<-EMCensArpT(cc=cc,y=y,x=x,z=z,ttc=tt,nj=nj,Arp=order,initial=initial,cens.type=cens.type,LL=LI,LU=LS,nu.fixed=nu.fixed,
iter.max=MaxIter,precision=error)
}
if(struc=="UNC"){
out<-EMCensArpT(cc=cc,y=y,x=x,z=z,ttc=tt,nj=nj,Arp=struc,initial=initial,cens.type=cens.type,LL=LI,LU=LS,nu.fixed=nu.fixed,
iter.max=MaxIter,precision=error)
}
if(struc=="DEC"|struc=="DEC(AR)"|struc=="SYM"){
out<-EMCensDECT(cc=cc,y=y,x=x,z=z,ttc=tt,nj=nj,struc=struc,initial=initial,cens.type=cens.type,LL=LI,LU=LS,nu.fixed=nu.fixed,
iter.max=MaxIter,precision=error)
}
}
if(struc=="ARp")
{
cat('\n')
cat('---------------------------------------------------\n')
cat('Autoregressive censored mixed-effects models \n')
cat('---------------------------------------------------\n')
cat('\n')
cat("Autoregressive order =",order)
cat('\n')
cat("Distribution =",typeModel)
cat('\n')
if(typeModel=="Student") cat("nu =",out$nu); cat('\n')
cat("Subjects =",length(nj),";",'Observations =',sum(nj))
cat('\n')
cat('\n')
cat('-----------\n')
cat('Estimates\n')
cat('-----------\n')
cat('\n')
cat('- Fixed effects \n')
cat('\n')
print(out$tableB)
cat('\n')
cat('\n')
cat('- Sigma^2 \n')
cat('\n')
print(out$tableS)
cat('\n')
cat('\n')
cat('- Autoregressives parameters\n')
cat('\n')
print(out$tableP)
cat('\n')
cat('\n')
cat('- Random effects \n')
cat('\n')
print(out$tableA)
cat('\n')
cat('\n')
cat('------------------------\n')
cat('Model selection criteria\n')
cat('------------------------\n')
cat('\n')
critFin <- c(out$loglik, out$AIC, out$BIC)
critFin <- round(t(as.matrix(critFin)),digits=3)
dimnames(critFin) <- list(c("Value"),c("Loglik", "AIC", "BIC"))
print(critFin)
cat('\n')
cat('-------\n')
cat('Details\n')
cat('-------\n')
cat('\n')
cat("Convergence reached? =",(out$iter < MaxIter))
cat('\n')
cat('Iterations =',out$iter,"/",MaxIter)
cat('\n')
cat("Processing time =",out$time,units(out$time))
cat('\n')
}
if(struc!="ARp")
{
cat('\n')
cat('---------------------------------------------------\n')
cat('DEC censored mixed-effects models \n')
cat('---------------------------------------------------\n')
cat('\n')
cat("Case =",struc)
cat('\n')
cat("Distribution =",typeModel)
cat('\n')
if(typeModel=="Student") cat("nu =",out$nu); cat('\n')
cat("Subjects =",length(nj),";",'Observations =',sum(nj))
cat('\n')
cat('\n')
cat('-----------\n')
cat('Estimates\n')
cat('-----------\n')
cat('\n')
cat('- Fixed effects \n')
cat('\n')
print(out$tableB)
cat('\n')
cat('\n')
cat('- Sigma^2 \n')
cat('\n')
print(out$tableS)
cat('\n')
cat('\n')
if(struc!="UNC"){
cat('- Autoregressives parameters\n')
cat('\n')
print(out$tableP)
cat('\n')
cat('\n')}
cat('- Random effects \n')
cat('\n')
print(out$tableA)
cat('\n')
cat('\n')
cat('------------------------\n')
cat('Model selection criteria\n')
cat('------------------------\n')
cat('\n')
critFin <- c(out$loglik, out$AIC, out$BIC)
critFin <- round(t(as.matrix(critFin)),digits=3)
dimnames(critFin) <- list(c("Value"),c("Loglik", "AIC", "BIC"))
print(critFin)
cat('\n')
cat('-------\n')
cat('Details\n')
cat('-------\n')
cat('\n')
cat("Convergence reached? =",(out$iter < MaxIter))
cat('\n')
cat('Iterations =',out$iter,"/",MaxIter)
cat('\n')
cat("Processing time =",out$time,units(out$time))
cat('\n')
}
if(typeModel=="Student") nu<-out$nu
if(typeModel=="Normal") nu<-NULL
obj.out <- list(FixEffect=out$tableB, Sigma2=out$tableS, Phi=out$tableP,RandEffect=out$tableA, nu=nu,
Est=c(out$beta1, sigma2=out$sigmae, phi=out$phi, RnEffect=out$dd), SE=out$SE,Residual=out$residual,
loglik=out$loglik, AIC=out$AIC, BIC=out$BIC, AICc=out$AICcorr, iter=out$iter, Yfit=out$yfit,
MI=out$MI, Prev=out$Prev, time=out$time, others=list(ubi = out$ubi, ubbi = out$ubbi, uybi = out$uybi, uyi = out$uyi, uyyi = out$uyyi,varbeta=out$varbeta,yog=out$yorg) )
class(obj.out) = "ARpMMEC"
return(obj.out)
}
Try the ARpLMEC package in your browser
Any scripts or data that you put into this service are public.
ARpLMEC documentation built on June 27, 2022, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions ARpMMEC.est
Documented in ARpMMEC.est
#' @title Censored Mixed-Effects Models with Autoregressive Correlation Structure and DEC for Normal and t-Student Errors
#' @import numDeriv
#' @import TruncatedNormal
#' @import LaplacesDemon
#' @import tcltk
#' @import MASS
#' @import stats
#' @import relliptical
#' @import expm
#' @description This functino fits left, right or intervalar censored mixed-effects linear model, with autoregressive errors of order \code{p}, using the EM algorithm. It returns estimates, standard errors and prediction of future observations.
#' @param y Vector \code{1 x n} of censored responses, where \code{n} is the sum of the number of observations of each individual
#' @param x Design matrix of the fixed effects of order \code{n x s}, corresponding to vector of fixed effects.
#' @param z Design matrix of the random effects of order\code{n x b}, corresponding to vector of random effects.
#' @param cc Vector of censoring indicators of length \code{n}, where \code{n} is the total of observations. For each observation: \code{0} if non-censored, \code{1} if censored.
#' @param nj Vector \code{1 x m} with the number of observations for each subject, where \code{m} is the total number of individuals.
#' @param tt Vector \code{1 x n} with the time the measurements were made, where \code{n} is the total number of measurements for all individuals. Default it's considered regular times.
#' @param struc \code{UNC},\code{ARp},\code{DEC},\code{SYM} or \code{DEC(AR)} for uncorrelated ,autoregressive, DEC(phi1,phi2), DEC(phi1,phi2=1), DEC(DEC(phi1,phi2=1)) structure, respectively
#' @param order Order of the autoregressive process. Must be a positive integer value.
#' @param initial List with the initial values in the next orden: betas,sigma2,alphas,phi and nu. If it is not indicated it will be provided automatically. Default is \code{NULL}
#' @param nu.fixed Logical. Should estimate the parameter "nu" for the t-student distribution?. If is False indicates the value in the list of initial values. Default is \code{FALSE}
#' @param typeModel \code{Normal} for Normal distribution and \code{Student} for t-Student distribution. Default is \code{Normal}
#' @param cens.type \code{left} for left censoring, \code{right} for right censoring and \code{interval} for intervalar censoring. Default is \code{left}
#' @param LI Vector censoring lower limit indicator of length \code{n}. For each observation: \code{0} if non-censored, \code{-inf} if censored. It is only indicated for when \code{cens.type} is \code{both}. Default is \code{NULL}
#' @param LS Vector censoring upper limit indicator of length \code{n}. For each observation: \code{0} if non-censored, \code{inf} if censored.It is only indicated for when \code{cens.type} is \code{both}. Default is \code{NULL}
#' @param MaxIter The maximum number of iterations of the EM algorithm. Default is \code{200}
#' @param error The convergence maximum error. Default is \code{0.0001}
#' @param Prev Indicator of the prediction process. Available at the moment only for the \code{typeModel=normal} case. Default is \code{FALSE}
#' @param isubj Vector indicator of subject included in the prediction process. Default is \code{NULL}
#' @param step Number of steps for prediction. Default is \code{NULL}
#' @param xpre Design matrix of the fixed effects to be predicted. Default is \code{NULL}.
#' @param zpre Design matrix of the random effects to be predicted. Default is \code{NULL}.
#' @return returns list of class \dQuote{ARpMMEC}:
#' \item{FixEffect}{Data frame with: estimate, standar errors and confidence intervals of the fixed effects.}
#' \item{Sigma2}{Data frame with: estimate, standar errors and confidence intervals of the variance of the white noise process.}
#' \item{Phi}{Data frame with: estimate, standar errors and confidence intervals of the autoregressive parameters.}
#' \item{RandEffect}{Data frame with: estimate, standar errors and confidence intervals of the random effects.}
#' \item{nu}{the parameter "nu" for the t-student distribution}
#' \item{Est}{Vector of parameters estimate (fixed Effects, sigma2, phi, random effects).}
#' \item{SE}{Vector of the standard errors of (fixed Effects, sigma2, phi, random effects).}
#' \item{Residual}{Vector of the marginal residuals.}
#' \item{loglik}{Log-likelihood value.}
#' \item{AIC}{Akaike information criterion.}
#' \item{BIC}{Bayesian information criterion.}
#' \item{AICc}{Corrected Akaike information criterion.}
#' \item{iter}{Number of iterations until convergence.}
#' \item{Yfit}{Vector "y" fitted}
#' \item{MI}{Information matrix}
#' \item{Prev}{Predicted values (if xpre and zpre is not \code{NULL}).}
#' \item{time}{Processing time.}
#' \item{others}{The first and second moments of the random effect and vector Y}
#' @references Olivari, R. C., Garay, A. M., Lachos, V. H., & Matos, L. A. (2021). Mixed-effects
#' models for censored data with autoregressive errors. Journal of Biopharmaceutical Statistics, 31(3), 273-294.
#' \doi{10.1080/10543406.2020.1852246}
#' @examples
#' \dontrun{
#'p.cens = 0.1
#'m = 10
#'D = matrix(c(0.049,0.001,0.001,0.002),2,2)
#'sigma2 = 0.30
#'phi = 0.6
#'beta = c(1,2,1)
#'nj=rep(4,10)
#'tt=rep(1:4,length(nj))
#'x<-matrix(runif(sum(nj)*length(beta),-1,1),sum(nj),length(beta))
#'z<-matrix(runif(sum(nj)*dim(D)[1],-1,1),sum(nj),dim(D)[1])
#'data=ARpMMEC.sim(m,x,z,tt,nj,beta,sigma2,D,phi,struc="ARp",typeModel="Normal",p.cens=p.cens)
#'
#'teste1=ARpMMEC.est(data$y_cc,x,z,tt,data$cc,nj,struc="ARp",order=1,typeModel="Normal",MaxIter = 2)
#'teste2=ARpMMEC.est(data$y_cc,x,z,tt,data$cc,nj,struc="ARp",order=1,typeModel="Student",MaxIter = 2)
#'
#'xx=matrix(runif(6*length(beta),-1,1),6,length(beta))
#'zz=matrix(runif(6*dim(D)[1],-1,1),6,dim(D)[1])
#'isubj=c(1,4,5)
#'teste3=ARpMMEC.est(data$y_cc,x,z,tt,data$cc,nj,struc="ARp",order=1,typeModel="Normal",
#' MaxIter = 2,Prev=TRUE,step=2,isubj=isubj,xpre=xx,zpre=zz)
#'teste3$Prev
#'
#' }
#'
#'
#' @export
#'
#'
#'
ARpMMEC.est=function(y,x,z,tt,cc,nj,struc="UNC",order=1, initial=NULL,nu.fixed=TRUE,
typeModel="Normal",cens.type="left", LI=NULL,LS=NULL, MaxIter=200,
error=0.0001, Prev=FALSE,step=NULL,isubj=NULL,xpre=NULL,zpre=NULL)
{
m<-length(y); N<-sum(nj); p<-dim(x)[2]; q1<-dim(z)[2]; m1<-m*p; m2<-m*q1
if(is.matrix(y)) y <- y[as.vector(!is.na(as.vector(t(y))))]
if(is.matrix(cc)) cc <- cc[as.vector(!is.na(as.vector(t(cc))))]
if (!is.matrix(x)) x=as.matrix(x)
if (!is.matrix(z)) x=as.matrix(z)
if( is.matrix(nj)) nj <- nj[as.vector(!is.na(as.vector(t(nj))))]
if(!is.numeric(y)) stop("y must be a numeric vector. Check documentation!")
if(sum(is.na(y))>0) stop("Vector y does not support NA values.")
if(!is.vector(y)) stop("y must be a vector.Check documentation!")
if(length(y)!=nrow(as.matrix(x))) stop("x does not have the same number of lines than y.")
if(length(y)!=length(cc)) stop("cc does not have the same length than y.")
if(length(y)!=nrow(as.matrix(z))) stop("x does not have the same number of lines than y.")
if(length(y)!=sum(nj)) stop("not compatible sizes between the response y and the repetited measures nj")
if(length(y)==0) stop("The parameter y must be provided.")
if(length(y)!=length(tt)) stop("not compatible sizes between the response y and the vector time tt")
if(!is.numeric(x)) stop("x must be a numeric matrix. Check documentation!")
if(sum(is.na(x))>0) stop("There are some NA values in x.")
if(!is.matrix(x)) stop("x must be a matrix. Check documentation!")
if(det(t(x)%*%x)==0) stop("the columns of x must be linearly independent.")
if(length(x)==0) stop("The parameter x must be provided.")
if(!is.numeric(z)) stop("z must be a numeric matrix. Check documentation!")
if(!is.matrix(z)) stop("z must be a matrix. Check documentation!")
if(sum(is.na(z))>0) stop("There are some NA values in z.")
if(length(z)==0) stop("The parameter z must be provided.")
if(!is.numeric(cc)) stop("cc must be a numeric vector. Check documentation!")
if(!is.vector(cc)) stop("cc must be a vector.Check documentation!")
if(sum(is.na(cc))>0) stop("There are some NA values in cc.")
if(sum(cc%in%c(0,1))<length(cc)) stop("The elements of the vector cc must be 0 or 1.")
if(length(cc)==0) stop("The parameter cc must be provided.")
if(!is.numeric(nj)) stop("nj must be a numeric vector. Check documentation!")
if(!is.vector(nj)) stop("nj must be a vector. Check documentation!")
if(sum(is.na(nj))>0) stop("There are some NA values in nj")
if(length(nj)==0) stop("The parameter nj must be provided.")
if(struc!="DEC"&struc!="DEC(AR)"&struc!="SYM"&struc!="ARp"&struc!="UNC") stop("Struc must be UNC, DEC, DEC(AR), SYM or ARp. Check documentation!")
if(struc=="ARp"){
if(!is.numeric(order) ) stop("Orde must be a number. Check documentation!")
if(length(order)!=1) stop("Order must be a value.")
if(is.numeric(order))
{ if(order!=round(order)|order<=0) stop("Order must be a positive integer value.")}}
if(!is.null(initial))
{ if(!is.null(initial$betas))
{if(!is.numeric(initial$betas)) stop("betas must be a numeric vector. Check documentation!")
if(!is.vector(initial$betas)) stop("betas must be a vector. Check documentation!")
if(length(initial$betas)!=ncol(x)) stop("not compatible sizes between the matrix x and parameter betas.")}
if(!is.null(initial$sigma2))
{if(!is.numeric(initial$sigma2)) stop("sigma2 must be a scalar.")
if(length(initial$sigma2)>1) stop("sigma2 must be a scalar.")}
if(!is.null(initial$alphas))
{if(!is.matrix(initial$alphas)) stop("alphas must be a matrix.")
if(initial$alphas[upper.tri(initial$alphas)]!=initial$alphas[lower.tri(initial$alphas)])stop("alphas must be a simetric matrix.")
if(dim(initial$alphas)[2]!=ncol(z)) stop("not compatible sizes between the matrix z and parameter alphas.")}
if(struc=="ARp"){
if(!is.null(initial$phi))
{if(!is.numeric(initial$phi)) stop("phi must be a numeric vector. Check documentation!")
if(length(initial$phi)!=order) stop("not compatible sizes between the value Arp and parameter phi. Check documentation!")}
}
}
if(typeModel!='Normal'& typeModel!='Student') stop('typeModel must be Normal or Student. Check documentation!')
if(cens.type!="left" & cens.type!="right" & cens.type!="interval")stop('cens.type must be left, right or interval. Check documentation!')
if(cens.type=="interval"&is.null(LI)) stop("The parameter LI must be provided.. Check documentation!")
if(cens.type=="interval"&is.null(LS)) stop("The parameter LS must be provided.. Check documentation!")
if(!is.null(LI)&!is.numeric(LI)) stop("LI must be a numeric vector. Check documentation!")
if(!is.null(LS)&!is.numeric(LS)) stop("LS must be a numeric vector. Check documentation!")
if(length(LS)!=length(LI)) stop("not compatible sizes between the vectors LI and LS. Check documentation!")
if(cens.type=="interval") {
if(length(y)!=length(LI)) stop("not compatible sizes between the vectors y and LI. Check documentation!")
if(length(y)!=length(LS)) stop("not compatible sizes between the vectors y and LS. Check documentation!")
}
if (!is.numeric(MaxIter)) stop("MaxIter must be a positive number. Check documentation!")
if (length(MaxIter) > 1) stop("MaxIter parameter must be a scalar")
if (MaxIter <0) stop("MaxIter parameter must be positive number")
if (!is.numeric(error)) stop("error must be a positive number. Check documentation!")
if (length(error) > 1) stop("error parameter must be a scalar")
if (error <0) stop("error parameter must be positive number")
if (Prev) {
if(is.null(step)|is.null(xpre)|is.null(zpre)|is.null(isubj)) stop("step, isubj, xpre, zpre needs to be provided. Check documentation!")
if (!is.numeric(isubj)) stop("isubj must be a numeric vector. Check documentation!")
if (!is.numeric(step)) stop("step must be a positive number. Check documentation!")
if (step <0) stop("step parameter must be positive number")
if (length(step) > 1) stop("step parameter must be a scalar")
if (ncol(xpre)!=ncol(as.matrix(x))) stop("xpre must have the same number of columns than x")
if (sum(is.na(xpre))>0) stop("There are some NA values in xpre")
if (!is.numeric(xpre)) stop("xpred must be a numeric matrix")
if (ncol(zpre)!=ncol(as.matrix(z))) stop("zpre must have the same number of columns than z")
if (sum(is.na(zpre))>0) stop("There are some NA values in zpre")
if (!is.numeric(zpre)) stop("zpred must be a numeric matrix")
if(nrow(xpre)!=length(isubj)*step) stop("not compatible sizes between xpre and isubj. Check documentation!")
if(nrow(zpre)!=length(isubj)*step) stop("not compatible sizes between zpre and isubj. Check documentation!")
}
if(typeModel=="Normal"){
if(struc=="ARp"){
out<-EMCensArpN(cc=cc,y=y,x=x,z=z,tt=tt,nj=nj, Arp=order, initial=initial, cens.type=cens.type, LI=LI,LS=LS,MaxIter=MaxIter,ee=error,
Prev=Prev,step=step,isubj=isubj ,xpre=xpre,zpre=zpre)}
if(struc=="UNC"){
out<-EMCensArpN(cc=cc,y=y,x=x,z=z,tt=tt,nj=nj, Arp=struc, initial=initial, cens.type=cens.type, LI=LI,LS=LS,MaxIter=MaxIter,ee=error,
Prev=Prev,step=step,isubj=isubj ,xpre=xpre,zpre=zpre)}
if(struc=="DEC"|struc=="DEC(AR)"|struc=="SYM"){
out<-EMCensDECN(cc=cc,y=y,x=x,z=z,tt=tt,nj=nj, struc=struc, initial=initial, cens.type=cens.type, LI=LI,LS=LS,MaxIter=MaxIter,ee=error,
Prev=Prev,step=step,isubj=isubj ,xpre=xpre,zpre=zpre)}
}
if(typeModel=="Student"){
if(struc=="ARp"){
out<-EMCensArpT(cc=cc,y=y,x=x,z=z,ttc=tt,nj=nj,Arp=order,initial=initial,cens.type=cens.type,LL=LI,LU=LS,nu.fixed=nu.fixed,
iter.max=MaxIter,precision=error)
}
if(struc=="UNC"){
out<-EMCensArpT(cc=cc,y=y,x=x,z=z,ttc=tt,nj=nj,Arp=struc,initial=initial,cens.type=cens.type,LL=LI,LU=LS,nu.fixed=nu.fixed,
iter.max=MaxIter,precision=error)
}
if(struc=="DEC"|struc=="DEC(AR)"|struc=="SYM"){
out<-EMCensDECT(cc=cc,y=y,x=x,z=z,ttc=tt,nj=nj,struc=struc,initial=initial,cens.type=cens.type,LL=LI,LU=LS,nu.fixed=nu.fixed,
iter.max=MaxIter,precision=error)
}
}
if(struc=="ARp")
{
cat('\n')
cat('---------------------------------------------------\n')
cat('Autoregressive censored mixed-effects models \n')
cat('---------------------------------------------------\n')
cat('\n')
cat("Autoregressive order =",order)
cat('\n')
cat("Distribution =",typeModel)
cat('\n')
if(typeModel=="Student") cat("nu =",out$nu); cat('\n')
cat("Subjects =",length(nj),";",'Observations =',sum(nj))
cat('\n')
cat('\n')
cat('-----------\n')
cat('Estimates\n')
cat('-----------\n')
cat('\n')
cat('- Fixed effects \n')
cat('\n')
print(out$tableB)
cat('\n')
cat('\n')
cat('- Sigma^2 \n')
cat('\n')
print(out$tableS)
cat('\n')
cat('\n')
cat('- Autoregressives parameters\n')
cat('\n')
print(out$tableP)
cat('\n')
cat('\n')
cat('- Random effects \n')
cat('\n')
print(out$tableA)
cat('\n')
cat('\n')
cat('------------------------\n')
cat('Model selection criteria\n')
cat('------------------------\n')
cat('\n')
critFin <- c(out$loglik, out$AIC, out$BIC)
critFin <- round(t(as.matrix(critFin)),digits=3)
dimnames(critFin) <- list(c("Value"),c("Loglik", "AIC", "BIC"))
print(critFin)
cat('\n')
cat('-------\n')
cat('Details\n')
cat('-------\n')
cat('\n')
cat("Convergence reached? =",(out$iter < MaxIter))
cat('\n')
cat('Iterations =',out$iter,"/",MaxIter)
cat('\n')
cat("Processing time =",out$time,units(out$time))
cat('\n')
}
if(struc!="ARp")
{
cat('\n')
cat('---------------------------------------------------\n')
cat('DEC censored mixed-effects models \n')
cat('---------------------------------------------------\n')
cat('\n')
cat("Case =",struc)
cat('\n')
cat("Distribution =",typeModel)
cat('\n')
if(typeModel=="Student") cat("nu =",out$nu); cat('\n')
cat("Subjects =",length(nj),";",'Observations =',sum(nj))
cat('\n')
cat('\n')
cat('-----------\n')
cat('Estimates\n')
cat('-----------\n')
cat('\n')
cat('- Fixed effects \n')
cat('\n')
print(out$tableB)
cat('\n')
cat('\n')
cat('- Sigma^2 \n')
cat('\n')
print(out$tableS)
cat('\n')
cat('\n')
if(struc!="UNC"){
cat('- Autoregressives parameters\n')
cat('\n')
print(out$tableP)
cat('\n')
cat('\n')}
cat('- Random effects \n')
cat('\n')
print(out$tableA)
cat('\n')
cat('\n')
cat('------------------------\n')
cat('Model selection criteria\n')
cat('------------------------\n')
cat('\n')
critFin <- c(out$loglik, out$AIC, out$BIC)
critFin <- round(t(as.matrix(critFin)),digits=3)
dimnames(critFin) <- list(c("Value"),c("Loglik", "AIC", "BIC"))
print(critFin)
cat('\n')
cat('-------\n')
cat('Details\n')
cat('-------\n')
cat('\n')
cat("Convergence reached? =",(out$iter < MaxIter))
cat('\n')
cat('Iterations =',out$iter,"/",MaxIter)
cat('\n')
cat("Processing time =",out$time,units(out$time))
cat('\n')
}
if(typeModel=="Student") nu<-out$nu
if(typeModel=="Normal") nu<-NULL
obj.out <- list(FixEffect=out$tableB, Sigma2=out$tableS, Phi=out$tableP,RandEffect=out$tableA, nu=nu,
Est=c(out$beta1, sigma2=out$sigmae, phi=out$phi, RnEffect=out$dd), SE=out$SE,Residual=out$residual,
loglik=out$loglik, AIC=out$AIC, BIC=out$BIC, AICc=out$AICcorr, iter=out$iter, Yfit=out$yfit,
MI=out$MI, Prev=out$Prev, time=out$time, others=list(ubi = out$ubi, ubbi = out$ubbi, uybi = out$uybi, uyi = out$uyi, uyyi = out$uyyi,varbeta=out$varbeta,yog=out$yorg) )
class(obj.out) = "ARpMMEC"
return(obj.out)
}
Try the ARpLMEC package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.