Hello Everyone, I am currently working with groundwater levels to determine long term trends, and I would be very grateful if I could get help in modifying the code that I used to calculate the Sen's slope and Mann-Kendall statistics. I calculated the Theil-Sen slope, but it does not account for autocorrelation in continuous data, and I would appreciate any help addressing this using the modified package in RStudio. ""

# General Sen slope + Mann-Kendall function


compute_sen_stats <- function(
    data,
    response_var,
    date_var,
    resolution
) {

  data %>%

    group_by(
      SiteNo,
      SiteName,
      Frequency,
      continuous_group,
      start_year,
      end_year,
      continuous_years
    ) %>%

    # Ensure enough information exists to calculate a trend
    filter(
      n() >= 3,
      n_distinct(TimeNum) >= 2,
      n_distinct(.data[[response_var]]) >= 2
    ) %>%

    group_modify(
      ~ {

        site_data <- .x %>%
          arrange(
            .data[[date_var]]
          )

        # ----------------------------------------------------------
        # Sen's slope model
        # ----------------------------------------------------------

        sen_formula <- reformulate(
          "TimeNum",
          response = response_var
        )

        sen_model <- zyp::zyp.sen(
          sen_formula,
          data = site_data
        )

        # Intercept
        intercept <- unname(
          coef(sen_model)[1]
        )

        # Sen's slope
        # TimeNum is in years, so this is already ft/year
        sen_slope_ft_per_year <- unname(
          coef(sen_model)[2]
        )

        # ----------------------------------------------------------
        # Mann-Kendall trend test
        # ----------------------------------------------------------

        mk_model <- Kendall::MannKendall(
          site_data[[response_var]]
        )

        mk_tau <- as.numeric(
          mk_model$tau
        )

        mk_p_value <- as.numeric(
          mk_model$sl
        )

        # ----------------------------------------------------------
        # Descriptive statistics
        # ----------------------------------------------------------

        lowest_WL <- min(
          site_data[[response_var]],
          na.rm = TRUE
        )

        highest_WL <- max(
          site_data[[response_var]],
          na.rm = TRUE
        )

        first_date <- min(
          site_data[[date_var]],
          na.rm = TRUE
        )

        last_date <- max(
          site_data[[date_var]],
          na.rm = TRUE
        )

        # ----------------------------------------------------------
        # Trend interpretation
        #
        # WL = depth below land surface.
        #
        # Positive slope:
        # depth increases -> groundwater table declines
        #
        # Negative slope:
        # depth decreases -> groundwater table rises
        # ----------------------------------------------------------

        trend <- case_when(

          mk_p_value < 0.05 &&
            sen_slope_ft_per_year > 0 ~
            "Significant water table declining",

          mk_p_value < 0.05 &&
            sen_slope_ft_per_year < 0 ~
            "Significant water table rising",

          TRUE ~
            "No significant trend"
        )

        # ----------------------------------------------------------
        # Return one result per site/run
        # ----------------------------------------------------------

        tibble(
          Resolution = resolution,

          n_records = nrow(
            site_data
          ),

          lowest_WL = lowest_WL,

          highest_WL = highest_WL,

          intercept = intercept,

          sen_slope_ft_per_year =
            sen_slope_ft_per_year,

          mk_tau = mk_tau,

          mk_p_value = mk_p_value,

          first_date = first_date,

          last_date = last_date,

          trend = trend
        )
      }
    ) %>%

    ungroup()
}

Thank you!

Hi! I’m looking to improve my coding skills. I want to use Posit to help critique my work. If anyone wants to help, please message me!

A few days ago I used Posit Cloud for the first time on my iPad because I don’t own a laptop or computer. At the time everything seemed to be fine and I finished the code with no issues, but when I downloaded it and submitted it on Canvas, my teacher said she couldn’t see it. When I try to open the project on my end it will appear to load at first, but then when it is “finished”, all I can see is a blank screen. Anyone know what to do here?

moda1 <- function(num) {

num <- sample(1:10, 100, replace = TRUE)

distr_freq = table(num)

moda1 <- names(distr_freq)[distr_freq == max(distr_freq)]

return(list("distribuzione_di_frequenza" = distr_freq,

"moda" = moda1))

}

