Skip to content
Snippets Groups Projects
Commit d66c2e84 authored by Boulakis Paradeisios Alexandros's avatar Boulakis Paradeisios Alexandros
Browse files

solutions

parent c14729dc
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 803 additions and 44 deletions
.Rhistory 0 → 100644
+ 512
− 0
View file @ d66c2e84
??pwr
sesoi.eff <- pwr.t.test(n=previous.sample,
sig.level=alpha.level,
power = 1/3,
type="two.sample",
alternative="two.sided" )
# Packages ---------------------------------------------------------------------
list.of.packages <- c("pwr")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(pwr)
# Find SESOI -------------------------------------------------------------------
previous.sample <- 25
alpha.level <- .01
beta.level <- .05
# find the SESOI
# ask for help using ??pwr
sesoi.eff <- pwr.t.test(n=previous.sample,
sig.level=alpha.level,
power = 1/3,
type="two.sample",
alternative="two.sided" )
sesoi.eff
sesoi.eff$n
pwr.t.test(sig.level=alpha.level,
power = beta.level,
d=sesoi.eff$d
type="two.sample",
pwr.t.test(sig.level=alpha.level,
power = beta.level,
d=sesoi.eff$d,
type="two.sample",
alternative="two.sided" )
pwr.t.test(sig.level=alpha.level,
power = 1- beta.level,
d=sesoi.eff$d,
type="two.sample",
alternative="two.sided" )
# Simulation parameters --------------------------------------------------------
n.sim <- 500
cohen.d <- .6
sd <- 2
alpha.error <- .01
beta.error <- .05
sample.sizes <- seq(10,200,5)
numeric(len(sample.sizes))
storage.power <- numeric(length(sample.sizes))
storage.power
?rnorm
?t.test
set.seed(2024)
# Simulation parameters --------------------------------------------------------
n.sim <- 500
cohen.d <- .6
sd <- 2
alpha.error <- .01
beta.error <- .05
sample.sizes <- seq(10,200,5)
# Simulation ------------------------------------------------------------------
storage.power <- numeric(length(sample.sizes))
#for every sample size
for (sample.size in sample.sizes){
# create a storage vector for pvals
storage.statistic <- numeric(n.sim)
# create a storage vector for statistic
storage.pvals <- numeric(n.sim)
# for every simulation
for (simulation in seq(n.sim)){
control.data <- rnorm(n=sample.size,mean = 0, sd= sd)
smoker.data <- rnorm(n=sample.size,mean = cohen.d, sd= sd)
results <- t.test(control.data,smoker.data)
}
}
results
control.data
smoker.data
results$statistic
results$p.value
# for every simulation
for (simulation in seq(n.sim)){
# simulate random normally distributed samples for group 1 and add noise
control.data <- rnorm(n=sample.size,mean = 0, sd= sd)
# simulate random normally distributed samples for group 2, add your effect and noise
smoker.data <- rnorm(n=sample.size,mean = cohen.d, sd= sd)
# do an independent t-test
results <- t.test(control.data,smoker.data)
# store statistic and pvals
storage.statistic[simulation] <- results$statistic
storage.pvals[simulation] <- results$p.value
}
storage.pvals
storage.pvals
storage.pvals < alpha.error
mean(storage.pvals < alpha.error)
set.seed(2024)
# Simulation parameters --------------------------------------------------------
n.sim <- 500
cohen.d <- .6
sd <- 2
alpha.error <- .01
beta.error <- .05
sample.sizes <- seq(10,200,5)
# Simulation ------------------------------------------------------------------
storage.power <- numeric(length(sample.sizes))
#for every sample size
for (ii in seq(length(sample.sizes)){
set.seed(2024)
# Simulation parameters --------------------------------------------------------
n.sim <- 500
cohen.d <- .6
sd <- 2
alpha.error <- .01
beta.error <- .05
sample.sizes <- seq(10,200,5)
# Simulation ------------------------------------------------------------------
storage.power <- numeric(length(sample.sizes))
#for every sample size
for (ii in seq(length(sample.sizes))){
# create a storage vector for pvals
storage.statistic <- numeric(n.sim)
# create a storage vector for statistic
storage.pvals <- numeric(n.sim)
# for every simulation
for (simulation in seq(n.sim)){
# simulate random normally distributed samples for group 1 and add noise
control.data <- rnorm(n=sample.sizes[ii],mean = 0, sd= sd)
# simulate random normally distributed samples for group 2, add your effect and noise
smoker.data <- rnorm(n=sample.sizes[ii],mean = cohen.d, sd= sd)
# do an independent t-test
results <- t.test(control.data,smoker.data)
# store statistic and pvals
storage.statistic[simulation] <- results$statistic
storage.pvals[simulation] <- results$p.value
}
# count how many results where significant
achieved.power <- mean(storage.pvals < alpha.error)
storage.power[ii] <- achieved.power
}
storage.power
set.seed(2024)
# Simulation parameters --------------------------------------------------------
n.sim <- 500
cohen.d <- .6
sd <- 2
alpha.error <- .01
beta.error <- .05
sample.sizes <- seq(10,200,5)
# Simulation ------------------------------------------------------------------
storage.power <- numeric(length(sample.sizes))
#for every sample size
for (ii in seq(length(sample.sizes))){
print(ii,sample.sizes[ii])
# create a storage vector for pvals
storage.statistic <- numeric(n.sim)
# create a storage vector for statistic
storage.pvals <- numeric(n.sim)
# for every simulation
for (simulation in seq(n.sim)){
# simulate random normally distributed samples for group 1 and add noise
control.data <- rnorm(n=sample.sizes[ii],mean = 0, sd= sd)
# simulate random normally distributed samples for group 2, add your effect and noise
smoker.data <- rnorm(n=sample.sizes[ii],mean = cohen.d, sd= sd)
# do an independent t-test
results <- t.test(control.data,smoker.data)
# store statistic and pvals
storage.statistic[simulation] <- results$statistic
storage.pvals[simulation] <- results$p.value
}
# count how many results where significant
achieved.power <- mean(storage.pvals < alpha.error)
storage.power[ii] <- achieved.power
}
sample.sizes
sample.sizes[ii]
set.seed(2024)
# Simulation parameters --------------------------------------------------------
n.sim <- 500
cohen.d <- .6
sd <- 2
alpha.error <- .01
beta.error <- .05
sample.sizes <- seq(10,200,5)
# Simulation ------------------------------------------------------------------
storage.power <- numeric(length(sample.sizes))
#for every sample size
for (ii in seq(length(sample.sizes))){
print(sample.sizes[ii])
# create a storage vector for pvals
storage.statistic <- numeric(n.sim)
# create a storage vector for statistic
storage.pvals <- numeric(n.sim)
# for every simulation
for (simulation in seq(n.sim)){
# simulate random normally distributed samples for group 1 and add noise
control.data <- rnorm(n=sample.sizes[ii],mean = 0, sd= sd)
# simulate random normally distributed samples for group 2, add your effect and noise
smoker.data <- rnorm(n=sample.sizes[ii],mean = cohen.d, sd= sd)
# do an independent t-test
results <- t.test(control.data,smoker.data)
# store statistic and pvals
storage.statistic[simulation] <- results$statistic
storage.pvals[simulation] <- results$p.value
}
# count how many results where significant
achieved.power <- mean(storage.pvals < alpha.error)
storage.power[ii] <- achieved.power
}
storage.power
set.seed(2024)
# Simulation parameters --------------------------------------------------------
n.sim <- 500
cohen.d <- .6
sd <- 1
alpha.error <- .01
beta.error <- .05
sample.sizes <- seq(10,200,5)
# Simulation ------------------------------------------------------------------
storage.power <- numeric(length(sample.sizes))
#for every sample size
for (ii in seq(length(sample.sizes))){
print(sample.sizes[ii])
# create a storage vector for pvals
storage.statistic <- numeric(n.sim)
# create a storage vector for statistic
storage.pvals <- numeric(n.sim)
# for every simulation
for (simulation in seq(n.sim)){
# simulate random normally distributed samples for group 1 and add noise
control.data <- rnorm(n=sample.sizes[ii],mean = 0, sd= sd)
# simulate random normally distributed samples for group 2, add your effect and noise
smoker.data <- rnorm(n=sample.sizes[ii],mean = cohen.d, sd= sd)
# do an independent t-test
results <- t.test(control.data,smoker.data)
# store statistic and pvals
storage.statistic[simulation] <- results$statistic
storage.pvals[simulation] <- results$p.value
}
# count how many results where significant
achieved.power <- mean(storage.pvals < alpha.error)
storage.power[ii] <- achieved.power
}
storage.power
# Packages ---------------------------------------------------------------------
list.of.packages <- c("ggplot2")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(ggplot2)
set.seed(2024)
# Simulation parameters --------------------------------------------------------
n.sim <- 500
cohen.d <- .6
sd <- 1
alpha.error <- .01
beta.error <- .05
sample.sizes <- seq(10,200,5)
# Simulation ------------------------------------------------------------------
storage.power <- numeric(length(sample.sizes))
#for every sample size
for (ii in seq(length(sample.sizes))){
print(sample.sizes[ii])
# create a storage vector for pvals
storage.statistic <- numeric(n.sim)
# create a storage vector for statistic
storage.pvals <- numeric(n.sim)
# for every simulation
for (simulation in seq(n.sim)){
# simulate random normally distributed samples for group 1 and add noise
control.data <- rnorm(n=sample.sizes[ii],mean = 0, sd= sd)
# simulate random normally distributed samples for group 2, add your effect and noise
smoker.data <- rnorm(n=sample.sizes[ii],mean = cohen.d, sd= sd)
# do an independent t-test
results <- t.test(control.data,smoker.data)
# store statistic and pvals
storage.statistic[simulation] <- results$statistic
storage.pvals[simulation] <- results$p.value
}
# count how many results where significant
achieved.power <- mean(storage.pvals < alpha.error)
storage.power[ii] <- achieved.power
}
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()+
geom_point()
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()+
geom_point()+
labs(
title = "Simulation for Power Analysis",
x = "Sample Size",
y = "Achieved Power",
) +
theme_minimal()
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()+
geom_point()+
labs(
title = "Simulation for Power Analysis",
x = "Sample Size",
y = "Achieved Power",
)
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()+
geom_point()+
geom_vline(xintercept=1-beta.error)
labs(
title = "Simulation for Power Analysis",
x = "Sample Size",
y = "Achieved Power",
)
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()+
geom_point()+
geom_hline(xintercept=1-beta.error)
labs(
title = "Simulation for Power Analysis",
x = "Sample Size",
y = "Achieved Power",
)
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()+
geom_point()+
geom_hline(yintercept=1-beta.error)
labs(
title = "Simulation for Power Analysis",
x = "Sample Size",
y = "Achieved Power",
)
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()+
geom_point()+
geom_hline(yintercept=1-beta.error,color = 'red', linetype = "dashed", size = 1)
labs(
title = "Simulation for Power Analysis",
x = "Sample Size",
y = "Achieved Power",
)
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()+
geom_point()+
geom_hline(yintercept=1-beta.error,color = 'red', linetype = "dashed", linewidth = 1) +
labs(
title = "Simulation for Power Analysis",
x = "Sample Size",
y = "Achieved Power",
)
target.sample <- sample.sizes[storage.power > 1-beta.error][1]
target.sample
# Packages ---------------------------------------------------------------------
list.of.packages <- c("ggplot2")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(ggplot2)
set.seed(2024)
# Simulation parameters --------------------------------------------------------
n.sim <- 500
cohen.d <- .6
sd <- 1
alpha.error <- .01
beta.error <- .05
sample.sizes <- seq(10,200,5)
# Simulation ------------------------------------------------------------------
storage.power <- numeric(length(sample.sizes))
#for every sample size
for (ii in seq(length(sample.sizes))){
print(sample.sizes[ii])
# create a storage vector for pvals
storage.statistic <- numeric(n.sim)
# create a storage vector for statistic
storage.pvals <- numeric(n.sim)
# for every simulation
for (simulation in seq(n.sim)){
# simulate random normally distributed samples for group 1 and add noise
control.data <- rnorm(n=sample.sizes[ii],mean = 0, sd= sd)
# simulate random normally distributed samples for group 2, add your effect and noise
smoker.data <- rnorm(n=sample.sizes[ii],mean = cohen.d, sd= sd)
# do an independent t-test
results <- t.test(control.data,smoker.data)
# store statistic and pvals
storage.statistic[simulation] <- results$statistic
storage.pvals[simulation] <- results$p.value
}
# count how many results where significant
achieved.power <- mean(storage.pvals < alpha.error)
storage.power[ii] <- achieved.power
}
# print when your results was above your threshold
target.sample <- sample.sizes[storage.power > 1-beta.error][1]
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()+
geom_point()+
geom_hline(yintercept=1-beta.error,color = 'red', linetype = "dashed", linewidth = 1) +
labs(
title = sprintf("Simulation for Power Analysis: Required sample size = %s",target.sample),
x = "Sample Size",
y = "Achieved Power",
)
# Packages ---------------------------------------------------------------------
list.of.packages <- c("ggplot2")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(ggplot2)
set.seed(2024)
# Simulation parameters --------------------------------------------------------
n.sim <- 5000
cohen.d <- .6
sd <- 1
alpha.error <- .01
beta.error <- .05
sample.sizes <- seq(10,200,5)
# Simulation ------------------------------------------------------------------
storage.power <- numeric(length(sample.sizes))
#for every sample size
for (ii in seq(length(sample.sizes))){
print(sample.sizes[ii])
# create a storage vector for pvals
storage.statistic <- numeric(n.sim)
# create a storage vector for statistic
storage.pvals <- numeric(n.sim)
# for every simulation
for (simulation in seq(n.sim)){
# simulate random normally distributed samples for group 1 and add noise
control.data <- rnorm(n=sample.sizes[ii],mean = 0, sd= sd)
# simulate random normally distributed samples for group 2, add your effect and noise
smoker.data <- rnorm(n=sample.sizes[ii],mean = cohen.d, sd= sd)
# do an independent t-test
results <- t.test(control.data,smoker.data)
# store statistic and pvals
storage.statistic[simulation] <- results$statistic
storage.pvals[simulation] <- results$p.value
}
# count how many results where significant
achieved.power <- mean(storage.pvals < alpha.error)
storage.power[ii] <- achieved.power
}
# print when your results was above your threshold
target.sample <- sample.sizes[storage.power > 1-beta.error][1]
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()+
geom_point()+
geom_hline(yintercept=1-beta.error,color = 'red', linetype = "dashed", linewidth = 1) +
labs(
title = sprintf("Simulation for Power Analysis: Required sample size = %s",target.sample),
x = "Sample Size",
y = "Achieved Power",
)
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()+
geom_point()+
geom_hline(yintercept=1-beta.error,color = 'red', linetype = "dashed", linewidth = 1) +
geom_hline(xintercept=target.sample,color = 'red', linetype = "dashed", linewidth = 1) +
labs(
title = sprintf("Simulation for Power Analysis: Required sample size = %s",target.sample),
x = "Sample Size",
y = "Achieved Power",
)
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()+
geom_point()+
geom_hline(yintercept=1-beta.error,color = 'red', linetype = "dashed", linewidth = 1) +
geom_vline(xintercept=target.sample,color = 'red', linetype = "dashed", linewidth = 1) +
labs(
title = sprintf("Simulation for Power Analysis: Required sample size = %s",target.sample),
x = "Sample Size",
y = "Achieved Power",
)
# Packages ---------------------------------------------------------------------
list.of.packages <- c("pwr","rpact")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(pwr)
library(rpact)
??pwr
??pwr
sesoi.eff <- pwr.t.test(n=previous.sample.cont + previous.sample.depr,
sig.level=alpha.level,
power = 1/3,
type="two.sample",
alternative="two.sided" )
previous.sample.cont <- 50
previous.sample.depr <- 44
alpha.error <- .01
beta.error <- .05
sesoi.eff <- pwr.t.test(n=previous.sample.cont + previous.sample.depr,
sig.level=alpha.level,
power = 1/3,
type="two.sample",
alternative="two.sided" )
sesoi.eff
sesoi.eff <- pwr.t.test(n=50,
sig.level=alpha.level,
power = 1/3,
type="two.sample",
alternative="two.sided" )
sesoi.eff
power <- pwr.t.test(sig.level=alpha.level,
power = 1- beta.level,
d=sesoi.eff$d,
type="two.sample",
alternative="two.sided" )
power
??rpact
design <- getDesignGroupSequential(
kMax = 4,
typeOfDesign = "P",
sided = 2,
alpha = alpha.error,
beta = beta.error
)
print(summary(design))
?getDesignGroupSequential
ex_01_params.R
+ 28
− 22
View file @ d66c2e84
# Packages ---------------------------------------------------------------------
list.of.packages <- c("pwrss")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(pwrss)
########################################################
# parameters for independent t-test
control.group.mean <- 130
control.group.sd <- 15
drug.group.mean <- 120
drug.group.sd <- 15
# Independent t-test -----------------------------------------------------------
# parameters
cont.mean <- 130
drug.mean <- 120
alpha.level <- .05
required.power <- .95
cont.sd <- 15
drug.sd <- 15
alternative <- "not equal" # also try "greater", "less"
alpha.error <- .05
beta.error <- .05
alternative <- "not equal" # also try "greater" or "less"
########################################################
# find the correct function
# assign the correct parameters
# you can look around the documentation of pwrss library using ??pwrss
t.test.power <-
t.test.power <-
plot(t.test.power)
#######################################################
# Paired t-test ----------------------------------------------------------------
# parameters for paired-sample t-test
drug.before.mean <- 130
drug.before.sd <- 15
# parameters
before.mean <- 130
after.mean <- 120
drug.after.mean <- 120
drug.after.sd <- 15
before.sd <- 15
after.sd <- 15
alpha.level <- .05
required.power <- .95
alpha.error <- .05
beta.error <- .05
alternative <- "not equal" # also try "greater", "less"
alternative <- "not equal" # also try "greater" or "less"
########################################################
# find the correct function
# assign the correct parameters
# you can look around the documentation of pwrss library using ??pwrss
paired.t.power <-
plot(paired.t.power)
\ No newline at end of file
plot(paired.t.power)
ex_01_params_solution.R 0 → 100644
+ 61
− 0
View file @ d66c2e84
# Packages ---------------------------------------------------------------------
list.of.packages <- c("pwrss")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(pwrss)
# Independent t-test -----------------------------------------------------------
# parameters
cont.mean <- 130
drug.mean <- 120
cont.sd <- 15
drug.sd <- 15
alpha.error <- .05
beta.error <- .05
alternative <- "not equal"
t.test.power <- pwrss.t.2means(mu1 = cont.mean,
mu2 = drug.mean,
sd1 = cont.sd,
sd2 = cont.sd,
alpha = alpha.error,
power = 1-beta.error,
alternative = alternative,
kappa = 1)
plot(t.test.power)
# Paired t-test ----------------------------------------------------------------
# parameters
before.mean <- 130
after.mean <- 120
before.sd <- 15
after.sd <- 15
alpha.error <- .05
beta.error <- .05
alternative <- "not equal"
paired.t.power <- pwrss.t.2means(mu1 = before.mean,
mu2 = after.mean,
sd1 = before.sd,
sd2 = after.sd,
alpha = alpha.error,
power = 1-beta.error,
alternative = alternative,
paired=TRUE)
plot(paired.t.power)
ex_02_sesoi.R
+ 10
− 6
View file @ d66c2e84
# Packages ---------------------------------------------------------------------
list.of.packages <- c("pwr")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
###############################################
previous.sample <-
alpha.level <-
beta.level <-
library(pwr)
# Find SESOI -------------------------------------------------------------------
previous.sample <- 25
alpha.level <- .01
beta.level <- .05
###############################################
# find the SESOI
# ask for help using ??pwr
sesoi.eff <-
# Find power -------------------------------------------------------------------
############################################
# run a power analysis using the SESOI
# ask for help using ??pwr
ex_02_sesoi_solution.R 0 → 100644
+ 31
− 0
View file @ d66c2e84
# Packages ---------------------------------------------------------------------
list.of.packages <- c("pwr")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(pwr)
# Find SESOI -------------------------------------------------------------------
previous.sample <- 25
alpha.level <- .01
beta.level <- .05
# find the SESOI
# ask for help using ??pwr
sesoi.eff <- pwr.t.test(n=previous.sample,
sig.level=alpha.level,
power = 1/3,
type="two.sample",
alternative="two.sided" )
# Find power -------------------------------------------------------------------
# run a power analysis using the SESOI
# ask for help using ??pwr
pwr.t.test(sig.level=alpha.level,
power = 1- beta.level,
d=sesoi.eff$d,
type="two.sample",
alternative="two.sided" )
\ No newline at end of file
ex_03_simulation.R
+ 27
− 15
View file @ d66c2e84
set.seed(2024)
# Packages ---------------------------------------------------------------------
list.of.packages <- c("ggplot2")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
# Specify simulation parameters
library(ggplot2)
# n_sim
# cohen_d
# random noise
# alpha_level
# beta_level
# sample_size_list
set.seed(2024)
# Simulation parameters --------------------------------------------------------
n.sim <-
cohen.d <-
sd <-
alpha.error <-
beta.error <-
sample.sizes <-
# Main loop
# Simulation ------------------------------------------------------------------
# create a storage vector for achieved power per sample size
#create storage so you track every sample size's power
# for every sample size
#for every sample size (try to iterate over the indices)
# create a storage vector for pvals
# create a storage vector for statistic
# for every simulation
# simulate random normally distributed samples for group 1 and add noise
# simulate random normally distributed samples for group 2, add your effect and noise
# do an independent t-test
# store statistic and pvals
# count how many results where significant
# plot power curve
# print when your results was above your threshold
\ No newline at end of file
# print when your results was above your threshold
# think of when you sample sizes become larger than the power you want to achieve
ex_03_simulation_solution.R 0 → 100644
+ 65
− 0
View file @ d66c2e84
# Packages ---------------------------------------------------------------------
list.of.packages <- c("ggplot2")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(ggplot2)
set.seed(2024)
# Simulation parameters --------------------------------------------------------
n.sim <- 5000
cohen.d <- .6
sd <- 1
alpha.error <- .01
beta.error <- .05
sample.sizes <- seq(10,200,5)
# Simulation ------------------------------------------------------------------
storage.power <- numeric(length(sample.sizes))
#for every sample size
for (ii in seq(length(sample.sizes))){
print(sample.sizes[ii])
# create a storage vector for pvals
storage.statistic <- numeric(n.sim)
# create a storage vector for statistic
storage.pvals <- numeric(n.sim)
# for every simulation
for (simulation in seq(n.sim)){
# simulate random normally distributed samples for group 1 and add noise
control.data <- rnorm(n=sample.sizes[ii],mean = 0, sd= sd)
# simulate random normally distributed samples for group 2, add your effect and noise
smoker.data <- rnorm(n=sample.sizes[ii],mean = cohen.d, sd= sd)
# do an independent t-test
results <- t.test(control.data,smoker.data)
# store statistic and pvals
storage.statistic[simulation] <- results$statistic
storage.pvals[simulation] <- results$p.value
}
# count how many results where significant
achieved.power <- mean(storage.pvals < alpha.error)
storage.power[ii] <- achieved.power
}
# print when your results was above your threshold
target.sample <- sample.sizes[storage.power > 1-beta.error][1]
ggplot(data=NULL,aes(x=sample.sizes,y=storage.power))+
geom_line()+
geom_point()+
geom_hline(yintercept=1-beta.error,color = 'red', linetype = "dashed", linewidth = 1) +
geom_vline(xintercept=target.sample,color = 'red', linetype = "dashed", linewidth = 1) +
labs(
title = sprintf("Simulation for Power Analysis: Required sample size = %s",target.sample),
x = "Sample Size",
y = "Achieved Power",
)
\ No newline at end of file
ex_04_seq_analysis.R 0 → 100644
+ 26
− 0
View file @ d66c2e84
# Packages ---------------------------------------------------------------------
list.of.packages <- c("pwr","rpact")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(pwr)
library(rpact)
# Find SESOI -------------------------------------------------------------------
previous.sample.cont <- 50
previous.sample.depr <- 50
alpha.error <- .01
beta.error <- .05
sesoi.eff <-
# Find power -------------------------------------------------------------------
power <-
# Find sequential analysis p values --------------------------------------------
# run a power analysis using the SESOI
# ask for help using ??rpact and especially the ?getDesignGroupSequential
ex_04_seq_analysis_solution.R 0 → 100644
+ 39
− 0
View file @ d66c2e84
# Packages ---------------------------------------------------------------------
list.of.packages <- c("pwr","rpact")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(pwr)
library(rpact)
# Find SESOI -------------------------------------------------------------------
previous.sample.cont <- 50
previous.sample.depr <- 50
alpha.error <- .01
beta.error <- .05
sesoi.eff <- pwr.t.test(n=50,
sig.level=alpha.level,
power = 1/3,
type="two.sample",
alternative="two.sided" )
# Find power -------------------------------------------------------------------
power <- pwr.t.test(sig.level=alpha.level,
power = 1- beta.level,
d=sesoi.eff$d,
type="two.sample",
alternative="two.sided" )
# Find sequential analysis p values --------------------------------------------
design <- getDesignGroupSequential(
kMax = 4,
typeOfDesign = "P",
sided = 2,
alpha = alpha.error,
beta = beta.error
)
print(summary(design))
plot_eff_against_power.R
+ 4
− 1
View file @ d66c2e84
library(ggplot2)
# make results replicable
set.seed(2024)
# specify how many simulations you want to run
nSimul <- 10000
# specify how many effects size you want to test
effects <- seq(.1,.9,.2)
effects <- seq(.1,.9,.1)
# specify how many sample sizes you want to test
samples <- seq(10,150,10)
......
power_per_effect.png 0 → 100644
+ 0
− 0
View file @ d66c2e84
power_per_effect.png

11 KiB

0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment