# R/utilite.R In jmdl: Joint Mean-Correlation Regression Approach for Discrete Longitudinal Data

#### Documented in angle.plotgetJMDLgetJMDL.jmdlModlrt.testun.angle.plot

#' @title Extract or Get Generalized Components from a Fitted Joint Mean
#' Correlation Model
#'
#' @description Extract (or "get") "components" - in a generalized sense - from
#' a fitted joint mean correlation model from an object of class "JmdlMod".
#'
#' @param object a fitted joint mean correlation model of class "JmdlMod", i.e.,
#' typically the result of jmdl().
#' @param name a character vector specifying the name(s) of the "component".
#'
#' possible values are:
#' \describe{
#'   \item{\code{"m"}}{a vector of number of measurement for each subject}
#'   \item{\code{"Y"}}{response matrix}
#'   \item{\code{"X"}}{model matrix for mean structure}
#'   \item{\code{"W"}}{model matrix for correlation structure (the lower
#'   triangular matrix)}
#'   \item{\code{"offset"}}{a vecter to be added to a linear predictor}
#'   \item{\code{"theta"}}{parameter estimates of joint mean correlation model}
#'   \item{\code{"beta"}}{parameter estimates for mean structure model}
#'   \item{\code{"delta"}}{parameter estimates for mean structure model (for
#'   Nbinom model)}
#'   \item{\code{"gamma"}}{parameter estimates for correlation structure (the
#'   lower triangular matrix)}
#'   \item{\code{"stdbeta"}}{standard error for parameter beta}
#'   \item{\code{"stddelta"}}{standard error for parameter delta}
#'   \item{\code{"stdgamma"}}{standard error for parameter gamma}
#'   \item{\code{"loglik"}}{log-likelihood, except for a constant}
#'   \item{\code{"family"}}{the marginal distributions of the discrete variables}
#'   \item{\code{"q"}}{degree of polynomial of the time lag to model the lower
#'   triangular matrix}
#'   \item{\code{"time"}}{a vector of time from the data}
#' }
#'
#'
#' @examples
#'
#' mydat <- toydata
#' fit <- jmdl(Y|id|time ~ X, data = mydat, q = 2, family ='Bernoulli')
#' beta <- getJMDL(fit, "beta")
#' beta
#' loglik  <- getJMDL(fit, "loglik")
#' loglik
#'
#' @export
getJMDL <- function(object, name) UseMethod("getJMDL")

#' @describeIn getJMDL Extract or Get Generalized Components from a Fitted Joint
#' Mean Correlation Model
#' @export
getJMDL.jmdlMod <- function(object,
name = c("m", "Y", "X", "W", "offset", "theta", "beta", "gamma", "delta",
"loglik", "family", "q", "time", "stdbeta", "stdgamma",
"stddelta"))
{
if(missing(name)) stop("'name' must not be missing")

opt     <- object@opt
args    <- object@args
devcomp <- object@devcomp
std     <- object@std
offset  <- object@offset

m = args$m Y = args$Y
X = args$X W = args$W
time = args$time theta = drop(opt$par)

if (devcomp$dims['Bernoulli']) family <- 'Bernoulli' if (devcomp$dims['Poisson'])   family <- 'Poisson'
if (devcomp$dims['Nbinom']) family <- 'Nbinom' switch(name, "m" = args$m,
"Y" = args$Y, "X" = args$X,
"W" = args$W, "offset" = offset, "time" = args$time,
"theta"  = drop(opt$theta), "beta" = drop(opt$beta),
"gamma"  = drop(opt$gamma), "delta" = drop(opt$delta),
"loglik" = opt$loglik, "q" = object$q,
"family" = family,
"stdbeta"  = std$stdbeta, "stdgamma" = std$stdgamma,
"stddelta" = std$stddelta) } ##-------------------------------------------------------------------------- lagseq <- function(time) { res <- NULL if(length(time) != 1) { for(i in 2:length(time)) { for(j in 1:(i-1)) res <- c(res, (time[i] - time[j])) } } res } #' @title Plot Fitted Results and Model Diagnostics #' #' @description plot (a) fitted angles; (b) fitted correlations vs time lag; #' (c) the empirical distribution function vs the fitted distribution function; #' (d) the empirical correlations vs the fitted correlations, when the descrete #' longitudinal dataset is balanced. #' #' @param object a fitted joint mean correlation model of class "JmdlMod", i.e., #' typically the result of jmdl(). #' @param time a vector of obeservation time points #' #' #' #' @export angle.plot <- function(object, time) { debug <- 0 op <- par(mfrow = c(2, 2)) opt <- object@opt gamma <- opt$gamma
beta   <- opt$beta lgma <- length(gamma) args <- object@args dims <- object@devcomp$dims
family = getJMDL(object,"family")