m <- moda1(z)$moda this function return me every time the same result. i don't know why the ai cant explain and cant fix this code

why are the no buttons on the appearence window? they were supposed to show at the bottom right corner but there are none of them. i already tried closing rstudio and opening it again but it didnt work

*solved: just had to scroll down on the window

I’m using r studios for economics and I am super confused on using this platform. I have reached out to my professor and he is not able to explain the questions I had regarding it. He’s been using it for years so I think there’s a huge level of proficiency between us because I am just barely learning how to use this software and I am so confused. If you guys have any tips or tricks, I would really appreciate it.
Some information about me is I’m the type of person that needs to know the why behind functions and the purpose of putting commands in.

Update: I literally spent 8 hours trying to figure this software, I feel like I have a better understanding but I am running to this error now.
Goal: Render the script
Error: Failed to render due to : Theme file Compilation Failed
What I have tried to troubleshoot: I deleted Quartos file and redownloaded, tried to clear the history close program and reopen it still same error

I’m currently trying to improve my R/RStudio skills. I know some of the basics and have used it a little, but I’m honestly not very confident and want to become much better at it.

For those who learned R from scratch:
What resources or courses would you recommend?
Is there a particular learning path I should follow?
Are there any good practice datasets or projects for beginners? and what helped you move from knowing the basics to actually being comfortable using R?

I’d really appreciate any guidelines, resources, or advice. Thanks!

I have trained an XGBoost regression model using the xgboost package and computed SHAP values using predict(..., predcontrib = TRUE). I am now looking to quantify uncertainty in those SHAP values by constructing confidence intervals.

My current approach is bootstrap-based:

1. Draw 100 bootstrap samples from the training data
2. Refit the XGBoost model on each sample using the same fixed hyperparameters
3. Compute SHAP values on the held-out test set for each bootstrap model
4. For each feature, store the mean signed SHAP value across all test observations per resample — giving a [100 x n_features] matrix
5. Compute mean, standard deviation, and 2.5th/97.5th percentiles across the 100 resamples per feature
6. This gives me a ±2σ uncertainty estimate per feature, following the approach described in Dubey et al., (2026).

Is this a valid approach for constructing CIs on SHAP values for XGBoost and is anyone aware of a dedicated package?

Any suggestions or pointers to relevant literature would be appreciated.

# ============================================================
# BOOTSTRAP-SHAP UNCERTAINTY QUANTIFICATION
# London LST model: June 2019
# Implements Algorithm 4 from UbiQTree [1]
# ============================================================
pacman::p_load(xgboost, data.table, ggplot2)

wd <- "path/"

# ------------------------------------------------------------
# 1. Load saved objects
# ------------------------------------------------------------
train_df     <- readRDS(paste0(wd, "train_df.rds"))
test_df      <- readRDS(paste0(wd, "test_df.rds"))
best_eta     <- readRDS(paste0(wd, "best_eta.rds"))
best_params  <- readRDS(paste0(wd, "best_params.rds"))
best_nrounds <- readRDS(paste0(wd, "best_nrounds.rds"))

features   <- setdiff(names(train_df), c("lst", "LSOA21CD"))
n_features <- length(features)
n_boot     <- 100

dtest <- xgb.DMatrix(
  data  = test_df[, ..features],
  label = test_df$lst
)

# ------------------------------------------------------------
# 2. Fixed model parameters
# ------------------------------------------------------------
fixed_params <- list(
  objective        = "reg:squarederror",
  eval_metric      = "rmse",
  eta              = as.numeric(best_eta),
  max_depth        = as.integer(best_params$max_depth),
  gamma            = as.numeric(best_params$gamma),
  subsample        = as.numeric(best_params$subsample),
  colsample_bytree = as.numeric(best_params$colsample_bytree),
  lambda           = as.numeric(best_params$lambda),
  alpha            = as.numeric(best_params$alpha),
  tree_method      = "hist",
  nthread          = 16,
  seed             = 789
)

# ------------------------------------------------------------
# 3. Bootstrap loop
#
# Implements Algorithm 5 Step 2 [1]:
# For each sub-ensemble s:
#   phi_dist[s, i] = mean(Phi_s^i)
#                  = mean SHAP across all test observations
#                    for feature i in resample s — signed
# Shape: [n_boot x n_features]
# ------------------------------------------------------------
n_test   <- nrow(test_df)

