inst/doc/simulation_functions.R

In mlpwr: A Power Analysis Toolbox to Find Cost-Efficient Study Designs

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  tidy = TRUE,
  tidy.opts=list(arrow=TRUE,width.cutoff = 50),
  eval=F
)

## -----------------------------------------------------------------------------
#  simfun <- function(N) {
#      # Generate a data set
#      # Test the hypothesis
#  }

## -----------------------------------------------------------------------------
#  library(mlpwr)

## -----------------------------------------------------------------------------
#  simfun_ttest <- function(N) {
#      # Generate a data set
#      dat <- rnorm(n = N, mean = 0.3)
#      # Test the hypothesis
#      res <- t.test(dat)
#      res$p.value < 0.01
#  }

## ----eval = F-----------------------------------------------------------------
#   res <- find.design(simfun = simfun_ttest,
#       boundaries = c(100,300), power = .95)

## -----------------------------------------------------------------------------
#  library(mlpwr)

## -----------------------------------------------------------------------------
#  simfun_anova <- function(n, n.groups) {
#  
#      # Generate a data set
#      groupmeans <- rnorm(n.groups, sd = 0.2)  # generate groupmeans using cohen's f=.2
#      dat <- sapply(groupmeans, function(x) rnorm(n,
#          mean = x, sd = 1))  # generate data
#      dat <- dat |>
#          as.data.frame() |>
#          gather()  # format
#  
#      # Test the hypothesis
#      res <- aov(value ~ key, data = dat)  # perform ANOVA
#      summary(res)[[1]][1, 5] < 0.01  # extract significance
#  }
#  

## ----eval = F-----------------------------------------------------------------
#  res <- find.design(simfun = simfun_anova,
#     costfun = function(n,n.groups) 5*n+20*n.groups,
#     boundaries = list(n = c(10, 150), n.groups = c(5, 30)),
#     power = .95)

## -----------------------------------------------------------------------------
#  library(mlpwr)

## -----------------------------------------------------------------------------
#  dat.original <- data.frame(
#    counts = c(18, 17, 15, 20,
#      10, 20, 25, 13, 12),
#    treatment = gl(3, 1, 9),
#    outcome = gl(3, 3))
#  mod.original <- glm(counts ~ outcome + treatment, data = dat.original,
#      family = poisson) # setting up the generalized linear model
#  summary(mod.original)

## -----------------------------------------------------------------------------
#  simfun_glm1 <- function(N) {
#  
#      # generate data
#      dat <- data.frame(outcome = gl(3, 1, ceiling(N/3)),
#          treatment = gl(3, ceiling(N/3)))[1:N, ] # predictors
#      a <- predict(mod.original, newdata = dat, type = "response") # criterion 'raw'
#      dat$counts <- rpois(N, a) # criterion applying poisson distribution
#  
#      # test hypothesis
#      mod <- glm(counts ~ outcome + treatment, data = dat,
#          family = poisson) # fit a glm
#      summary(mod)$coefficients["treatment2", "Pr(>|z|)"] <
#          0.01 # test the coefficient of the treatment
#  }

## ----eval = F-----------------------------------------------------------------
#  res <- find.design(simfun = simfun_glm1,
#       boundaries = c(20,100), power = .95)

## -----------------------------------------------------------------------------
#  logistic <- function(x) 1/(1 + exp(-x)) # logistic function

## -----------------------------------------------------------------------------
#  simfun_glm2 <- function(N) {
#  
#      # generate data
#      dat <- data.frame(pred1 = rnorm(N), pred2 = rnorm(N))
#      beta <- c(1.2, 0.8)  # parameter weights
#      prob <- logistic(as.matrix(dat) %*% beta)  # get probability
#      dat$criterion <- runif(N) < prob  # draw according to probability
#  
#      # test hypothesis
#      mod <- glm(criterion ~ pred1 + pred2, data = dat,
#          family = binomial) # fit a glm
#      summary(mod)$coefficients["pred2", "Pr(>|z|)"] <
#          0.01 # test the coefficient of the predictor
#  }

## ----eval = F-----------------------------------------------------------------
#  res <- find.design(simfun = simfun_glm2,
#       boundaries = c(90,200), power = .95)

## -----------------------------------------------------------------------------
#  library(mlpwr)
#  library(mirt)

## -----------------------------------------------------------------------------
#  simfun_irt1 <- function(N) {
#  
#      # generate data
#      dat <- simdata(a = c(1.04, 1.2, 1.19, 0.61, 1.31,
#          0.83, 1.46, 1.27, 0.51, 0.81), d = c(0.06,
#          -1.79, -1.15, 0.88, -0.2, -1.87, 1.23, -0.08,
#          -0.71, 0.6), N = N, itemtype = "2PL") # uses a 2PL model with a and d parameters
#  
#      # test hypothesis
#      mod <- mirt(dat)  # Fit 2PL Model
#      constrained <- "F = 1-4
#            CONSTRAIN = (1-4, a1)" # specifying that the slopes should be kept equal for items 1 to 4
#      mod_constrained <- mirt(dat, constrained)  # Fit 2PL with equal slopes
#  
#      res <- anova(mod_constrained, mod)  # perform model comparison
#      res$p[2] < 0.01  # extract significance
#  }

## ----eval = F-----------------------------------------------------------------
#  res <- find.design(simfun = simfun_irt1,
#       boundaries = c(100,500), power = .95,evaluations =500)

