source("generExampData.r")
family <- "binomial"
m <- 100 ; nVec <- rep(10,m)
betaATrue <- c(0.35,0.96)
betaBTrue <- c(-0.47,1.06,-1.31)
dA <- length(betaATrue) ; dB <- length(betaBTrue)
SigmaTrue <- matrix(c(0.56,-0.34,-0.34,0.89),2,2)
XAdistn <- "unif" ; XBdistn <- "unif"
allData <- generExampData(m,nVec,XAdistn,XBdistn,betaATrue,betaBTrue,SigmaTrue,family,seed=currSeed)
allData <- generExampData(m,nVec,XAdistn,XBdistn,betaATrue,betaBTrue,SigmaTrue,family)
setwd("~/Library/CloudStorage/Dropbox/UTS/GLMMasymptotics/Code/MBWcode")
source("generExampData.r")
allData <- generExampData(m,nVec,XAdistn,XBdistn,betaATrue,betaBTrue,SigmaTrue,family,seed=currSeed)
allData <- generExampData(m,nVec,XAdistn,XBdistn,betaATrue,betaBTrue,SigmaTrue,family)
library(MASS)
allData <- generExampData(m,nVec,XAdistn,XBdistn,betaATrue,betaBTrue,SigmaTrue,family)
aaa <- glmmTMB(y ~ 1 + XA[,-1] + XB + (1 + XA[,-1]|idNum),data=allData,family=family)
library(glmmTMB)
aaa <- glmmTMB(y ~ 1 + XA[,-1] + XB + (1 + XA[,-1]|idNum),data=allData,family=family)
aaa <- glmmTMB::glmmTMB(y ~ 1 + XA[,-1] + XB + (1 + XA[,-1]|idNum),data=allData,family=family)
aaa <- glmmTMB:::glmmTMB(y ~ 1 + XA[,-1] + XB + (1 + XA[,-1]|idNum),data=allData,family=family)
fitMLE <- glmmTMB(y ~ 1 + XA[,-1] + XB + (1 + XA[,-1]|idNum),data=allData,family=family)
betaMLEtable <- coefficients(summary(fitMLE))
betaMLE <- fitMLE$fit$par[1:(dA+dB)]
betaAMLE <- betaMLE[1:dA]
betaBMLE <- betaMLE[(dA+1):(dA+dB)]
SDbetaXCT <- as.numeric(betaMLEtable$cond[,2])
SDbetaAXCT <- SDbetaXCT[1:dA]
SDbetaBXCT <- SDbetaXCT[(dA+1):(dA+dB)]
SigmaMLE <- matrix((summary(fitMLE)$varcor$cond$idNum)[1:(dA^2)],dA,dA)
# source("vecInverse.r")
# source("vechInverse.r")
# source("starProdCpp.r")
# source("OmegaAAAdCpp.r")
# source("OmegaAABdCpp.r")
source("OneTermAsyCov.r")
source("TwoTermAsyCov.r")
JWBres <- OneTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB)
library(glmmTMB)
library(MASS)
library(matrixcalc)
library(cubature)
library(Rcpp)
library(RcppArmadillo)
library(parallel)
library(doParallel)
JWBres <- OneTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB)
cores <- detectCores()
clus <- makeCluster(cores[1]-1)
registerDoParallel(clus)
MBWres <- TwoTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB,m,nVec)
MBWres <- TwoTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB,m,nVec)
source("vecInverse.r")
MBWres <- TwoTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB,m,nVec)
remove(vecInverse)
source("vechInverse.r")
source("vechInverse.r")
MBWres <- TwoTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB,m,nVec)
source("TwoTermAsyCov.r")
source("vecInverse.r")
cores <- detectCores()
clus <- makeCluster(cores[1]-1)
registerDoParallel(clus)
# Obtain the two-term asymptotic covariance for beta and Sigma based on the
# Maestrini, Bhaskaran & Wand asymptotic results:
MBWres <- TwoTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB,m,nVec)
stopCluster(clus)
1-0.05/2
warnings()
lwdVal <- 3
source("graphicCompareCIs.r")
betaATrue
betaBTrue
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(fitMLE,JWBRes,MBWRes,betaATrue,betaBTrue,SigmaTrue)
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(fitMLE,JWBres,MBWres,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
confint(fitMLE)
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
betaTrue
betaATrue
betaBTrue
sqrt(diag(SigmaTrue))
SigmaTrue[2,1]
SigmaTrue
SigmaTrue[2,1]/(sqrt(SigmaTrue[1,1]*SigmaTrue[2,2]))
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
setwd("~/Library/CloudStorage/Dropbox/UTS/GLMMasymptotics/Code/MBWcode")
confint(fitMLE,method="Wald")
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
legend("topleft",legend=c("M.B.W. theory","J.W.B. theory","glmmTMB","true value"),
col=c(colVec,"black"),lwd=c(rep(lwdVal,3),0.5*lwdVal),bty="n",pch = "|")
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
source("graphicCompareCIs.r")
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
setwd("~/Library/CloudStorage/Dropbox/UTS/GLMMasymptotics/Code/MBWcode")
######################## R script: MBWinferenceExample #########################
# For computing confidence intervals for fixed effects parameters and standard
# deviation and correlation parameters from the random effects covariance matrix
# in logistic generalized linear mixed models with bivariate random effects.
# Confidence intervals based on the theory of Maestrini, Bhaskaran & Wand are
# obtained through the function TwoTermAsyCov() and compared with those from
# Jiang, Wand & Bhaskaran (2022, JRSSB) produced using OneTermAsyCov() and those
# from the 'glmmTMB' package. Graphical comparison of the confidence intervals
# is performed via the function graphicCompareCIs().
# Last changed: 26 JUN 2023
# Load required package and functions:
library(glmmTMB)
library(MASS)
library(matrixcalc)
library(cubature)
library(Rcpp)
library(RcppArmadillo)
library(parallel)
library(doParallel)
library(stipkg)
source("generExampData.r")
source("OneTermAsyCov.r")
source("TwoTermAsyCov.r")
source("graphicCompareCIs.r")
# Set values corresponding to the specified example:
family <- "binomial"
m <- 100 ; nVec <- rep(10,m)
betaATrue <- c(0.35,0.96)
betaBTrue <- c(-0.47,1.06,-1.31)
dA <- length(betaATrue) ; dB <- length(betaBTrue)
SigmaTrue <- matrix(c(0.56,-0.34,-0.34,0.89),2,2)
XAdistn <- "unif" ; XBdistn <- "unif"
# Obtain data:
allData <- generExampData(m,nVec,XAdistn,XBdistn,betaATrue,betaBTrue,SigmaTrue,family)
# Obtain the maximimum likelihood fit using the glmmTMB::glmmTMB() function:
fitMLE <- glmmTMB(y ~ 1 + XA[,-1] + XB + (1 + XA[,-1]|idNum),data=allData,family=family)
betaMLEtable <- coefficients(summary(fitMLE))
betaMLE <- fitMLE$fit$par[1:(dA+dB)]
betaAMLE <- betaMLE[1:dA]
betaBMLE <- betaMLE[(dA+1):(dA+dB)]
# Obtain the one-term asymptotic covariance for beta and Sigma based on the
# Jiang, Wand & Bhaskaran (2022) asymptotic normality theorem:
JWBres <- OneTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB)
SigmaMLE <- matrix((summary(fitMLE)$varcor$cond$idNum)[1:(dA^2)],dA,dA)
######################## R script: MBWinferenceExample #########################
# For computing confidence intervals for fixed effects parameters and standard
# deviation and correlation parameters from the random effects covariance matrix
# in logistic generalized linear mixed models with bivariate random effects.
# Confidence intervals based on the theory of Maestrini, Bhaskaran & Wand are
# obtained through the function TwoTermAsyCov() and compared with those from
# Jiang, Wand & Bhaskaran (2022, JRSSB) produced using OneTermAsyCov() and those
# from the 'glmmTMB' package. Graphical comparison of the confidence intervals
# is performed via the function graphicCompareCIs().
# Last changed: 26 JUN 2023
# Load required package and functions:
library(glmmTMB)
library(MASS)
library(matrixcalc)
library(cubature)
library(Rcpp)
library(RcppArmadillo)
library(parallel)
library(doParallel)
library(stipkg)
source("generExampData.r")
source("OneTermAsyCov.r")
source("TwoTermAsyCov.r")
source("graphicCompareCIs.r")
# Set values corresponding to the specified example:
family <- "binomial"
m <- 100 ; nVec <- rep(10,m)
betaATrue <- c(0.35,0.96)
betaBTrue <- c(-0.47,1.06,-1.31)
dA <- length(betaATrue) ; dB <- length(betaBTrue)
SigmaTrue <- matrix(c(0.56,-0.34,-0.34,0.89),2,2)
XAdistn <- "unif" ; XBdistn <- "unif"
# Obtain data:
allData <- generExampData(m,nVec,XAdistn,XBdistn,betaATrue,betaBTrue,SigmaTrue,family)
# Obtain the maximimum likelihood fit using the glmmTMB::glmmTMB() function:
fitMLE <- glmmTMB(y ~ 1 + XA[,-1] + XB + (1 + XA[,-1]|idNum),data=allData,family=family)
betaMLEtable <- coefficients(summary(fitMLE))
betaMLE <- fitMLE$fit$par[1:(dA+dB)]
betaAMLE <- betaMLE[1:dA]
betaBMLE <- betaMLE[(dA+1):(dA+dB)]
SigmaMLE <- matrix((summary(fitMLE)$varcor$cond$idNum)[1:(dA^2)],dA,dA)
# Obtain the one-term asymptotic covariance for beta and Sigma based on the
# Jiang, Wand & Bhaskaran (2022) asymptotic normality theorem:
JWBres <- OneTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB)
# Obtain the two-term asymptotic covariance for beta and Sigma based on the
# Maestrini, Bhaskaran & Wand asymptotic results using a parallel backend:
cores <- detectCores()
clus <- makeCluster(cores[1]-1)
registerDoParallel(clus)
MBWres <- TwoTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB,m,nVec)
stopCluster(clus)
# Obtain the approximate standard deviations based on the
# Jiang, Wand & Bhaskaran (2022) asymptotic normality theorem:
JWBres <- OneTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB)
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
######################### End of MBWinferenceExample ###########################
setwd("~/Library/CloudStorage/Dropbox/UTS/GLMMasymptotics/Code/MBWcode")
######################## R script: MBWinferenceExample #########################
# For computing confidence intervals for fixed effects parameters and standard
# deviation and correlation parameters from the random effects covariance matrix
# in logistic generalized linear mixed models with bivariate random effects.
# Confidence intervals based on the theory of Maestrini, Bhaskaran & Wand are
# obtained through the function TwoTermAsyCov() and compared with those from
# Jiang, Wand & Bhaskaran (2022, JRSSB) produced using OneTermAsyCov() and those
# from the 'glmmTMB' package. Graphical comparison of the confidence intervals
# is performed via the function graphicCompareCIs().
# Last changed: 26 JUN 2023
# Load required package and functions:
library(glmmTMB)
library(MASS)
library(matrixcalc)
library(cubature)
library(Rcpp)
library(RcppArmadillo)
library(parallel)
library(doParallel)
library(stipkg)
source("generExampData.r")
source("OneTermAsyCov.r")
source("TwoTermAsyCov.r")
source("graphicCompareCIs.r")
# Set values corresponding to the specified example:
family <- "binomial"
m <- 150 ; nVec <- rep(15,m)
betaATrue <- c(0.35,0.96)
betaBTrue <- c(-0.47,1.06,-1.31)
dA <- length(betaATrue) ; dB <- length(betaBTrue)
SigmaTrue <- matrix(c(0.56,-0.34,-0.34,0.89),2,2)
XAdistn <- "unif" ; XBdistn <- "unif"
# Obtain data:
allData <- generExampData(m,nVec,XAdistn,XBdistn,betaATrue,betaBTrue,SigmaTrue,family)
# Obtain the maximimum likelihood fit using the glmmTMB::glmmTMB() function:
fitMLE <- glmmTMB(y ~ 1 + XA[,-1] + XB + (1 + XA[,-1]|idNum),data=allData,family=family)
betaMLEtable <- coefficients(summary(fitMLE))
betaMLE <- fitMLE$fit$par[1:(dA+dB)]
betaAMLE <- betaMLE[1:dA]
betaBMLE <- betaMLE[(dA+1):(dA+dB)]
SigmaMLE <- matrix((summary(fitMLE)$varcor$cond$idNum)[1:(dA^2)],dA,dA)
# Obtain the one-term asymptotic covariance for beta and Sigma based on the
# Jiang, Wand & Bhaskaran (2022) asymptotic normality theorem:
JWBres <- OneTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB)
# Obtain the two-term asymptotic covariance for beta and Sigma based on the
# Maestrini, Bhaskaran & Wand asymptotic results using a parallel backend:
cores <- detectCores()
clus <- makeCluster(cores[1]-1)
registerDoParallel(clus)
MBWres <- TwoTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB,m,nVec)
stopCluster(clus)
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
######################### End of MBWinferenceExample ###########################
setwd("~/Library/CloudStorage/Dropbox/UTS/GLMMasymptotics/Code/MBWcode")
system("R CMD build stipkg")
system("R CMD INSTALL stipkg_1.0.tar.gz")
setwd("~/Library/CloudStorage/Dropbox/UTS/GLMMasymptotics/Code/MBWcode")
######################## R script: MBWinferenceExample #########################
# For computing confidence intervals for fixed effects parameters and standard
# deviation and correlation parameters from the random effects covariance matrix
# in logistic generalized linear mixed models with bivariate random effects.
# Confidence intervals based on the theory of Maestrini, Bhaskaran & Wand are
# obtained through the function TwoTermAsyCov() and compared with those from
# Jiang, Wand & Bhaskaran (2022, JRSSB) produced using OneTermAsyCov() and those
# from the 'glmmTMB' package. Graphical comparison of the confidence intervals
# is performed via the function graphicCompareCIs().
# Last changed: 26 JUN 2023
# Load required package and functions:
library(glmmTMB)
library(MASS)
library(matrixcalc)
library(cubature)
library(Rcpp)
library(RcppArmadillo)
library(parallel)
library(doParallel)
library(stipkg)
source("generExampData.r")
source("OneTermAsyCov.r")
source("TwoTermAsyCov.r")
source("graphicCompareCIs.r")
# Set values corresponding to the specified example:
family <- "binomial"
m <- 150 ; nVec <- rep(15,m)
betaATrue <- c(0.35,0.96)
betaBTrue <- c(-0.47,1.06,-1.31)
dA <- length(betaATrue) ; dB <- length(betaBTrue)
SigmaTrue <- matrix(c(0.56,-0.34,-0.34,0.89),2,2)
XAdistn <- "unif" ; XBdistn <- "unif"
# Obtain data:
allData <- generExampData(m,nVec,XAdistn,XBdistn,betaATrue,betaBTrue,SigmaTrue,family)
# Obtain the maximimum likelihood fit using the glmmTMB::glmmTMB() function:
fitMLE <- glmmTMB(y ~ 1 + XA[,-1] + XB + (1 + XA[,-1]|idNum),data=allData,family=family)
betaMLEtable <- coefficients(summary(fitMLE))
betaMLE <- fitMLE$fit$par[1:(dA+dB)]
betaAMLE <- betaMLE[1:dA]
betaBMLE <- betaMLE[(dA+1):(dA+dB)]
SigmaMLE <- matrix((summary(fitMLE)$varcor$cond$idNum)[1:(dA^2)],dA,dA)
# Obtain the one-term asymptotic covariance for beta and Sigma based on the
# Jiang, Wand & Bhaskaran (2022) asymptotic normality theorem:
JWBres <- OneTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB)
# Obtain the two-term asymptotic covariance for beta and Sigma based on the
# Maestrini, Bhaskaran & Wand asymptotic results using a parallel backend:
cores <- detectCores()
clus <- makeCluster(cores[1]-1)
registerDoParallel(clus)
MBWres <- TwoTermAsyCov(betaAMLE,betaBMLE,SigmaMLE,allData$XA,allData$XB,m,nVec)
stopCluster(clus)
# Graphically compare confidence intervals from different approaches:
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
######################### End of MBWinferenceExample ###########################
sqrt(diag(SigmaTrue))
SigmaTrue[2,1]/sqrt(SigmaTrue[1,1]*SigmaTrue[2,2])
source("graphicCompareCIs.r")
graphicCompareCIs(MBWres,JWBres,fitMLE,betaATrue,betaBTrue,SigmaTrue)