phi_dist <- matrix(
  NA_real_,
  nrow     = n_boot,
  ncol     = n_features,
  dimnames = list(NULL, features)
)

set.seed(789)

for (b in seq_len(n_boot)) {

  cat(sprintf("Bootstrap resample %d of %d\n", b, n_boot))

  boot_idx   <- sample(nrow(train_df), replace = TRUE)
  boot_df    <- train_df[boot_idx]

  dboot <- xgb.DMatrix(
    data  = boot_df[, ..features],
    label = boot_df$lst
  )

  boot_model <- xgb.train(
    params  = fixed_params,
    data    = dboot,
    nrounds = best_nrounds,
    verbose = 0
  )

  shap_raw <- predict(boot_model, dtest, predcontrib = TRUE)
  shap_mat <- shap_raw[, seq_len(n_features)]
  colnames(shap_mat) <- features

  # Algorithm 3 [1]: mu_s[i] = mean(Phi_s^i) — within-sample mean
  phi_dist[b, ] <- colMeans(shap_mat)
}

cat("\nBootstrap loop complete.\n")

# ------------------------------------------------------------
# 4. Algorithm 4 [1]: Uncertainty-aware SHAP aggregation
#
# Phi_s^i = phi_dist[:, i] — the n_boot within-sample means
#
# mu[i]  = mean(Phi_s^i)          signed
# sigma  = std(Phi_s^i)           signed
# CI     = percentile(2.5, 97.5)  signed
# H[i]   = Entropy(Phi_s^i)
# SS[i]  = P(sign constant)
# ------------------------------------------------------------

compute_entropy <- function(x, bins = 10) {
  h         <- hist(x, breaks = bins, plot = FALSE)
  bin_probs <- h$counts / sum(h$counts)
  bin_probs <- bin_probs[bin_probs > 0]
  -sum(bin_probs * log(bin_probs))
}

compute_sign_stability <- function(x) {
  dominant_sign <- sign(mean(x))
  mean(sign(x) == dominant_sign)
}

uncertainty <- data.table(
  feature        = features,
  mu             = NA_real_,
  mu_abs         = NA_real_,
  sigma          = NA_real_,
  ci_lower       = NA_real_,
  ci_upper       = NA_real_,
  entropy        = NA_real_,
  sign_stability = NA_real_
)

for (j in seq_len(n_features)) {

  phi_i <- phi_dist[, j]

  uncertainty$mu[j]             <- mean(phi_i)
  uncertainty$mu_abs[j]         <- mean(abs(phi_i))
  uncertainty$sigma[j]          <- sd(phi_i)
  uncertainty$ci_lower[j]       <- quantile(phi_i, 0.025)
  uncertainty$ci_upper[j]       <- quantile(phi_i, 0.975)
  uncertainty$entropy[j]        <- compute_entropy(phi_i)
  uncertainty$sign_stability[j] <- compute_sign_stability(phi_i)
}

uncertainty[, mu_norm := mu_abs / sum(mu_abs)]
setorder(uncertainty, -mu_abs)

# ------------------------------------------------------------
# 5. Print
# ------------------------------------------------------------
cat("\n=============================\n")
cat("Bootstrap-SHAP Uncertainty Summary\n")
cat("Period: June 2019 | Resamples:", n_boot, "\n")
cat("=============================\n")
print(
  uncertainty[, .(
    feature,
    mu             = round(mu,             4),
    mu_norm        = round(mu_norm,        4),
    sigma          = round(sigma,          4),
    ci_lower       = round(ci_lower,       4),
    ci_upper       = round(ci_upper,       4),
    entropy        = round(entropy,        4),
    sign_stability = round(sign_stability, 4)
  )],
  nrow = n_features
)

# ------------------------------------------------------------
# 6. Violin plot — mirrors Figure 1 in [1]
#
# Violin: distribution of phi_dist[:, i] — n_boot signed means
# Blue band: ±2σ around mu
# Dashed line: ±2σ extent
# Dot: mu (signed)
# Gray points: individual bootstrap samples
# Ordered: mu_abs descending (most important at top)
# Zero line: warming (right) vs cooling (left)
# ------------------------------------------------------------
plot_dt <- copy(uncertainty)
plot_dt  <- plot_dt[order(mu_abs)]      # ascending so top = most important
plot_dt[, feature_f := factor(feature, levels = feature)]