m <- args[["m"]]
Y <- args[["Y"]]
X <- args[["X"]]
W <- args[["W"]]
W1 <- W[1,,,]
offset <- getJMDL(object,"offset")

if (length(unique(m)) != 1)
stop("Unbalanced longitudinal dataset.")

if(is.null(offset)) offset <- c(rep(0,length(m)*m[1]))

DataMat <- Y

tlag  <- lagseq(time)
tslag <- seq(min(tlag), max(tlag), length.out = 100)

W.tslag <- NULL
for(i in 0:(lgma-1)) W.tslag <- cbind(W.tslag, tslag^i)

Wgma <- W.tslag %*% gamma

n <- length(m)
zscore<-matrix(0,n,m)

for(i in 1:n){
if(i==1){
xi<-X[1:m[1],,drop=F]
osi<-offset[1:m[1]]
}else{
xi<-X[(sum(m[1:(i-1)])+1):(sum(m[1:(i-1)])+m[i]),,drop=F]
osi<-offset[(sum(m[1:(i-1)])+1):(sum(m[1:(i-1)])+m[i])]
}
yi = Y[i,]

if(family == 'Nbinom'){
delta = opt$delta mui<-as.vector(exp(osi+xi%*%beta)) sizei<-delta probi<-sizei/(sizei+mui) zscore[i,]<-qnorm(pnbinom(yi,sizei,probi)) } if(family == 'Bernoulli'){ pij<-as.vector(boot::inv.logit(osi+xi%*%beta)) zscore[i,]<-qnorm(pbinom(yi,size=1,prob=pij)) } if(family == 'Poisson'){ phi.theta<-as.vector(exp(osi+xi%*%beta)) zscore[i,]<-qnorm(ppois(yi,phi.theta)) } } R <- cor(zscore) # first plot B <- t(chol(R)) PhiMat <- matrix(0, dim(B)[1], dim(B)[2]) for(j in 2:dim(B)[1]) { for(k in 1:(j-1)) { tmp <- 1 if (k != 1) { tmp <- prod(sin(PhiMat[j, 1:(k-1)])) } # if PhiMat[j,k] <- acos(B[j, k]/tmp) } # for k } # for j PhiMatt <- t(PhiMat) phi <- PhiMatt[upper.tri(PhiMatt, diag=FALSE)] plot(tlag, tan(pi/2-phi), xlab="Lag", ylab=expression(tan(pi/2-w)), pch=16,main="(a)") lines(tslag, Wgma) # second plot fit.Tmat<-fit.phimat<-diag(1, dim(B)[1], dim(B)[2]) for(j in 2:dim(B)[1]){ for(k in 1:(j-1)){ fit.phimat[j,k]<-pi/2-atan(W1[j,k,]%*%gamma) } fit.Tmat[j,j]<-prod(sin(fit.phimat[j,1:(j-1)])) fit.Tmat[j,1]<-cos(fit.phimat[j,1]) if(j>2){ for(l in 2:(j-1)) fit.Tmat[j,l]<-cos(fit.phimat[j,l])*prod(sin(fit.phimat[j,1:(l-1)])) } } fit.cormat<-fit.Tmat%*%t(fit.Tmat) cor <- fit.cormat[upper.tri(fit.cormat, diag=FALSE)] plot(tlag, cor, xlab="Lag", ylab="Correlation",pch=16,main="(b)") # third plot plot(R[lower.tri(cor(zscore))], fit.cormat[lower.tri(fit.cormat)], xlim=c(0.2,0.8),ylim=c(0.2,0.8), xlab="Sample-based Correlation", ylab="Fitted Correlation",pch=16,main="(c)") abline(0,1) # forth plot Ecdf = c(0) for(i in 1:n){ sum = 0 for(j in 1:n){ temp = prod(I(Y[j,]<Y[i,])) sum = sum + temp } Ecdf[i] = sum/n } f.fit <- c(0) for(i in 1:n){ f.fit[i] <- mvtnorm::pmvnorm(c(rep(-Inf, m[i])),zscore[i,],corr=fit.cormat) } plot(sort(f.fit),sort(Ecdf), xlab="Fitted CDF", ylab="Empirical CDF",pch=16,main="(d)") abline(0,1) } #' @title Plot Fitted Results #' #' @description plot (a) fitted angles; (b) fitted correlations vs time lag, when #' the descrete longitudinal dataset is unbalanced. #' #' @param object a fitted joint mean correlation model of class "JmdlMod", i.e., #' typically the result of jmdl(). #' #' #' @export un.angle.plot <- function(object) { debug <- 0 op <- par(mfrow = c(1, 2)) opt <- object@opt gamma <- opt$gamma
beta   <- opt$beta lgma <- length(gamma) args <- object@args dims <- object@devcomp$dims
family = getJMDL(object,"family")
stdgma = getJMDL(object,"stdgamma")

