From 7ce6af2bc9c1b710e29f8458cec0efdfcd3ae414 Mon Sep 17 00:00:00 2001 From: ctbus Date: Thu, 21 Sep 2023 14:53:25 -0400 Subject: [PATCH] Update formatting --- R/simulate_beta_data.R | 41 +++++++--------- R/simulate_zero_inflated_beta_data.R | 48 +++++++++---------- R/zibr.R | 4 +- ...e_zero_inflated_beta_random_effect_data.Rd | 16 +++---- man/zibr.Rd | 4 +- tests/testthat/test-zibr.R | 4 +- vignettes/zibr.Rmd | 5 +- 7 files changed, 58 insertions(+), 64 deletions(-) diff --git a/R/simulate_beta_data.R b/R/simulate_beta_data.R index 1c3b719..bc7d32b 100755 --- a/R/simulate_beta_data.R +++ b/R/simulate_beta_data.R @@ -1,15 +1,13 @@ #' @importFrom stats rnorm rbeta -simulate_beta_random_effect_data <- function(subject.n = 50, time.n = 5, v = 2, +simulate_beta_random_effect_data <- function(subject_n = 50, time_n = 5, v = 2, beta = as.matrix(c(-0.5, -0.5, 0.5)), - Z = NA, s2 = 1, sim.seed = 100) { - ########################## - ########################## + Z = NA, s2 = 1, sim_seed = 100) { #### Beta regression if (length(NA) == 1 & any(is.na(Z))) { - set.seed(sim.seed * 5000 + 10) + set.seed(sim_seed * 5000 + 10) Z <- as.matrix(data.frame( - log.Time = as.matrix(log(rep(seq(1, time.n), subject.n))), - Treatment = as.matrix(c(rep(0, subject.n * time.n / 2), rep(1, subject.n * time.n / 2))) + log.Time = as.matrix(log(rep(seq(1, time_n), subject_n))), + Treatment = as.matrix(c(rep(0, subject_n * time_n / 2), rep(1, subject_n * time_n / 2))) )) } if (is.null(colnames(Z))) { @@ -17,23 +15,23 @@ simulate_beta_random_effect_data <- function(subject.n = 50, time.n = 5, v = 2, colnames(Z) <- paste("var", seq(1, ncol(Z)), sep = "") } - set.seed(sim.seed * 5000 + 1) - c <- as.matrix(rnorm(subject.n, mean = 0, sd = s2)) - c.rep <- as.matrix(as.vector(matrix(c, nrow = time.n, ncol = length(c), byrow = TRUE))) + set.seed(sim_seed * 5000 + 1) + c <- as.matrix(rnorm(subject_n, mean = 0, sd = s2)) + c_rep <- as.matrix(as.vector(matrix(c, nrow = time_n, ncol = length(c), byrow = TRUE))) ##### - subject.ind <- as.vector(matrix(paste("Subject_", seq(1, subject.n), sep = ""), nrow = time.n, ncol = subject.n, byrow = TRUE)) - time.ind <- rep(seq(1, time.n), subject.n) + subject_ind <- as.vector(matrix(paste("Subject_", seq(1, subject_n), sep = ""), nrow = time_n, ncol = subject_n, byrow = TRUE)) + time_ind <- rep(seq(1, time_n), subject_n) ###### - Z.aug <- cbind(intersept = 1, Z) + z_aug <- cbind(intersept = 1, Z) # browser() - logit.u <- Z.aug %*% beta + c.rep + logit_u <- z_aug %*% beta + c_rep ###### - u <- 1 / (1 + exp(-logit.u)) + u <- 1 / (1 + exp(-logit_u)) ## v is the phi ###### - set.seed(sim.seed * 5000 + 4) - Y <- rbeta(subject.n * time.n, shape1 = u * v, shape2 = (1 - u) * v) + set.seed(sim_seed * 5000 + 4) + Y <- rbeta(subject_n * time_n, shape1 = u * v, shape2 = (1 - u) * v) if (any(Y > 1 - 10^(-6))) { Y[Y > 1 - 10^(-6)] <- 1 - 10^(-6) } @@ -42,7 +40,7 @@ simulate_beta_random_effect_data <- function(subject.n = 50, time.n = 5, v = 2, #### betareg can not fit random effect model #### so set the s2 to a small value (small random effect) # library(betareg) - # tdata <- data.frame(Y=Y,Z,SID=subject.ind) + # tdata <- data.frame(Y=Y,Z,SID=subject_ind) # gy <- betareg(Y ~ log.Time + as.factor(Treatment) , data = tdata,type='ML') # summary(gy) # gy <- betareg(Y ~ log.Time + as.factor(Treatment) , data = tdata,type='BC') @@ -51,11 +49,8 @@ simulate_beta_random_effect_data <- function(subject.n = 50, time.n = 5, v = 2, # summary(gy) # # - # fit_beta_random_effect(Z=Z,Y=Y,subject.ind=subject.ind,time.ind=time.ind, + # fit_beta_random_effect(Z=Z,Y=Y,subject_ind=subject_ind,time_ind=time_ind, # quad.n=30,verbose=FALSE) - return(list( - Y = Y, Z = Z, c = c, u = u, v = v, beta = beta, s2 = s2, - subject.ind = subject.ind, time.ind = time.ind - )) + list(Y = Y, Z = Z, c = c, u = u, v = v, beta = beta, s2 = s2, subject_ind = subject_ind, time_ind = time_ind) } diff --git a/R/simulate_zero_inflated_beta_data.R b/R/simulate_zero_inflated_beta_data.R index 2a02830..e9cf378 100755 --- a/R/simulate_zero_inflated_beta_data.R +++ b/R/simulate_zero_inflated_beta_data.R @@ -1,7 +1,7 @@ #' Simulate data according to zero-inflated beta random effects model #' -#' @param subject.n number of subjects -#' @param time.n number of time points for each subject +#' @param subject_n number of subjects +#' @param time_n number of time points for each subject #' @param v the dispersion parameter in beta component #' @param alpha the coefficients in logistic component #' @param beta the coefficients in beta component @@ -9,7 +9,7 @@ #' @param Z the covariates in beta component #' @param s1 the stardard deviation of random effect in logistic component #' @param s2 the stardard deviation of random effect in beta component -#' @param sim.seed the random seed +#' @param sim_seed the random seed #' @return Y the bacterial abundance generated from the model #' @return X the covariates in logistic component #' @return Z the covariates in beta component @@ -24,29 +24,29 @@ #' @examples #' \dontrun{ #' simulate_zero_inflated_beta_random_effect_data( -#' subject.n = 100, time.n = 5, +#' subject_n = 100, time_n = 5, #' X = as.matrix(c(rep(0, 50 * 5), rep(1, 50 * 5))), #' alpha = as.matrix(c(-0.5, 1)), #' beta = as.matrix(c(-0.5, 0.5)), #' s1 = 1, s2 = 0.8, #' v = 5, -#' sim.seed = 100 +#' sim_seed = 100 #' ) #' } -simulate_zero_inflated_beta_random_effect_data <- function(subject.n = 50, time.n = 5, v = 2, +simulate_zero_inflated_beta_random_effect_data <- function(subject_n = 50, time_n = 5, v = 2, alpha = as.matrix(c(0, 0.5, -1)), beta = as.matrix(c(-0.5, -0.5, 0.5)), X = NA, Z = NA, - s1 = 0.2, s2 = 0.2, sim.seed = 100) { + s1 = 0.2, s2 = 0.2, sim_seed = 100) { ###### if (length(NA) == 1 & any(is.na(X))) { - set.seed(sim.seed * 5000 + 10) + set.seed(sim_seed * 5000 + 10) X <- as.matrix(data.frame( - log.Time = as.matrix(log(rep(seq(1, time.n), subject.n))), - Treatment = as.matrix(c(rep(0, subject.n * time.n / 2), rep(1, subject.n * time.n / 2))) + log.Time = as.matrix(log(rep(seq(1, time_n), subject_n))), + Treatment = as.matrix(c(rep(0, subject_n * time_n / 2), rep(1, subject_n * time_n / 2))) )) - # X <- as.matrix(runif(subject.n*time.n,-1,1)) - # X <- as.matrix(c(rep(0,N*time.n/2),rep(1,N*time.n/2))) + # X <- as.matrix(runif(subject_n*time_n,-1,1)) + # X <- as.matrix(c(rep(0,N*time_n/2),rep(1,N*time_n/2))) } if (is.null(colnames(X))) { X <- as.matrix(X) @@ -57,15 +57,15 @@ simulate_zero_inflated_beta_random_effect_data <- function(subject.n = 50, time. } - set.seed(sim.seed * 5000 + 1) - b <- as.matrix(rnorm(subject.n, mean = 0, sd = s1)) - b.rep <- as.matrix(as.vector(matrix(b, nrow = time.n, ncol = length(b), byrow = TRUE))) - set.seed(sim.seed * 5000 + 2) - c <- as.matrix(rnorm(subject.n, mean = 0, sd = s2)) - c.rep <- as.matrix(as.vector(matrix(c, nrow = time.n, ncol = length(c), byrow = TRUE))) + set.seed(sim_seed * 5000 + 1) + b <- as.matrix(rnorm(subject_n, mean = 0, sd = s1)) + b.rep <- as.matrix(as.vector(matrix(b, nrow = time_n, ncol = length(b), byrow = TRUE))) + set.seed(sim_seed * 5000 + 2) + c <- as.matrix(rnorm(subject_n, mean = 0, sd = s2)) + c.rep <- as.matrix(as.vector(matrix(c, nrow = time_n, ncol = length(c), byrow = TRUE))) ##### - subject.ind <- as.vector(matrix(paste("Subject_", seq(1, subject.n), sep = ""), nrow = time.n, ncol = subject.n, byrow = TRUE)) - time.ind <- rep(seq(1, time.n), subject.n) + subject.ind <- as.vector(matrix(paste("Subject_", seq(1, subject_n), sep = ""), nrow = time_n, ncol = subject_n, byrow = TRUE)) + time.ind <- rep(seq(1, time_n), subject_n) ###### X.aug <- cbind(interespt = 1, X) Z.aug <- cbind(intersept = 1, Z) @@ -78,7 +78,7 @@ simulate_zero_inflated_beta_random_effect_data <- function(subject.n = 50, time. ###### p <- 1 / (1 + exp(-logit.p)) u <- 1 / (1 + exp(-logit.u)) - # set.seed(sim.seed+2) + # set.seed(sim_seed+2) # v <- runif(1,0,5) ## v is the phi if (is.na(v)) { v <- 2 @@ -86,8 +86,8 @@ simulate_zero_inflated_beta_random_effect_data <- function(subject.n = 50, time. ###### - Y <- rep(NA, subject.n * time.n) - set.seed(sim.seed * 5000 + 3) + Y <- rep(NA, subject_n * time_n) + set.seed(sim_seed * 5000 + 3) ind.mix <- rbinom(length(Y), 1, p) for (i in 1:length(Y)) { if (ind.mix[i] == 0) { @@ -95,7 +95,7 @@ simulate_zero_inflated_beta_random_effect_data <- function(subject.n = 50, time. } if (ind.mix[i] == 1) { # rbeta(n, shape1, shape2, ncp = 0) - set.seed(sim.seed * 5000 + i) + set.seed(sim_seed * 5000 + i) Y[i] <- rbeta(1, shape1 = u[i] * v, shape2 = (1 - u[i]) * v) if (Y[i] > 1 - 10^(-6)) { Y[i] <- 1 - 10^(-6) diff --git a/R/zibr.R b/R/zibr.R index 6bed496..52e799e 100755 --- a/R/zibr.R +++ b/R/zibr.R @@ -20,14 +20,14 @@ #' \dontrun{ #' ## simulate some data #' sim <- simulate_zero_inflated_beta_random_effect_data( -#' subject.n = 100, time.n = 5, +#' subject_n = 100, time_n = 5, #' X = as.matrix(c(rep(0, 50 * 5), rep(1, 50 * 5))), #' Z = as.matrix(c(rep(0, 50 * 5), rep(1, 50 * 5))), #' alpha = as.matrix(c(-0.5, 1)), #' beta = as.matrix(c(-0.5, 0.5)), #' s1 = 1, s2 = 0.8, #' v = 5, -#' sim.seed = 100 +#' sim_seed = 100 #' ) #' ## run zibr on the simulated data #' zibr.fit <- zibr( diff --git a/man/simulate_zero_inflated_beta_random_effect_data.Rd b/man/simulate_zero_inflated_beta_random_effect_data.Rd index c910aa7..5d76c67 100755 --- a/man/simulate_zero_inflated_beta_random_effect_data.Rd +++ b/man/simulate_zero_inflated_beta_random_effect_data.Rd @@ -5,8 +5,8 @@ \title{Simulate data according to zero-inflated beta random effects model} \usage{ simulate_zero_inflated_beta_random_effect_data( - subject.n = 50, - time.n = 5, + subject_n = 50, + time_n = 5, v = 2, alpha = as.matrix(c(0, 0.5, -1)), beta = as.matrix(c(-0.5, -0.5, 0.5)), @@ -14,13 +14,13 @@ simulate_zero_inflated_beta_random_effect_data( Z = NA, s1 = 0.2, s2 = 0.2, - sim.seed = 100 + sim_seed = 100 ) } \arguments{ -\item{subject.n}{number of subjects} +\item{subject_n}{number of subjects} -\item{time.n}{number of time points for each subject} +\item{time_n}{number of time points for each subject} \item{v}{the dispersion parameter in beta component} @@ -36,7 +36,7 @@ simulate_zero_inflated_beta_random_effect_data( \item{s2}{the stardard deviation of random effect in beta component} -\item{sim.seed}{the random seed} +\item{sim_seed}{the random seed} } \value{ Y the bacterial abundance generated from the model @@ -63,13 +63,13 @@ Simulate data according to zero-inflated beta random effects model \examples{ \dontrun{ simulate_zero_inflated_beta_random_effect_data( - subject.n = 100, time.n = 5, + subject_n = 100, time_n = 5, X = as.matrix(c(rep(0, 50 * 5), rep(1, 50 * 5))), alpha = as.matrix(c(-0.5, 1)), beta = as.matrix(c(-0.5, 0.5)), s1 = 1, s2 = 0.8, v = 5, - sim.seed = 100 + sim_seed = 100 ) } } diff --git a/man/zibr.Rd b/man/zibr.Rd index f47533c..e8d04f9 100755 --- a/man/zibr.Rd +++ b/man/zibr.Rd @@ -54,14 +54,14 @@ Fit zero-inflated beta regression with random effects \dontrun{ ## simulate some data sim <- simulate_zero_inflated_beta_random_effect_data( - subject.n = 100, time.n = 5, + subject_n = 100, time_n = 5, X = as.matrix(c(rep(0, 50 * 5), rep(1, 50 * 5))), Z = as.matrix(c(rep(0, 50 * 5), rep(1, 50 * 5))), alpha = as.matrix(c(-0.5, 1)), beta = as.matrix(c(-0.5, 0.5)), s1 = 1, s2 = 0.8, v = 5, - sim.seed = 100 + sim_seed = 100 ) ## run zibr on the simulated data zibr.fit <- zibr( diff --git a/tests/testthat/test-zibr.R b/tests/testthat/test-zibr.R index c36f4cb..7ba7a82 100755 --- a/tests/testthat/test-zibr.R +++ b/tests/testthat/test-zibr.R @@ -15,14 +15,14 @@ test_that("zibr main function works on package data", { set.seed(19683) sim <- simulate_zero_inflated_beta_random_effect_data( - subject.n = 100, time.n = 5, + subject_n = 100, time_n = 5, X = as.matrix(c(rep(0, 50 * 5), rep(1, 50 * 5))), Z = as.matrix(c(rep(0, 50 * 5), rep(1, 50 * 5))), alpha = as.matrix(c(-0.5, 1)), beta = as.matrix(c(-0.5, 0.5)), s1 = 1, s2 = 0.8, v = 5, - sim.seed = 100 + sim_seed = 100 ) sim_expected <- list(logistic.est.table = structure(c(-0.571900346260645, 0.828072571379781, 0.0069577928092327, 0.00547006539016526), .Dim = c(2L, 2L), .Dimnames = list(c("intersept", "var1"), c("Estimate", "Pvalue"))), logistic.s1.est = 1.06801391077711, beta.est.table = structure(c(-0.593090577296489, 0.591745582740125, 4.18162741963046e-05, 0.00169945318503251), .Dim = c(2L, 2L), .Dimnames = list(c("intersept", "var1"), c("Estimate", "Pvalue"))), beta.s2.est = 0.630386186272789, beta.v.est = 4.73406414312415, loglikelihood = NA, joint.p = NA) diff --git a/vignettes/zibr.Rmd b/vignettes/zibr.Rmd index dd19b17..abf1bbf 100755 --- a/vignettes/zibr.Rmd +++ b/vignettes/zibr.Rmd @@ -17,7 +17,6 @@ knitr::opts_chunk$set( ```{r setup} library(ZIBR) library(dplyr) -library(nlme) set.seed(19683) ``` @@ -41,14 +40,14 @@ The following function will simulate some data according to the zero-inflated be ```{r} sim <- simulate_zero_inflated_beta_random_effect_data( - subject.n = 100, time.n = 5, + subject_n = 100, time_n = 5, X = as.matrix(c(rep(0, 50 * 5), rep(1, 50 * 5))), Z = as.matrix(c(rep(0, 50 * 5), rep(1, 50 * 5))), alpha = as.matrix(c(-0.5, 1)), beta = as.matrix(c(-0.5, 0.5)), s1 = 1, s2 = 0.8, v = 5, - sim.seed = 100 + sim_seed = 100 ) names(sim) ```