Nothing
R/covest.R
In JointModel: Semiparametric Joint Models for Longitudinal and Counting Processes
Defines functions
covest
covest <- function(par_Y, par_T, par_b, K, dl1, dl2, npar, GQl){
if (K==0){
H <- function(x) (x)
Hinv <- function(x) {x}
H1d <- function(x) 1
H2d <- function(x) 0
H3d <- function(x) 0
} else if (K > 0){
H <- function(x) {log(1+K*x)/K}
Hinv <- function(x) {(exp(K*x)-1)/K}
H1d <- function(x) {1/(1+K*x)}
H2d <- function(x) {-K/(1+K*x)^2}
H3d <- function(x) {2*K^2/(1+K*x)^3}
}
N=npar$N; nb=length(par_b);
thetay=par_Y[-npar$nY]; sigma2e= par_Y[npar$nY];
gamma=par_T[1:npar$nG]; phi = par_T[npar$nG+(1:npar$nP)];
lambda= par_T[(npar$nG+npar$nP)+1:npar$mt];
r1= par_b[1]; r2= par_b[2]; r3= par_b[3];
bbv <- Tbb(par_Y,par_T,par_b, K, dl1,dl2, npar, GQl)
b1= bbv[1:N, ]; b2=bbv[N+(1:N), ];
wwN <- GetWeight(bbv, par_Y, par_T, par_b, K, dl1, dl2, npar, GQl)
YAXX1B <- c(dl1$yt- dl1$Xy%*%thetay)
L_Vi <- c(dl2$ind_L %*% lambda)
exptp <- exp(c(dl2$X2%*%gamma) + phi*b1)
expL <- exptp*L_Vi
H1=H1d(expL); H2=H2d(expL); H3=H3d(expL);
gH2H1 = H2/H1*exptp
gdH2H1= H3/H1*exptp*exptp-(H2/H1*exptp)^2
gH1 = H1*exptp
gH2 = H2*exptp*exptp
Eb1 <- rowSums(b1*wwN)
Eb2 <- rowSums(b2*wwN)
Eb1sq <- rowSums(b1^2*wwN)
Eb2sq <- rowSums(b2^2*wwN)
Eb1b2 <- rowSums(b1*b2*wwN)
EH2H1 <- rowSums(gH2H1*wwN)
EdH2H1 <- rowSums(gdH2H1*wwN)
EH1 <- rowSums(gH1*wwN)
EH2 <- rowSums(gH2*wwN)
EH1b1 <- rowSums(gH1*b1*wwN)
EH2b1<- rowSums(gH2*b1*wwN)
EH1b1sq <- rowSums(gH1*b1^2*wwN)
EH2b1sq <- rowSums(gH2*b1^2*wwN)
EH2H1b1 <- rowSums(gH2H1*b1*wwN)
EdH2H1b1 <- rowSums(gdH2H1*b1*wwN)
EH2H1b1sq <- rowSums(gH2H1*b1^2*wwN)
EdH2H1b1sq <- rowSums(gdH2H1*b1^2*wwN)
ED2LY <- matrix(0, npar$nY, npar$nY)
ED2LT <- matrix(0, npar$nTl, npar$nTl)
ED2Lb <- matrix(0, nb, nb)
E_DLDL <- matrix(0, (npar$nY+npar$nTl+nb), (npar$nY+npar$nTl+nb))
EDLEDL <- matrix(0, (npar$nY+npar$nTl+nb), (npar$nY+npar$nTl+nb))
for (ii in 1:N){
IDXY <- which(dl1$id.vec.y==dl2$id.vec[ii])
D1thetay <- matrix(0, length(thetay), GQl$Nnode)
D1sigma2e <- matrix(0, 1, GQl$Nnode)
D2thetay <- matrix(0, length(thetay), length(thetay))
D2sigma2e <- 0
D2thetaysigma2e <-matrix(0, length(thetay), 1)
for (k in IDXY){
### Y
D1thetay=D1thetay + dl1$Xy[k, ]%*%t((YAXX1B[k]-dl1$ZZ1[k,1]*b1[ii,]-dl1$ZZ1[k,2]*b2[ii,])/sigma2e)
D1sigma2e=D1sigma2e + t((YAXX1B[k]-dl1$ZZ1[k,1]*b1[ii,]-dl1$ZZ1[k,2]*b2[ii,])^2/(2*sigma2e^2)-1/(2*sigma2e))
D2thetay <- D2thetay -dl1$Xy[k, ]%*%t(dl1$Xy[k, ])/sigma2e
D2thetaysigma2e=D2thetaysigma2e -(YAXX1B[k]-dl1$ZZ1[k,1]*Eb1[ii]-dl1$ZZ1[k,2]*Eb2[ii])/sigma2e^2*dl1$Xy[k, ]
temp=(YAXX1B[k]-dl1$ZZ1[k,1]*b1[ii,]-dl1$ZZ1[k,2]*b2[ii,])^2
D2sigma2e=D2sigma2e + 1/2/sigma2e^2-sum(temp*wwN[ii,])/sigma2e^3
}
RISK <- which(dl2$wt <= dl2$EventTime[ii])
indX <- ifelse(dl2$wt <= dl2$EventTime[ii], 1, 0)
coef2 <- dl2$DELTA[ii]*gH2H1[ii, ]-gH1[ii, ]
coef4 <- dl2$DELTA[ii]*(EH2H1[ii]+EdH2H1[ii]*L_Vi[ii])-(EH1[ii]+EH2[ii]*L_Vi[ii])
coef5b1=dl2$DELTA[ii]*(EH2H1b1[ii]+EdH2H1b1[ii]*L_Vi[ii])-(EH1b1[ii]+EH2b1[ii]*L_Vi[ii])
D1G <- dl2$DELTA[ii]*dl2$X2[ii, ] + dl2$X2[ii, ]%*%t(coef2*L_Vi[ii])
D1phi1 <- t(dl2$DELTA[ii]*b1[ii,] + coef2*b1[ii,]*L_Vi[ii])
D1lambda <- dl2$DELTA[ii]*ifelse(dl2$wt==dl2$EventTime[ii],1,0)/lambda + indX%*%t(coef2)
D2G <- coef4*L_Vi[ii]*dl2$X2[ii, ]%*%t(dl2$X2[ii, ])
D2Gphi1 <- matrix(coef5b1*L_Vi[ii]*dl2$X2[ii, ])
D2Glambda <- matrix(0, npar$nG, npar$mt)
D2Glambda[ ,RISK] <- coef4*dl2$X2[ii, ]
D2phi1 <- -(EH1b1sq[ii]+EH2b1sq[ii]*L_Vi[ii])*L_Vi[ii]
if (dl2$DELTA[ii]==1) {D2phi1 =D2phi1+(EH2H1b1sq[ii]+EdH2H1b1sq[ii]*L_Vi[ii])*L_Vi[ii]}
D2phi1lambda <- matrix(0, 1, npar$mt)
D2phi1lambda[ ,RISK] <- coef5b1
D2lambda <- indX%*%t(indX)*(dl2$DELTA[ii]*EdH2H1[ii]-EH2[ii])
if (dl2$DELTA[ii]==1) {D2lambda = D2lambda-diag(npar$mt)*ifelse(dl2$wt==dl2$EventTime[ii],1,0)/lambda^2}
# random effect
cb1 <- t((b1[ii,]/r1)^2)
cb2 <- t((b2[ii,]/r2)^2)
cb12 <- t(r3*b1[ii,]*b2[ii,]/(r1*r2)^0.5)
Ecb1 <- Eb1sq[ii]/r1^2
Ecb2 <- Eb2sq[ii]/r2^2
Ecb1b2 <- r3*Eb1b2[ii]/(r1*r2)^0.5
Dr1 <- -1/2/r1 -1/2/(1-r3^2)*(-cb1 + cb12/r1)
Dr2 <- -1/2/r2 -1/2/(1-r3^2)*(-cb2 + cb12/r2)
Dr3 <- (r3-r3/(1-r3^2)*(r1*cb1-2*cb12+r2*cb2)+cb12/r3)/(1-r3^2)
D2r1<- 1/2/r1^2 - (2*Ecb1/r1-1.5*Ecb1b2/r1^2)/(2*(1-r3^2))
D2r2<- 1/2/r2^2 - (2*Ecb2/r2-1.5*Ecb1b2/r2^2)/(2*(1-r3^2))
D2r3<- ((1+r3^2) - (1+3*r3^2)/(1-r3^2)*(Ecb1*r1-2*Ecb1b2+Ecb2*r2) + 4*Ecb1b2)/(1-r3^2)^2
D2r1r2<- Ecb1b2/(r1*r2*4*(1-r3^2))
D2r1r3<- r3/(1-r3^2)^2*(Ecb1-Ecb1b2/r1)-Ecb1b2/r3/(2*(1-r3^2)*r1)
D2r2r3<- r3/(1-r3^2)^2*(Ecb2-Ecb1b2/r2)-Ecb1b2/r3/(2*(1-r3^2)*r2)
ED2LY <- ED2LY + rbind(cbind(D2thetay, D2thetaysigma2e), cbind(t(D2thetaysigma2e), D2sigma2e))
ED2LT <- ED2LT + rbind(cbind(D2G, D2Gphi1, D2Glambda),
cbind(t(D2Gphi1), D2phi1, D2phi1lambda),
cbind(t(D2Glambda), t(D2phi1lambda), D2lambda))
ED2Lb <- ED2Lb + cbind(c(D2r1,D2r1r2,D2r1r3), c(D2r1r2,D2r2,D2r2r3), c(D2r1r3,D2r2r3,D2r3))
DL <- rbind(D1thetay,D1sigma2e, D1G,D1phi1,D1lambda, Dr1,Dr2,Dr3)
EDL <- t(DL)*wwN[ii,]
E_DLDL <- E_DLDL + DL%*%(EDL)
EDL <- colSums(EDL)
EDLEDL <- EDLEDL + EDL%*%t(EDL)
}
ED2L <- rbind(cbind(ED2LY, matrix(0,npar$nY,npar$nTl), matrix(0,npar$nY,nb)),
cbind(matrix(0,npar$nTl,npar$nY), ED2LT, matrix(0,npar$nTl,nb)),
cbind(matrix(0,nb,npar$nY), matrix(0,nb,npar$nTl), ED2Lb))
HS <- -ED2L - (E_DLDL-EDLEDL)
COV <- solve(HS); COV;
errCOV=0;
if (length(which(diag(COV)<0))>0){
out.eigen<-eigen(HS)
values<-out.eigen$values
temp1<-which(values < 0)
if (length(temp1) > 0){
values[-temp1]<-1/values[-temp1]
values[temp1]<-0
errCOV=1;
} else {values <- 1/values}
COV <- (out.eigen$vectors) %*% diag(values) %*% t(out.eigen$vectors);
}
outCOV <- list(COV=COV, errCOV=errCOV); outCOV;
}
Try the JointModel package in your browser
Any scripts or data that you put into this service are public.
JointModel documentation built on May 2, 2019, 12:40 p.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 covest
covest <- function(par_Y, par_T, par_b, K, dl1, dl2, npar, GQl){
if (K==0){
H <- function(x) (x)
Hinv <- function(x) {x}
H1d <- function(x) 1
H2d <- function(x) 0
H3d <- function(x) 0
} else if (K > 0){
H <- function(x) {log(1+K*x)/K}
Hinv <- function(x) {(exp(K*x)-1)/K}
H1d <- function(x) {1/(1+K*x)}
H2d <- function(x) {-K/(1+K*x)^2}
H3d <- function(x) {2*K^2/(1+K*x)^3}
}
N=npar$N; nb=length(par_b);
thetay=par_Y[-npar$nY]; sigma2e= par_Y[npar$nY];
gamma=par_T[1:npar$nG]; phi = par_T[npar$nG+(1:npar$nP)];
lambda= par_T[(npar$nG+npar$nP)+1:npar$mt];
r1= par_b[1]; r2= par_b[2]; r3= par_b[3];
bbv <- Tbb(par_Y,par_T,par_b, K, dl1,dl2, npar, GQl)
b1= bbv[1:N, ]; b2=bbv[N+(1:N), ];
wwN <- GetWeight(bbv, par_Y, par_T, par_b, K, dl1, dl2, npar, GQl)
YAXX1B <- c(dl1$yt- dl1$Xy%*%thetay)
L_Vi <- c(dl2$ind_L %*% lambda)
exptp <- exp(c(dl2$X2%*%gamma) + phi*b1)
expL <- exptp*L_Vi
H1=H1d(expL); H2=H2d(expL); H3=H3d(expL);
gH2H1 = H2/H1*exptp
gdH2H1= H3/H1*exptp*exptp-(H2/H1*exptp)^2
gH1 = H1*exptp
gH2 = H2*exptp*exptp
Eb1 <- rowSums(b1*wwN)
Eb2 <- rowSums(b2*wwN)
Eb1sq <- rowSums(b1^2*wwN)
Eb2sq <- rowSums(b2^2*wwN)
Eb1b2 <- rowSums(b1*b2*wwN)
EH2H1 <- rowSums(gH2H1*wwN)
EdH2H1 <- rowSums(gdH2H1*wwN)
EH1 <- rowSums(gH1*wwN)
EH2 <- rowSums(gH2*wwN)
EH1b1 <- rowSums(gH1*b1*wwN)
EH2b1<- rowSums(gH2*b1*wwN)
EH1b1sq <- rowSums(gH1*b1^2*wwN)
EH2b1sq <- rowSums(gH2*b1^2*wwN)
EH2H1b1 <- rowSums(gH2H1*b1*wwN)
EdH2H1b1 <- rowSums(gdH2H1*b1*wwN)
EH2H1b1sq <- rowSums(gH2H1*b1^2*wwN)
EdH2H1b1sq <- rowSums(gdH2H1*b1^2*wwN)
ED2LY <- matrix(0, npar$nY, npar$nY)
ED2LT <- matrix(0, npar$nTl, npar$nTl)
ED2Lb <- matrix(0, nb, nb)
E_DLDL <- matrix(0, (npar$nY+npar$nTl+nb), (npar$nY+npar$nTl+nb))
EDLEDL <- matrix(0, (npar$nY+npar$nTl+nb), (npar$nY+npar$nTl+nb))
for (ii in 1:N){
IDXY <- which(dl1$id.vec.y==dl2$id.vec[ii])
D1thetay <- matrix(0, length(thetay), GQl$Nnode)
D1sigma2e <- matrix(0, 1, GQl$Nnode)
D2thetay <- matrix(0, length(thetay), length(thetay))
D2sigma2e <- 0
D2thetaysigma2e <-matrix(0, length(thetay), 1)
for (k in IDXY){
### Y
D1thetay=D1thetay + dl1$Xy[k, ]%*%t((YAXX1B[k]-dl1$ZZ1[k,1]*b1[ii,]-dl1$ZZ1[k,2]*b2[ii,])/sigma2e)
D1sigma2e=D1sigma2e + t((YAXX1B[k]-dl1$ZZ1[k,1]*b1[ii,]-dl1$ZZ1[k,2]*b2[ii,])^2/(2*sigma2e^2)-1/(2*sigma2e))
D2thetay <- D2thetay -dl1$Xy[k, ]%*%t(dl1$Xy[k, ])/sigma2e
D2thetaysigma2e=D2thetaysigma2e -(YAXX1B[k]-dl1$ZZ1[k,1]*Eb1[ii]-dl1$ZZ1[k,2]*Eb2[ii])/sigma2e^2*dl1$Xy[k, ]
temp=(YAXX1B[k]-dl1$ZZ1[k,1]*b1[ii,]-dl1$ZZ1[k,2]*b2[ii,])^2
D2sigma2e=D2sigma2e + 1/2/sigma2e^2-sum(temp*wwN[ii,])/sigma2e^3
}
RISK <- which(dl2$wt <= dl2$EventTime[ii])
indX <- ifelse(dl2$wt <= dl2$EventTime[ii], 1, 0)
coef2 <- dl2$DELTA[ii]*gH2H1[ii, ]-gH1[ii, ]
coef4 <- dl2$DELTA[ii]*(EH2H1[ii]+EdH2H1[ii]*L_Vi[ii])-(EH1[ii]+EH2[ii]*L_Vi[ii])
coef5b1=dl2$DELTA[ii]*(EH2H1b1[ii]+EdH2H1b1[ii]*L_Vi[ii])-(EH1b1[ii]+EH2b1[ii]*L_Vi[ii])
D1G <- dl2$DELTA[ii]*dl2$X2[ii, ] + dl2$X2[ii, ]%*%t(coef2*L_Vi[ii])
D1phi1 <- t(dl2$DELTA[ii]*b1[ii,] + coef2*b1[ii,]*L_Vi[ii])
D1lambda <- dl2$DELTA[ii]*ifelse(dl2$wt==dl2$EventTime[ii],1,0)/lambda + indX%*%t(coef2)
D2G <- coef4*L_Vi[ii]*dl2$X2[ii, ]%*%t(dl2$X2[ii, ])
D2Gphi1 <- matrix(coef5b1*L_Vi[ii]*dl2$X2[ii, ])
D2Glambda <- matrix(0, npar$nG, npar$mt)
D2Glambda[ ,RISK] <- coef4*dl2$X2[ii, ]
D2phi1 <- -(EH1b1sq[ii]+EH2b1sq[ii]*L_Vi[ii])*L_Vi[ii]
if (dl2$DELTA[ii]==1) {D2phi1 =D2phi1+(EH2H1b1sq[ii]+EdH2H1b1sq[ii]*L_Vi[ii])*L_Vi[ii]}
D2phi1lambda <- matrix(0, 1, npar$mt)
D2phi1lambda[ ,RISK] <- coef5b1
D2lambda <- indX%*%t(indX)*(dl2$DELTA[ii]*EdH2H1[ii]-EH2[ii])
if (dl2$DELTA[ii]==1) {D2lambda = D2lambda-diag(npar$mt)*ifelse(dl2$wt==dl2$EventTime[ii],1,0)/lambda^2}
# random effect
cb1 <- t((b1[ii,]/r1)^2)
cb2 <- t((b2[ii,]/r2)^2)
cb12 <- t(r3*b1[ii,]*b2[ii,]/(r1*r2)^0.5)
Ecb1 <- Eb1sq[ii]/r1^2
Ecb2 <- Eb2sq[ii]/r2^2
Ecb1b2 <- r3*Eb1b2[ii]/(r1*r2)^0.5
Dr1 <- -1/2/r1 -1/2/(1-r3^2)*(-cb1 + cb12/r1)
Dr2 <- -1/2/r2 -1/2/(1-r3^2)*(-cb2 + cb12/r2)
Dr3 <- (r3-r3/(1-r3^2)*(r1*cb1-2*cb12+r2*cb2)+cb12/r3)/(1-r3^2)
D2r1<- 1/2/r1^2 - (2*Ecb1/r1-1.5*Ecb1b2/r1^2)/(2*(1-r3^2))
D2r2<- 1/2/r2^2 - (2*Ecb2/r2-1.5*Ecb1b2/r2^2)/(2*(1-r3^2))
D2r3<- ((1+r3^2) - (1+3*r3^2)/(1-r3^2)*(Ecb1*r1-2*Ecb1b2+Ecb2*r2) + 4*Ecb1b2)/(1-r3^2)^2
D2r1r2<- Ecb1b2/(r1*r2*4*(1-r3^2))
D2r1r3<- r3/(1-r3^2)^2*(Ecb1-Ecb1b2/r1)-Ecb1b2/r3/(2*(1-r3^2)*r1)
D2r2r3<- r3/(1-r3^2)^2*(Ecb2-Ecb1b2/r2)-Ecb1b2/r3/(2*(1-r3^2)*r2)
ED2LY <- ED2LY + rbind(cbind(D2thetay, D2thetaysigma2e), cbind(t(D2thetaysigma2e), D2sigma2e))
ED2LT <- ED2LT + rbind(cbind(D2G, D2Gphi1, D2Glambda),
cbind(t(D2Gphi1), D2phi1, D2phi1lambda),
cbind(t(D2Glambda), t(D2phi1lambda), D2lambda))
ED2Lb <- ED2Lb + cbind(c(D2r1,D2r1r2,D2r1r3), c(D2r1r2,D2r2,D2r2r3), c(D2r1r3,D2r2r3,D2r3))
DL <- rbind(D1thetay,D1sigma2e, D1G,D1phi1,D1lambda, Dr1,Dr2,Dr3)
EDL <- t(DL)*wwN[ii,]
E_DLDL <- E_DLDL + DL%*%(EDL)
EDL <- colSums(EDL)
EDLEDL <- EDLEDL + EDL%*%t(EDL)
}
ED2L <- rbind(cbind(ED2LY, matrix(0,npar$nY,npar$nTl), matrix(0,npar$nY,nb)),
cbind(matrix(0,npar$nTl,npar$nY), ED2LT, matrix(0,npar$nTl,nb)),
cbind(matrix(0,nb,npar$nY), matrix(0,nb,npar$nTl), ED2Lb))
HS <- -ED2L - (E_DLDL-EDLEDL)
COV <- solve(HS); COV;
errCOV=0;
if (length(which(diag(COV)<0))>0){
out.eigen<-eigen(HS)
values<-out.eigen$values
temp1<-which(values < 0)
if (length(temp1) > 0){
values[-temp1]<-1/values[-temp1]
values[temp1]<-0
errCOV=1;
} else {values <- 1/values}
COV <- (out.eigen$vectors) %*% diag(values) %*% t(out.eigen$vectors);
}
outCOV <- list(COV=COV, errCOV=errCOV); outCOV;
}
Try the JointModel 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.