Nothing
R/isnimgm.R
In isni: Index of Local Sensitivity to Nonignorability
Defines functions
isnimgm
Documented in
isnimgm
## Updates
## Nov 16, 2018
## Add codes to permit ISNI calculation when there are subjects with all outcomes are missing.
## The added codes are in isnimgm(), fun.mgmsubi(),
#' A data set for Psychiatric Drug Treatment
#'
#' The variables are as follows:
#' \itemize{
#' \item y. number of days
#' \item time. in weeks
#' \item sub. patients id
#' \item group. treatment group with placebo (PBO) and desipramine (DMI)
#' \item basey. baseline days
#' \item g1. dummy variables for treatment group
#' \item g2. dummy variables for treatment group
#' \item time1. dummy variables for time
#' \item time2. dummy variables for time
#' \item grptime. group by time
#' }
#'
#' @docType data
#' @keywords datasets
#' @name coc
#' @usage data(coc)
#' @format A data frame with 869 rows and 10 variables
NULL
#' A data set for Quality of Life Emontional Functioning outcome.
#'
#' The variables are as follows:
#' \itemize{
#' \item id: patients id
#' \item y: EF score
#' \item time: time in months since randomization.
#' \item group: placebo (0) or flutamide (1)
#' \item perf: baseline performance score
#' \item sever: baseline disease severity
## \item t1: indicator variable for month 1
## \item t3: indicator variable for month 3
## \item t6: indicator variable for month 6
# \item yp: most recently observed prior outcome
#' \item basey: EF at baseline
#' \item g: missingness status ("O"=observed, "D"=dropout, "I"=intermittent missingness)
#' \item gp: missingness status in the prior visit ("O"=observed, "D"=dropout, "I"=intermittent missingness, "U"=undefined)
#' }
#'
#' @docType data
#' @keywords datasets
#' @name qolef
#' @usage data(qolef)
#' @format A data frame with 2860 rows and 10 variables
NULL
#' Function for ISNI computation when the outcome follows marginal multivariate Gaussian Models.
#'
#' Calculate the ISNI when the regression outcome is subject to missingness and follows marginal multivaraite Gaussian models.
#' @param formula an object of class "Formula": a symbolic description of the models to be fitted for the outcome and missingness status variable.
#' The details of model specification are given under "Details".
#' @param data the name of data frame containing all the variables in the model and all the observations including those intended to be collected
#' but became missing.
## @param gmodel an object of class "formula": a symbolic description of the selection model to be fitted for the missingness indictor g
#' @param cortype the type of within-subject correlation structure.
#' @param id the name of variable for the level-2 clustering variable.
#' @param subset an optional vector specifying a subset of observations to be used in the fitting process for the outcome model and the missingness mechanism model.
#' @param weights frequency weights to be assigned to each \code{id}. when supplied, indicates differential weights are used; otherwise each \code{id} is weighted equally.
#' @param predprobobs Null if using buil-in multinomial transitional logistic model to obtain predicted probabilities of being observed;
#' otherwise user supply the name of the variable in \code{data} that gives these probabilities for all the observations in the data.
#' @param misni FALSE if using the default approach to computing ISNI with a scalar nonignorability parameter;
#' TRUE when computing ISNI with multiple nonignorability parameters.
#' @details The ISNI analysis is based on a joint selection model and requires specifying two model equations: the complete-data model and the missing data mechanism model.
#' To specify the variables in the models that are required for computing the isni measures, we make use of the \code{R} package "Formula" designed for handling model
#' equations with multiple responses and multiple sets of predictors . At a minimum, the user should supply a single-equation
#' in the typical form: \code{response} ~ \code{Xterms} where \code{response} is the (numeric or factor) vector for the outcome of interest and \code{Xterms}
#' is a series of terms, separated by + operators, which specify a linear predictor for response. With the signle-equation specification, the \code{isniglm} function
#' will by default use the utility function \code{definemissingstatus} provided in the package to generate the
#' missingness status variables at the current and prior visits and then use \code{Xterms} as the observed missingness predictors for fitting a first-order transitional missing data model.
#' It is important to sort within-\code{id} observations by time so that the missingness status variables can be defined correctly in this default setting. The \code{isnimgm} then computes the MAR estimates and conducts ISNI computation
#' to evaluate the rate of change of model estimates in the neighborhood of the MAR model where the missingness probability is allowed to depend on the unobserved value
#' of \code{response}, even after conditioning on the other missingness predictors.
#'
#' The above single-equation formula specification uses the same set of predictors for the outcome model and the missing data mechanism model for computing ISNI.
#' To use different sets of predictors, one can explicitly specifiy a two-equation formula as: \code{response} | \code{miss + missprior} ~ \code{Xterms} | \code{Sterms},
#' which specifies the formula for the outcome model as \code{response} ~ \code{Xterms} and that for the missing data mechanism model as \code{miss + missprior} ~ \code{Sterms},
#' where \code{Xterms} and the observed missingness predictors \code{Sterms} can be different, \code{miss} and \code{missprior} are the variable names in \code{data} denoting the missingness status at the
#' current and prior visits, respectively.
#' @references Ma G, Troxel AB, Heitjan DF. An Index of Local Sensitivity to Nonignorable Dropout in Longitudinal Modeling. Stat Med. 2005;24:2129-2150. \cr
#' Xie H. Analyzing Longitudinal Clinical Trial Data with Nonignorable Missingness and Unknown Missingness Reasons. Comput Stat Data Anal. 2012;56:1287-1300. \cr
#' Xie H, Qian Y. Measuring the Impact of Nonignorability in Panel Data with Non-Monotone Nonresponse. Journal of Applied Econometrics. 2012;27:129-159.
#' @name isnimgm
#' @aliases isnimgm
#' @import matrixcalc mvtnorm nlme nnet stats Formula
#' @importFrom nlme corCompSymm
#' @export isnimgm
#' @examples
#' models= y | g+gp ~ perf + sever+ as.factor(time) + group +as.factor(time):group |
#' as.factor(time) * group + yp+ perf + sever
## gmodel= as.factor(g)~ t3+t6+group+yp+perf+sever
## qolef$t12<-qolef$t1*qolef$group
## qolef$t32<-qolef$t3*qolef$group
## qolef$t62<-qolef$t6*qolef$group
##
#' qolef.isni=isnimgm(models, data=qolef, id=id)
#' summary(qolef.isni)
#'
## ymodel= y | g+gp ~ t1+t3+t6 + group + t12 + t32+t62 + perf + sever | t3+t6+group+yp+perf+sever
## summary(isnimgm(ymodel, data=qolef, id=id))
#'
isnimgm = function(formula, data, cortype="CS", id, subset, weights, predprobobs, misni=FALSE) {
## (1) process the call and set up model frame.
cl <- match.call()
if (missing(data)) data <- environment(formula)
if (missing(id)) stop("Users need to supply the name of the level-2 clusting variable in the data via the argument id.")
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "weights", "subset", "id","predprobobs"), names(mf), 0L)
mf <- mf[c(1L, m)]
f <- Formula(formula)
if (any(length(f)>2)) stop("Cannot have more than two models")
mf$formula <- f
options(na.action=na.pass)
mf[[1L]] <- as.name("get_all_vars")
mfvar<- eval(mf, parent.frame())
if (!missing(subset)) {
SubSet <-mfvar$subset
mfvar <-mfvar[SubSet,]
}
mf[[1L]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
ID= model.extract(mf, "id")
if (missing(weights)) WTs <- NULL else
WTs <- model.extract(mf, "weights")
if (!is.null(WTs) && !is.numeric(WTs))
stop("'weights' must be a numeric vector")
if (!is.null(WTs) && any(WTs < 0))
stop("negative weights not allowed")
if (is.null(WTs)){
IWT <- 0
WTs <- rep(1, nrow(mf))
}
isni_WTs <- isni_id <- isni_outcome_ <- NULL
mf= cbind(mf, mfvar[!(names(mfvar) %in% names(mf))], isni_id=ID, isni_WTs=WTs)
if (! missing(predprobobs)) mf$fitprobs_= model.extract(mf, "predprobobs")
##(2) Extract responses and predictors from the model frame
ymodel=formula(f, lhs=1, rhs=1)
x<- as.matrix(model.matrix(f, data=mf, rhs=1))
y<- model.part(f,data=mf, lhs=1)[,1] ## will be a vector
if (nrow(x)<1) print("No predictor variables specified for the outcome")
options(na.action=na.omit)
## check if the missing status indicator g and gp are specified, if not add g and gp in the lhs of formula
if (length(f)[1]==1) {
cat("\n NOTE: The model formula specification does not include the missingness status variables at the current and the prior visits.
The program will assume that data has already been sorted by the time variable within each level-2 ID variable and proceed to create the missingness status variables and then fit a transitional missingness mechanism model. \n")
f <- update(f, .| g_ +gp_ ~.)
mf <- cbind(mf, isni_outcome_=y)
mf <- definemissingstatus(data=mf, id=isni_id, y=isni_outcome_)
if (ncol(model.part(f,data=mf, lhs=2)) != 2) stop("The missing data model should be in the form of g+gp ~ Sterm")
}
gmodel <- formula(f, lhs=2,rhs=length(f)[2])
## ## Drop all observations with missing values in predictor matrix X to avoid NAs in ISNI computation.
missX<- apply(x, 2, FUN=function (u) any(is.na(u)))
if (any(missX)) cat(paste("\n Note: There are missing values in fixed-effect outcome predictors-- ", paste(dimnames(x)[[2]][missX], collapse=" "), ". Observations with missing values in these predictors will be dropped out from computing ISNI,
which may reduce ISNI values if the outcome is missing concurrently. Alternatively one can impute these missing predictors values and re-compute ISNI. ", sep=""))
missX <- apply(x, 1, FUN=function (u) any(is.na(u)))
WTs=WTs[!missX]
mf=mf[!missX,,drop=F]
options(na.action=na.pass)
x<- as.matrix(model.matrix(f, data=mf, rhs=1))
y<- model.part(f,data=mf, lhs=1)[,1]
xomit <- as.matrix(model.matrix(f, data=mf[!is.na(y), ], rhs=1))
if (ncol(x) != ncol(xomit)) {
cat(paste("\n All variable names in fixed-effects design matrix of ISNI analysis including observations with missing outcomes: \n", paste (dimnames(x)[[2]], collapse=" ")))
cat(paste("\n All variable names in fixed-effects design matrix of the MAR model excluding observations with missing outcomes: \n", paste (dimnames(xomit)[[2]], collapse=" ")))
stop("\n The design matrix for the MAR model and ISNI analysis are different as shown above. Please modify your model formula specification, e.g. avoiding using as.factor function
in model formula. ")
}
options(na.action=na.omit)
maxT = max(table(mf$isni_id))
sdy <- sd(y,na.rm=T)
## fit the mgm with the observed data.
if (cortype=="CS") iggls= gls(ymodel, data=mf,correlation = corCompSymm(form = ~ 1 | isni_id), weights=~1/isni_WTs, method="ML", na.action=na.exclude) else
if (cortype=="AR1") iggls= gls(ymodel, data=mf,correlation = corAR1(form = ~ 1 | isni_id), weights=~1/isni_WTs, method="ML",na.action=na.exclude) else
if (cortype=="UN") iggls= gls(ymodel, data=mf,correlation = corSymm(form = ~ 1 | isni_id), weights=~1/isni_WTs, method="ML",na.action=na.exclude) else
stop ("This Type of correlation is not implmented yet for ISNI computation")
## fit a missing data model to obtain fitted probabilities for being observed.
if (missing(predprobobs)) {
mf.mdm= tmdm(gmodel, data=mf)
mf$fitprobs=mf.mdm$obsprob
mf$A10= mf.mdm$A10; mf$A20= mf.mdm$A20; mf$A11= mf.mdm$A11
} else if ( misni==T)
stop("predprobobs needs to to a matrix to obtain ISNI vector for multiple nonignorability parameter")
D=list(sigma=iggls$sigma, rho=intervals(iggls)$corStruct[,"est."])
b<- iggls$coef
nb = length(b)
nD=length(D$sigma)+length(D$rho)
npar= nb + nD
uid = unique(mf$isni_id)
nabla12=matrix(0, nrow=npar, ncol=3)
## compute nabla11 and nabla12
nabla11=matrix(0, nrow=npar, ncol=npar)
for (i in 1:length(uid)) {
xi= as.matrix(x[mf$isni_id==uid[i],,drop=F])
yi= y[mf$isni_id==uid[i]]
## Increase nabla11 only if there is at least one observed outcome
if (!all(is.na(yi))) nabla11=nabla11+ fun.mgmsubi(yi=yi,xi=xi,maxT=maxT, b=b,D=D, cortype=cortype, transform=FALSE, case=1)
gfiti= NULL; Afiti=NULL
if (! missing(predprobobs)) gfiti= mf$fitprobs_[mf$isni_id==uid[i]]
else Afiti= cbind(mf$A10[mf$isni_id==uid[i]], mf$A20[mf$isni_id==uid[i]], mf$A11[mf$isni_id==uid[i]])
if (any(is.na(yi))) nabla12 = nabla12 + fun.mgmsubi(yi=yi,xi=xi,maxT=maxT, b=b,D=D, cortype=cortype, transform=FALSE, gfiti=gfiti,Afiti=Afiti, case=2)
}
invnabla11= solve(nabla11)
invnabla11[1:nb,1:nb]= -iggls$varBeta
isni=-invnabla11%*%nabla12
se=sqrt(diag((-invnabla11)))
bD= c(b,unlist(D))
if (cortype=="CS" | cortype=="AR1" ) {
names(bD)=c(names(b), names(D))
} else
if (cortype=="UN" ) {
names(bD)=c(names(b), 'sigma', row.names(intervals(iggls)$corStruct)) ##names(D[[2]]))
}
isnivec=isni
dimnames(isnivec)<-list(names(bD),c("ISNI_IO", "ISNI_DO", "ISNI_II"))
if (misni==FALSE) {
isni=apply(isni,1, sum); senstran<-abs((sdy*se)/isni)
} else
{
isni=apply(abs(isni),1, sum); senstran<-abs((sdy*se)/isni)
}
res=list(coefficients=bD,se=se,isni=c(isni),c=c(senstran), call=cl, isnivec=isnivec, misni=misni,
logLik=iggls$logLik, aic=summary(iggls)$AIC, bic=summary(iggls)$BIC)
class(res) = c(res$class, "isnimgm")
return(res)
}
#' Function to print out a summary of isnimgm object in a matrix form.
#'
#' @param object the isnimgm object obtained from the isnimgm function
#' @param digits the number of significant digits to use when printing
#' @param ... additional arguements
#' @export summary.isnimgm
#' @export
#' @name summary.isnimgm
#' @aliases summary.isnimgm
summary.isnimgm<-function(object, digits = max(3, getOption("digits") - 2),
...) {
if (class(object) != "isnimgm") stop('Invalid object class')
cat("\nCall:\n", paste(deparse(object$call), sep = "\n",
collapse = "\n"), "\n\n", sep = "")
## Name the columns
if (object$misni==T) isniname<-c('MAR Est.','Std. Err','MISNI','c') else
isniname<-c('MAR Est.','Std. Err','ISNI','c')
## Set up matrix to hold result
res<-matrix(0,length(object$coef),length(isniname))
dimnames(res)<-list(names(object$coef),isniname)
for (i in 1:length(object$coef))
res[i,]<- c(object$coef[i],object$se[i],object$isni[i],object$c[i])
printCoefmat(res, digits = digits, cs.ind = 1:2)
}
#' Function to print the isnimgm object.
#'
#' \code{print} method for class isnimgm
#' @param x the isnimgm object obtained from the isnimgm function
#' @param digits the number of significant digits to use when printing
#' @param ... further arguments passed to or from other methods.
#'
#' @return The function print.isnimgm prints the model call, isni and c statistics from the isnimgm object.
#' @name print.isnimgm
#' @aliases print.isnimgm
#' @export print.isnimgm
#' @export
print.isnimgm<-function(x, digits = max(3, getOption("digits") - 2), ...) {
if (class(x) != "isnimgm") stop('Invalid object class')
cat("\nCall:\n", paste(deparse(x$call), sep = "\n",
collapse = "\n"), "\n\n", sep = "")
if (length(x$coef)>0) {
cat("ISNIs:\n")
print.default(format(x$isni, digits = digits), print.gap = 2L,
quote = FALSE)
cat("\n")
cat("c statistics:\n")
print.default(format(x$c, digits = digits), print.gap = 2L,
quote = FALSE)
}
else cat("No coefficients\n")
cat("\nlogLik of the MAR model: ", paste(format(x$logLik,digits=digits), sep = "\n",
collapse = "\n"), "\n", sep = "")
cat("\nAIC of the MAR model: ", paste(format(x$aic, digits=digits), sep = "\n",
collapse = "\n"), "\n", sep = "")
cat("\nBIC of the MAR model: ", paste(format(x$bic, digits=digits), sep = "\n",
collapse = "\n"), "\n", sep = "")
cat("\n")
invisible(x)
}
## Internal Function for various subject-level computation required for ISNIGLS.
##
## Calculate subject-level quantities when the regression outcome is subject to missingness and follows marginal multivaraite Gaussian Models.
## @param yi vector of the response for the ith subject
## @param xi matrix of the covariates for the ith subject
## @param maxT maximum number of visits
## @param b the mean parameter vector beta
## @param D the vector of unique parameters in the variance-covariance matrix for the error term in the GLS model for Y
## @param correlation the form of within-subject correlation structure in the GLS model for Y
## @param transform logical indicating wether or not the parameter in D is transformed.
## @param gfiti vector of predicted probabilities of being observed for all the observations from the ith subject
## @param Afiti matrix of 3 columns of predicted transitional probabilities for the missing observations from the ith subject.
## @param case 1: calculated nabla11_i; 2: calculate nabla12_i
## @aliases fun.glssubi
## @export
fun.mgmsubi <- function ( yi, xi, maxT=maxT, b, D, cortype,transform=FALSE, gfiti=NULL, Afiti=NULL, case=1) {
## output
## res -- if case=1, res=hessian matrix
## if case=2, res=\nabla_12 for subject i.
## if case=3, res=negative loglikelihood, retired and unused in the package.
## if case=4, res=negative score value, retired and unused in the pacakge
## if case=5, res=der.meanymi, retired and unused in the package.
obsi = !is.na(yi)
xiobs = as.matrix(xi[obsi,,drop=F])
yiobs = as.vector(yi[obsi])
niobs = length(yiobs)
muiobs = as.vector(xiobs %*% as.matrix(b))
riobs= yiobs - muiobs
nb = length(b)
nD=length(D$sigma)+length(D$rho)
npar= nb + nD
sigma= D$sigma; rho=D$rho
## obtain the variance-covariance matrix
if (cortype=="CS") {
covyi= fun.csmat(sigma=sigma, rho=rho, obsi=length(yi),transform=transform, case=1)
}
if (cortype=="AR1") {
covyi= fun.ar1mat(sigma=sigma, rho=rho, obsi=length(yi),transform=transform, case=1)
}
if (cortype=="UN") {
col=matrix(0, nrow=maxT, ncol=maxT)
for (i in 1:(maxT-1)){
col[(i+1):maxT,i]=rho[((i-1)*maxT-i*(i-1)/2+1) : (i*maxT-i*(i+1)/2)]
}
rho2=col+t(col)+diag(maxT)
covyi= fun.generalmat(sigma=sigma, rho=rho2, obsi=length(yi),transform=transform, case=1)
}
covyiobs= covyi[obsi,obsi]
## only invert if there is yiobs is not empty
if (niobs>0) icovyiobs= solve(covyiobs)
##print(icovyiobs)
## case=1, compute the nabla11 for subject i
if (case==1) {
nabla11= matrix(0, npar, npar)
## return 0 if there is no observed outcome.
if (niobs==0) return(nabla11)
nabla11[1:nb, 1:nb]=-t(xiobs)%*%icovyiobs%*%xiobs
if (cortype=="CS") {
covid1 = fun.csmat(sigma=sigma,rho=rho,obsi=niobs,transform=transform,case=2)
covid2 = fun.csmat(sigma=sigma,rho=rho,obsi=niobs,transform=transform,case=3)
Ai= covid1
derAi= matrix(NA, nrow=niobs, ncol=niobs)
for (j in 1:length(Ai)) Ai[[j]]= icovyiobs%*%covid1[[j]]%*%icovyiobs
for (i in 1:nD) {
for (j in 1:nD){
derAi= -icovyiobs%*%(covid1[[i]]%*%icovyiobs%*%covid1[[j]]-covid2[i,j,,]+covid1[[j]]%*%icovyiobs%*%covid1[[i]])%*%icovyiobs
nabla11[(nb+i),(nb+j)]= -0.5*sum(diag(-Ai[[i]]%*%covid1[[j]]+icovyiobs%*%covid2[i,j,,]))+0.5*t(as.matrix(riobs))%*%derAi%*%riobs
}
}
for (j in 1:nD) {nabla11[1:nb,(nb+j)]= -0.5*t(xiobs)%*%(Ai[[j]]+t(Ai[[j]]))%*%riobs}
nabla11[(nb+1):npar, 1:nb]=t(nabla11[1:nb, (nb+1):(npar)])
}
if (cortype=="AR1") {
covid1 = fun.ar1mat(sigma=sigma,rho=rho,obsi=niobs,transform=transform,case=2)
covid2 = fun.ar1mat(sigma=sigma,rho=rho,obsi=niobs,transform=transform,case=3)
Ai= covid1
derAi= matrix(NA, nrow=niobs, ncol=niobs)
for (j in 1:length(Ai)) Ai[[j]]= icovyiobs%*%covid1[[j]]%*%icovyiobs
for (i in 1:nD) {
for (j in 1:nD){
derAi= -icovyiobs%*%(covid1[[i]]%*%icovyiobs%*%covid1[[j]]-covid2[i,j,,]+covid1[[j]]%*%icovyiobs%*%covid1[[i]])%*%icovyiobs
nabla11[(nb+i),(nb+j)]= -0.5*sum(diag(-Ai[[i]]%*%covid1[[j]]+icovyiobs%*%covid2[i,j,,]))+0.5*t(as.matrix(riobs))%*%derAi%*%riobs
}
}
for (j in 1:nD) nabla11[1:nb,(nb+j)]= -0.5*t(xiobs)%*%(Ai[[j]]+t(Ai[[j]]))%*%riobs
nabla11[(nb+1):npar, 1:nb] = t(nabla11[1:nb, (nb+1):(npar)])
}
if (cortype=="UN") {
covid1 = fun.generalmat(sigma=sigma,rho=rho2,obsi=obsi,transform=transform,case=2, case2=1)
covid2 = fun.generalmat(sigma=sigma,rho=rho2,obsi=obsi,transform=transform,case=3)
Ai= covid1
derAi= matrix(NA, nrow=niobs, ncol=niobs)
for (j in 1:length(Ai)) Ai[[j]]= icovyiobs%*%covid1[[j]]%*%icovyiobs
for (i in 1:nD) {
for (j in 1:nD){
derAi= -icovyiobs%*%(covid1[[i]]%*%icovyiobs%*%covid1[[j]]-covid2[i,j,,]+covid1[[j]]%*%icovyiobs%*%covid1[[i]])%*%icovyiobs
nabla11[(nb+i),(nb+j)]= -0.5*sum(diag(-Ai[[i]]%*%covid1[[j]]+icovyiobs%*%covid2[i,j,,]))+0.5*t(as.matrix(riobs))%*%derAi%*%riobs
}
}
for (j in 1:nD) nabla11[1:nb,(nb+j)]= -0.5*t(xiobs)%*%(Ai[[j]]+t(Ai[[j]]))%*%riobs
nabla11[(nb+1):npar, 1:nb]=t(nabla11[1:nb, (nb+1):(npar)])
}
return(nabla11)
}## end if case=1
## case=2, calculate the \nabla_12 for subject i.
if (case==2) {
if (any(is.na(yi))){
misi = is.na(yi)
ximis = as.matrix(xi[misi,,drop=F])
yimis = as.vector(yi[misi])
##fitp0i = gfiti[misi]
Afit0i=matrix(0, nrow=sum(misi==T), ncol=3)
if (is.null(gfiti)) Afit0i =Afiti[misi, ,drop=F]
else Afit0i[,1]= gfiti[misi]
cit = covyi[misi, obsi,drop=F]
dmean.D= matrix(0, nrow=nD, ncol=ncol(Afit0i))
if (niobs==0) {dmean.beta=t(ximis)%*% Afit0i; return(rbind(dmean.beta, dmean.D)) }
dmean.beta = t(ximis - cit %*% icovyiobs %*% xiobs) %*% Afit0i ## as.matrix(fitp0i)
##numeric(nD)
if (cortype=="CS") {
derD= fun.csmat(sigma=sigma,rho=rho,obsi=length(yi),transform=transform,case=2)
} else
if (cortype=="AR1") {
derD= fun.ar1mat(sigma=sigma,rho=rho,obsi=length(yi),transform=transform,case=2)
}
if (cortype=="UN") {
derD= fun.generalmat(sigma=sigma,rho=rho2,obsi=length(yi),transform=transform,case=2, case2=2)
}
for (k in 1:nD) {
Aij<-icovyiobs%*%derD[[k]][obsi,obsi]%*%icovyiobs
dercitDk<- derD[[k]][misi,obsi,drop=F]
dmean.D[k,]<-t(dercitDk%*%icovyiobs%*%as.matrix(riobs) -cit%*%Aij%*%as.matrix(riobs)) %*% Afit0i ##as.matrix(fitp0i)
}
## dercov.sigma<-fun.csmat(sigma=sigma,rho=rho,nobs=length(yi),transform=transform,case=2)
## Aij<-icovyiobs%*%dercov.sigma[obsi,obsi]%*%icovyiobs
## dercit.sigma<- dercov.sigma[misi,obsi,drop=F]
## dmean.sigma<-t(dercit.sigma%*%icovyiobs%*%as.matrix(riobs) -cit%*%Aij%*%as.matrix(riobs)) %*% Afit0i ##as.matrix(fitp0i)
## dercov.rho<-fun.csmat(sigma=sigma,rho=rho,nobs=length(yi),transform=transform,case=3)
## Aij<-icovyiobs%*%dercov.rho[obsi,obsi]%*%icovyiobs
## dercit.rho<- dercov.rho[misi,obsi,drop=F]
## dmean.rho<-t(dercit.rho%*%icovyiobs%*%as.matrix(riobs) -cit%*%Aij%*%as.matrix(riobs))%*% Afit0i ##as.matrix(fitp0i)
res = rbind(dmean.beta, dmean.D) ##c(dmean.beta, dmean.D)
return(res)
} else
res = numeric(dim(xiobs)[2]+2)
} ## end if case=2
res
}
## Internal Function for generate compound symmetry matrix and its first and second derivatives wrt D
##
## Calculate the compound symmetry matrix and its first and second derivatives wrt D.
## @param sigma standard deviation
## @param rho correlation coefficient
## @param obsi dimension of the square matrix
## @param transform logical indicating wether or not the parameter in D is transformed.
## @param case 1: return nobs x nobs CS matrix; 2: return a list for first derivatives of the compound symmetry matrix wrt sigma and rho;
## 3: return an array containing 2nd derivatives of the compound symmetry matrix wrt sigma and rho.
## @aliases fun.csmat
## @export
fun.csmat<-function(sigma,rho,obsi,transform,case=1) {
## calculate the the related quantity for compound symmetry vcov matrix.
## input
## sigma -- the standard deviation
## rho --- the correlation coefficient
## nobs --- the number of observations, nobs>=1
## transform --- logical number indicator if the parameters are transformed
## output
## case=1 --- return nobs x nobs CS matrix
## case=2 --- return a list for first derivatives of the compound symmetry matrix wrt sigma and rho
## case=3 --- return an array containing 2nd derivatives of the compound symmetry matrix wrt sigma and rho
if (transform==FALSE) {
if (case==1) {
nobs = obsi
res<-sigma^2*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs))
}
else if (case==2) {
nobs = obsi
derSigma<- 2*sigma*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs))
derRho<- sigma^2*(matrix(1,nrow=nobs,ncol=nobs)-diag(1,nrow=nobs))
res <- list(derSigma=derSigma, derRho=derRho)
res
}
## derivative wrt to rho
else if (case==3) {
nobs = obsi
res = array(NA, c(2,2, nobs,nobs))
res[1,1,,]=2*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs))
res[1,2,,]=res[2,1,,]=2*sigma*(matrix(1,nrow=nobs,ncol=nobs)- diag(1, nrow=nobs))
res[2,2,,]=matrix(0, nrow=nobs, ncol=nobs)
res
}
} else
{
sigma<-fun.trans2pos(sigma)
rho <- fun.trans2rho(rho)
if (case==1) {
nobs = obsi
res<-sigma^2*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs))
}
else if (case==2) {
nobs = obsi
derSigma<- 2*sigma*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs))
derRho<- sigma^2*(matrix(1,nrow=nobs,ncol=nobs)-diag(1,nrow=nobs))
res <- list(derSigma=derSigma, derRho=derRho)
res
}
## derivative wrt to rho
else if (case==3) {
nobs = obsi
res = array(NA, c(2,2, nobs,nobs))
res[1,1,,]=2*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs))
res[1,2,,]=res[2,1,,]=2*sigma*(matrix(1,nrow=nobs,ncol=nobs)- diag(1, nrow=nobs))
res[2,2,,]=matrix(0, nrow=nobs, ncol=nobs)
res
}
}
res
}
## Internal Function for generate AR1 matrix and its first and second derivatives wrt D
##
## Calculate the autoregressive matrix and its first and second derivatives wrt D.
## @param sigma standard deviation
## @param rho autoregressive parameter
## @param obsi dimision of the square matrix
## @param transform logical indicating wether or not the parameter in D is transformed.
## @param case 1: return nobs x nobs AR1 matrix; 2: return a list for first derivatives of the compound symmetry matrix wrt sigma and rho;
## 3: return an array containing 2nd derivatives of the compound symmetry matrix wrt sigma and rho.
## @aliases fun.ar1mat
## @export
fun.ar1mat<-function(sigma,rho,obsi,transform,case=1) {
## calculate the the related quantity for AR1 (autoregressive process of order 1) vcov matrix.
## input
## sigma -- the standard deviation
## rho --- the autoregressive parameter
## nobs --- the number of observations, nobs>=1
## transform --- logical number indicator if the parameters are transformed
## output
## case=1 --- return nobs x nobs AR1 matrix
## case=2 --- return a list for first derivatives of the AR1 matrix wrt sigma and rho
## case=3 --- return an array containing 2nd derivatives of the AR1 matrix wrt sigma and rho
if (transform==FALSE) {
if (case==1) {nobs = obsi
res<-matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
res[i,j]<-sigma^2*rho^abs(i-j)
}
}
} else
if (case==2) {nobs = obsi
derSigma<-matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
derSigma[i,j] <- 2*sigma*rho^abs(i-j)
}}
derRho<-matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
derRho[i,j] <- sigma^2*abs(i-j)*rho^(abs(i-j)-1)
}}
res <- list(derSigma=derSigma, derRho=derRho)
} else
## derivative wrt to rho
if (case==3) {nobs = obsi
res = array(NA, c(2,2, nobs,nobs))
res11 <- matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
res11[i,j]<-2*rho^abs(i-j)
}}
res[1,1,,]=res11
res12 <- matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs){
for (j in 1:nobs){
res12[i,j]<-2*sigma*abs(i-j)*rho^(abs(i-j)-1)
}}
res[1,2,,]=res[2,1,,]=res12
res22 <- matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs){
for (j in 1:nobs){
if(abs(i-j)<=1) res22[i,j] <- 0 else
res22[i,j]<-sigma^2*abs(i-j)*(abs(i-j)-1)*rho^(abs(i-j)-2)
}}
res[2,2,,]=res22
}
} else
{
sigma<-fun.trans2pos(sigma)
rho <- fun.trans2rho(rho)
if (case==1) {nobs = obsi
res<-matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
res[i,j]<-sigma^2*rho^abs(i-j)
}}
} else
if (case==2) {nobs = obsi
derSigma<-matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
derSigma[i,j] <- 2*sigma*rho^abs(i-j)
}}
derRho<-matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
derRho[i,j] <- sigma^2*abs(i-j)*rho^(abs(i-j)-1)
}}
res <- list(derSigma=derSigma, derRho=derRho)
} else
## derivative wrt to rho
if (case==3) {nobs = obsi
res = array(NA, c(2,2, nobs,nobs))
res11 <- matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
res11[i,j]<-2*rho^abs(i-j)
}}
res[1,1,,]=res11
res12 <- matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs){
for (j in 1:nobs){
res12[i,j]<-2*sigma*abs(i-j)*rho^(abs(i-j)-1)
}}
res[1,2,,]=res[2,1,,]=res12
res22 <- matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs){
for (j in 1:nobs){
if(abs(i-j)<=1) res22[i,j] <- 0 else
res22[i,j]<-abs(i-j)*(abs(i-j)-1)*rho^(abs(i-j)-2)
}}
res[2,2,,]=res22
}
}
res
}
## Internal Function for generate a general symmetric matrix and its first and second derivatives wrt D
## Calculate the general symmetric matrix and its first and second derivatives wrt D.
## @param sigma standard deviation
## @param rho correlation coefficient
## @param obsi dimension of the square matrix
## @param transform logical indicating wether or not the parameter in D is transformed.
## @param case 1: return nobs x nobs general symmetric matrix; 2: return a list for first derivatives of the general symmetric matrix wrt D;
## 3: return an array containing 2nd derivatives of the general symmetric matrix wrt D.
## @param case2 indicator for variance and covariance parameters.
## @aliases fun.generalmat
## @export
fun.generalmat<-function(sigma,rho,obsi,transform,case=1, case2=NULL) {
## calculate the the related quantity for compound symmetry vcov matrix.
## input
## sigma -- the standard deviation
## rho --- the correlation coefficient
## nobs --- the number of observations, nobs>=1
## transform --- logical number indicator if the parameters are transformed
## output
## case=1 --- return nobs x nobs CS matrix
## case=2 --- return a list for first derivatives of the compound symmetry matrix wrt sigma and rho
## case=3 --- return an array containing 2nd derivatives of the compound symmetry matrix wrt sigma and rho
if (transform==FALSE) {
if (case==1) {rhoi=rho[1:obsi, 1:obsi]
res<-sigma^2*rhoi
}
else if (case==2) {if (case2==1) {
rhoi=rho[obsi, obsi]
nobs=sum(obsi)
nful=length(obsi)
nrho=sum(lower.tri(rho))
res=vector("list", 1+length(nrho))
row<-c()
for (i in 1:(nful-1)){
row<- append(row,rep(i, (nful-i)))
}
col<-c()
for (i in 1:(nful-1)){
col<- append(col,seq(i+1, nful))
}
w<- which(obsi==FALSE)
for (i in 1:(1+nrho)){
if (i==1){
res[[i]]=2*sigma*rhoi
} else
{
if (nobs==nful){
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=sigma^2
a[col[i-1],row[i-1]]=sigma^2
res[[i]]=a
}else
{
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=2*sigma
a[col[i-1],row[i-1]]=2*sigma
res[[i]]=a[-w,-w]
}
}
}
}
else if (case2==2) {
rhoi=rho[1:obsi, 1:obsi]
nrho=sum(lower.tri(rho))
res=vector("list", 1+length(nrho))
row<-c()
for (i in 1:(obsi-1)){
row<- append(row,rep(i, (obsi-i)))
}
col<-c()
for (i in 1:(obsi-1)){
col<- append(col,seq(i+1, obsi))
}
for (i in 1:(1+nrho)){
if (i==1){
res[[i]]=2*sigma*rhoi
} else
{
a=matrix(0, nrow=obsi, ncol=obsi)
a[row[i-1],col[i-1]]=sigma^2
a[col[i-1],row[i-1]]=sigma^2
res[[i]]=a
}
}
}
}
else if (case==3) {rhoi=rho[obsi, obsi]
nobs=sum(obsi)
nful=length(obsi)
nrho=sum(lower.tri(rho))
res = array(NA, c(1+nrho,1+nrho, nobs,nobs))
row<-c()
for (i in 1:(nful-1)){
row<- append(row,rep(i, (nful-i)))
}
col<-c()
for (i in 1:(nful-1)){
col<- append(col,seq(i+1, nful))
}
w<- which(obsi==FALSE)
for (i in 1:(1+nrho)){
for (j in 1:(1+nrho)){
if (i==1 & j==1){
res[i,j,,]=2*rhoi
} else
if (i==1 & j>1){
if (nobs==nful){
a=matrix(0, nrow=nful, ncol=nful)
a[row[j-1],col[j-1]]=2*sigma
a[col[j-1],row[j-1]]=2*sigma
res[i,j,,]=a
}else
{
a=matrix(0, nrow=nful, ncol=nful)
a[row[j-1],col[j-1]]=2*sigma
a[col[j-1],row[j-1]]=2*sigma
res[i,j,,]=a[-w,-w]
}
} else
if (j==1 & i>1){
if (nobs==nful){
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=2*sigma
a[col[i-1],row[i-1]]=2*sigma
res[i,j,,]=a
}else
{
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=2*sigma
a[col[i-1],row[i-1]]=2*sigma
res[i,j,,]=a[-w,-w]
}
} else
if(i>1 & j>1) {
res[i,j,,]=matrix(0, nobs, nobs)
}
}
}
}
} else
{sigma<-fun.trans2pos(sigma)
rho <- fun.trans2rho(rho)
if (case==1) {rhoi=rho[1:obsi, 1:obsi]
res<-sigma^2*rhoi
}
else if (case==2) {if (case2==1) {
rhoi=rho[obsi, obsi]
nobs=sum(obsi)
nful=length(obsi)
nrho=sum(lower.tri(rho))
res=vector("list", 1+length(nrho))
row<-c()
for (i in 1:(nful-1)){
row<- append(row,rep(i, (nful-i)))
}
col<-c()
for (i in 1:(nful-1)){
col<- append(col,seq(i+1, nful))
}
w<- which(obsi==FALSE)
for (i in 1:(1+nrho)){
if (i==1){
res[[i]]=2*sigma*rhoi
} else
{
if (nobs==nful){
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=sigma^2
a[col[i-1],row[i-1]]=sigma^2
res[[i]]=a
}else
{
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=2*sigma
a[col[i-1],row[i-1]]=2*sigma
res[[i]]=a[-w,-w]
}
}
}
}
else if (case2==2) {
rhoi=rho[1:obsi, 1:obsi]
nrho=sum(lower.tri(rho))
res=vector("list", 1+length(nrho))
row<-c()
for (i in 1:(nful-1)){
row<- append(row,rep(i, (nful-i)))
}
col<-c()
for (i in 1:(nful-1)){
col<- append(col,seq(i+1, nful))
}
for (i in 1:(1+nrho)){
if (i==1){
res[[i]]=2*sigma*rhoi
} else
{
a=matrix(0, nrow=obsi, ncol=obsi)
a[row[i-1],col[i-1]]=sigma^2
a[col[i-1],row[i-1]]=sigma^2
res[[i]]=a
}
}
}
}
else if (case==3) {rhoi=rho[obsi, obsi]
nobs=sum(obsi)
nful=length(obsi)
nrho=sum(lower.tri(rho))
res = array(NA, c(1+nrho,1+nrho, nobs,nobs))
row<-c()
for (i in 1:(nful-1)){
row<- append(row,rep(i, (nful-i)))
}
col<-c()
for (i in 1:(nful-1)){
col<- append(col,seq(i+1, nful))
}
w<- which(obsi==FALSE)
#to.upper<-function(X) t(X)[lower.tri(X,diag=FALSE)]
#v<-to.upper(temp)
#w<-which(is.na(v))
#if (row[j-1] %in% w | col[j-1] %in% w)
for (i in 1:(1+nrho)){
for (j in 1:(1+nrho)){
if (i==1 & j==1){
res[i,j,,]=2*rhoi
} else
if (i==1 & j>1){
if (nobs==nful){
a=matrix(0, nrow=nful, ncol=nful)
a[row[j-1],col[j-1]]=2*sigma
a[col[j-1],row[j-1]]=2*sigma
res[i,j,,]=a
}else
{
a=matrix(0, nrow=nful, ncol=nful)
a[row[j-1],col[j-1]]=2*sigma
a[col[j-1],row[j-1]]=2*sigma
res[i,j,,]=a[-w,-w]
}
} else
if (j==1 & i>1){
if (nobs==nful){
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=2*sigma
a[col[i-1],row[i-1]]=2*sigma
res[i,j,,]=a
}else
{
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=2*sigma
a[col[i-1],row[i-1]]=2*sigma
res[i,j,,]=a[-w,-w]
}
} else
if(i>1 & j>1) {
res[i,j,,]=matrix(0, nobs, nobs)
}
}
}
}
}
res
}
#' Utility function to generate missing status variables in longitudinal data with dropout and/or intermittent missingness.
#'
#' @param data the name of the panel dataset in the long format with each row denoting a subject-visit observation
#' for ALL the planned visits, regardless of being missed or not. When a subject is lost to follow up, the data set
#' must include the observation at the first time of being lost to follow up.
#' @param id the name of the level-2 clustering variable.
#' @param time the name of the variable denoting the time of the visit. Can set time=NULL if data is already sorted by id and time within
#' id.
#' @param y the name of the outcome variable of the interest that is subject to missingness.
#' @return a dataset with the following three new variables added:
#' \itemize{
#' \item g_ : missingness indicator, "O"-observed, "I"-intermittent missing, "D"-dropout
#' \item gp_: missingness indicator in the previous visit, "O"-observed, "I"-intermittent missing, "D"-dropout, "U"-undefined.
#' \item yp_: the immediately observed prior outcome.
#' }
#' @export definemissingstatus
#' @name definemissingstatus
#' @examples
#' qolefnew <- definemissingstatus(qolef, id=id, time=time, y=y)
definemissingstatus<-function(data, id, time, y) {
if (missing(data)) stop("Need to specify the name of the dataset for defining missingness status")
if (missing(id)) stop("Need to specify the name of the level-2 clustering variable for defining missingness status")
if (missing(y)) stop("Need to specify the name of the outcome variable for defining missingness status")
id=as.character(substitute(id))
if (!missing(time)) time=as.character(substitute(time))
y=as.character(substitute(y))
## first sort the data by id and time
if (!missing(time)) data = data[order(data[,id], data[,time]),]
uid = unique(data[,id])
lastobsT=numeric(length(uid))
g=gp=res=NULL
for (i in 1:length(uid)) {
datai <- data[data[,id]==uid[i],,drop=F]
ypi <-gi<-gpi<-numeric(dim(datai)[1])
maxT=nrow(datai)
if (any(!is.na(datai[,y]))) lastobsT[i] <- max((1:nrow(datai))[!is.na(datai[,y])]) else
lastobsT[i] <- 0
for (j in 1:nrow(datai)) {
## assign values to yp, g, gp
if (j==1) {
ypi[1]=NA;
if (!is.na(datai[j,y])) gi[j]= "O" else
if (lastobsT[i]<1) gi[j]="D" else
gi[j]="I"
gpi[j]="U"
if (j < maxT) ypi[j+1] <- datai[j,y]
} else
{
if (j < maxT) ypi[j+1] = ifelse(is.na(datai[j,y]), ypi[j], datai[j, y])
gpi[j] = gi[j-1]
if (j < lastobsT[i]) gi[j] = ifelse( is.na(datai[j,y]), "I", "O") else
if (j == lastobsT[i]) gi[j] = ifelse( is.na(datai[j,y]), "D","O") else
if (j > lastobsT[i]) gi[j] ="D"
}
} ## end of for loop for index j
res = rbind(res, data.frame(datai, yp_=ypi,g_=gi,gp_=gpi))
} ## end of the outer for loop for index i.
res
}
#' Function to fit the transitional missing data model and obtain the predicted probabilities of being observed for all observations.
#'
#' Fit missing data model and obtain predicted probabilities of being observed for all observations.
#' @param formula a formula to specify a multinomial transitional missing data model in the form of g+gp ~Sterms.
#' @param data the name of the dataset for fitting missing data mechansim model
#' @param subset an optional vector specifying a subset of observations to be used in the fitting process for the outcome model and the missingness mechanism model.
#' @param weights an optional vector of "prior weights" to be used in the fitting process for the missingness mechanism model.
#' Should be NULL or a numeric vector.
#' @aliases tmdm
#' @export tmdm
#' @examples
#' qolefitg <- tmdm(g+gp~as.factor(time)+group+perf+sever,data=qolef)
tmdm<- function(formula, data, weights, subset) {
## (1) process the call and set up model frame.
cl <- match.call()
if (missing(data)) data <- environment(formula)
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "weights"), names(mf), 0L)
mf <- mf[c(1L, m)]
f <- Formula(formula)
if (any(length(f) !=1 )) stop ("Model specification for the missing data model is wrong. Please check if missing data model is specified as g +gp ~Sterms.")
mf$formula <- f
mf[[1L]] <- as.name("model.frame")
options(na.action=na.pass)
mf <- eval(mf, parent.frame())
## This is the function to fit the multinomial logit transitional model for the missingness
if (missing(weights)) WTs <- NULL else
WTs <- model.extract(mf, "weights")
if (!is.null(WTs) && !is.numeric(WTs))
stop("'weights' must be a numeric vector")
if (!is.null(WTs) && any(WTs < 0))
stop("negative weights not allowed")
if (is.null(WTs))
WTs <- rep(1, nrow(mf))
WTs_ <- NULL
mf$WTs_= WTs
mf$obsprob_= NA ## initialize the hhat column.
mf$A10_= NA
mf$A20_= NA
mf$A11_= NA
##(2) Extract responses and predictors from the model frame
if (ncol(model.part(f,data=mf, lhs=1)) != 2) stop("Need to specify both g_ and gp_ as the responses for the transitional missing data model")
mf$g_<- model.part(f,data=mf, lhs=1)[,1]
mf$gp_<- model.part(f,data=mf, lhs=1)[,2]
gmodel=formula(f, lhs=1, rhs=1)
gmodel <- as.formula(paste("g_ ~ ", gmodel[3], sep=" " ))
s<- as.data.frame(model.matrix(gmodel, data=mf))
if (nrow(s)<1) print("No predictor variables specified for the outcome")
options(na.action=na.omit)
if (any(!(mf$g_ %in% c("O","I","D")))) stop("Data error: missing status variable must take D, I or O")
if (any(!(mf$gp_ %in% c("O","I","D","U")))) stop("Data error: missing status variable
for the prior visit must take D, I, O or U")
if (any(mf$gp_=="I" & mf$g_=="D")) {
if (!("D" %in% levels(mf$gp_))) levels(mf$gp_) <- c(levels(mf$gp_), "D")
mf$gp_[mf$gp_=="I" & mf$g_=="D"]= "D"
}
## set reference level of the response for logistic regression.
mf$g_=as.factor(mf$g_)
mf$g_=relevel(mf$g_, "O")
## (1) gp="D: set hhat=0 for all observations after dropout time
mf$obsprob_[mf$gp_=="D"]=0; mf$A10_[mf$gp_=="D"]=0; mf$A20_[mf$gp_=="D"]=0; mf$A11_[data$gp=="D"]=0
## (2) gp="O": fit multinomial logistic regression to obtain hhat
mf.v0 = mf[mf$gp_=="O",,drop=F]
s.v0= s[mf$gp_=="O",,drop=F]
if (! any(mf.v0$g_=="O")) { ## If the set of g values is all on "I", "D" or "I, D" in this subsample.This is essentially dropouts.
mf$obsprob_[mf$gp_=="O"]=0; mf$A10_[mf$gp_=="O"]=0; mf$A20_[mf$gp_=="O"]=0; mf$A11_[mf$gp_=="O"]=0;
} else if (all(mf.v0$g_=="O")) { ## If all data are observed in this subsample
mf$obsprob_[mf$gp_=="O"]=1; mf$A10_[mf$gp_=="O"]=1; mf$A20_[mf$gp_=="O"]=0; mf$A11_[mf$gp_=="O"]=0;
} else if (any(mf.v0$g_!="O")) {
missS <- apply(s.v0, 2, FUN=function (u) any(is.na(u)))
s.v0$WTs_=mf.v0$WTs_
s.v0$g_=mf.v0$g_
if (length(unique(mf.v0$g_))>2) { ## g includes "O, I, D" in this subsample.
mdm.v0 =multinom(g_ ~ .-1- WTs_, data=s.v0[,!missS], weights=WTs_, Hess=TRUE,maxit=500)
mdm.v0.fit=as.data.frame(mdm.v0$fit)
mf$obsprob_[mf$gp_=="O"]=mdm.v0.fit$O; mf$A10_[mf$gp_=="O"]=as.numeric(mf.v0$g_=="I")-mdm.v0.fit$I;
mf$A20_[mf$gp_=="O"]=as.numeric(mf.v0$g_=="D")-mdm.v0.fit$D; mf$A11_[mf$gp_=="O"]=0;
} else if (length(unique(mf.v0$g_))==2) { ## g includes "O, I", or "O,D" in this subsample.
mdm.v0 =multinom(g_ ~ .-1- WTs_, data=droplevels(s.v0[,!missS]), weights=WTs_, Hess=TRUE,maxit=500)
## mdm.v0 =glm(g_ ~ .-1- WTs_, data=s.v0[,!missS], weights=WTs_, family=binomial(link=logit))
mf$obsprob_[mf$gp_=="O"]=1-mdm.v0$fit;
mf$A10_[mf$gp_=="O"]=as.numeric(mf.v0$g_=="I")*(1-mdm.v0$fit);
mf$A20_[mf$gp_=="O"]=as.numeric(mf.v0$g_=="D")*(1-mdm.v0$fit);
mf$A11_[mf$gp_=="O"]=0;
}
} else print("Warning: check the values of g_ for the subsample gp=O")
## (3) gp="I":
if (any(mf$gp_=="I")) {
## The missingness status at prior visit gp_ is essentially a dropout.
##mf$gp_[mf$gp_=="I" & mf$g_=="D"]= "D"
mf.v1 = mf[mf$gp_=="I",,drop=F]
s.v1= s[mf$gp_=="I",,drop=F]
missS <- apply(s.v1, 2, FUN=function (u) any(is.na(u)))
s.v1$WTs_=mf.v1$WTs_
s.v1$g_=mf.v1$g_
if (any(mf.v1$g_=="D")) { ## if the set of g values is "D", "O, D", "I, D", "O,I,D"
stop("Data error: There should no dropout immediately after intermittent missingness. Change intermittent missingness to dropout ")
} else if (! any(mf.v1$g_=="O")) {## if the set of g values is "I" only. This means all are dropouts.
mf$obsprob_[mf$gp_=="I"]=0; mf$A10_[mf$gp_=="I"]=0; mf$A20_[mf$gp_=="I"]=0; mf$A11_[mf$gp_=="I"]=0;
} else if (all(mf.v1$g_=="O")) { ## If all data are observed in this subsample
mf$obsprob_[mf$gp_=="I"]=1; mf$A10_[mf$gp_=="I"]=1; mf$A20_[mf$gp_=="I"]=0; mf$A11_[mf$gp_=="I"]=0;
} else if (any(mf.v1$g_!="O")) { ## the set of g values is " O, I"
mdm.v1 =multinom(g_ ~ .-1- WTs_, data=droplevels(s.v1[,!missS]), weights=WTs_, Hess=TRUE,maxit=500)
## mdm.v1 =glm(g_ ~ .-1- WTs_, data=s.v1[,!missS], family=binomial(link=logit), weights=WTs_)
mf$obsprob_[mf$gp_=="I"]=mf$A11_[mf$gp_=="I"]=1-mdm.v1$fit;
mf$A10_[mf$gp_=="I"]=mf$A20_[mf$gp_=="I"]=0
}
}
## (4) gp="U": baseline observations
##mf$A10_[mf$gp_=="U"]=0; mf$A20_[mf$gp_=="U"]=0; mf$A11_[mf$gp_=="U"]=0
##mf$obsprob_[mf$gp_=="U"]=mf$A10_[mf$gp_=="U"]+ mf$A20_[mf$gp_=="U"]+ mf$A11_[mf$gp_=="U"];
mf.vU = mf[mf$gp_=="U",,drop=F]
s.vU= s[mf$gp_=="U",,drop=F]
if (! any(mf.vU$g_=="O")) { ## If the set of g values is all on "I", "D" or "I, D" in this subsample.Then the probability of being observed is zero.
mf$obsprob_[mf$gp_=="U"]=0; mf$A10_[mf$gp_=="U"]=0; mf$A20_[mf$gp_=="U"]=0; mf$A11_[mf$gp_=="U"]=0;
} else if (all(mf.vU$g_=="O")) { ## If all data are observed in this subsample
mf$obsprob_[mf$gp_=="U"]=1; mf$A10_[mf$gp_=="U"]=1; mf$A20_[mf$gp_=="U"]=0; mf$A11_[mf$gp_=="U"]=0;
} else if (any(mf.vU$g_!="O")) {
missS <- apply(s.vU, 2, FUN=function (u) any(is.na(u)))
s.vU$WTs_=mf.vU$WTs_
s.vU$g_=mf.vU$g_
if (length(unique(mf.vU$g_))>2) { ## g includes "O, I, D" in this subsample.
mdm.vU =multinom(g_ ~ .-1- WTs_, data=s.vU[,!missS], weights=WTs_, Hess=TRUE,maxit=500)
mdm.vU.fit=as.data.frame(mdm.vU$fit)
mf$obsprob_[mf$gp_=="U"]=mdm.vU.fit$O; mf$A10_[mf$gp_=="U"]=as.numeric(mf.vU$g_=="I")-mdm.vU.fit$I;
mf$A20_[mf$gp_=="U"]=as.numeric(mf.vU$g_=="D")-mdm.vU.fit$D; mf$A11_[mf$gp_=="U"]=0;
} else if (length(unique(mf.vU$g_))==2) { ## g includes "O, I", or "O,D" in this subsample.
mdm.vU =multinom(g_ ~ .-1- WTs_, data=droplevels(s.vU[,!missS]), weights=WTs_, Hess=TRUE,maxit=500)
## mdm.v0 =glm(g_ ~ .-1- WTs_, data=s.v0[,!missS], weights=WTs_, family=binomial(link=logit))
mf$obsprob_[mf$gp_=="U"]=1-mdm.vU$fit;
mf$A10_[mf$gp_=="U"]=as.numeric(mf.vU$g_=="I")*(1-mdm.vU$fit);
mf$A20_[mf$gp_=="U"]=as.numeric(mf.vU$g_=="D")*(1-mdm.vU$fit);
mf$A11_[mf$gp_=="U"]=0;
}
} ##else print("Warning: check the values of g_ for the subsample gp=U")
return(mf)
}
fun.trans2pos <- function(u) {
res<-exp(u)
res
}
fun.trans2rho <- function(u) {
tol=10^308
if (u==Inf) u= tol
if (u== -Inf) u= -tol
rho<- 2*exp(u)/(1+exp(u))-1
rho
}
fun.postrans <- function(u) {
if (u) stop("Argument has to be positive")
res <- log(u)
res
}
fun.rhotrans <- function(rho) {
if (abs(rho) >1) stop("correlation coefficient must be between 1 and -1")
tol= 1e-10
u = (1+rho)/2
if (u==0) u= u + tol
if (u==1) u= 1- tol
res<- log(u/(1-u))
res
}
Try the isni package in your browser
Any scripts or data that you put into this service are public.
isni documentation built on Aug. 23, 2021, 9: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 isnimgm
Documented in isnimgm
## Updates
## Nov 16, 2018
## Add codes to permit ISNI calculation when there are subjects with all outcomes are missing.
## The added codes are in isnimgm(), fun.mgmsubi(),
#' A data set for Psychiatric Drug Treatment
#'
#' The variables are as follows:
#' \itemize{
#' \item y. number of days
#' \item time. in weeks
#' \item sub. patients id
#' \item group. treatment group with placebo (PBO) and desipramine (DMI)
#' \item basey. baseline days
#' \item g1. dummy variables for treatment group
#' \item g2. dummy variables for treatment group
#' \item time1. dummy variables for time
#' \item time2. dummy variables for time
#' \item grptime. group by time
#' }
#'
#' @docType data
#' @keywords datasets
#' @name coc
#' @usage data(coc)
#' @format A data frame with 869 rows and 10 variables
NULL
#' A data set for Quality of Life Emontional Functioning outcome.
#'
#' The variables are as follows:
#' \itemize{
#' \item id: patients id
#' \item y: EF score
#' \item time: time in months since randomization.
#' \item group: placebo (0) or flutamide (1)
#' \item perf: baseline performance score
#' \item sever: baseline disease severity
## \item t1: indicator variable for month 1
## \item t3: indicator variable for month 3
## \item t6: indicator variable for month 6
# \item yp: most recently observed prior outcome
#' \item basey: EF at baseline
#' \item g: missingness status ("O"=observed, "D"=dropout, "I"=intermittent missingness)
#' \item gp: missingness status in the prior visit ("O"=observed, "D"=dropout, "I"=intermittent missingness, "U"=undefined)
#' }
#'
#' @docType data
#' @keywords datasets
#' @name qolef
#' @usage data(qolef)
#' @format A data frame with 2860 rows and 10 variables
NULL
#' Function for ISNI computation when the outcome follows marginal multivariate Gaussian Models.
#'
#' Calculate the ISNI when the regression outcome is subject to missingness and follows marginal multivaraite Gaussian models.
#' @param formula an object of class "Formula": a symbolic description of the models to be fitted for the outcome and missingness status variable.
#' The details of model specification are given under "Details".
#' @param data the name of data frame containing all the variables in the model and all the observations including those intended to be collected
#' but became missing.
## @param gmodel an object of class "formula": a symbolic description of the selection model to be fitted for the missingness indictor g
#' @param cortype the type of within-subject correlation structure.
#' @param id the name of variable for the level-2 clustering variable.
#' @param subset an optional vector specifying a subset of observations to be used in the fitting process for the outcome model and the missingness mechanism model.
#' @param weights frequency weights to be assigned to each \code{id}. when supplied, indicates differential weights are used; otherwise each \code{id} is weighted equally.
#' @param predprobobs Null if using buil-in multinomial transitional logistic model to obtain predicted probabilities of being observed;
#' otherwise user supply the name of the variable in \code{data} that gives these probabilities for all the observations in the data.
#' @param misni FALSE if using the default approach to computing ISNI with a scalar nonignorability parameter;
#' TRUE when computing ISNI with multiple nonignorability parameters.
#' @details The ISNI analysis is based on a joint selection model and requires specifying two model equations: the complete-data model and the missing data mechanism model.
#' To specify the variables in the models that are required for computing the isni measures, we make use of the \code{R} package "Formula" designed for handling model
#' equations with multiple responses and multiple sets of predictors . At a minimum, the user should supply a single-equation
#' in the typical form: \code{response} ~ \code{Xterms} where \code{response} is the (numeric or factor) vector for the outcome of interest and \code{Xterms}
#' is a series of terms, separated by + operators, which specify a linear predictor for response. With the signle-equation specification, the \code{isniglm} function
#' will by default use the utility function \code{definemissingstatus} provided in the package to generate the
#' missingness status variables at the current and prior visits and then use \code{Xterms} as the observed missingness predictors for fitting a first-order transitional missing data model.
#' It is important to sort within-\code{id} observations by time so that the missingness status variables can be defined correctly in this default setting. The \code{isnimgm} then computes the MAR estimates and conducts ISNI computation
#' to evaluate the rate of change of model estimates in the neighborhood of the MAR model where the missingness probability is allowed to depend on the unobserved value
#' of \code{response}, even after conditioning on the other missingness predictors.
#'
#' The above single-equation formula specification uses the same set of predictors for the outcome model and the missing data mechanism model for computing ISNI.
#' To use different sets of predictors, one can explicitly specifiy a two-equation formula as: \code{response} | \code{miss + missprior} ~ \code{Xterms} | \code{Sterms},
#' which specifies the formula for the outcome model as \code{response} ~ \code{Xterms} and that for the missing data mechanism model as \code{miss + missprior} ~ \code{Sterms},
#' where \code{Xterms} and the observed missingness predictors \code{Sterms} can be different, \code{miss} and \code{missprior} are the variable names in \code{data} denoting the missingness status at the
#' current and prior visits, respectively.
#' @references Ma G, Troxel AB, Heitjan DF. An Index of Local Sensitivity to Nonignorable Dropout in Longitudinal Modeling. Stat Med. 2005;24:2129-2150. \cr
#' Xie H. Analyzing Longitudinal Clinical Trial Data with Nonignorable Missingness and Unknown Missingness Reasons. Comput Stat Data Anal. 2012;56:1287-1300. \cr
#' Xie H, Qian Y. Measuring the Impact of Nonignorability in Panel Data with Non-Monotone Nonresponse. Journal of Applied Econometrics. 2012;27:129-159.
#' @name isnimgm
#' @aliases isnimgm
#' @import matrixcalc mvtnorm nlme nnet stats Formula
#' @importFrom nlme corCompSymm
#' @export isnimgm
#' @examples
#' models= y | g+gp ~ perf + sever+ as.factor(time) + group +as.factor(time):group |
#' as.factor(time) * group + yp+ perf + sever
## gmodel= as.factor(g)~ t3+t6+group+yp+perf+sever
## qolef$t12<-qolef$t1*qolef$group
## qolef$t32<-qolef$t3*qolef$group
## qolef$t62<-qolef$t6*qolef$group
##
#' qolef.isni=isnimgm(models, data=qolef, id=id)
#' summary(qolef.isni)
#'
## ymodel= y | g+gp ~ t1+t3+t6 + group + t12 + t32+t62 + perf + sever | t3+t6+group+yp+perf+sever
## summary(isnimgm(ymodel, data=qolef, id=id))
#'
isnimgm = function(formula, data, cortype="CS", id, subset, weights, predprobobs, misni=FALSE) {
## (1) process the call and set up model frame.
cl <- match.call()
if (missing(data)) data <- environment(formula)
if (missing(id)) stop("Users need to supply the name of the level-2 clusting variable in the data via the argument id.")
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "weights", "subset", "id","predprobobs"), names(mf), 0L)
mf <- mf[c(1L, m)]
f <- Formula(formula)
if (any(length(f)>2)) stop("Cannot have more than two models")
mf$formula <- f
options(na.action=na.pass)
mf[[1L]] <- as.name("get_all_vars")
mfvar<- eval(mf, parent.frame())
if (!missing(subset)) {
SubSet <-mfvar$subset
mfvar <-mfvar[SubSet,]
}
mf[[1L]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
ID= model.extract(mf, "id")
if (missing(weights)) WTs <- NULL else
WTs <- model.extract(mf, "weights")
if (!is.null(WTs) && !is.numeric(WTs))
stop("'weights' must be a numeric vector")
if (!is.null(WTs) && any(WTs < 0))
stop("negative weights not allowed")
if (is.null(WTs)){
IWT <- 0
WTs <- rep(1, nrow(mf))
}
isni_WTs <- isni_id <- isni_outcome_ <- NULL
mf= cbind(mf, mfvar[!(names(mfvar) %in% names(mf))], isni_id=ID, isni_WTs=WTs)
if (! missing(predprobobs)) mf$fitprobs_= model.extract(mf, "predprobobs")
##(2) Extract responses and predictors from the model frame
ymodel=formula(f, lhs=1, rhs=1)
x<- as.matrix(model.matrix(f, data=mf, rhs=1))
y<- model.part(f,data=mf, lhs=1)[,1] ## will be a vector
if (nrow(x)<1) print("No predictor variables specified for the outcome")
options(na.action=na.omit)
## check if the missing status indicator g and gp are specified, if not add g and gp in the lhs of formula
if (length(f)[1]==1) {
cat("\n NOTE: The model formula specification does not include the missingness status variables at the current and the prior visits.
The program will assume that data has already been sorted by the time variable within each level-2 ID variable and proceed to create the missingness status variables and then fit a transitional missingness mechanism model. \n")
f <- update(f, .| g_ +gp_ ~.)
mf <- cbind(mf, isni_outcome_=y)
mf <- definemissingstatus(data=mf, id=isni_id, y=isni_outcome_)
if (ncol(model.part(f,data=mf, lhs=2)) != 2) stop("The missing data model should be in the form of g+gp ~ Sterm")
}
gmodel <- formula(f, lhs=2,rhs=length(f)[2])
## ## Drop all observations with missing values in predictor matrix X to avoid NAs in ISNI computation.
missX<- apply(x, 2, FUN=function (u) any(is.na(u)))
if (any(missX)) cat(paste("\n Note: There are missing values in fixed-effect outcome predictors-- ", paste(dimnames(x)[[2]][missX], collapse=" "), ". Observations with missing values in these predictors will be dropped out from computing ISNI,
which may reduce ISNI values if the outcome is missing concurrently. Alternatively one can impute these missing predictors values and re-compute ISNI. ", sep=""))
missX <- apply(x, 1, FUN=function (u) any(is.na(u)))
WTs=WTs[!missX]
mf=mf[!missX,,drop=F]
options(na.action=na.pass)
x<- as.matrix(model.matrix(f, data=mf, rhs=1))
y<- model.part(f,data=mf, lhs=1)[,1]
xomit <- as.matrix(model.matrix(f, data=mf[!is.na(y), ], rhs=1))
if (ncol(x) != ncol(xomit)) {
cat(paste("\n All variable names in fixed-effects design matrix of ISNI analysis including observations with missing outcomes: \n", paste (dimnames(x)[[2]], collapse=" ")))
cat(paste("\n All variable names in fixed-effects design matrix of the MAR model excluding observations with missing outcomes: \n", paste (dimnames(xomit)[[2]], collapse=" ")))
stop("\n The design matrix for the MAR model and ISNI analysis are different as shown above. Please modify your model formula specification, e.g. avoiding using as.factor function
in model formula. ")
}
options(na.action=na.omit)
maxT = max(table(mf$isni_id))
sdy <- sd(y,na.rm=T)
## fit the mgm with the observed data.
if (cortype=="CS") iggls= gls(ymodel, data=mf,correlation = corCompSymm(form = ~ 1 | isni_id), weights=~1/isni_WTs, method="ML", na.action=na.exclude) else
if (cortype=="AR1") iggls= gls(ymodel, data=mf,correlation = corAR1(form = ~ 1 | isni_id), weights=~1/isni_WTs, method="ML",na.action=na.exclude) else
if (cortype=="UN") iggls= gls(ymodel, data=mf,correlation = corSymm(form = ~ 1 | isni_id), weights=~1/isni_WTs, method="ML",na.action=na.exclude) else
stop ("This Type of correlation is not implmented yet for ISNI computation")
## fit a missing data model to obtain fitted probabilities for being observed.
if (missing(predprobobs)) {
mf.mdm= tmdm(gmodel, data=mf)
mf$fitprobs=mf.mdm$obsprob
mf$A10= mf.mdm$A10; mf$A20= mf.mdm$A20; mf$A11= mf.mdm$A11
} else if ( misni==T)
stop("predprobobs needs to to a matrix to obtain ISNI vector for multiple nonignorability parameter")
D=list(sigma=iggls$sigma, rho=intervals(iggls)$corStruct[,"est."])
b<- iggls$coef
nb = length(b)
nD=length(D$sigma)+length(D$rho)
npar= nb + nD
uid = unique(mf$isni_id)
nabla12=matrix(0, nrow=npar, ncol=3)
## compute nabla11 and nabla12
nabla11=matrix(0, nrow=npar, ncol=npar)
for (i in 1:length(uid)) {
xi= as.matrix(x[mf$isni_id==uid[i],,drop=F])
yi= y[mf$isni_id==uid[i]]
## Increase nabla11 only if there is at least one observed outcome
if (!all(is.na(yi))) nabla11=nabla11+ fun.mgmsubi(yi=yi,xi=xi,maxT=maxT, b=b,D=D, cortype=cortype, transform=FALSE, case=1)
gfiti= NULL; Afiti=NULL
if (! missing(predprobobs)) gfiti= mf$fitprobs_[mf$isni_id==uid[i]]
else Afiti= cbind(mf$A10[mf$isni_id==uid[i]], mf$A20[mf$isni_id==uid[i]], mf$A11[mf$isni_id==uid[i]])
if (any(is.na(yi))) nabla12 = nabla12 + fun.mgmsubi(yi=yi,xi=xi,maxT=maxT, b=b,D=D, cortype=cortype, transform=FALSE, gfiti=gfiti,Afiti=Afiti, case=2)
}
invnabla11= solve(nabla11)
invnabla11[1:nb,1:nb]= -iggls$varBeta
isni=-invnabla11%*%nabla12
se=sqrt(diag((-invnabla11)))
bD= c(b,unlist(D))
if (cortype=="CS" | cortype=="AR1" ) {
names(bD)=c(names(b), names(D))
} else
if (cortype=="UN" ) {
names(bD)=c(names(b), 'sigma', row.names(intervals(iggls)$corStruct)) ##names(D[[2]]))
}
isnivec=isni
dimnames(isnivec)<-list(names(bD),c("ISNI_IO", "ISNI_DO", "ISNI_II"))
if (misni==FALSE) {
isni=apply(isni,1, sum); senstran<-abs((sdy*se)/isni)
} else
{
isni=apply(abs(isni),1, sum); senstran<-abs((sdy*se)/isni)
}
res=list(coefficients=bD,se=se,isni=c(isni),c=c(senstran), call=cl, isnivec=isnivec, misni=misni,
logLik=iggls$logLik, aic=summary(iggls)$AIC, bic=summary(iggls)$BIC)
class(res) = c(res$class, "isnimgm")
return(res)
}
#' Function to print out a summary of isnimgm object in a matrix form.
#'
#' @param object the isnimgm object obtained from the isnimgm function
#' @param digits the number of significant digits to use when printing
#' @param ... additional arguements
#' @export summary.isnimgm
#' @export
#' @name summary.isnimgm
#' @aliases summary.isnimgm
summary.isnimgm<-function(object, digits = max(3, getOption("digits") - 2),
...) {
if (class(object) != "isnimgm") stop('Invalid object class')
cat("\nCall:\n", paste(deparse(object$call), sep = "\n",
collapse = "\n"), "\n\n", sep = "")
## Name the columns
if (object$misni==T) isniname<-c('MAR Est.','Std. Err','MISNI','c') else
isniname<-c('MAR Est.','Std. Err','ISNI','c')
## Set up matrix to hold result
res<-matrix(0,length(object$coef),length(isniname))
dimnames(res)<-list(names(object$coef),isniname)
for (i in 1:length(object$coef))
res[i,]<- c(object$coef[i],object$se[i],object$isni[i],object$c[i])
printCoefmat(res, digits = digits, cs.ind = 1:2)
}
#' Function to print the isnimgm object.
#'
#' \code{print} method for class isnimgm
#' @param x the isnimgm object obtained from the isnimgm function
#' @param digits the number of significant digits to use when printing
#' @param ... further arguments passed to or from other methods.
#'
#' @return The function print.isnimgm prints the model call, isni and c statistics from the isnimgm object.
#' @name print.isnimgm
#' @aliases print.isnimgm
#' @export print.isnimgm
#' @export
print.isnimgm<-function(x, digits = max(3, getOption("digits") - 2), ...) {
if (class(x) != "isnimgm") stop('Invalid object class')
cat("\nCall:\n", paste(deparse(x$call), sep = "\n",
collapse = "\n"), "\n\n", sep = "")
if (length(x$coef)>0) {
cat("ISNIs:\n")
print.default(format(x$isni, digits = digits), print.gap = 2L,
quote = FALSE)
cat("\n")
cat("c statistics:\n")
print.default(format(x$c, digits = digits), print.gap = 2L,
quote = FALSE)
}
else cat("No coefficients\n")
cat("\nlogLik of the MAR model: ", paste(format(x$logLik,digits=digits), sep = "\n",
collapse = "\n"), "\n", sep = "")
cat("\nAIC of the MAR model: ", paste(format(x$aic, digits=digits), sep = "\n",
collapse = "\n"), "\n", sep = "")
cat("\nBIC of the MAR model: ", paste(format(x$bic, digits=digits), sep = "\n",
collapse = "\n"), "\n", sep = "")
cat("\n")
invisible(x)
}
## Internal Function for various subject-level computation required for ISNIGLS.
##
## Calculate subject-level quantities when the regression outcome is subject to missingness and follows marginal multivaraite Gaussian Models.
## @param yi vector of the response for the ith subject
## @param xi matrix of the covariates for the ith subject
## @param maxT maximum number of visits
## @param b the mean parameter vector beta
## @param D the vector of unique parameters in the variance-covariance matrix for the error term in the GLS model for Y
## @param correlation the form of within-subject correlation structure in the GLS model for Y
## @param transform logical indicating wether or not the parameter in D is transformed.
## @param gfiti vector of predicted probabilities of being observed for all the observations from the ith subject
## @param Afiti matrix of 3 columns of predicted transitional probabilities for the missing observations from the ith subject.
## @param case 1: calculated nabla11_i; 2: calculate nabla12_i
## @aliases fun.glssubi
## @export
fun.mgmsubi <- function ( yi, xi, maxT=maxT, b, D, cortype,transform=FALSE, gfiti=NULL, Afiti=NULL, case=1) {
## output
## res -- if case=1, res=hessian matrix
## if case=2, res=\nabla_12 for subject i.
## if case=3, res=negative loglikelihood, retired and unused in the package.
## if case=4, res=negative score value, retired and unused in the pacakge
## if case=5, res=der.meanymi, retired and unused in the package.
obsi = !is.na(yi)
xiobs = as.matrix(xi[obsi,,drop=F])
yiobs = as.vector(yi[obsi])
niobs = length(yiobs)
muiobs = as.vector(xiobs %*% as.matrix(b))
riobs= yiobs - muiobs
nb = length(b)
nD=length(D$sigma)+length(D$rho)
npar= nb + nD
sigma= D$sigma; rho=D$rho
## obtain the variance-covariance matrix
if (cortype=="CS") {
covyi= fun.csmat(sigma=sigma, rho=rho, obsi=length(yi),transform=transform, case=1)
}
if (cortype=="AR1") {
covyi= fun.ar1mat(sigma=sigma, rho=rho, obsi=length(yi),transform=transform, case=1)
}
if (cortype=="UN") {
col=matrix(0, nrow=maxT, ncol=maxT)
for (i in 1:(maxT-1)){
col[(i+1):maxT,i]=rho[((i-1)*maxT-i*(i-1)/2+1) : (i*maxT-i*(i+1)/2)]
}
rho2=col+t(col)+diag(maxT)
covyi= fun.generalmat(sigma=sigma, rho=rho2, obsi=length(yi),transform=transform, case=1)
}
covyiobs= covyi[obsi,obsi]
## only invert if there is yiobs is not empty
if (niobs>0) icovyiobs= solve(covyiobs)
##print(icovyiobs)
## case=1, compute the nabla11 for subject i
if (case==1) {
nabla11= matrix(0, npar, npar)
## return 0 if there is no observed outcome.
if (niobs==0) return(nabla11)
nabla11[1:nb, 1:nb]=-t(xiobs)%*%icovyiobs%*%xiobs
if (cortype=="CS") {
covid1 = fun.csmat(sigma=sigma,rho=rho,obsi=niobs,transform=transform,case=2)
covid2 = fun.csmat(sigma=sigma,rho=rho,obsi=niobs,transform=transform,case=3)
Ai= covid1
derAi= matrix(NA, nrow=niobs, ncol=niobs)
for (j in 1:length(Ai)) Ai[[j]]= icovyiobs%*%covid1[[j]]%*%icovyiobs
for (i in 1:nD) {
for (j in 1:nD){
derAi= -icovyiobs%*%(covid1[[i]]%*%icovyiobs%*%covid1[[j]]-covid2[i,j,,]+covid1[[j]]%*%icovyiobs%*%covid1[[i]])%*%icovyiobs
nabla11[(nb+i),(nb+j)]= -0.5*sum(diag(-Ai[[i]]%*%covid1[[j]]+icovyiobs%*%covid2[i,j,,]))+0.5*t(as.matrix(riobs))%*%derAi%*%riobs
}
}
for (j in 1:nD) {nabla11[1:nb,(nb+j)]= -0.5*t(xiobs)%*%(Ai[[j]]+t(Ai[[j]]))%*%riobs}
nabla11[(nb+1):npar, 1:nb]=t(nabla11[1:nb, (nb+1):(npar)])
}
if (cortype=="AR1") {
covid1 = fun.ar1mat(sigma=sigma,rho=rho,obsi=niobs,transform=transform,case=2)
covid2 = fun.ar1mat(sigma=sigma,rho=rho,obsi=niobs,transform=transform,case=3)
Ai= covid1
derAi= matrix(NA, nrow=niobs, ncol=niobs)
for (j in 1:length(Ai)) Ai[[j]]= icovyiobs%*%covid1[[j]]%*%icovyiobs
for (i in 1:nD) {
for (j in 1:nD){
derAi= -icovyiobs%*%(covid1[[i]]%*%icovyiobs%*%covid1[[j]]-covid2[i,j,,]+covid1[[j]]%*%icovyiobs%*%covid1[[i]])%*%icovyiobs
nabla11[(nb+i),(nb+j)]= -0.5*sum(diag(-Ai[[i]]%*%covid1[[j]]+icovyiobs%*%covid2[i,j,,]))+0.5*t(as.matrix(riobs))%*%derAi%*%riobs
}
}
for (j in 1:nD) nabla11[1:nb,(nb+j)]= -0.5*t(xiobs)%*%(Ai[[j]]+t(Ai[[j]]))%*%riobs
nabla11[(nb+1):npar, 1:nb] = t(nabla11[1:nb, (nb+1):(npar)])
}
if (cortype=="UN") {
covid1 = fun.generalmat(sigma=sigma,rho=rho2,obsi=obsi,transform=transform,case=2, case2=1)
covid2 = fun.generalmat(sigma=sigma,rho=rho2,obsi=obsi,transform=transform,case=3)
Ai= covid1
derAi= matrix(NA, nrow=niobs, ncol=niobs)
for (j in 1:length(Ai)) Ai[[j]]= icovyiobs%*%covid1[[j]]%*%icovyiobs
for (i in 1:nD) {
for (j in 1:nD){
derAi= -icovyiobs%*%(covid1[[i]]%*%icovyiobs%*%covid1[[j]]-covid2[i,j,,]+covid1[[j]]%*%icovyiobs%*%covid1[[i]])%*%icovyiobs
nabla11[(nb+i),(nb+j)]= -0.5*sum(diag(-Ai[[i]]%*%covid1[[j]]+icovyiobs%*%covid2[i,j,,]))+0.5*t(as.matrix(riobs))%*%derAi%*%riobs
}
}
for (j in 1:nD) nabla11[1:nb,(nb+j)]= -0.5*t(xiobs)%*%(Ai[[j]]+t(Ai[[j]]))%*%riobs
nabla11[(nb+1):npar, 1:nb]=t(nabla11[1:nb, (nb+1):(npar)])
}
return(nabla11)
}## end if case=1
## case=2, calculate the \nabla_12 for subject i.
if (case==2) {
if (any(is.na(yi))){
misi = is.na(yi)
ximis = as.matrix(xi[misi,,drop=F])
yimis = as.vector(yi[misi])
##fitp0i = gfiti[misi]
Afit0i=matrix(0, nrow=sum(misi==T), ncol=3)
if (is.null(gfiti)) Afit0i =Afiti[misi, ,drop=F]
else Afit0i[,1]= gfiti[misi]
cit = covyi[misi, obsi,drop=F]
dmean.D= matrix(0, nrow=nD, ncol=ncol(Afit0i))
if (niobs==0) {dmean.beta=t(ximis)%*% Afit0i; return(rbind(dmean.beta, dmean.D)) }
dmean.beta = t(ximis - cit %*% icovyiobs %*% xiobs) %*% Afit0i ## as.matrix(fitp0i)
##numeric(nD)
if (cortype=="CS") {
derD= fun.csmat(sigma=sigma,rho=rho,obsi=length(yi),transform=transform,case=2)
} else
if (cortype=="AR1") {
derD= fun.ar1mat(sigma=sigma,rho=rho,obsi=length(yi),transform=transform,case=2)
}
if (cortype=="UN") {
derD= fun.generalmat(sigma=sigma,rho=rho2,obsi=length(yi),transform=transform,case=2, case2=2)
}
for (k in 1:nD) {
Aij<-icovyiobs%*%derD[[k]][obsi,obsi]%*%icovyiobs
dercitDk<- derD[[k]][misi,obsi,drop=F]
dmean.D[k,]<-t(dercitDk%*%icovyiobs%*%as.matrix(riobs) -cit%*%Aij%*%as.matrix(riobs)) %*% Afit0i ##as.matrix(fitp0i)
}
## dercov.sigma<-fun.csmat(sigma=sigma,rho=rho,nobs=length(yi),transform=transform,case=2)
## Aij<-icovyiobs%*%dercov.sigma[obsi,obsi]%*%icovyiobs
## dercit.sigma<- dercov.sigma[misi,obsi,drop=F]
## dmean.sigma<-t(dercit.sigma%*%icovyiobs%*%as.matrix(riobs) -cit%*%Aij%*%as.matrix(riobs)) %*% Afit0i ##as.matrix(fitp0i)
## dercov.rho<-fun.csmat(sigma=sigma,rho=rho,nobs=length(yi),transform=transform,case=3)
## Aij<-icovyiobs%*%dercov.rho[obsi,obsi]%*%icovyiobs
## dercit.rho<- dercov.rho[misi,obsi,drop=F]
## dmean.rho<-t(dercit.rho%*%icovyiobs%*%as.matrix(riobs) -cit%*%Aij%*%as.matrix(riobs))%*% Afit0i ##as.matrix(fitp0i)
res = rbind(dmean.beta, dmean.D) ##c(dmean.beta, dmean.D)
return(res)
} else
res = numeric(dim(xiobs)[2]+2)
} ## end if case=2
res
}
## Internal Function for generate compound symmetry matrix and its first and second derivatives wrt D
##
## Calculate the compound symmetry matrix and its first and second derivatives wrt D.
## @param sigma standard deviation
## @param rho correlation coefficient
## @param obsi dimension of the square matrix
## @param transform logical indicating wether or not the parameter in D is transformed.
## @param case 1: return nobs x nobs CS matrix; 2: return a list for first derivatives of the compound symmetry matrix wrt sigma and rho;
## 3: return an array containing 2nd derivatives of the compound symmetry matrix wrt sigma and rho.
## @aliases fun.csmat
## @export
fun.csmat<-function(sigma,rho,obsi,transform,case=1) {
## calculate the the related quantity for compound symmetry vcov matrix.
## input
## sigma -- the standard deviation
## rho --- the correlation coefficient
## nobs --- the number of observations, nobs>=1
## transform --- logical number indicator if the parameters are transformed
## output
## case=1 --- return nobs x nobs CS matrix
## case=2 --- return a list for first derivatives of the compound symmetry matrix wrt sigma and rho
## case=3 --- return an array containing 2nd derivatives of the compound symmetry matrix wrt sigma and rho
if (transform==FALSE) {
if (case==1) {
nobs = obsi
res<-sigma^2*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs))
}
else if (case==2) {
nobs = obsi
derSigma<- 2*sigma*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs))
derRho<- sigma^2*(matrix(1,nrow=nobs,ncol=nobs)-diag(1,nrow=nobs))
res <- list(derSigma=derSigma, derRho=derRho)
res
}
## derivative wrt to rho
else if (case==3) {
nobs = obsi
res = array(NA, c(2,2, nobs,nobs))
res[1,1,,]=2*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs))
res[1,2,,]=res[2,1,,]=2*sigma*(matrix(1,nrow=nobs,ncol=nobs)- diag(1, nrow=nobs))
res[2,2,,]=matrix(0, nrow=nobs, ncol=nobs)
res
}
} else
{
sigma<-fun.trans2pos(sigma)
rho <- fun.trans2rho(rho)
if (case==1) {
nobs = obsi
res<-sigma^2*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs))
}
else if (case==2) {
nobs = obsi
derSigma<- 2*sigma*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs))
derRho<- sigma^2*(matrix(1,nrow=nobs,ncol=nobs)-diag(1,nrow=nobs))
res <- list(derSigma=derSigma, derRho=derRho)
res
}
## derivative wrt to rho
else if (case==3) {
nobs = obsi
res = array(NA, c(2,2, nobs,nobs))
res[1,1,,]=2*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs))
res[1,2,,]=res[2,1,,]=2*sigma*(matrix(1,nrow=nobs,ncol=nobs)- diag(1, nrow=nobs))
res[2,2,,]=matrix(0, nrow=nobs, ncol=nobs)
res
}
}
res
}
## Internal Function for generate AR1 matrix and its first and second derivatives wrt D
##
## Calculate the autoregressive matrix and its first and second derivatives wrt D.
## @param sigma standard deviation
## @param rho autoregressive parameter
## @param obsi dimision of the square matrix
## @param transform logical indicating wether or not the parameter in D is transformed.
## @param case 1: return nobs x nobs AR1 matrix; 2: return a list for first derivatives of the compound symmetry matrix wrt sigma and rho;
## 3: return an array containing 2nd derivatives of the compound symmetry matrix wrt sigma and rho.
## @aliases fun.ar1mat
## @export
fun.ar1mat<-function(sigma,rho,obsi,transform,case=1) {
## calculate the the related quantity for AR1 (autoregressive process of order 1) vcov matrix.
## input
## sigma -- the standard deviation
## rho --- the autoregressive parameter
## nobs --- the number of observations, nobs>=1
## transform --- logical number indicator if the parameters are transformed
## output
## case=1 --- return nobs x nobs AR1 matrix
## case=2 --- return a list for first derivatives of the AR1 matrix wrt sigma and rho
## case=3 --- return an array containing 2nd derivatives of the AR1 matrix wrt sigma and rho
if (transform==FALSE) {
if (case==1) {nobs = obsi
res<-matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
res[i,j]<-sigma^2*rho^abs(i-j)
}
}
} else
if (case==2) {nobs = obsi
derSigma<-matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
derSigma[i,j] <- 2*sigma*rho^abs(i-j)
}}
derRho<-matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
derRho[i,j] <- sigma^2*abs(i-j)*rho^(abs(i-j)-1)
}}
res <- list(derSigma=derSigma, derRho=derRho)
} else
## derivative wrt to rho
if (case==3) {nobs = obsi
res = array(NA, c(2,2, nobs,nobs))
res11 <- matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
res11[i,j]<-2*rho^abs(i-j)
}}
res[1,1,,]=res11
res12 <- matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs){
for (j in 1:nobs){
res12[i,j]<-2*sigma*abs(i-j)*rho^(abs(i-j)-1)
}}
res[1,2,,]=res[2,1,,]=res12
res22 <- matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs){
for (j in 1:nobs){
if(abs(i-j)<=1) res22[i,j] <- 0 else
res22[i,j]<-sigma^2*abs(i-j)*(abs(i-j)-1)*rho^(abs(i-j)-2)
}}
res[2,2,,]=res22
}
} else
{
sigma<-fun.trans2pos(sigma)
rho <- fun.trans2rho(rho)
if (case==1) {nobs = obsi
res<-matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
res[i,j]<-sigma^2*rho^abs(i-j)
}}
} else
if (case==2) {nobs = obsi
derSigma<-matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
derSigma[i,j] <- 2*sigma*rho^abs(i-j)
}}
derRho<-matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
derRho[i,j] <- sigma^2*abs(i-j)*rho^(abs(i-j)-1)
}}
res <- list(derSigma=derSigma, derRho=derRho)
} else
## derivative wrt to rho
if (case==3) {nobs = obsi
res = array(NA, c(2,2, nobs,nobs))
res11 <- matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs) {
for (j in 1:nobs) {
res11[i,j]<-2*rho^abs(i-j)
}}
res[1,1,,]=res11
res12 <- matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs){
for (j in 1:nobs){
res12[i,j]<-2*sigma*abs(i-j)*rho^(abs(i-j)-1)
}}
res[1,2,,]=res[2,1,,]=res12
res22 <- matrix(0, nrow=nobs, ncol=nobs)
for (i in 1:nobs){
for (j in 1:nobs){
if(abs(i-j)<=1) res22[i,j] <- 0 else
res22[i,j]<-abs(i-j)*(abs(i-j)-1)*rho^(abs(i-j)-2)
}}
res[2,2,,]=res22
}
}
res
}
## Internal Function for generate a general symmetric matrix and its first and second derivatives wrt D
## Calculate the general symmetric matrix and its first and second derivatives wrt D.
## @param sigma standard deviation
## @param rho correlation coefficient
## @param obsi dimension of the square matrix
## @param transform logical indicating wether or not the parameter in D is transformed.
## @param case 1: return nobs x nobs general symmetric matrix; 2: return a list for first derivatives of the general symmetric matrix wrt D;
## 3: return an array containing 2nd derivatives of the general symmetric matrix wrt D.
## @param case2 indicator for variance and covariance parameters.
## @aliases fun.generalmat
## @export
fun.generalmat<-function(sigma,rho,obsi,transform,case=1, case2=NULL) {
## calculate the the related quantity for compound symmetry vcov matrix.
## input
## sigma -- the standard deviation
## rho --- the correlation coefficient
## nobs --- the number of observations, nobs>=1
## transform --- logical number indicator if the parameters are transformed
## output
## case=1 --- return nobs x nobs CS matrix
## case=2 --- return a list for first derivatives of the compound symmetry matrix wrt sigma and rho
## case=3 --- return an array containing 2nd derivatives of the compound symmetry matrix wrt sigma and rho
if (transform==FALSE) {
if (case==1) {rhoi=rho[1:obsi, 1:obsi]
res<-sigma^2*rhoi
}
else if (case==2) {if (case2==1) {
rhoi=rho[obsi, obsi]
nobs=sum(obsi)
nful=length(obsi)
nrho=sum(lower.tri(rho))
res=vector("list", 1+length(nrho))
row<-c()
for (i in 1:(nful-1)){
row<- append(row,rep(i, (nful-i)))
}
col<-c()
for (i in 1:(nful-1)){
col<- append(col,seq(i+1, nful))
}
w<- which(obsi==FALSE)
for (i in 1:(1+nrho)){
if (i==1){
res[[i]]=2*sigma*rhoi
} else
{
if (nobs==nful){
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=sigma^2
a[col[i-1],row[i-1]]=sigma^2
res[[i]]=a
}else
{
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=2*sigma
a[col[i-1],row[i-1]]=2*sigma
res[[i]]=a[-w,-w]
}
}
}
}
else if (case2==2) {
rhoi=rho[1:obsi, 1:obsi]
nrho=sum(lower.tri(rho))
res=vector("list", 1+length(nrho))
row<-c()
for (i in 1:(obsi-1)){
row<- append(row,rep(i, (obsi-i)))
}
col<-c()
for (i in 1:(obsi-1)){
col<- append(col,seq(i+1, obsi))
}
for (i in 1:(1+nrho)){
if (i==1){
res[[i]]=2*sigma*rhoi
} else
{
a=matrix(0, nrow=obsi, ncol=obsi)
a[row[i-1],col[i-1]]=sigma^2
a[col[i-1],row[i-1]]=sigma^2
res[[i]]=a
}
}
}
}
else if (case==3) {rhoi=rho[obsi, obsi]
nobs=sum(obsi)
nful=length(obsi)
nrho=sum(lower.tri(rho))
res = array(NA, c(1+nrho,1+nrho, nobs,nobs))
row<-c()
for (i in 1:(nful-1)){
row<- append(row,rep(i, (nful-i)))
}
col<-c()
for (i in 1:(nful-1)){
col<- append(col,seq(i+1, nful))
}
w<- which(obsi==FALSE)
for (i in 1:(1+nrho)){
for (j in 1:(1+nrho)){
if (i==1 & j==1){
res[i,j,,]=2*rhoi
} else
if (i==1 & j>1){
if (nobs==nful){
a=matrix(0, nrow=nful, ncol=nful)
a[row[j-1],col[j-1]]=2*sigma
a[col[j-1],row[j-1]]=2*sigma
res[i,j,,]=a
}else
{
a=matrix(0, nrow=nful, ncol=nful)
a[row[j-1],col[j-1]]=2*sigma
a[col[j-1],row[j-1]]=2*sigma
res[i,j,,]=a[-w,-w]
}
} else
if (j==1 & i>1){
if (nobs==nful){
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=2*sigma
a[col[i-1],row[i-1]]=2*sigma
res[i,j,,]=a
}else
{
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=2*sigma
a[col[i-1],row[i-1]]=2*sigma
res[i,j,,]=a[-w,-w]
}
} else
if(i>1 & j>1) {
res[i,j,,]=matrix(0, nobs, nobs)
}
}
}
}
} else
{sigma<-fun.trans2pos(sigma)
rho <- fun.trans2rho(rho)
if (case==1) {rhoi=rho[1:obsi, 1:obsi]
res<-sigma^2*rhoi
}
else if (case==2) {if (case2==1) {
rhoi=rho[obsi, obsi]
nobs=sum(obsi)
nful=length(obsi)
nrho=sum(lower.tri(rho))
res=vector("list", 1+length(nrho))
row<-c()
for (i in 1:(nful-1)){
row<- append(row,rep(i, (nful-i)))
}
col<-c()
for (i in 1:(nful-1)){
col<- append(col,seq(i+1, nful))
}
w<- which(obsi==FALSE)
for (i in 1:(1+nrho)){
if (i==1){
res[[i]]=2*sigma*rhoi
} else
{
if (nobs==nful){
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=sigma^2
a[col[i-1],row[i-1]]=sigma^2
res[[i]]=a
}else
{
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=2*sigma
a[col[i-1],row[i-1]]=2*sigma
res[[i]]=a[-w,-w]
}
}
}
}
else if (case2==2) {
rhoi=rho[1:obsi, 1:obsi]
nrho=sum(lower.tri(rho))
res=vector("list", 1+length(nrho))
row<-c()
for (i in 1:(nful-1)){
row<- append(row,rep(i, (nful-i)))
}
col<-c()
for (i in 1:(nful-1)){
col<- append(col,seq(i+1, nful))
}
for (i in 1:(1+nrho)){
if (i==1){
res[[i]]=2*sigma*rhoi
} else
{
a=matrix(0, nrow=obsi, ncol=obsi)
a[row[i-1],col[i-1]]=sigma^2
a[col[i-1],row[i-1]]=sigma^2
res[[i]]=a
}
}
}
}
else if (case==3) {rhoi=rho[obsi, obsi]
nobs=sum(obsi)
nful=length(obsi)
nrho=sum(lower.tri(rho))
res = array(NA, c(1+nrho,1+nrho, nobs,nobs))
row<-c()
for (i in 1:(nful-1)){
row<- append(row,rep(i, (nful-i)))
}
col<-c()
for (i in 1:(nful-1)){
col<- append(col,seq(i+1, nful))
}
w<- which(obsi==FALSE)
#to.upper<-function(X) t(X)[lower.tri(X,diag=FALSE)]
#v<-to.upper(temp)
#w<-which(is.na(v))
#if (row[j-1] %in% w | col[j-1] %in% w)
for (i in 1:(1+nrho)){
for (j in 1:(1+nrho)){
if (i==1 & j==1){
res[i,j,,]=2*rhoi
} else
if (i==1 & j>1){
if (nobs==nful){
a=matrix(0, nrow=nful, ncol=nful)
a[row[j-1],col[j-1]]=2*sigma
a[col[j-1],row[j-1]]=2*sigma
res[i,j,,]=a
}else
{
a=matrix(0, nrow=nful, ncol=nful)
a[row[j-1],col[j-1]]=2*sigma
a[col[j-1],row[j-1]]=2*sigma
res[i,j,,]=a[-w,-w]
}
} else
if (j==1 & i>1){
if (nobs==nful){
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=2*sigma
a[col[i-1],row[i-1]]=2*sigma
res[i,j,,]=a
}else
{
a=matrix(0, nrow=nful, ncol=nful)
a[row[i-1],col[i-1]]=2*sigma
a[col[i-1],row[i-1]]=2*sigma
res[i,j,,]=a[-w,-w]
}
} else
if(i>1 & j>1) {
res[i,j,,]=matrix(0, nobs, nobs)
}
}
}
}
}
res
}
#' Utility function to generate missing status variables in longitudinal data with dropout and/or intermittent missingness.
#'
#' @param data the name of the panel dataset in the long format with each row denoting a subject-visit observation
#' for ALL the planned visits, regardless of being missed or not. When a subject is lost to follow up, the data set
#' must include the observation at the first time of being lost to follow up.
#' @param id the name of the level-2 clustering variable.
#' @param time the name of the variable denoting the time of the visit. Can set time=NULL if data is already sorted by id and time within
#' id.
#' @param y the name of the outcome variable of the interest that is subject to missingness.
#' @return a dataset with the following three new variables added:
#' \itemize{
#' \item g_ : missingness indicator, "O"-observed, "I"-intermittent missing, "D"-dropout
#' \item gp_: missingness indicator in the previous visit, "O"-observed, "I"-intermittent missing, "D"-dropout, "U"-undefined.
#' \item yp_: the immediately observed prior outcome.
#' }
#' @export definemissingstatus
#' @name definemissingstatus
#' @examples
#' qolefnew <- definemissingstatus(qolef, id=id, time=time, y=y)
definemissingstatus<-function(data, id, time, y) {
if (missing(data)) stop("Need to specify the name of the dataset for defining missingness status")
if (missing(id)) stop("Need to specify the name of the level-2 clustering variable for defining missingness status")
if (missing(y)) stop("Need to specify the name of the outcome variable for defining missingness status")
id=as.character(substitute(id))
if (!missing(time)) time=as.character(substitute(time))
y=as.character(substitute(y))
## first sort the data by id and time
if (!missing(time)) data = data[order(data[,id], data[,time]),]
uid = unique(data[,id])
lastobsT=numeric(length(uid))
g=gp=res=NULL
for (i in 1:length(uid)) {
datai <- data[data[,id]==uid[i],,drop=F]
ypi <-gi<-gpi<-numeric(dim(datai)[1])
maxT=nrow(datai)
if (any(!is.na(datai[,y]))) lastobsT[i] <- max((1:nrow(datai))[!is.na(datai[,y])]) else
lastobsT[i] <- 0
for (j in 1:nrow(datai)) {
## assign values to yp, g, gp
if (j==1) {
ypi[1]=NA;
if (!is.na(datai[j,y])) gi[j]= "O" else
if (lastobsT[i]<1) gi[j]="D" else
gi[j]="I"
gpi[j]="U"
if (j < maxT) ypi[j+1] <- datai[j,y]
} else
{
if (j < maxT) ypi[j+1] = ifelse(is.na(datai[j,y]), ypi[j], datai[j, y])
gpi[j] = gi[j-1]
if (j < lastobsT[i]) gi[j] = ifelse( is.na(datai[j,y]), "I", "O") else
if (j == lastobsT[i]) gi[j] = ifelse( is.na(datai[j,y]), "D","O") else
if (j > lastobsT[i]) gi[j] ="D"
}
} ## end of for loop for index j
res = rbind(res, data.frame(datai, yp_=ypi,g_=gi,gp_=gpi))
} ## end of the outer for loop for index i.
res
}
#' Function to fit the transitional missing data model and obtain the predicted probabilities of being observed for all observations.
#'
#' Fit missing data model and obtain predicted probabilities of being observed for all observations.
#' @param formula a formula to specify a multinomial transitional missing data model in the form of g+gp ~Sterms.
#' @param data the name of the dataset for fitting missing data mechansim model
#' @param subset an optional vector specifying a subset of observations to be used in the fitting process for the outcome model and the missingness mechanism model.
#' @param weights an optional vector of "prior weights" to be used in the fitting process for the missingness mechanism model.
#' Should be NULL or a numeric vector.
#' @aliases tmdm
#' @export tmdm
#' @examples
#' qolefitg <- tmdm(g+gp~as.factor(time)+group+perf+sever,data=qolef)
tmdm<- function(formula, data, weights, subset) {
## (1) process the call and set up model frame.
cl <- match.call()
if (missing(data)) data <- environment(formula)
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "weights"), names(mf), 0L)
mf <- mf[c(1L, m)]
f <- Formula(formula)
if (any(length(f) !=1 )) stop ("Model specification for the missing data model is wrong. Please check if missing data model is specified as g +gp ~Sterms.")
mf$formula <- f
mf[[1L]] <- as.name("model.frame")
options(na.action=na.pass)
mf <- eval(mf, parent.frame())
## This is the function to fit the multinomial logit transitional model for the missingness
if (missing(weights)) WTs <- NULL else
WTs <- model.extract(mf, "weights")
if (!is.null(WTs) && !is.numeric(WTs))
stop("'weights' must be a numeric vector")
if (!is.null(WTs) && any(WTs < 0))
stop("negative weights not allowed")
if (is.null(WTs))
WTs <- rep(1, nrow(mf))
WTs_ <- NULL
mf$WTs_= WTs
mf$obsprob_= NA ## initialize the hhat column.
mf$A10_= NA
mf$A20_= NA
mf$A11_= NA
##(2) Extract responses and predictors from the model frame
if (ncol(model.part(f,data=mf, lhs=1)) != 2) stop("Need to specify both g_ and gp_ as the responses for the transitional missing data model")
mf$g_<- model.part(f,data=mf, lhs=1)[,1]
mf$gp_<- model.part(f,data=mf, lhs=1)[,2]
gmodel=formula(f, lhs=1, rhs=1)
gmodel <- as.formula(paste("g_ ~ ", gmodel[3], sep=" " ))
s<- as.data.frame(model.matrix(gmodel, data=mf))
if (nrow(s)<1) print("No predictor variables specified for the outcome")
options(na.action=na.omit)
if (any(!(mf$g_ %in% c("O","I","D")))) stop("Data error: missing status variable must take D, I or O")
if (any(!(mf$gp_ %in% c("O","I","D","U")))) stop("Data error: missing status variable
for the prior visit must take D, I, O or U")
if (any(mf$gp_=="I" & mf$g_=="D")) {
if (!("D" %in% levels(mf$gp_))) levels(mf$gp_) <- c(levels(mf$gp_), "D")
mf$gp_[mf$gp_=="I" & mf$g_=="D"]= "D"
}
## set reference level of the response for logistic regression.
mf$g_=as.factor(mf$g_)
mf$g_=relevel(mf$g_, "O")
## (1) gp="D: set hhat=0 for all observations after dropout time
mf$obsprob_[mf$gp_=="D"]=0; mf$A10_[mf$gp_=="D"]=0; mf$A20_[mf$gp_=="D"]=0; mf$A11_[data$gp=="D"]=0
## (2) gp="O": fit multinomial logistic regression to obtain hhat
mf.v0 = mf[mf$gp_=="O",,drop=F]
s.v0= s[mf$gp_=="O",,drop=F]
if (! any(mf.v0$g_=="O")) { ## If the set of g values is all on "I", "D" or "I, D" in this subsample.This is essentially dropouts.
mf$obsprob_[mf$gp_=="O"]=0; mf$A10_[mf$gp_=="O"]=0; mf$A20_[mf$gp_=="O"]=0; mf$A11_[mf$gp_=="O"]=0;
} else if (all(mf.v0$g_=="O")) { ## If all data are observed in this subsample
mf$obsprob_[mf$gp_=="O"]=1; mf$A10_[mf$gp_=="O"]=1; mf$A20_[mf$gp_=="O"]=0; mf$A11_[mf$gp_=="O"]=0;
} else if (any(mf.v0$g_!="O")) {
missS <- apply(s.v0, 2, FUN=function (u) any(is.na(u)))
s.v0$WTs_=mf.v0$WTs_
s.v0$g_=mf.v0$g_
if (length(unique(mf.v0$g_))>2) { ## g includes "O, I, D" in this subsample.
mdm.v0 =multinom(g_ ~ .-1- WTs_, data=s.v0[,!missS], weights=WTs_, Hess=TRUE,maxit=500)
mdm.v0.fit=as.data.frame(mdm.v0$fit)
mf$obsprob_[mf$gp_=="O"]=mdm.v0.fit$O; mf$A10_[mf$gp_=="O"]=as.numeric(mf.v0$g_=="I")-mdm.v0.fit$I;
mf$A20_[mf$gp_=="O"]=as.numeric(mf.v0$g_=="D")-mdm.v0.fit$D; mf$A11_[mf$gp_=="O"]=0;
} else if (length(unique(mf.v0$g_))==2) { ## g includes "O, I", or "O,D" in this subsample.
mdm.v0 =multinom(g_ ~ .-1- WTs_, data=droplevels(s.v0[,!missS]), weights=WTs_, Hess=TRUE,maxit=500)
## mdm.v0 =glm(g_ ~ .-1- WTs_, data=s.v0[,!missS], weights=WTs_, family=binomial(link=logit))
mf$obsprob_[mf$gp_=="O"]=1-mdm.v0$fit;
mf$A10_[mf$gp_=="O"]=as.numeric(mf.v0$g_=="I")*(1-mdm.v0$fit);
mf$A20_[mf$gp_=="O"]=as.numeric(mf.v0$g_=="D")*(1-mdm.v0$fit);
mf$A11_[mf$gp_=="O"]=0;
}
} else print("Warning: check the values of g_ for the subsample gp=O")
## (3) gp="I":
if (any(mf$gp_=="I")) {
## The missingness status at prior visit gp_ is essentially a dropout.
##mf$gp_[mf$gp_=="I" & mf$g_=="D"]= "D"
mf.v1 = mf[mf$gp_=="I",,drop=F]
s.v1= s[mf$gp_=="I",,drop=F]
missS <- apply(s.v1, 2, FUN=function (u) any(is.na(u)))
s.v1$WTs_=mf.v1$WTs_
s.v1$g_=mf.v1$g_
if (any(mf.v1$g_=="D")) { ## if the set of g values is "D", "O, D", "I, D", "O,I,D"
stop("Data error: There should no dropout immediately after intermittent missingness. Change intermittent missingness to dropout ")
} else if (! any(mf.v1$g_=="O")) {## if the set of g values is "I" only. This means all are dropouts.
mf$obsprob_[mf$gp_=="I"]=0; mf$A10_[mf$gp_=="I"]=0; mf$A20_[mf$gp_=="I"]=0; mf$A11_[mf$gp_=="I"]=0;
} else if (all(mf.v1$g_=="O")) { ## If all data are observed in this subsample
mf$obsprob_[mf$gp_=="I"]=1; mf$A10_[mf$gp_=="I"]=1; mf$A20_[mf$gp_=="I"]=0; mf$A11_[mf$gp_=="I"]=0;
} else if (any(mf.v1$g_!="O")) { ## the set of g values is " O, I"
mdm.v1 =multinom(g_ ~ .-1- WTs_, data=droplevels(s.v1[,!missS]), weights=WTs_, Hess=TRUE,maxit=500)
## mdm.v1 =glm(g_ ~ .-1- WTs_, data=s.v1[,!missS], family=binomial(link=logit), weights=WTs_)
mf$obsprob_[mf$gp_=="I"]=mf$A11_[mf$gp_=="I"]=1-mdm.v1$fit;
mf$A10_[mf$gp_=="I"]=mf$A20_[mf$gp_=="I"]=0
}
}
## (4) gp="U": baseline observations
##mf$A10_[mf$gp_=="U"]=0; mf$A20_[mf$gp_=="U"]=0; mf$A11_[mf$gp_=="U"]=0
##mf$obsprob_[mf$gp_=="U"]=mf$A10_[mf$gp_=="U"]+ mf$A20_[mf$gp_=="U"]+ mf$A11_[mf$gp_=="U"];
mf.vU = mf[mf$gp_=="U",,drop=F]
s.vU= s[mf$gp_=="U",,drop=F]
if (! any(mf.vU$g_=="O")) { ## If the set of g values is all on "I", "D" or "I, D" in this subsample.Then the probability of being observed is zero.
mf$obsprob_[mf$gp_=="U"]=0; mf$A10_[mf$gp_=="U"]=0; mf$A20_[mf$gp_=="U"]=0; mf$A11_[mf$gp_=="U"]=0;
} else if (all(mf.vU$g_=="O")) { ## If all data are observed in this subsample
mf$obsprob_[mf$gp_=="U"]=1; mf$A10_[mf$gp_=="U"]=1; mf$A20_[mf$gp_=="U"]=0; mf$A11_[mf$gp_=="U"]=0;
} else if (any(mf.vU$g_!="O")) {
missS <- apply(s.vU, 2, FUN=function (u) any(is.na(u)))
s.vU$WTs_=mf.vU$WTs_
s.vU$g_=mf.vU$g_
if (length(unique(mf.vU$g_))>2) { ## g includes "O, I, D" in this subsample.
mdm.vU =multinom(g_ ~ .-1- WTs_, data=s.vU[,!missS], weights=WTs_, Hess=TRUE,maxit=500)
mdm.vU.fit=as.data.frame(mdm.vU$fit)
mf$obsprob_[mf$gp_=="U"]=mdm.vU.fit$O; mf$A10_[mf$gp_=="U"]=as.numeric(mf.vU$g_=="I")-mdm.vU.fit$I;
mf$A20_[mf$gp_=="U"]=as.numeric(mf.vU$g_=="D")-mdm.vU.fit$D; mf$A11_[mf$gp_=="U"]=0;
} else if (length(unique(mf.vU$g_))==2) { ## g includes "O, I", or "O,D" in this subsample.
mdm.vU =multinom(g_ ~ .-1- WTs_, data=droplevels(s.vU[,!missS]), weights=WTs_, Hess=TRUE,maxit=500)
## mdm.v0 =glm(g_ ~ .-1- WTs_, data=s.v0[,!missS], weights=WTs_, family=binomial(link=logit))
mf$obsprob_[mf$gp_=="U"]=1-mdm.vU$fit;
mf$A10_[mf$gp_=="U"]=as.numeric(mf.vU$g_=="I")*(1-mdm.vU$fit);
mf$A20_[mf$gp_=="U"]=as.numeric(mf.vU$g_=="D")*(1-mdm.vU$fit);
mf$A11_[mf$gp_=="U"]=0;
}
} ##else print("Warning: check the values of g_ for the subsample gp=U")
return(mf)
}
fun.trans2pos <- function(u) {
res<-exp(u)
res
}
fun.trans2rho <- function(u) {
tol=10^308
if (u==Inf) u= tol
if (u== -Inf) u= -tol
rho<- 2*exp(u)/(1+exp(u))-1
rho
}
fun.postrans <- function(u) {
if (u) stop("Argument has to be positive")
res <- log(u)
res
}
fun.rhotrans <- function(rho) {
if (abs(rho) >1) stop("correlation coefficient must be between 1 and -1")
tol= 1e-10
u = (1+rho)/2
if (u==0) u= u + tol
if (u==1) u= 1- tol
res<- log(u/(1-u))
res
}
Try the isni 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.