R research package development with devtools, roxygen2 documentation, testthat testing, CRAN submission, and vignette creation for statistical methods.
Instalar com Codex ou Claude Copie este prompt, cole no Codex, Claude ou outro assistente e deixe que ele revise a página da skill e instale para você.
Um comando direto ignora o prompt de revisão. Verifique a origem antes de executá-lo.
R research package development with devtools, roxygen2 documentation, testthat testing, CRAN submission, and vignette creation for statistical methods.
S4: Formal. setClass(), setMethod(), slot access via @
R6: Reference semantics (mutable). Ideal for stateful objects.
DESCRIPTION File
Key fields:
Package, Title, Description, Version
: ORCID-linked author metadata
Authors@R
Depends, Imports, Suggests
License, Encoding, LazyData
Environment Setup
# Install development tools (run in R)# install.packages(c("devtools", "roxygen2", "testthat", "usethis", "pkgdown"))# Python interface (optional)
pip install rpy2>=3.5 pandas>=2.0 numpy>=1.24
# Verify R tools
Rscript -e "packageVersion('devtools')"
Rscript -e "packageVersion('testthat')"
Core Workflow
Step 1: Create Package Scaffolding
# Run in R console or Rscript
library(usethis)
library(devtools)# Create new package
usethis::create_package("~/research/mystatpack")
setwd("~/research/mystatpack")# Set up testing
usethis::use_testthat()# Set up version control
usethis::use_git()# Set up CI (GitHub Actions)
usethis::use_github_action("check-standard")# Add a license
usethis::use_mit_license()# Create vignette template
usethis::use_vignette("introduction")
# File: DESCRIPTION
Package: mystatpack
Title: Statistical Methods for Research Analysis
Version: 0.1.0
Authors@R:
person("Research", "Team",
email = "research@university.edu",
role = c("aut", "cre"),
comment = c(ORCID = "0000-0001-2345-6789"))
Description: Provides functions for common statistical analyses in
academic research, including effect size computation, power
analysis, and reproducible reporting utilities.
License: MIT + file LICENSE
Encoding: UTF-8
LazyData: true
RoxygenNote: 7.2.3
Imports:
stats,
utils
Suggests:
testthat (>= 3.0.0),
knitr,
rmarkdown
Config/testthat/edition: 3
Step 2: Write Documented R Functions
# File: R/effect_size.R#' Compute Effect Size Between Two Groups#'#' @description#' Calculates Cohen's d, Glass's delta, or Hedges' g effect size for#' comparing means between two independent groups. Includes bootstrap#' confidence intervals when \code{ci = TRUE}.#'#' @param group1 Numeric vector. Observations from the first group.#' @param group2 Numeric vector. Observations from the second group.#' @param method Character string. Effect size measure:#' \code{"cohens_d"} (default), \code{"glass_delta"}, or \code{"hedges_g"}.#' @param ci Logical. Compute bootstrap confidence interval? Default \code{FALSE}.#' @param n_boot Integer. Number of bootstrap resamples if \code{ci = TRUE}. Default 1000.#' @param alpha Numeric. Significance level for CI. Default 0.05.#'#' @return An object of class \code{"effect_size"} with components:#' \describe{#' \item{estimate}{Point estimate of the effect size.}#' \item{method}{Effect size method used.}#' \item{ci_lower}{Lower CI bound (if \code{ci = TRUE}, else \code{NA}).}#' \item{ci_upper}{Upper CI bound (if \code{ci = TRUE}, else \code{NA}).}#' \item{n1}{Sample size of group 1.}#' \item{n2}{Sample size of group 2.}#' \item{interpretation}{Qualitative interpretation following Cohen (1988).}#' }#'#' @references#' Cohen, J. (1988). \emph{Statistical Power Analysis for the Behavioral Sciences}#' (2nd ed.). Lawrence Erlbaum Associates.#'#' @seealso\code{\link{power_analysis}}, \code{\link{cohens_d_ci}}#'#' @examples
#' set.seed(42)
#' g1 <- rnorm(50, mean = 10, sd = 2)
#' g2 <- rnorm(50, mean = 8, sd = 2)
#'
#' # Cohen's d
#' result <- compute_effect_size(g1, g2)
#' print(result)
#'
#' # Hedges' g with bootstrap CI
#' result_g <- compute_effect_size(g1, g2, method = "hedges_g", ci = TRUE)
#' summary(result_g)
#'#' @export
compute_effect_size <-function(
group1,
group2,
method =c("cohens_d","glass_delta","hedges_g"),
ci =FALSE,
n_boot =1000L,
alpha =0.05){
method <- match.arg(method)# Input validationif(!is.numeric(group1)||!is.numeric(group2)){
stop("'group1' and 'group2' must be numeric vectors.", call. =FALSE)}if(length(group1)<2||length(group2)<2){
stop("Both groups must have at least 2 observations.", call. =FALSE)}if(alpha <=0|| alpha >=1){
stop("'alpha' must be between 0 and 1 exclusive.", call. =FALSE)}# Compute effect size
estimate <- .es_compute(group1, group2, method)# Bootstrap CI
ci_lower <-NA_real_
ci_upper <-NA_real_if(ci){
boot_estimates <- vapply(seq_len(n_boot),function(i){
b1 <- sample(group1, replace =TRUE)
b2 <- sample(group2, replace =TRUE)
.es_compute(b1, b2, method)},
numeric(1L))
ci_lower <- stats::quantile(boot_estimates, alpha /2,names=FALSE)
ci_upper <- stats::quantile(boot_estimates,1- alpha /2,names=FALSE)}# Interpret
abs_d <-abs(estimate)
interpretation <- dplyr_like_case(
abs_d <0.2~"negligible",
abs_d <0.5~"small",
abs_d <0.8~"medium",TRUE~"large")
structure(list(
estimate = estimate,
method = method,
ci_lower = ci_lower,
ci_upper = ci_upper,
n1 =length(group1),
n2 =length(group2),
interpretation = interpretation
),class="effect_size")}# Internal helper (not exported)
.es_compute <-function(g1, g2, method){
n1 <-length(g1); n2 <-length(g2)
mean_diff <- mean(g1)- mean(g2)switch(method,
cohens_d ={
pooled_sd <-sqrt(((n1 -1)* var(g1)+(n2 -1)* var(g2))/(n1 + n2 -2))
mean_diff / pooled_sd
},
glass_delta = mean_diff / sd(g2),
hedges_g ={
d <- .es_compute(g1, g2,"cohens_d")
j <- 1 -3/(4*(n1 + n2 -2)-1)
d * j
})}# Workaround for dplyr::case_when without dependency
dplyr_like_case <-function(...){
cases <-list(...)for(formula in cases){
cond <- eval(formula[[2]])if(isTRUE(cond))return(eval(formula[[3]]))}NA_character_}#' @export
print.effect_size <-function(x, ...){
cat(sprintf("Effect Size (%s)\n Estimate: %.3f [%s]\n n1 = %d, n2 = %d\n",
x$method, x$estimate, x$interpretation, x$n1, x$n2
))if(!is.na(x$ci_lower)){
cat(sprintf(" 95%% CI: [%.3f, %.3f]\n", x$ci_lower, x$ci_upper))}invisible(x)}#' @export
summary.effect_size <-function(object, ...){
print(object, ...)
cat(sprintf(" Interpretation: %s effect (%s, 1988)\n",
object$interpretation,"Cohen"))invisible(object)}
# File: R/power_analysis.R#' Sample Size Determination for Two-Sample t-Test#'#' @description#' Computes the required sample size per group to achieve specified#' statistical power for detecting a given effect size.#'#' @param d Numeric. Expected Cohen's d effect size.#' @param power Numeric. Desired statistical power (default 0.80).#' @param alpha Numeric. Significance level (default 0.05, two-tailed).#' @param ratio Numeric. Ratio n2/n1 for unequal group sizes (default 1).#'#' @return Named numeric vector with \code{n1}, \code{n2}, and \code{total}.#'#' @examples
#' # Power for medium effect, 80% power, alpha = 0.05
#' power_analysis(d = 0.5)
#'
#' # Small effect, 90% power
#' power_analysis(d = 0.2, power = 0.90)
#'#' @export
power_analysis <-function(d, power =0.80, alpha =0.05, ratio =1){if(d <=0) stop("'d' must be positive.")if(power <=0|| power >=1) stop("'power' must be in (0, 1).")
z_alpha <- stats::qnorm(1- alpha /2)
z_beta <- stats::qnorm(power)# Derivation for unequal samples:# n1 = (z_alpha + z_beta)^2 * (1 + 1/ratio) / d^2
n1 <-ceiling((z_alpha + z_beta)^2*(1+1/ ratio)/ d^2)
n2 <-ceiling(n1 * ratio)c(n1 = n1, n2 = n2, total = n1 + n2)}
Step 3: testthat Unit Tests
# File: tests/testthat/test-effect_size.R
library(testthat)
test_that("compute_effect_size returns correct class",{
set.seed(42)
g1 <- rnorm(50, mean =10, sd =2)
g2 <- rnorm(50, mean =8, sd =2)
result <- compute_effect_size(g1, g2)
expect_s3_class(result,"effect_size")})
test_that("Cohen's d is approximately 1.0 for known groups",{
set.seed(42)
n <- 1000
g1 <- rnorm(n, mean =2, sd =1)# true d = (2-0)/1 = 2
g2 <- rnorm(n, mean =0, sd =1)
result <- compute_effect_size(g1, g2, method ="cohens_d")
expect_equal(result$estimate,2, tolerance =0.1)})
test_that("Reversing groups negates effect size",{
set.seed(42)
g1 <- rnorm(50, mean =10, sd =2)
g2 <- rnorm(50, mean =8, sd =2)
d_fwd <- compute_effect_size(g1, g2)$estimate
d_rev <- compute_effect_size(g2, g1)$estimate
expect_equal(d_fwd,-d_rev, tolerance =1e-10)})
test_that("Hedges' g is smaller than Cohen's d",{
set.seed(42)
g1 <- rnorm(20, mean =10, sd =2)
g2 <- rnorm(20, mean =8, sd =2)
d <-abs(compute_effect_size(g1, g2, method ="cohens_d")$estimate)
g <-abs(compute_effect_size(g1, g2, method ="hedges_g")$estimate)
expect_lt(g, d)})
test_that("Bootstrap CI has correct coverage",{
set.seed(42)
g1 <- rnorm(100, mean =10, sd =2)
g2 <- rnorm(100, mean =8, sd =2)
result <- compute_effect_size(g1, g2, ci =TRUE, n_boot =200)
expect_true(!is.na(result$ci_lower))
expect_true(result$ci_lower < result$estimate)
expect_true(result$ci_upper > result$estimate)})
test_that("Invalid inputs raise errors",{
expect_error(compute_effect_size("a",c(1,2,3)),"numeric")
expect_error(compute_effect_size(c(1),c(1,2,3)),"at least 2")
expect_error(compute_effect_size(c(1,2),c(1,2), alpha =1.5),"alpha")})
test_that("power_analysis gives sensible n for d=0.5",{
result <- power_analysis(d =0.5, power =0.80, alpha =0.05)
expect_true(result["n1"]>=60&& result["n1"]<=70)})
test_that("All effect size methods return numeric",{
g1 <-c(8,9,10,11,12)
g2 <-c(6,7,8,9,10)for(m inc("cohens_d","glass_delta","hedges_g")){
result <- compute_effect_size(g1, g2, method = m)
expect_type(result$estimate,"double")}})
# Run all tests
Rscript -e "devtools::test()"# Run R CMD check (full CRAN check)
Rscript -e "devtools::check()"# Install package
Rscript -e "devtools::install()"# Build documentation
Rscript -e "devtools::document()"# Build pkgdown website
Rscript -e "pkgdown::build_site()"
Advanced Usage
rpy2: Call R Package from Python
import os
import numpy as np
import pandas as pd
# Set up rpy2 environment
os.environ["R_HOME"] = os.popen("R RHOME").read().strip()
try:
import rpy2.robjects as ro
from rpy2.robjects import numpy2ri, pandas2ri
from rpy2.robjects.packages import importr
numpy2ri.activate()
pandas2ri.activate()
# Install package if needed (run once)# utils = importr("utils")# utils.install_packages("mystatpack")# Import the R package# mystatpack = importr("mystatpack")# Call R function from Python# g1_py = np.random.normal(10, 2, 50)# g2_py = np.random.normal(8, 2, 50)# result_r = mystatpack.compute_effect_size(g1_py, g2_py)# print(f"Cohen's d (via rpy2): {result_r.rx2('estimate')[0]:.3f}")# Use base R stats functions
stats = importr("stats")
g1 = np.random.normal(10, 2, 50)
g2 = np.random.normal(8, 2, 50)
t_test = stats.t_test(g1, g2, alternative="two.sided")
print(f"t-test via rpy2: t = {t_test.rx2('statistic')[0]:.3f}, "f"p = {t_test.rx2('p.value')[0]:.4f}")
except ImportError:
print("rpy2 not available — using Python statsmodels instead")
import scipy.stats as stats
np.random.seed(42)
g1 = np.random.normal(10, 2, 50)
g2 = np.random.normal(8, 2, 50)
t_stat, p_val = stats.ttest_ind(g1, g2)
print(f"t-test (scipy): t = {t_stat:.3f}, p = {p_val:.4f}")
Package Vignette Template
# File: vignettes/introduction.Rmd---
title:"Introduction to mystatpack"
output: rmarkdown::html_vignette
vignette:>%\VignetteIndexEntry{Introduction to mystatpack}
%\VignetteEngine{knitr::rmarkdown}%\VignetteEncoding{UTF-8}
---
\`\`\`{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 7, fig.height = 4)
library(mystatpack)
set.seed(42)
\`\`\`
## Basic Usage
### Computing Effect Sizes
\`\`\`{r effect-size-example}
# Generate two groups
group1 <- rnorm(50, mean = 10, sd = 2)
group2 <- rnorm(50, mean = 8, sd = 2)
# Cohen's d
result <- compute_effect_size(group1, group2)
print(result)
\`\`\`
### Power Analysis
\`\`\`{r power-example}
# Required sample size for medium effect (d=0.5), 80% power
n_required <- power_analysis(d =0.5)
cat("Required n per group:", n_required["n1"],"\n")
cat("Total sample size: ", n_required["total"],"\n")\`\`\`
## Advanced Usage
### Bootstrap Confidence Intervals
\`\`\`{r bootstrap-ci}
result_ci <- compute_effect_size(
group1, group2,
method = "hedges_g",
ci = TRUE,
n_boot = 1000
)
summary(result_ci)
\`\`\`
