R/GENMETA.summary.R

Defines functions GENMETA.summary

Documented in GENMETA.summary

#' Summarizing Generalized Meta Analysis
#'
#' This function prints the summary of GENMETA results.
#' @param object an object of class "GENMETA"
#' @param signi_digits an optional numeric indicating the number of significant digits to be shown in the summary. Default is 3.
#' @import MASS
#' @importFrom Matrix rankMatrix
#' @examples
#' # This example shows how to obtain the summary of GENMETA object.
#' #####
#' ### Basic setting
#' #####
#'d.X = 3 # number of covariates.
#' mu = matrix(rep(0,d.X), nrow=d.X) # mean vector of the covariates.

#' r1 = 0.3 # correlation coefficient of the covariates.
#' r2 = 0.6
#' r3 = 0.1
#' Sigma = matrix( 
#'   c(1, r1, r2,  
#'     r1, 1, r3, 
#'     r2, r3, 1), 
#'   nrow=d.X, 
#'   ncol=d.X) # covariance matrix of the covariates.


#' beta.star = matrix(c(-1.2, log(1.3), log(1.3), log(1.3)),nrow = d.X+1) # beta.star
#' #beta.star = matrix(c(-1.2, 0.26, 0.26, 0.26),nrow = d.X+1) # beta.star
#' #beta.star = matrix(c(-3, 1, 2, 3),nrow = d.X+1) # beta.star

#' n1 = 300 # sample size of the 1st data set.
#' n2 = 500 # 2nd
#' n3 = 1000 # 3rd

#' n = 50
#' sim=1
#' set.seed(sim)
# Generate the reference data set
#' X.rf = MASS::mvrnorm(n = n, mu, Sigma)


# Generate data set 1. m1 means model 1.
#' X.m1 = MASS::mvrnorm(n = n1, mu, Sigma) # Generate the covariates.
#' X.m1.1 = cbind(rep(1, n1), X.m1) # Add a column of 1's to X.m1.
#' p.m1 = 1/(1+exp(-X.m1.1%*%beta.star)) # the vector of probabilities
#' Y.m1 = rbinom(n1, size=1, p.m1) # the Bernoulli responses
#' # print(p.m1[1])
#' # print(mean(Y.m1))
#' # print(mean(p.m1))

#' # Generate data set 2. m1 means model 2.
#' X.m2 = MASS::mvrnorm(n = n2, mu, Sigma)
#' X.m2.1 = cbind(rep(1, n2), X.m2)
#' p.m2 = 1/(1+exp(-X.m2.1%*%beta.star))
#' Y.m2 = rbinom(n2, size=1, p.m2)

# Generate data set 3. m1 means model 3.
#' X.m3 = MASS::mvrnorm(n = n3, mu, Sigma)
#' X.m3.1 = cbind(rep(1, n3), X.m3)
#' p.m3 = 1/(1+exp(-X.m3.1%*%beta.star))
#' Y.m3 = rbinom(n3, size=1, p.m3)

#' #####
#' ### Create data sets in the format of data frame.
#' #####
#' data.m1 = data.frame(Y=Y.m1, X.m1)
#' data.m2 = data.frame(Y=Y.m2, X.m2)
#' data.m3 = data.frame(Y=Y.m3, X.m3)
# str(data.m1)



#' #####
#' ### Apply logistic regression with reduced models to the data sets 
#' #####
#' logit.m1 <- glm(Y ~ X1 + X2, data = data.m1, family = "binomial")
#' # print(logit.m1)
#' if(logit.m1$converged == FALSE)
#' {
#'   print("glm for logit.m1 is not convergent.")
#'   next
#' }

#' logit.m2 <- glm(Y ~ X2 + X3, data = data.m2, family = "binomial")
#' # print(logit.m2)
#' if(logit.m2$converged == FALSE)
#' {
#'   print("glm for logit.m2 is not convergent.")
#'   next
#' }

#' logit.m3 <- glm(Y ~ X1 + X3, data = data.m3, family = "binomial")
#' # print(logit.m3)
#' if(logit.m3$converged == FALSE)
#' {
#'   print("glm for logit.m3 is not convergent.")
#'   next
#' }


#' #####
#' ### Obtain the estimators of the parameters in the reduced models. 
#' #####
#' theta.m1 = logit.m1$coefficients
#' theta.m2 = logit.m2$coefficients
#' theta.m3 = logit.m3$coefficients



#' #####
#' ### Find the covariance matrix estimators for the reduced models
#' #####



#' #####
#' # Basic notations for inputs
#' #####

#' K = 3 # Number of data sets

#' A1 = c(1, 2) # index set A1, the indexes of the covariates of data set 1.
#' A2 = c(2, 3) # index set A2
#' A3 = c(1, 3) # index set A3


#' X.m1.used = cbind(rep(1, n1), X.m1[, A1, drop=FALSE])
#' X.m2.used = cbind(rep(1, n2), X.m2[, A2, drop=FALSE])
#' X.m3.used = cbind(rep(1, n3), X.m3[, A3, drop=FALSE])
#' # str(X.m1.used)
#' # str(X.m2.used)
#' # str(X.m3.used)

#' ##### Find Sigma.m1

#' T.1 = matrix(rep(0, (length(A1)+1)^2), nrow=length(A1)+1)
#' T.2 = T.1

