Skip to content

Commit

Permalink
tabs2spaces
Browse files Browse the repository at this point in the history
  • Loading branch information
@vnijs
vnijs committed Jul 1, 2017
1 parent 3fed01b commit eb88ed3
Show file tree
Hide file tree
Showing 9 changed files with 187 additions and 180 deletions.
74 changes: 37 additions & 37 deletions R/doe.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,50 +54,50 @@ doe <- function(factors, int = "", trials = NA, seed = NA) {

part_fac <- function(df, model = ~ ., int = 0, trials = NA, seed = 172110) {

full <- expand.grid(df)
full <- expand.grid(df)

###############################################
# eliminate combinations from full
# by removing then from the variable _experiment_
# http://stackoverflow.com/questions/18459311/creating-a-fractional-factorial-design-in-r-without-prohibited-pairs?rq=1
###############################################
###############################################
# eliminate combinations from full
# by removing then from the variable _experiment_
# http://stackoverflow.com/questions/18459311/creating-a-fractional-factorial-design-in-r-without-prohibited-pairs?rq=1
###############################################

levs <- sapply(df, length)
nr_levels <- sum(levs)
min_trials <- nr_levels - length(df) + 1
max_trials <- nrow(full)
levs <- sapply(df, length)
nr_levels <- sum(levs)
min_trials <- nr_levels - length(df) + 1
max_trials <- nrow(full)

## make sure the number of trials set by the user is within an appropriate range
if (!is_empty(trials))
max_trials <- min_trials <- max(min(trials, max_trials), min_trials)

## define a data.frame that will store design spec
eff <-
data.frame(
Trials = min_trials:max_trials,
"D-efficiency" = NA,
"Balanced" = NA,
check.names = FALSE
)

for (i in min_trials:max_trials) {
eff <-
data.frame(
Trials = min_trials:max_trials,
"D-efficiency" = NA,
"Balanced" = NA,
check.names = FALSE
)

for (i in min_trials:max_trials) {
seed %>% gsub("[^0-9]","",.) %>% { if (!is_empty(.)) set.seed(seed) }
design <- try(AlgDesign::optFederov(model, data = full, nRepeats = 50,
design <- try(AlgDesign::optFederov(model, data = full, nRepeats = 50,
nTrials = i, maxIteration=1000,
approximate = FALSE), silent = TRUE)

if (is(design, "try-error")) next
ind <- which(eff$Trials %in% i)
eff[ind,"D-efficiency"] <- design$Dea
eff[ind,"Balanced"] <- all(i %% levs == 0)
if (is(design, "try-error")) next
ind <- which(eff$Trials %in% i)
eff[ind, "D-efficiency"] <- design$Dea
eff[ind, "Balanced"] <- all(i %% levs == 0)

if (design$Dea == 1) break
}
if (design$Dea == 1) break
}

if (!is(design, "try-error"))
cor_mat <- sshhr(polycor::hetcor(design$design, std.err = FALSE)$correlations)

if (exists("cor_mat")) {
if (exists("cor_mat")) {
detcm <- det(cor_mat)

full <- arrange_(full, .dots = names(df)) %>%
Expand All @@ -106,17 +106,17 @@ doe <- function(factors, int = "", trials = NA, seed = NA) {
part <- arrange_(design$design, .dots = names(df)) %>%
{suppressMessages(dplyr::right_join(full, .))}

list(df = df, cor_mat = cor_mat, detcm = detcm, Dea = design$Dea,
list(df = df, cor_mat = cor_mat, detcm = detcm, Dea = design$Dea,
part = part,
full = full,
eff = na.omit(eff),
seed = seed)
} else if (!is.na(trials)) {
"No solution exists for the selected number of trials"
} else {
"No solution found"
}
}
full = full,
eff = na.omit(eff),
seed = seed)
} else if (!is.na(trials)) {
"No solution exists for the selected number of trials"
} else {
"No solution found"
}
}

part_fac(df_list, model = as.formula(model), int = nInt, trials = trials, seed = seed) %>%
add_class("doe")
Expand Down
80 changes: 42 additions & 38 deletions R/sample_size.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,23 +31,27 @@ sample_size <- function(type,
pop_correction = "no",
pop_size = 1000000) {

if (pop_correction == "yes" && is_not(pop_size)) pop_size <- 1000000
if (pop_correction == "yes" && is_not(pop_size)) pop_size <- 1000000

if (is_not(conf_lev) || conf_lev < 0) conf_lev <- 1.96
if (is_not(conf_lev) || conf_lev < 0) conf_lev <- 1.96

zval <- conf_lev
zval <- conf_lev

if (type == "mean") {
if (is_not(err_mean)) return("Please select an acceptable error greater than 0" %>% add_class("sample_size"))
n <- (zval^2 * sd_mean^2) / err_mean^2
} else {
if (is_not(err_prop)) return("Please select an acceptable error greater than 0" %>% add_class("sample_size"))
n <- (zval^2 * p_prop * (1 - p_prop)) / err_prop^2
}
if (type == "mean") {
if (is_not(err_mean))
return("Please select an acceptable error greater than 0" %>%
add_class("sample_size"))
n <- (zval^2 * sd_mean^2) / err_mean^2
} else {
if (is_not(err_prop))
return("Please select an acceptable error greater than 0" %>%
add_class("sample_size"))
n <- (zval^2 * p_prop * (1 - p_prop)) / err_prop^2
}

if (pop_correction == 'yes') n <- n * pop_size / ((n - 1) + pop_size)
if (pop_correction == 'yes') n <- n * pop_size / ((n - 1) + pop_size)

n <- ceiling(n)
n <- ceiling(n)

as.list(environment()) %>% add_class("sample_size")
}
Expand All @@ -70,32 +74,32 @@ summary.sample_size <- function(object, ...) {

if (is.character(object)) return(object)

cat("Sample size calculation\n")

if (object$type == "mean") {
cat("Calculation type : Mean\n")
cat("Acceptable Error :", object$err_mean, "\n")
cat("Sample std. deviation:", object$sd_mean, "\n")
} else {
cat("Type: Proportion\n")
cat("Acceptable Error :", object$err_prop, "\n")
cat("Sample proportion :", object$p_prop, "\n")
}

cat("Confidence level :", object$conf_lev, "\n")
cat("Incidence rate :", object$incidence, "\n")
cat("Response rate :", object$response, "\n")

if (object$pop_correction == "no") {
cat("Population correction: None\n")
} else {
cat("Population correction: Yes\n")
cat("Population size :", formatnr(object$pop_size, dec = 0), "\n")
}

cat("\nRequired sample size :", formatnr(object$n, dec = 0))
cat("\nRequired contact attempts:", formatnr(ceiling(object$n / object$incidence / object$response), dec = 0))
cat("\n\nChoose a z.value:\n")
cat("Sample size calculation\n")

if (object$type == "mean") {
cat("Calculation type : Mean\n")
cat("Acceptable Error :", object$err_mean, "\n")
cat("Sample std. deviation:", object$sd_mean, "\n")
} else {
cat("Type: Proportion\n")
cat("Acceptable Error :", object$err_prop, "\n")
cat("Sample proportion :", object$p_prop, "\n")
}

cat("Confidence level :", object$conf_lev, "\n")
cat("Incidence rate :", object$incidence, "\n")
cat("Response rate :", object$response, "\n")

if (object$pop_correction == "no") {
cat("Population correction: None\n")
} else {
cat("Population correction: Yes\n")
cat("Population size :", formatnr(object$pop_size, dec = 0), "\n")
}

cat("\nRequired sample size :", formatnr(object$n, dec = 0))
cat("\nRequired contact attempts:", formatnr(ceiling(object$n / object$incidence / object$response), dec = 0))
cat("\n\nChoose a z.value:\n")

for (z in c(.80, .85, .90, .95, .99))
cat(paste0(100*z,"%\t"),-qnorm((1 - z)/2) %>% round(2),"\n")
Expand Down
68 changes: 34 additions & 34 deletions R/sample_size_comp.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,58 +35,58 @@ sample_size_comp <- function(type,
sig.level <- if (is.null(conf_lev)) NULL else 1 - conf_lev
adj <- ifelse (alternative == "two.sided", 2, 1)

if (type == 'mean') {
if (type == 'mean') {
if (!is.null(delta) && is.na(delta)) delta <- NULL
if (!is.null(delta)) delta <- abs(delta)
if (!is.null(sd) && is.na(sd)) sd <- NULL

nr_null <- is.null(n)+is.null(power)+is.null(delta)+is.null(sd)+is.null(conf_lev)
if (nr_null == 0 || nr_null > 1)
return("Exactly one of 'Sample size', 'Delta', 'Std. deviation',\n'Confidence level', and 'Power' must be blank or NULL" %>% add_class("sample_size_comp"))
return("Exactly one of 'Sample size', 'Delta', 'Std. deviation',\n'Confidence level', and 'Power' must be blank or NULL" %>% add_class("sample_size_comp"))

res <-
power.t.test(n = n, delta = delta, sd = sd, sig.level = sig.level, power = power, alternative = alternative) %>%
tidy

## adjustment based on http://powerandsamplesize.com/Calculators/Compare-2-Means/2-Sample-Equality
if (is.null(ratio) || is.na(ratio)) ratio <- 1
## adjustment based on http://powerandsamplesize.com/Calculators/Compare-2-Means/2-Sample-Equality
if (is.null(ratio) || is.na(ratio)) ratio <- 1
n2 <-
(1 + 1 / ratio) *
(res$sd * (qnorm(1 - res$sig.level/adj) + qnorm(res$power)) /
(res$delta))^2
n2 <- ceiling(n2)
n1 <- ceiling(ratio * n2)

} else {
} else {
if (!is.null(p1) && is.na(p1)) p1 <- NULL
if (!is.null(p2) && is.na(p2)) p2 <- NULL

if (!is.null(p1) && !is.null(p2)) {
if (p1 == p2)
return("Proportion 1 and 2 should not be set equal. Please change the proportion values" %>% add_class("sample_size_comp"))
if (p1 == p2)
return("Proportion 1 and 2 should not be set equal. Please change the proportion values" %>% add_class("sample_size_comp"))
}

nr_null <- is.null(n)+is.null(power)+is.null(p1)+is.null(p2)+is.null(conf_lev)
if (nr_null == 0 || nr_null > 1)
return("Exactly one of 'Sample size', 'Proportion 1', 'Proportion 2',\n'Confidence level', and 'Power' must be blank or NULL" %>% add_class("sample_size_comp"))
return("Exactly one of 'Sample size', 'Proportion 1', 'Proportion 2',\n'Confidence level', and 'Power' must be blank or NULL" %>% add_class("sample_size_comp"))

res <-
power.prop.test(n = n, p1 = p1, p2 = p2, sig.level = sig.level, power = power, alternative = alternative) %>%
tidy

## adjustment based on http://powerandsamplesize.com/Calculators/Compare-2-Proportions/2-Sample-Equality
if (is.null(ratio) || is.na(ratio)) ratio <- 1
adj <- ifelse (alternative == "two.sided", 2, 1)
## adjustment based on http://powerandsamplesize.com/Calculators/Compare-2-Proportions/2-Sample-Equality
if (is.null(ratio) || is.na(ratio)) ratio <- 1
adj <- ifelse (alternative == "two.sided", 2, 1)
n2 <-
(res$p1 * (1 - res$p1)/ratio + res$p2 * (1 - res$p2)) *
((qnorm(1 - res$sig.level/adj) + qnorm(res$power))/(res$p1 - res$p2))^2
n2 <- ceiling(n2)
n1 <- ceiling(ratio * n2)
}
}

res$n <- formatnr(ceiling(res$n), dec = 0)
res$n1 <- formatnr(n1, dec = 0)
res$n2 <- formatnr(n2, dec = 0)
res$n <- formatnr(ceiling(res$n), dec = 0)
res$n1 <- formatnr(n1, dec = 0)
res$n2 <- formatnr(n2, dec = 0)

as.list(environment()) %>% add_class("sample_size_comp")
}
Expand All @@ -105,25 +105,25 @@ summary.sample_size_comp <- function(object, ...) {

if (is.character(object)) return(object)

cat("Sample size calculation for comparisons\n")
cat("Sample size calculation for comparisons\n")

cat("Type :", object$type, "\n")
if (object$ratio == 1) {
cat("Type :", object$type, "\n")
if (object$ratio == 1) {
cat("Sample size :", object$res$n, "\n")
} else {
cat(paste0("Sample size 1 : ", object$res$n1, " (", object$res$n2, " x ", object$ratio, ")\n"))
cat(paste0("Sample size 2 : ", object$res$n2, "\n"))
}

if (object$type == "mean") {
cat("Delta :", object$res$delta, "\n")
cat("Std. deviation :", object$res$sd, "\n")
} else {
cat("Proportion 1 :", object$res$p1, "\n")
cat("Proportion 2 :", object$res$p2, "\n")
}
cat("Confidence level:", 1 - object$res$sig.level, "\n")
cat("Power :", object$res$power, "\n")
cat("Ratio (n1 / n2) :", object$ratio, "\n")
cat("Alternative :", object$alternative, "\n\n")
} else {
cat(paste0("Sample size 1 : ", object$res$n1, " (", object$res$n2, " x ", object$ratio, ")\n"))
cat(paste0("Sample size 2 : ", object$res$n2, "\n"))
}

if (object$type == "mean") {
cat("Delta :", object$res$delta, "\n")
cat("Std. deviation :", object$res$sd, "\n")
} else {
cat("Proportion 1 :", object$res$p1, "\n")
cat("Proportion 2 :", object$res$p2, "\n")
}
cat("Confidence level:", 1 - object$res$sig.level, "\n")
cat("Power :", object$res$power, "\n")
cat("Ratio (n1 / n2) :", object$ratio, "\n")
cat("Alternative :", object$alternative, "\n\n")
}
4 changes: 2 additions & 2 deletions R/sampling.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ sampling <- function(dataset, var, sample_size,
## use seed if provided
seed %>% gsub("[^0-9]","",.) %>% { if (!is_empty(.)) set.seed(seed) }

## example list of names obtained from http://listofrandomnames.com
## example list of names obtained from http://listofrandomnames.com
# dat$rnd_number <- runif(nrow(dat), min = 0, max = 1) %>% round(3)
dat$rnd_number <- runif(nrow(dat), min = 0, max = 1)
seldat <- dat %>%
Expand Down Expand Up @@ -64,7 +64,7 @@ summary.sampling <- function(object, prn = TRUE, ...) {
cat("Random seed:", object$seed,"\n")
cat("Sample size:", object$sample_size, "\n\n")
cat("Selected:\n")
print(formatdf(object$seldat, dec = 3), row.names = FALSE)
print(formatdf(object$seldat, dec = 3), row.names = FALSE)
if (prn) {
cat("\nSampling frame:\n")
print(formatdf(object$dat, dec = 3), row.names = FALSE)
Expand Down
44 changes: 22 additions & 22 deletions build/build_radiant.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,9 +24,9 @@ if (!file.exists(dirwin)) dir.create(dirwin, recursive = TRUE)

## delete older version of radiant
rem_old <- function(app) {
unlink(paste0(dirsrc, "/", app, "*"))
unlink(paste0(dirmac, "/", app, "*"))
unlink(paste0(dirwin, "/", app, "*"))
unlink(paste0(dirsrc, "/", app, "*"))
unlink(paste0(dirmac, "/", app, "*"))
unlink(paste0(dirwin, "/", app, "*"))
}

sapply("radiant.design", rem_old)
Expand All @@ -37,24 +37,24 @@ system(paste0(Sys.which("R"), " -e \"source('~/gh/radiant.design/build/build_mac
win <- readline(prompt = "Did you build on Windows? y/n: ")
if (grepl("[yY]", win)) {

## move packages to radiant_miniCRAN. must package in Windows first
setwd("~/gh/")
sapply(list.files(".",pattern = "*.tar.gz"), file.copy, dirsrc)
unlink("*.tar.gz")
sapply(list.files(".",pattern = "*.tgz"), file.copy, dirmac)
unlink("*.tgz")
sapply(list.files(".",pattern = "*.zip"), file.copy, dirwin)
unlink("*.zip")

tools::write_PACKAGES(dirmac,type = 'mac.binary')
tools::write_PACKAGES(dirwin,type = 'win.binary')
tools::write_PACKAGES(dirsrc,type = 'source')

# commit to repo
setwd("~/gh/minicran")
system("git add --all .")
mess <- paste0("radiant package update: ", format(Sys.Date(), format="%m-%d-%Y"))
system(paste0("git commit -m '",mess,"'"))
system("git push")
## move packages to radiant_miniCRAN. must package in Windows first
setwd("~/gh/")
sapply(list.files(".",pattern = "*.tar.gz"), file.copy, dirsrc)
unlink("*.tar.gz")
sapply(list.files(".",pattern = "*.tgz"), file.copy, dirmac)
unlink("*.tgz")
sapply(list.files(".",pattern = "*.zip"), file.copy, dirwin)
unlink("*.zip")

tools::write_PACKAGES(dirmac,type = 'mac.binary')
tools::write_PACKAGES(dirwin,type = 'win.binary')
tools::write_PACKAGES(dirsrc,type = 'source')

# commit to repo
setwd("~/gh/minicran")
system("git add --all .")
mess <- paste0("radiant package update: ", format(Sys.Date(), format="%m-%d-%Y"))
system(paste0("git commit -m '",mess,"'"))
system("git push")

}
Loading

0 comments on commit eb88ed3

Please sign in to comment.