m <- args[["m"]]
Y <- args[["Y"]]
X <- args[["X"]]
W <- args[["W"]]

DataMat <- Y

time <- args[["time"]]

n <- length(m)
ncor <- NULL
ntlag <- NULL

# first plot
for(i in 1:n){
tlag = NULL
cor = NULL
fit.Tmat<-fit.phimat<-diag(1, m[i], m[i])
if(m[i]>1){
for(j in 2:m[i]){
for(k in 1:(j-1)){
fit.phimat[j,k]<-pi/2-atan(W[i,j,k,]%*%gamma)
}
fit.Tmat[j,j]<-prod(sin(fit.phimat[j,1:(j-1)]))
fit.Tmat[j,1]<-cos(fit.phimat[j,1])
if(j>2){
for(l in 2:(j-1))
fit.Tmat[j,l]<-cos(fit.phimat[j,l])*prod(sin(fit.phimat[j,1:(l-1)]))
}
}
if(i==1){
timei<-time[1:m[1],drop=F]
}else{
timei<-time[(sum(m[1:(i-1)])+1):(sum(m[1:(i-1)])+m[i]),drop=F]
}
tlag <- lagseq(timei)
}
fit.cormat <- fit.Tmat%*%t(fit.Tmat)
cor <- fit.cormat[upper.tri(fit.cormat, diag=FALSE)]
ncor <- c(ncor, cor)
ntlag <- c(ntlag, tlag)
}
plot(ntlag, ncor, xlab="Lag", ylab="Correlation",pch=20,cex=0.5,main="(a)")

# second plot
tslag <- seq(0, max(time) - min(time), length.out = n)
W.tslag <- NULL

for(i in 0:(lgma-1)) W.tslag <- cbind(W.tslag, tslag^i)
Wgma <- drop(W.tslag %*% gamma)

ft.sd<-(W.tslag%*%stdgma)

plot(tslag, Wgma, type = 'l', xlab = "Lag", ylab = expression(tan(pi/2-w)), ylim=c(-1,2),main="(b)")
lines(tslag, Wgma+1.96*ft.sd, lty = 2, lwd = 2)
lines(tslag, Wgma-1.96*ft.sd, lty = 2, lwd = 2)
}

#' @title Pairwise Likelihood Ratio Statistic Test
#'
#' @description Conducts a pairwise likelihood ratio test for joint mean-correlation regression.
#'
#' @param fit a fitted joint mean correlation model of class "JmdlMod", i.e.,
#' typically the result of jmdl().
#' @param id the id of paremeter to test
#'
#'
#' @export
lrt.test <- function(fit, id){
#theta0 = 0
m <- getJMDL(fit, "m")
Y <- getJMDL(fit, "Y")
X <- getJMDL(fit, "X")
W <- getJMDL(fit, "W")
family <- getJMDL(fit, "family")
offset <- getJMDL(fit, "offset")
lrt.full <- fit@opt
id <- c(id)
theta0 <- c(rep(0,length(id)))

lgma<-dim(W)[4]
Yv<-na.exclude(as.vector(t(Y)))

##Bernoulli
if(family == 'Bernoulli'){
#null comlik estimation
y.fit0<-glm(Yv~offset(offset)+X[,-id]-1,family='binomial')
theta<-c(coef(y.fit0),rep(0,lgma))
lrt0<-minqa::bobyqa(par=theta,fn=cond.comlik, theta0=theta0,id=id,Y=Y,X=X,W=W,m=m,
family=family, offset=offset)

lm.obj <- glm(Yv ~ offset(offset) + X-1 ,family='binomial')
bta1 <- coef(lm.obj)
gma1 <- rep(0, lgma)
theta1 <- c(bta1, gma1)
lrt1<-minqa::bobyqa(par=theta1,fn=comlik,m=m,Y=Y,X=X,W=W,family=family,offset=offset)
}

##Poisson
if(family == 'Poisson'){
#null comlik estimation
y.fit0<-glm(Yv~offset(offset)+X[,-id]-1,family='poisson')
theta<-c(coef(y.fit0),rep(0,lgma))
lrt0<-minqa::bobyqa(par=theta,fn=cond.comlik, theta0=theta0, id=id,Y=Y,X=X,W=W,m=m,
family=family, offset=offset)

lm.obj <- glm(Yv ~ offset(offset) + X-1 ,family='poisson')
bta1 <- coef(lm.obj)
gma1 <- rep(0, lgma)
theta1 <- c(bta1, gma1)
lrt1<-minqa::bobyqa(par=theta1,fn=comlik,m=m,Y=Y,X=X,W=W,family=family,offset=offset)
}

##Nbinom
if(family == 'Nbinom'){
lbta = ncol(X)
if(sum(id==(lbta+1))) return(list(lrt=NA, pval=NA, critval=NA))
#null comlik estimation
y.fit0<-MASS::glm.nb(Yv~offset(offset)+X[,-id]-1)
theta<-c(coef(y.fit0),y.fit0$theta,rep(0,lgma)) lrt0<-minqa::bobyqa(par=theta,fn=cond.comlik, theta0=theta0, id=id,Y=Y,X=X,W=W,m=m, family=family,offset=offset) lm.obj <- MASS::glm.nb(Yv ~ offset(offset) + X-1) bta1 <- coef(lm.obj) gma1 <- rep(0, lgma) theta1 <- c(bta1, lm.obj$theta, gma1)
lrt1<-minqa::bobyqa(par=theta1,fn=comlik,m=m,Y=Y,X=X,W=W,family=family,offset=offset)
}

lrt<-2*(lrt0$fval-lrt1$fval)

theta<-rep(0,length(lrt1$par)) theta[id]<-theta0 theta[-id]<-lrt1$par[-id]

#pvalue
#set.seed(1234)
cov.fit0<-asycov(m=m,Y=Y,X=X,W=W,family=family,theta=theta,offset=offset)
chimat<-matrix(rchisq(length(id)*10000, df=1),ncol=length(id))
if(length(id)==1){
lmds<-cov.fit0$Gmat.inv[id,id]/cov.fit0$Kmat.inv[id,id]
a<-lmds*chimat
}else{
lmds<-eigen(solve(cov.fit0$Kmat.inv[id,id])%*%cov.fit0$Gmat.inv[id,id])
a<-chimat%*%lmds\$values
}
a<-as.vector(a)
pval<-sum(a>lrt)/10000
return(list(lrt.statistic=lrt, pval=pval, critval=round(quantile(a,c(0.025, 0.975)), 4))) #up 95% quantile
}


## Try the jmdl package in your browser

Any scripts or data that you put into this service are public.

jmdl documentation built on May 2, 2019, 11:04 a.m.