Skip to content

Commit

Permalink
Update formatting
Browse files Browse the repository at this point in the history
  • Loading branch information
@Ulthran
Ulthran committed Sep 21, 2023
1 parent bdc4fc2 commit 7ce6af2
Show file tree
Hide file tree
Showing 7 changed files with 58 additions and 64 deletions.
41 changes: 18 additions & 23 deletions R/simulate_beta_data.R
View file
Original file line number Diff line number Diff line change
@@ -1,39 +1,37 @@
#' @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))) {
Z <- as.matrix(Z)
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)
}
Expand All @@ -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')
Expand All @@ -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)
}
48 changes: 24 additions & 24 deletions R/simulate_zero_inflated_beta_data.R
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,15 @@
#' 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
#' @param X the covariates in logistic component
#' @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
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -78,24 +78,24 @@ 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
}


######
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) {
Y[i] <- 0
}
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)
Expand Down
4 changes: 2 additions & 2 deletions R/zibr.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
16 changes: 8 additions & 8 deletions man/simulate_zero_inflated_beta_random_effect_data.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

4 changes: 2 additions & 2 deletions man/zibr.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

4 changes: 2 additions & 2 deletions tests/testthat/test-zibr.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
5 changes: 2 additions & 3 deletions vignettes/zibr.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,6 @@ knitr::opts_chunk$set(
```{r setup}
library(ZIBR)
library(dplyr)
library(nlme)
set.seed(19683)
```
Expand All @@ -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)
```
Expand Down

0 comments on commit 7ce6af2

Please sign in to comment.