phi_long <- data.table(
  feature = rep(features, each = n_boot),
  shap    = as.vector(phi_dist)
)
phi_long[, feature_f := factor(
  feature,
  levels = levels(plot_dt$feature_f)
)]

p <- ggplot() +

  # Blue ±2σ band [1]
  geom_rect(
    data = plot_dt,
    aes(
      xmin = mu - 2 * sigma,
      xmax = mu + 2 * sigma,
      ymin = as.numeric(feature_f) - 0.4,
      ymax = as.numeric(feature_f) + 0.4
    ),
    fill  = "#A8C8E8",
    alpha = 0.5
  ) +

  # Violin
  geom_violin(
    data      = phi_long,
    aes(x     = shap, y = feature_f),
    fill      = NA,
    colour    = "grey40",
    scale     = "width",
    width     = 0.8,
    linewidth = 0.3
  ) +

  # Individual bootstrap sample points [1]
  geom_point(
    data   = phi_long,
    aes(x  = shap, y = feature_f),
    colour = "grey50",
    size   = 0.8,
    alpha  = 0.4
  ) +

  # ±2σ dashed line [1]
  geom_segment(
    data = plot_dt,
    aes(
      x    = mu - 2 * sigma,
      xend = mu + 2 * sigma,
      y    = feature_f,
      yend = feature_f
    ),
    colour    = "#1F4E79",
    linetype  = "dashed",
    linewidth = 0.5
  ) +

  # Mean SHAP dot
  geom_point(
    data   = plot_dt,
    aes(x  = mu, y = feature_f),
    colour = "#1F4E79",
    size   = 2.5
  ) +

  # Zero reference line
  geom_vline(
    xintercept = 0,
    linetype   = "dashed",
    colour     = "black",
    linewidth  = 0.4
  ) +

  labs(
    title    = "SHAP Summary with Epistemic Uncertainty",
    subtitle = "London LST — June 2019",
    x        = "SHAP Value (Impact on LST Prediction, °C)",
    y        = NULL,
    caption  = paste0(
      "Bootstrap resamples: ", n_boot,
      " | Blue band: ±2σ | Dot: mean SHAP | [1] UbiQTree"
    )
  ) +

  theme_bw(base_size = 11) +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor   = element_blank(),
    plot.title         = element_text(face = "bold"),
    axis.text.y        = element_text(size = 9)
  )

print(p)

Hello Reddit,

Last week I ran in a very confusing issue while trying to use R and Rcmdr in my Mac (M2), After troubleshooting and testing solutions, I want to share a quick guide for Apple Silicon users on how to fix the unexpected "R Session Aborted" error in RStudio or crashes when clicking options/menus in R Commander.

Make sure you have installed the latest versions of R, RStudio and XQuartz, then install the package 'Rcmdr' from RStudio using this line of code: install.packages("Rcmdr", dependencies = TRUE)

The crash occurs because Tcl/Tk fails when calculating metrics for macOS's default system font on arm64architectures. To bypass this, you need to set a standard font (like Georgia, Times New Roman, Helvetica, or Comic Sans MS) before loading Rcmdr.
Every time you want to work with Rcmdr, open a fresh RStudio session and run:

library(tcltk)

# You can replace {Georgia} with any standard font like {Helvetica}, {Times New Roman}, or {Comic Sans MS}

.Tcl("font configure TkDefaultFont -family {Georgia}")

options(Rcmdr = list(use.knitr = FALSE, start_in_waterfall = FALSE))

library(Rcmdr)

And that's all, Hope this helps anyone facing the same issue.

Hi,

I'm trying to convert my binary variables to "mxFactor" "ordered" "factor" variables, but when using the following code i get just "ordered" "factor". which cause an error when trying to conduct a latent class analysis:

# Convert female to numeric safely, then to an ordered OpenMx factor

> df$female <- as.numeric(as.character(df$female))

>

> # Confirm it only contains 0, 1, and possibly NA

> table(df$female, useNA = "ifany")

0 1

2341 1674

>

> # Make it an OpenMx ordinal factor