## -----------------------------------------------------------------------------
#  costfun_irt2 <- function(N1, N2) 5 * N1 + 7 * N2 # specifying a cost function

## -----------------------------------------------------------------------------
#  simfun_irt2 <- function(N1, N2) {
#  
#      # generate data
#      a1 <- a2 <- c(1.04, 1.2, 1.19, 0.61, 1.31, 0.83,
#          1.46, 1.27, 0.51, 0.81) # specifying the slope
#      d1 <- d2 <- c(0.06, -1.79, -1.15, 0.88, -0.2, -1.87,
#          1.23, -0.08, -0.71, 0.6) # specifying the difficulty
#      a2[1] <- a2[1] + 0.3 # the slope is different for the first item in group2
#      d2[1] <- d2[1] + 0.5 # the difficulty is different for the first item in group2
#      dat1 <- simdata(a = a1, d = d1, N = N1, itemtype = "2PL") # creating artificial data for both groups
#      dat2 <- simdata(a = a2, d = d2, N = N2, itemtype = "2PL")
#      dat <- as.data.frame(rbind(dat1, dat2)) # combining the data sets into one object
#      group <- c(rep("1", N1), rep("2", N2)) # create a variable that indicates the group membership
#  
#      # fit models
#      mod1 <- multipleGroup(dat, 1, group = group) # fit model with different parameters for each group
#      mod2 <- multipleGroup(dat, 1, group = group, invariance = c("slopes",
#          "intercepts")) # fit model with the same parameters for each group
#  
#      # test hypothesis
#      res <- anova(mod2, mod1)
#  
#      # extract significance
#      res$p[2] < 0.01
#  }

## ----eval = F-----------------------------------------------------------------
#  res <- find.design(simfun = simfun_irt2,
#       boundaries = list(N1 = c(100,700), N2 = c(100,700)),
#       costfun = costfun_irt2,
#       power = .95)

## -----------------------------------------------------------------------------
#  library(mlpwr)
#  library(lme4)
#  library(lmerTest)

## -----------------------------------------------------------------------------
#  simfun_multi1 <- function(N) {
#  
#      # generate data
#      params <- list(theta = 0.5, beta = c(2, -0.2, -0.4,
#          -0.6)) # specifying the standard deviation of the random effects and parameter weights
#      dat <- expand.grid(herd = 1:ceiling(N/4), period = factor(1:4))[1:N,
#          ] # creating predictors
#      dat$x <- simulate(~period + (1 | herd), newdata = dat,
#          family = poisson, newparams = params)[[1]] # creating criterion
#  
#      # test hypothesis
#      mod <- glmer(x ~ period + (1 | herd), data = dat,
#          family = poisson)  # fit model
#      pvalues <- summary(mod)[["coefficients"]][2:4,
#          "Pr(>|z|)"] # extract p-values
#      any(pvalues < 0.01) # test hypothes that any is significant
#  }

## ----eval = F-----------------------------------------------------------------
#  res <- find.design(simfun = simfun_multi1,
#       boundaries = c(100, 500),
#       power = .95)

## -----------------------------------------------------------------------------
#  logistic <- function(x) 1/(1 + exp(-x))
#  
#  N.original <- 300
#  n.countries.original <- 20
#  
#  # generate original data
#  dat.original <- data.frame(country = rep(1:n.countries.original,
#      length.out = N.original), pred1 = rnorm(N.original),
#      pred2 = rnorm(N.original)) # creating predictors
#  country.intercepts <- rnorm(n.countries.original, sd = 0.5) # creating random intercepts
#  dat.original$intercepts <- country.intercepts[dat.original$country] # add interecepts to data
#  beta <- c(1, 0.4, -0.3)  # parameter weights
#  prob <- logistic(as.matrix(dat.original[c("intercepts",
#      "pred1", "pred2")]) %*% as.matrix(beta))  # get probability
#  dat.original$criterion <- runif(N.original) < prob  # draw according to probability
#  
#  # fit original model to obtain parameters
#  mod.original <- glmer(criterion ~ pred1 + pred2 + 0 +
#      (1 | country), data = dat.original, family = binomial)

## -----------------------------------------------------------------------------
#  simfun_multi2 <- function(n, n.countries) {
#  
#      # generate data
#      dat <- data.frame(country = rep(1:n.countries,
#          length.out = n * n.countries), pred1 = rnorm(n *
#          n.countries), pred2 = rnorm(n * n.countries))
#      dat$criterion <- simulate(mod.original, nsim = 1,
#          newdata = dat, allow.new.levels = TRUE, use.u = FALSE) |>
#          unlist()  # criterion data from the fitted model
#  
#      # test hypothesis
#      mod <- glmer(criterion ~ pred1 + pred2 + 0 + (1 |
#          country), data = dat, family = binomial)
#      summary(mod)[["coefficients"]]["pred2", "Pr(>|z|)"] <
#          0.01 # check if significant
#  }

## -----------------------------------------------------------------------------
#  costfun_multi2 <- function(n, n.countries) 5 * m +
#      100 * n.countries
#  

## ----eval = F-----------------------------------------------------------------
#  res <- find.design(simfun = simfun_multi2,
#       boundaries = list(n=c(10,40),n.countries=c(5,20)),
#       costfun = costfun_multi2,
#       power = .95)