## ----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)