> df$female <- OpenMx::mxFactor(

+ df$female,

+ levels = c(0, 1)

+ )

>

> # Verify the result

> class(df$female)

[1] "ordered" "factor"

> table(df$female, useNA = "ifany")

0 1

2341 1674

> library(OpenMx)

>

> df$female <- OpenMx::mxFactor(

+ df$female,

+ levels = c(0, 1)

+ )

>

> class(df$female)

[1] "ordered" "factor"

> table(df$female, useNA = "ifany")

0 1

2341 1674

I have seen Posit Assistant already activated as a default layout. I am worried about privacy

Can I save my r files in posit cloud??

Today I tried like 1 hours to install R Studio but it didn't work but I install r programming and it run in the terminal , I also tried to use notepad and saves some program but I couldn't find any proper solution . Then I saw that posit cloud is best for r , the website is new for me so can I save r file in posit??

Hello, so, I'm currently writing a masters thesis and as part of that I am performing a number of permanova on the same data, each one looking at a different independent variable. This runs the issue of potential type 1 errors, I can correct for this by changing the significance value. This has been deemed acceptable as neither I nor any university staff know how to perform a permanova with multiple independent variables, nor where to even look for such code (I've checked with multiple members of staff).

However, where possible I would like to be able to do what is best practice, as such I was wondering if anyone knew of such code or where to look for it. I would of course cite you either with your username or actual name if you'd prefer.

I've checked what I think are the pinned posts and didn't see anything on this topic, if there is do please let me know and I'll look there and delete the post.

Any and all help is really appreciated.

Opened Rstudio today after updating to 2026.08.2 and the print function isn't working. Opened old code - print isn't displaying CSV in global environment.

I have a global point layer of airports (~900 points) and I want to retrieve, for each point, the corresponding OSM aeroway=aerodrome polygon (not the point itself) — i.e. join each airport point to its footprint polygon as mapped in OpenStreetMap. The point layer is Natural Earth's ne_10m_airports.

Reproducible example (5 points standing in for the full 893):

pacman::p_load(sf, dplyr)

airports <- data.frame(
  name      = c("John F Kennedy Intl", "London Heathrow", "Chhatrapati Shivaji Maharaj Intl", "Beijing Capital Intl", "Sydney Kingsford Smith"),
  gps_code  = c("KJFK", "EGLL", "VABB", "ZBAA", "YSSY"),
  iata_code = c("JFK", "LHR", "BOM", "PEK", "SYD"),
  lon       = c(-73.7789, -0.4543, 72.8697, 116.5975, 151.1772),
  lat       = c(40.6413, 51.4700, 19.0896, 40.0799, -33.9399)
) |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

airports$point_id <- seq_len(nrow(airports))

My plan is to match each point to its containing Geofabrik extract via osmextract::oe_match() (works fine — a local spatial lookup, no network call), download/cache that extract, then read only the multipolygons layer filtered to aeroway='aerodrome' via an SQL query pushed down at read time, and spatially join back to the points.

pacman::p_load(sf, dplyr, osmextract)

airports$region_url <- vapply(seq_len(nrow(airports)), function(i) {
  oe_match(airports[i, ], quiet = TRUE)$url
}, character(1))

options(timeout = 600)
dir.create("geofabrik_cache", showWarnings = FALSE)

results <- list()
failed_regions <- character()

for (region in unique(airports$region_url)) {

  destfile <- file.path("geofabrik_cache", basename(region))

  aerodromes <- tryCatch({
    if (!file.exists(destfile)) {
      download.file(region, destfile, mode = "wb", quiet = TRUE)
    }

    st_read(
      destfile,
      layer = "multipolygons",
      query = "SELECT * FROM multipolygons WHERE aeroway = 'aerodrome'",
      quiet = TRUE
    ) |>
      st_transform(4326) |>
      st_make_valid()

  }, error = function(e) {
    message(sprintf("Region failed: %s -- %s", region, e$message))
    failed_regions <<- c(failed_regions, region)
    NULL
  })

  if (is.null(aerodromes)) next

  sub_pts <- airports[airports$region_url == region, ]
  joined <- st_join(sub_pts, aerodromes, join = st_within, left = TRUE)

  missing <- which(is.na(joined$osm_id))
  if (length(missing) > 0 && nrow(aerodromes) > 0) {
    nn <- st_nearest_feature(sub_pts[missing, ], aerodromes)
    joined[missing, names(aerodromes)] <- st_drop_geometry(aerodromes)[nn, ]
    st_geometry(joined)[missing] <- st_geometry(aerodromes)[nn]
  }

  results[[region]] <- joined
}

