Upload 3 files

Dockerfile

@@ -0,0 +1,23 @@
+FROM rocker/r-ver:4.3.2
+
+# Install system dependencies
+RUN apt-get update && apt-get install -y \
+libcurl4-openssl-dev \
+libssl-dev \
+libxml2-dev \
+&& rm -rf /var/lib/apt/lists/*
+  
+  # Install R packages
+  RUN R -e "install.packages('plumber', repos='https://cloud.r-project.org/')"
+
+# Set working directory
+WORKDIR /app
+
+# Copy application files
+COPY app.R /app/
+  
+  # Expose port (Hugging Face uses 7860)
+  EXPOSE 7860
+
+# Start the API
+CMD ["R", "-e", "pr <- plumber::plumb('app.R'); pr$run(host='0.0.0.0', port=7860)"]

README.md

@@ -1,12 +1,39 @@
----
-title: D.quantile
-emoji: 🔥
-colorFrom: gray
-colorTo: red
-sdk: docker
-pinned: false
-license: mit
-short_description: Distribution-free effect size
----
-
-Check out the configuration reference at https://huggingface.co/docs/hub/spaces-config-reference
+---
+title: Effect Size Calculator
+emoji: 📊
+colorFrom: blue
+colorTo: green
+sdk: docker
+pinned: false
+license: mit
+---
+
+# Distribution-Free Effect Size Calculator
+
+This API calculates distribution-free effect sizes using quantile function modeling.
+
+## Usage
+
+Send GET or POST request to `/calculate` with:
+- `group1`: comma-separated numeric values
+- `group2`: comma-separated numeric values
+- `degree`: polynomial degree (optional, default=5)
+
+Example:
+```
+/calculate?group1=1,2,3,4,5&group2=6,7,8,9,10
+```
+
+## API Endpoints
+
+- `GET /` - Health check
+- `GET /calculate` - Calculate effect size
+- `POST /calculate` - Calculate effect size
+
+## Author
+
+Wolfgang Lenhard and Alexandra Lenhard
+
+## Citation
+
+Lenhard, W. (submitted). Distribution-Free Effect Size Estimation via Quantile Function Modeling.

app.R

@@ -0,0 +1,261 @@
+library(plumber)
+
+#* @apiTitle Effect Size Calculator API
+
+# Your functions here
+d.quantile <- function(x1, x2, degree = 5, silent = TRUE) {
+  
+  # Input validation
+  if (!is.numeric(x1) || !is.numeric(x2)) {
+    stop("Both x1 and x2 must be numeric vectors.")
+  }
+  
+  n1 <- length(x1)
+  n2 <- length(x2)
+  
+  if (n1 < degree + 1) {
+    stop("Group 1 has insufficient data: need at least ", degree + 1, 
+         " observations for degree ", degree, " polynomial (got ", n1, ").")
+  }
+  
+  if (n2 < degree + 1) {
+    stop("Group 2 has insufficient data: need at least ", degree + 1, 
+         " observations for degree ", degree, " polynomial (got ", n2, ").")
+  }
+  
+  if (any(is.na(x1)) || any(is.na(x2))) {
+    warning("Missing values detected and will be removed.")
+    x1 <- x1[!is.na(x1)]
+    x2 <- x2[!is.na(x2)]
+    n1 <- length(x1)
+    n2 <- length(x2)
+  }
+  
+  if (n1 == 0 || n2 == 0) {
+    stop("Cannot compute effect size with empty groups after removing NAs.")
+  }
+  
+  model1 <- fit_quantile_function(x1, degree)
+  model2 <- fit_quantile_function(x2, degree)
+  
+  tie1 <- attr(model1, "tie_proportion")
+  tie2 <- attr(model2, "tie_proportion")
+  
+  if (!silent && (tie1 > 0.3 || tie2 > 0.3)) {
+    message(sprintf(
+      "Note: Substantial ties detected (Group 1: %.1f%%, Group 2: %.1f%%).",
+      tie1 * 100, tie2 * 100
+    ))
+  }
+  
+  moments1 <- get_moments(model1, group_label = "Group 1")
+  moments2 <- get_moments(model2, group_label = "Group 2")
+  
+  weighted_pooled_variance <- ((n1 -1) * moments1$variance + (n2 - 1) * moments2$variance) / (n1 + n2 - 2)
+  pooled_sd <- sqrt(weighted_pooled_variance)
+  
+  mean_diff <- moments2$mean - moments1$mean
+  
+  if (pooled_sd == 0) {
+    if (mean_diff == 0) {
+      d_q <- 0
+    } else {
+      d_q <- sign(mean_diff) * Inf
+      warning("Pooled SD is zero but means differ. Returning Inf with appropriate sign.")
+    }
+  } else {
+    d_q <- mean_diff / pooled_sd
+  }
+  
+  result <- list(
+    d_q = d_q,
+    group1_mean = moments1$mean,
+    group1_variance = moments1$variance,
+    group1_sd = sqrt(moments1$variance),
+    group2_mean = moments2$mean,
+    group2_variance = moments2$variance,
+    group2_sd = sqrt(moments2$variance),
+    pooled_sd = pooled_sd,
+    n1 = n1,
+    n2 = n2,
+    degree = degree
+  )
+  
+  result$tie_proportion_1 <- tie1
+  result$tie_proportion_2 <- tie2
+  
+  return(result)
+}
+
+fit_quantile_function <- function(x, poly_degree, 
+                                  check_monotonicity = FALSE,
+                                  min_degree = 2) {
+  
+  n <- length(x)
+  
+  if (n < 3) {
+    stop("Need at least 3 observations to fit a polynomial quantile function.")
+  }
+  
+  n_unique <- length(unique(x))
+  tie_proportion <- 1 - (n_unique / n)
+  
+  max_possible_degree <- n_unique - 1
+  
+  if (poly_degree > max_possible_degree) {
+    poly_degree <- max_possible_degree
+  }
+  
+  if (tie_proportion > 0.3 && poly_degree > 3) {
+    recommended_degree <- min(poly_degree, max(3, floor(n_unique / 2)))
+    if (recommended_degree < poly_degree) {
+      poly_degree <- recommended_degree
+    }
+  }
+  
+  if (poly_degree < min_degree) {
+    stop(sprintf(
+      "Insufficient unique values (%d) to fit minimum polynomial degree (%d).",
+      n_unique, min_degree
+    ))
+  }
+  
+  avg_ranks <- rank(x, ties.method = "average")
+  p <- (avg_ranks - 0.5) / n
+  z <- qnorm(p)
+  
+  current_degree <- poly_degree
+  model <- lm(x ~ poly(z, current_degree, raw = TRUE))
+  
+  attr(model, "sample_size") <- n
+  attr(model, "n_unique") <- n_unique
+  attr(model, "tie_proportion") <- tie_proportion
+  attr(model, "poly_degree") <- current_degree
+  
+  return(model)
+}
+
+get_moments <- function(model, group_label = "Unknown") {
+  coeffs <- coef(model)
+  poly_degree <- length(coeffs) - 1
+  
+  f <- function(z) {
+    val <- coeffs[poly_degree + 1]
+    for (i in poly_degree:1) {
+      val <- val * z + coeffs[i]
+    }
+    return(val)
+  }
+  
+  mean_integrand <- function(z) { 
+    f(z) * dnorm(z) 
+  }
+  
+  mean_result <- integrate(
+    mean_integrand, 
+    lower = -Inf, 
+    upper = Inf, 
+    subdivisions = 2000L,
+    rel.tol = 1e-8,
+    abs.tol = 1e-10,
+    stop.on.error = FALSE
+  )
+  
+  mu <- mean_result$value
+  
+  variance_integrand <- function(z) { 
+    deviation <- f(z) - mu
+    deviation^2 * dnorm(z) 
+  }
+  
+  variance_result <- integrate(
+    variance_integrand, 
+    lower = -Inf, 
+    upper = Inf, 
+    subdivisions = 2000L,
+    rel.tol = 1e-8,
+    abs.tol = 1e-10,
+    stop.on.error = FALSE
+  )
+  
+  var <- variance_result$value
+  
+  if (var < 0 && abs(var) < 1e-10) {
+    var <- 0
+  }
+  
+  return(list(
+    mean = mu, 
+    variance = var
+  ))
+}
+
+# API endpoint
+#* Calculate effect size from two groups
+#* @param group1 Comma-separated numeric values for group 1
+#* @param group2 Comma-separated numeric values for group 2
+#* @param degree Polynomial degree (default 5)
+#* @post /calculate
+#* @get /calculate
+function(group1, group2, degree = 5) {
+  
+  tryCatch({
+    # Parse input
+    x1 <- as.numeric(unlist(strsplit(group1, ",")))
+    x2 <- as.numeric(unlist(strsplit(group2, ",")))
+    deg <- as.numeric(degree)
+    
+    # Remove any NA values from parsing
+    x1 <- x1[!is.na(x1)]
+    x2 <- x2[!is.na(x2)]
+    
+    if (length(x1) == 0 || length(x2) == 0) {
+      return(list(
+        error = "Invalid input: Could not parse numeric values from input strings"
+      ))
+    }
+    
+    # Calculate effect size
+    result <- d.quantile(x1, x2, degree = deg, silent = TRUE)
+    
+    # Return clean result
+    return(list(
+      success = TRUE,
+      d_q = result$d_q,
+      group1 = list(
+        n = result$n1,
+        mean = result$group1_mean,
+        sd = result$group1_sd
+      ),
+      group2 = list(
+        n = result$n2,
+        mean = result$group2_mean,
+        sd = result$group2_sd
+      ),
+      pooled_sd = result$pooled_sd,
+      degree = result$degree,
+      tie_info = list(
+        group1_ties = result$tie_proportion_1,
+        group2_ties = result$tie_proportion_2
+      )
+    ))
+    
+  }, error = function(e) {
+    return(list(
+      success = FALSE,
+      error = as.character(e$message)
+    ))
+  })
+}
+
+#* Health check endpoint
+#* @get /
+function() {
+  return(list(
+    status = "running",
+    message = "Effect Size Calculator API is ready",
+    endpoints = list(
+      calculate = "/calculate?group1=1,2,3&group2=4,5,6"
+    )
+  ))
+}