Nothing
QICW.gee <-
function(object) {
#########################################################################
# Arguments: mu.R weight, scale and rho are extracted from macro GEE
# data data frame with subject variable from 1:size
# model specify the model of interest
# logit_model specify the model for dropout
# corstr Working correlation structure: "independence", "AR-M", "exchangeable", "unstructured".
# family type of the outcomes: gaussian, binomial or poisson
#### continous case: x*beta;
#### binomial case: 1/exp(-x*beta)
#### poisson case: exp(x*beta)
#########################################################################
# input from wgee()
# para_est the estimates of the parameters
# alpha_drop the parameter estimates for the dropout model
# scale the estimate of variance for gaussian outcomes only based on weighted GEE
# rho the estimate of correlation coefficient based on weighted GEE
# mu.R the fitted values from the candidate model based on weighted GEE
# weight the esimtate of weights based on weighted GEE
#####weight######
# browser()
model=object$model
data=object$data
id=object$id
corstr=object$corstr
family=object$family
if (length(id) != nrow(data)){
stop("variable lengths differ (found for '(id)')")}
m.temp <- model.frame(model, data, na.action='na.pass')
if (sum(is.na(m.temp)) == 0){
stop("This data is complete. QICW.gee() is only for data with missingness")}
cluster <-cluster.size(id)$n #number of observations for each subject;
size <-cluster.size(id)$m # sample size;
subject= rep(1:size, cluster)
data$subject <- subject #1:size; create subject variable#
data$mu.R=object$mu_fit
data$weight=object$weight
model=object$model
model_drop=object$mis_fit$formula
para_est=as.vector(object$beta)
alpha_drop=as.vector(coef(object$mis_fit))
scale=object$scale
rho=object$rho
init <- model.frame(model, data)
init$num <- 1:length(init[,1])
### Get the design matrix;
m <- model.frame(model, data, na.action='na.pass')
data$response <- model.response(m, "numeric")
mat <- as.data.frame(model.matrix(model, m))
mat$subj <- rep(unique(data$subject), cluster)
beta <- para_est
px <- length(beta) # number non-redunant columns in design matrix
# Quasi Likelihood;
index.nmiss <- which(!is.na(data$response))
switch(family,
gaussian={model.R <- geeglm(model, id=subject, data=data, corstr=corstr, family=gaussian)
y <- model.R$y
quasi.R <- -1/2*sum((y-data$mu.R[index.nmiss])^2*data$weight[index.nmiss])
},
binomial={model.R <- geeglm(model, id=subject, data=data, corstr=corstr, family=binomial)
y <- model.R$y
quasi.R <- sum((y*log(data$mu.R[index.nmiss]/(1-data$mu.R[index.nmiss]))+log(1-data$mu.R[index.nmiss]))*data$weight[index.nmiss])
},
poisson={model.R <- geeglm(model, id=subject, data=data, corstr=corstr, family=poisson)
y <- model.R$y
quasi.R <- sum((y*log(data$mu.R[index.nmiss])-data$mu.R[index.nmiss])*data$weight[index.nmiss])
},
stop("Warnings: Invalid type of outcomes!")
)
##### get the trace term: the penalty for model complexity;
h.part<-matrix(0, nrow=length(mat[1,])-1, ncol=length(mat[1,])-1)
UF<-matrix(0,ncol=length(alpha_drop), nrow=length(mat[1,])-1)
FF<-matrix(0,ncol=length(alpha_drop), nrow=length(alpha_drop))
j.list <- k.list <- list()
for (i in 1:size){
cluster_alt <- cluster.size(data[!is.na(data$response),]$subject)$n
ncluster <- cluster_alt[i]
y<-as.matrix(data[data$subject==i,]$response)
mu.R<-as.matrix(data[data$subject==i,]$mu.R)
covariate<-as.matrix(subset(mat[,-length(mat[1,])], mat$subj==unique(data$subject)[i])) ### specify the covariate matrix
switch(family,
gaussian={
index <- which(!is.na(data[data$subject==i,]$response))
var <- switch(corstr,
independence=scale*cormax.ind(ncluster),
exchangeable=scale*cormax.exch(ncluster, rho),
ar1=scale*cormax.ar1(ncluster, rho)
)
if (is.matrix(covariate[index,])==FALSE) {cal.m=t(t(covariate[index,]))
}else{cal.m=t(covariate[index,])}
h.ind<-cal.m%*%ginv(var)%*%diag(c(data[data$subject==i,]$weight[index]))%*%covariate[index,]
#h.ind<-t(covariate[index,])%*%ginv(var)%*%covariate[index,]
h.part<-h.part+h.ind
###Calculate the G_i in the formula;
x<-diag(c(data[data$subject==i,]$weight[index]))%*%matrix(y[index]-mu.R[index])
U<-cal.m%*%ginv(var)%*%x
### Get the design matrix for dropout model;
m_drop <- model.frame(model_drop, data[data$subject==i,],na.action='na.pass')
mat_drop <- as.data.frame(model.matrix(model_drop, m_drop))
exp.x <- exp(as.matrix(mat_drop)%*%as.matrix(alpha_drop))
lam_d <- exp.x/(1+exp.x)
row.names(m_drop)=seq(1:cluster[i])
index_d=as.numeric(row.names(m_drop))[-1]
index_d2=seq(1:cluster[i])[-length(seq(1:cluster[i]))]
F= as.matrix(colSums( na.omit( (m_drop[index_d,][,1]* mat_drop[index_d,]- c(lam_d[index_d])*mat_drop[index_d,] )*
mat_drop[index_d2,]) ) )
UF<-UF+U%*%t(F)
FF<-FF+F%*%t(F)
},
binomial={
index <- which(!is.na(data[data$subject==i,]$response))
var <- switch(corstr,
independence=scale*diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster)%*%cormax.ind(ncluster)%*%diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster),
exchangeable=scale*diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster)%*%cormax.exch(ncluster, rho)%*%diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster),
ar1=scale*diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster)%*%cormax.ar1(ncluster, rho)%*%diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster)
)
D<-mat.prod(covariate[index,], exp(covariate[index,]%*%beta)/((1+exp(covariate[index,]%*%beta))^2))
h.ind<-t(D)%*%ginv(var)%*%diag(c(data[data$subject==i,]$weight[index]))%*%D
h.part<-h.part+h.ind
###Calculate the G_i in the formula;
x<-diag(c(data[data$subject==i,]$weight[index]))%*%matrix(y[index]-mu.R[index])
U<-t(D)%*%ginv(var)%*%x
### Get the design matrix for dropout model;
m_drop <- model.frame(model_drop, data[data$subject==i,], na.action='na.pass')
mat_drop <- as.data.frame(model.matrix(model_drop, m_drop))
exp.x <- exp(as.matrix(mat_drop)%*%as.matrix(alpha_drop))
lam_d <- exp.x/(1+exp.x)
row.names(m_drop)=seq(1:cluster[i])
index_d=as.numeric(row.names(m_drop))[-1]
index_d2=seq(1:cluster[i])[-length(seq(1:cluster[i]))]
F= as.matrix(colSums( na.omit( (m_drop[index_d,][,1]* mat_drop[index_d,]- c(lam_d[index_d])*mat_drop[index_d,] )*
mat_drop[index_d2,]) ) )
UF<-UF+U%*%t(F)
FF<-FF+F%*%t(F)
},
poisson={
index <- which(!is.na(data[data$subject==i,]$response))
var <- switch(corstr,
independence=scale*diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster)%*%cormax.ind(ncluster)%*%diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster),
exchangeable=scale*diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster)%*%cormax.exch(ncluster, rho)%*%diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster),
ar1=scale*diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster)%*%cormax.ar1(ncluster, rho)%*%diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster)
)
D<-mat.prod(covariate[index,], exp(covariate[index,]%*%beta))
h.ind<-t(D)%*%ginv(var)%*%diag(c(data[data$subject==i,]$weight[index]))%*%D
h.part<-h.part+h.ind
###Calculate the G_i in the formula;
x<-diag(c(data[data$subject==i,]$weight[index]))%*%matrix(y[index]-mu.R[index])
U<-t(D)%*%ginv(var)%*%x
### Get the design matrix for dropout model;
m_drop <- model.frame(model_drop, data[data$subject==i,], na.action='na.pass')
mat_drop <- as.data.frame(model.matrix(model_drop, m_drop))
exp.x <- exp(as.matrix(mat_drop)%*%as.matrix(alpha_drop))
lam_d <- exp.x/(1+exp.x)
row.names(m_drop)=seq(1:cluster[i])
index_d=as.numeric(row.names(m_drop))[-1]
index_d2=seq(1:cluster[i])[-length(seq(1:cluster[i]))]
F= as.matrix(colSums( na.omit( (m_drop[index_d,][,1]* mat_drop[index_d,]- c(lam_d[index_d])*mat_drop[index_d,] )*
mat_drop[index_d2,]) ) )
UF<-UF+U%*%t(F)
FF<-FF+F%*%t(F)
},
stop("Warnings: Invalid type of outcomes!")
)
}
G_pre<-UF%*%ginv(FF) ####the first two terms of Gi;
j.part<-matrix(0,nrow=length(mat[1,])-1,ncol=length(mat[1,])-1)
k.part<-matrix(0,nrow=length(mat[1,])-1,ncol=length(mat[1,])-1)
for (i in 1:size){
cluster_alt <- cluster.size(data[!is.na(data$response),]$subject)$n
ncluster <- cluster_alt[i]
y<-as.matrix(data[data$subject==i,]$response)
mu.R<-as.matrix(data[data$subject==i,]$mu.R)
covariate<-as.matrix(subset(mat[,-length(mat[1,])], mat$subj==unique(data$subject)[i])) ### specify the covariate matrix
switch(family,
gaussian={
index <- which(!is.na(data[data$subject==i,]$response))
var <- switch(corstr,
independence=scale*cormax.ind(ncluster),
exchangeable=scale*cormax.exch(ncluster, rho),
ar1=scale*cormax.ar1(ncluster, rho)
)
x<-diag(c(data[data$subject==i,]$weight[index]))%*%matrix(y[index]-mu.R[index])
### Get the design matrix for dropout model;
m_drop <- model.frame(model_drop, data[data$subject==i,],na.action='na.pass')
mat_drop <- as.data.frame(model.matrix(model_drop, m_drop))
exp.x <- exp(as.matrix(mat_drop)%*%as.matrix(alpha_drop))
lam_d <- exp.x/(1+exp.x)
row.names(m_drop)=seq(1:cluster[i])
index_d=as.numeric(row.names(m_drop))[-1]
index_d2=seq(1:cluster[i])[-length(seq(1:cluster[i]))]
F= as.matrix(colSums( na.omit( (m_drop[index_d,][,1]* mat_drop[index_d,]- c(lam_d[index_d])*mat_drop[index_d,] )*
mat_drop[index_d2,]) ) )
G<-G_pre%*%F
if (is.matrix(covariate[index,])==FALSE) {cal.m=t(t(covariate[index,]))
}else{cal.m=t(covariate[index,])}
E<-cal.m%*%ginv(var)%*%x
j.part<-j.part+(E-G)%*%t(E-G)
varI <- scale*cormax.ind(ncluster)
#k.part<-k.part+t(covariate[index,])%*%ginv(varI)%*%diag(c(data[data$subject==i,]$weight[index]))%*%covariate[index,]
k.part<-k.part+cal.m%*%ginv(varI)%*%covariate[index,]
},
binomial={
index <- which(!is.na(data[data$subject==i,]$response))
var <- switch(corstr,
independence=scale*diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster)%*%cormax.ind(ncluster)%*%diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster),
exchangeable=scale*diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster)%*%cormax.exch(ncluster, rho)%*%diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster),
ar1=scale*diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster)%*%cormax.ar1(ncluster, rho)%*%diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster)
)
D<-mat.prod(covariate[index,], exp(covariate[index,]%*%beta)/((1+exp(covariate[index,]%*%beta))^2))
x<-diag(c(data[data$subject==i,]$weight[index]))%*%matrix(y[index]-mu.R[index])
### Get the design matrix for dropout model;
m_drop <- model.frame(model_drop, data[data$subject==i,],na.action='na.pass')
mat_drop <- as.data.frame(model.matrix(model_drop, m_drop))
exp.x <- exp(as.matrix(mat_drop)%*%as.matrix(alpha_drop))
lam_d <- exp.x/(1+exp.x)
row.names(m_drop)=seq(1:cluster[i])
index_d=as.numeric(row.names(m_drop))[-1]
index_d2=seq(1:cluster[i])[-length(seq(1:cluster[i]))]
F= as.matrix(colSums( na.omit( (m_drop[index_d,][,1]* mat_drop[index_d,]- c(lam_d[index_d])*mat_drop[index_d,] )*
mat_drop[index_d2,]) ) )
G<-G_pre%*%F
E<-t(D)%*%ginv(var)%*%x
j.part<-j.part+(E-G)%*%t(E-G)
varI <- scale*diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster)%*%cormax.ind(ncluster)%*%diag(sqrt(c(exp(covariate[index,]%*%beta)/(1+exp(covariate[index,]%*%beta))^2)),ncluster)
#k.part<-k.part+t(D)%*%ginv(varI)%*%diag(c(data[data$subject==i,]$weight[index]))%*%D
k.part<-k.part+t(D)%*%ginv(varI)%*%D
},
poisson={
index <- which(!is.na(data[data$subject==i,]$response))
var <- switch(corstr,
independence=scale*diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster)%*%cormax.ind(ncluster)%*%diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster),
exchangeable=scale*diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster)%*%cormax.exch(ncluster, rho)%*%diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster),
ar1=scale*diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster)%*%cormax.ar1(ncluster, rho)%*%diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster)
)
D<-mat.prod(covariate[index,], exp(covariate[index,]%*%beta))
x<-diag(c(data[data$subject==i,]$weight[index]))%*%matrix(y[index]-mu.R[index])
### Get the design matrix for dropout model;
m_drop <- model.frame(model_drop, data[data$subject==i,],na.action='na.pass')
mat_drop <- as.data.frame(model.matrix(model_drop, m_drop))
exp.x <- exp(as.matrix(mat_drop)%*%as.matrix(alpha_drop))
index <- which(!is.na(data[data$subject==i,]$response))
lam_d <- exp.x/(1+exp.x)
row.names(m_drop)=seq(1:cluster[i])
index_d=as.numeric(row.names(m_drop))[-1]
index_d2=seq(1:cluster[i])[-length(seq(1:cluster[i]))]
F= as.matrix(colSums( na.omit( (m_drop[index_d,][,1]* mat_drop[index_d,]- c(lam_d[index_d])*mat_drop[index_d,] )*
mat_drop[index_d2,]) ) )
G<-G_pre%*%F
E<-t(D)%*%ginv(var)%*%x
j.part<-j.part+(E-G)%*%t(E-G)
varI <- scale*diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster)%*%cormax.ind(ncluster)%*%diag(sqrt(c(exp(covariate[index,]%*%beta))),ncluster)
#k.part<-k.part+t(D)%*%ginv(varI)%*%diag(c(data[data$subject==i,]$weight[index]))%*%D
k.part<-k.part+t(D)%*%ginv(varI)%*%D
},
stop("Warnings: Invalid type of outcomes!")
)
}
trace.R <- sum(diag(ginv(h.part)%*%j.part%*%ginv(h.part)%*%k.part))
#trace.R <- sum(diag(k.part%*%ginv(h.part)%*%j.part%*%ginv(h.part)))
#WQIC;
WQIC <- (-2)*quasi.R + 2*trace.R
WQICu <- (-2)*quasi.R + 2*px
# output the results;
out <- data.frame(QICWr=WQIC, QICWp=WQICu, Wquasi_lik=quasi.R)
return(round(out,1))
#return(list(QICWr=WQIC, QICWp=WQICu, Wquasi_lik=quasi.R))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.