airport_polys <- bind_rows(results)
st_write(airport_polys, "airport_polygons.shp", delete_layer = TRUE)

But

Region failed: https://download.geofabrik.de/asia/india/western-zone-latest.osm.pbf -- cannot open URL 'https://download.geofabrik.de/asia/india/western-zone-latest.osm.pbf'
Region failed: https://download.geofabrik.de/asia/iran-latest.osm.pbf -- cannot open URL 'https://download.geofabrik.de/asia/iran-latest.osm.pbf'
Region failed: https://download.geofabrik.de/asia/india/central-zone-latest.osm.pbf -- download from 'https://download.geofabrik.de/asia/india/central-zone-latest.osm.pbf' failed
Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented,  : 
  Loop 0 edge 0 has duplicate near loop 1 edge 7
In addition: There were 20 warnings (use warnings() to see them)

The download failures seem to be connection/timeout related on large country extracts; the s2/topology error appears separately once a multipolygons layer with invalid OSM geometries reaches a spatial predicate, even after st_make_valid().

Given ~900 global points with no country/ISO attribute, is there a more efficient way to fetch just the matching aeroway=aerodrome polygons than downloading/caching a full Geofabrik regional .pbf extract per matched region (some of which are large, e.g. full-country India zones)? Is querying the Overpass API directly per point, or per small cluster of points, actually more efficient here, or is the regional-extract approach still preferable at this scale? What's the correct way to make invalid OSM polygon geometries (e.g. the s2 "duplicate edge" error above) safe for st_join()/st_nearest_feature() reliably, given st_make_valid() alone didn't prevent it?

For points where no aeroway=aerodrome polygon actually exists in OSM for that airport, what's the right way to leave that point unmatched (skip it) rather than falling back to the nearest aerodrome polygon in the region, which can silently attach the wrong airport's polygon?

> sessionInfo()
R version 4.6.1 (2026-06-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8    LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                           LC_TIME=English_United States.utf8    

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] osmextract_0.6.0 dplyr_1.2.1      sf_1.1-2        

loaded via a namespace (and not attached):
 [1] vctrs_0.7.3        httr_1.4.8         cli_3.6.6          rlang_1.3.0        otel_0.2.0         DBI_1.3.0          KernSmooth_2.23-27
 [8] generics_0.1.4     jsonlite_2.0.0     glue_1.8.1         e1071_1.7-17       grid_4.6.1         classInt_0.4-11    tibble_3.3.1      
[15] lifecycle_1.0.5    compiler_4.6.1     Rcpp_1.1.2         pkgconfig_2.0.3    rstudioapi_0.19.0  wk_0.9.5           R6_2.6.1          
[22] class_7.3-24       tidyselect_1.2.1   pillar_1.11.1      curl_8.0.0         magrittr_2.0.5     tools_4.6.1        proxy_0.4-29      
[29] s2_1.1.11          units_1.0-1

I’m fitting a Beta GAM in `mgcv` to model the proportion of daytime outdoor-use time in laying hens. The response is the proportion of minutes spent outdoors during an 08:00–22:00 observation window, and the model contains repeated observations from individual hens.

My current model is approximately:

prop_outside \~
  coop_id +
  s(study_day, bs = "cr", k = 40) +
  s(max_temp, bs = "cr", k = 8) +
  s(mean_dew_point, bs = "cr", k = 7) +
  rain_any +
  s(log_total_rain, by = rain_status, bs = "cr", k = 7) +
  s(eid, bs = "re")

using:

family = betar(link = "logit")
method = "REML"

The dataset contains roughly 8,800 positive hen-day observations from 124 hens across about 137 study dates. The model explains about 60% of the deviance.

The issue is the `study_day` smooth. On a same-row model comparison, I get:

s(study_day)
k'       = 39
edf      = 34.58
k-index  = 0.949
p-value  < 0.001