#' for (i in 1:n1)
#' {
#' a = as.vector(exp(-X.m1.used[i, , drop=FALSE]%*%theta.m1))
#'   T.1 = T.1 + (a/(1+a)^2) * (t(X.m1.used[i, , drop=FALSE])%*%X.m1.used[i, , drop=FALSE])
#' }

#' for (i in 1:n1)
#' {
#'   a = as.vector(1/( 1 + exp(-X.m1.used[i, , drop=FALSE]%*%theta.m1)))
#'   T.2 = T.2 + (Y.m1[i]-a)^2 * (t(X.m1.used[i, , drop=FALSE])%*%X.m1.used[i, , drop=FALSE])
#' }

#' Sigma.m1 = solve(T.1)%*%T.2%*%solve(T.1) # This is actually Sigma.m1.n1. 

#' ##### Find Sigma.m2

#' T.1 = matrix(rep(0, (length(A2)+1)^2), nrow=length(A2)+1)
#' T.2 = T.1

#' for (i in 1:n2)
#' {
#'   a = as.vector(exp(-X.m2.used[i, , drop=FALSE]%*%theta.m2))
#'   T.1 = T.1 + (a/(1+a)^2) * (t(X.m2.used[i, , drop=FALSE])%*%X.m2.used[i, , drop=FALSE])
#' }

#' for (i in 1:n2)
#' {
#'   a = as.vector(1/( 1 + exp(-X.m2.used[i, , drop=FALSE]%*%theta.m2)))
#'   T.2 = T.2 + (Y.m2[i]-a)^2 * (t(X.m2.used[i, , drop=FALSE])%*%X.m2.used[i, , drop=FALSE])
#' }

#' Sigma.m2 = solve(T.1)%*%T.2%*%solve(T.1)


#' ##### Find Sigma.m3

#' T.1 = matrix(rep(0, (length(A3)+1)^2), nrow=length(A3)+1)
#' T.2 = T.1

#' for (i in 1:n3)
#' {
#'   a = as.vector(exp(-X.m3.used[i, , drop=FALSE]%*%theta.m3))
#'   T.1 = T.1 + (a/(1+a)^2) * (t(X.m3.used[i, , drop=FALSE])%*%X.m3.used[i, , drop=FALSE])
#' }

#' for (i in 1:n3)
#' {
#'   a = as.vector(1/( 1 + exp(-X.m3.used[i, , drop=FALSE]%*%theta.m3)))
#'   T.2 = T.2 + (Y.m3[i]-a)^2 * (t(X.m3.used[i, , drop=FALSE])%*%X.m3.used[i, , drop=FALSE])
#' }

#' Sigma.m3 = solve(T.1)%*%T.2%*%solve(T.1)







#' names(theta.m1)=c("(Intercept)","Age","Height")
#' names(theta.m2)=c("(Intercept)","Height", "Weight")
#' names(theta.m3)=c("(Intercept)","Age", "Weight")

###now put in the GENMETA example

#' study1 = list(Coeff=theta.m1,Covariance=Sigma.m1,Sample_size=n1)
#' study2 = list(Coeff=theta.m2,Covariance=Sigma.m2,Sample_size=n2)
#' study3 = list(Coeff=theta.m3,Covariance=Sigma.m3,Sample_size=n3)

#' studies = list(study1,study2,study3)
#' model = "logistic"

#' reference = cbind(rep(1,n), X.rf)
#' colnames(reference) = c("(Intercept)","Age","Height", "Weight")
#' result.same = GENMETA(studies, reference, model, 
#' initial_val = c(-1.2, log(1.3), log(1.3), log(1.3)))
#' GENMETA.summary(result.same)
#' @export
# summary.GENMETA <-function(object, ...){
#   UseMethod("summary")
#   NextMethod("generic = NULL, object = NULL", ...)
# }
GENMETA.summary <- function(object, signi_digits = 3)
{
  x <- object
  GMeta_opt_estimate <- as.vector(x[[1]])
  GMeta_opt_std_error <- sqrt(as.vector(diag(x[[2]])))
  Var_opt_GMeta <- diag(GMeta_opt_std_error)
  z_stat_opt_GMeta <- GMeta_opt_estimate/GMeta_opt_std_error
  p_val_opt_GMeta <- 1 - pnorm(abs(z_stat_opt_GMeta))
  summary_data_frame_opt <- data.frame(cbind(GMeta_opt_estimate, GMeta_opt_std_error, z_stat_opt_GMeta, p_val_opt_GMeta))
  colnames(summary_data_frame_opt) <- c("Estimate", "Std.Error", "z value", "Pr(>|z|)")
  rownames(summary_data_frame_opt) <- names(x[[1]])
  summary_data_frame_opt <- signif(summary_data_frame_opt, signi_digits)
  signif_column <- factor(noquote(sapply(summary_data_frame_opt[, 4], sign.star)))
  summary_data_frame_opt <- cbind(summary_data_frame_opt, signif_column)
  colnames(summary_data_frame_opt)[5] <- paste0(' ')
  cat("Call:\n")
  print(x[[5]])
  cat("\n")
  cat("Coefficients: \n")
  print.data.frame(summary_data_frame_opt, print.gap = 2)
  cat("\n---\n")
  cat("Significant codes:\n")
  cat("0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 \n")
  cat("\n")
  cat("Total number of iterations: ")
  cat(x[[4]])

}
28pro92/packages-GENMETA documentation built on Nov. 13, 2022, 1:49 p.m.