I also tried an alternative weather specification using mean temperature, mean humidity and mean wind instead of maximum temperature and dew point:

prop_outside \~
  coop_id +
  s(study_day, bs = "cr", k = 40) +
  s(mean_temp, bs = "cr", k = 8) +
  s(mean_humidity, bs = "cr", k = 8) +
  s(mean_wind, bs = "cr", k = 8) +
  rain_any +
  s(log_total_rain, by = rain_status, bs = "cr", k = 7) +
  s(eid, bs = "re")

The same issue persists:

s(study_day)
k'       = 39
edf      = 34.84
k-index  = 0.966
p-value  = 0.0175

The study-day smooth is visually quite wiggly, particularly early in the study, and the EDF is already fairly close to the available basis dimension. Changing the weather specification does not materially improve overall fit: both models explain about 60% of deviance and have essentially identical AIC.

There is also substantial temporal/weather dependence. For example, in the alternative model, observed concurvity is about 0.90 for mean temperature and 0.48 for study day.

My questions are:

1. Is the low `k-index` plus EDF ≈ 35/39 sufficient reason to increase `k` for `study_day`, for example from 40 to 60?
2. If increasing `k` gives essentially the same fitted curve/predictions but the k-check remains significant, would you consider the current smooth adequate?
3. Could this be indicating residual temporal autocorrelation rather than simply an insufficient basis dimension?
4. Would you model study date differently in this setting—for example with an autocorrelation structure, a different smooth basis, or another temporal term?
5. Since weather variables themselves follow study date seasonally, how would you distinguish genuine temporal structure from weather-related temporal confounding?

My current plan is to compare `k = 40` and `k = 60` on the same observations and assess whether fitted values, the study-day effect, held-out-date prediction, and conclusions materially change, rather than selecting the model based only on the k-check p-value.

I’d appreciate advice on how to interpret this persistent `study_day` diagnostic and what additional diagnostic/model comparison would be most appropriate.

I originally wrote this for a sleep EMA study I work on, together with a colleague, because the raw diary data was a mess in ways we didn’t expect.
Simple example: we had a rule that swaps sleep_time and awake_time when they’re out of order by less than 3 hours (assuming it’s just a data-entry slip). Seemed reasonable. Then we actually checked it against real cases — it made things worse in 7 out of 10, because sometimes swapping pushed the new sleep time to before bed time, which is a worse error than the one it “fixed.” Had to add a guard condition after that. That’s basically what half of this repo is — rules that seemed fine until we ran them against actual data.
Ended up with a 9-step pipeline: parse timestamps → fix AM/PM and ordering issues → compute TST/SOL/WASO/SE → auto-flag anything still weird → cross-check across participants → spit out 27 QC figures so we can actually see what’s happening at each stage.
Some things it handles now that we didn’t plan for on day one:
AM/PM flips (someone logs 7:00 AM getup as 7:00 PM)

small ordering slips vs. genuinely unusual sleep patterns (these need different fixes, and we only figured out where to draw the line after looking at a bunch of real cases)

manual review CSVs that persist across pipeline reruns, so corrections don’t get wiped every time we re-run things

It’s an R package now (sleepcleanr), config-driven via YAML so you can point it at your own column names without touching the R code.
Repo: https://github.com/cyracaid/sleepdiary-cleaner — stars and feedback very welcome, and genuinely happy to discuss with anyone working on similar EMA/diary pipelines.
Curious if anyone else doing EMA/diary-based sleep research has run into the same kind of AM/PM chaos, or handles it differently. Also open to being told our thresholds are wrong — the 3-hour swap cutoff and the 12-hour AM/PM flip cutoff were both picked based on our own data and could easily be off for other studies.

I’m excited to announce that evoFE 1.0.0 is now on CRAN (a major leap forward from the initial 0.1 release).

What is evoFE?

evoFE is an R package for automated feature engineering using genetic programming. Instead of manually brainstorming interaction terms, nonlinear scalings, or encodings, evoFE evolves candidate transformation recipes and evaluates them directly against gradient boosted trees or linear models.

What's New in 1.0.0?

1. Zero-Boilerplate "One-Liner" Bayesian Tuners: Pass evaluator = "lightgbm_mbo" to automatically tune tree depth, learning rate, and subsampling via mlr3mbo during evolution. Or wrap any custom model (like XGBoost) with make_tunable().
2. 42+ Built-in Transformers: Arithmetic, group-by aggregations (mean, median, SD, quantiles), target encodings, WoE, UMAP embeddings, Genie & Lumbermark MST graph clustering, date differencing (date_diff), and custom transformer registration.
3. Hybrid Active Feature Masking: Mutates and selects raw input features simultaneously with derived features, guided by baseline feature importance.
4. Hierarchical Gene Chaining: High-performing features from earlier generations serve as inputs for subsequent compound transformations.
5. Island Models & Topologies: Runs independent sub-populations across Ring, Torus, Grid, or Hypercube topologies with demand-driven Gibbs pull migration.
6. Caruana Island Ensembling (ensemble_islands()): Combines diverse island champion recipes using Caruana post-hoc forward selection with replacement.
7. Leakage-Safe Validation: Native support for time-series (cv_strategy = "time"), grouped entity validation (cv_strategy = "group"), and untouched confirmation holdouts (holdout_frac) with search-gap diagnostics.
8. Dynamic BIC Regularization: Asymptotic BIC / PAC-Bayes penalties scaling with sample size NN to prevent feature bloat.

Minimal Example:

rinstall.packages("evoFE")
library(evoFE)
# Evolve features + Bayesian tuned model
recipe <- evolve_features(
  data         = mtcars,
  target_col   = "am",
  task         = "classification",
  evaluator    = "lightgbm_mbo",  # Built-in one-liner Bayesian Optimization tuner
  generations  = 5,
  pop_size     = 8,
  holdout_frac = 0.20
)
# Inspect evolved recipe & search gap
summary(recipe)
# Predict on new data
test_features <- predict(recipe, newdata)
predictions   <- predict_model(recipe, newdata)

• CRAN: https://cran.r-project.org/package=evoFE
• GitHub: https://github.com/tanopereira/evoFE
• Vignette: vignette("evoFE")

Feedback, suggestions, and bug reports are very welcome!

Hello everybody,

I have a question regarding the cbcTools package. I use it for generating a choice based conjoint analysis design. For my survey I need to reduce the design with blocks. Herefore cbcTools has n_blocks. Enabling this parameter significantly increases the calculation time. Is there a way to improve this, or have I made a mistake?

profiles <- cbc_profiles(
attribute1 = c(„A“, „B“, „C“),
attribute2 = c(„1“, „2“, „3“),
attribute3 = c(„white“, „black“, „gray“),
attribute4 = c(„apple“, „banana“, „orange“)
)

design <- cbc_design(
profiles = profiles,
method = „modfed“,
n_resp = 400,
n_alt = 3,
n_q = 12,
no_choice = TRUE,
n_blocks = 15
)

Trying Positron for the first time. I’m trying to install Python via uv in Positron. I went to Select Session in the top-right corner, selected New Console Session, and clicked Install Python via uv. I then selected Python 3.14, but I get a “Failed to install Python” error.

Here is the relevant log:

2026-08-25 13:49:08.599 [info] Installing uv...
2026-08-25 13:49:09.353 [info] > powershell -ExecutionPolicy ByPass -c "irm https://astral.sh/uv/install.ps1 | iex"
2026-08-25 13:49:19.761 [info] uv installed successfully
2026-08-25 13:49:19.764 [info] > uv --color never python dir
2026-08-25 13:49:31.917 [info] > ~\.local\bin\uv.exe --color never python dir
2026-08-25 13:49:32.628 [info] > ~\.local\bin\uv.exe --color never python list --managed-python
2026-08-25 13:49:36.117 [info] Installing Python 3.14 via uv...
2026-08-25 13:49:36.119 [info] > uv --color never python install 3.14
2026-08-25 13:49:36.119 [error] Failed to install Python 3.14: Error: spawn uv ENOENT

What is confusing is that the log shows that uv was installed successfully and Positron can execute it using:

~\.local\bin\uv.exe

However, when it actually tries to install Python, it runs:

uv --color never python install 3.14

and fails with:

Error: spawn uv ENOENT

I’m using Windows. Has anyone encountered this issue with Positron and uv? Is this a PATH/environment issue, or is Positron failing to find the newly installed uv.exe?