From fde1963b79a159dff12cfd7580842699b6c61a07 Mon Sep 17 00:00:00 2001 From: Theertha Kariyathan Date: Mon, 12 Oct 2026 12:02:08 +0100 Subject: [PATCH 1/2] dslagsarlm --- R/ds.lagsarlm.R | 148 +++++++++++++++++++++++++++++++ vignettes/Opal_functioncheck.Rmd | 6 +- 2 files changed, 151 insertions(+), 3 deletions(-) create mode 100644 R/ds.lagsarlm.R diff --git a/R/ds.lagsarlm.R b/R/ds.lagsarlm.R new file mode 100644 index 0000000..7ce698e --- /dev/null +++ b/R/ds.lagsarlm.R @@ -0,0 +1,148 @@ +# library(spatialreg) +# +# model <- lagsarlm( +# obesity_qof_2024_25 ~ dist_food_fastfood + dist_food_supermarket + dist_hc_emergency, +# data = cells, +# listw = W[["band800"]], +# zero.policy = TRUE, +# method = "Matrix" +# ) +# +#' Calculate Moran's I +#' +#' Calculates Moran's I for a numeric variable using a spatial weights +#' object across one or more DataSHIELD studies. +#' +#' @param x Character string specifying the name of the numeric vector on the +#' server to be analysed. +#' @param listw Character string specifying the name of the spatial weights +#' object of class \code{listw} on the server. +#' @param type Character string specifying the output format. Currently only +#' \code{"split"} is supported. +#' @param vars Optional variable specification. Currently unused. +#' @param do.checks Logical. If \code{TRUE}, checks that the input objects are +#' defined and have consistent classes across the studies. +#' @param nsim Numeric value specifying the number of simulations used for +#' calculating the Monte Carlo p-value. Default is 199. +#' @param datasources A list of \code{DSConnection} objects. If \code{NULL}, +#' active DataSHIELD connections are used. +#' +#' @return A named list with one element per study. Each study contains a +#' named numeric vector with: +#' \itemize{ +#' \item \code{moran_I}: Moran's I statistic. +#' \item \code{z}: standardised Moran's I statistic. +#' \item \code{p_mc}: Monte Carlo p-value. +#' } +#' +#' @examples +#' \dontrun{ +#' DSOpal::Opal() +#' builder <- DSI::newDSLoginBuilder() +#' builder$append( +#' server = "mersey-demo", +#' url = "http://localhost:8880", +#' user = "administrator", +#' password = "password", +#' driver = "OpalDriver", +#' profile = "ds-geospatial" +#' ) +#' +#' logindata <- builder$build() +#' conns <- DSI::datashield.login(logins = logindata) +#' DSI::datashield.assign.resource( +#' conns, +#' symbol = "cells", +#' resource = paste0(project, ".", "cells") +#' ) +#' DSI::datashield.assign.resource( +#' conns, +#' symbol = "W", +#' resource = paste0(project, ".", "W_D_5km") +#' ) +#' ## 2.8 Materialise arbitrary R objects ---- +#' DSI::datashield.assign.resource( +#' conns, +#' symbol = "cells_resource", +#' resource = paste0(project, ".cells") +#' ) +#' +#' DSI::datashield.assign.expr( +#' conns, +#' symbol = "cells_object", +#' expr = quote(as.resource.object(cells_resource)) +#' ) +#' ## 2.9 Map specific columns ---- +#' DSI::datashield.assign.expr( +#' conns, +#' symbol = "obesity_qof_2024_25", +#' expr = quote(cells_object$obesity_qof_2024_25) +#' ) +#' DSI::datashield.assign.expr( +#' conns, +#' symbol = "W_object", +#' expr = quote(as.resource.object(W)) +#' ) +#' DSI::datashield.assign.expr( +#' conns, +#' symbol = "W_band800", +#' expr = quote(W_object$band800) +#' ) +#' +#' ds.moransI( +#' x = "obesity_qof_2024_25", +#' listw = "W_band800", +#' nsim = 199, +#' datasources = conns +#' ) +#' +#' +#' @export + +ds.lagsarlm <- function(formula = NULL, data = NULL, listw = NULL, type='split', + do.checks = TRUE, + datasources=NULL){ + + # look for DS connections + if(is.null(datasources)){ + datasources <- datashield.connections_find() + } + + # ensure datasources is a list of DSConnection-class + if(!(is.list(datasources) && all(unlist(lapply(datasources, function(d) {methods::is(d,"DSConnection")}))))){ + stop("The 'datasources' were expected to be a list of DSConnection-class objects", call.=FALSE) + } + + if(is.null(formula)){ + stop("Please provide the input formula!", call.=FALSE) + } + + if(is.null(listw)){ + stop("Please provide the input vector 'listw'!, a list object", call.=FALSE) + } + + if(do.checks){ + + # check if the input objects are defined in all the studies + #isDefined(datasources, x) + isDefined(datasources, listw) + + # call the internal function that checks the input object is of the same class in all studies. + #typ1 <- checkClass(datasources, x) + typ2 <- checkClass(datasources, listw) + } + browser() + + calltext <- paste0("lagsarlmDS(", deparse(formula), ", ", listw, ", ", data, ")", collapse = "") + + output <- DSI::datashield.aggregate(datasources, as.symbol(calltext)) + + ## Only type split supported now + if(type == "split"){ + + result <- lapply(output, function(x) unlist(x)) + return(result) + } +} + + diff --git a/vignettes/Opal_functioncheck.Rmd b/vignettes/Opal_functioncheck.Rmd index 533984e..16b2ed4 100644 --- a/vignettes/Opal_functioncheck.Rmd +++ b/vignettes/Opal_functioncheck.Rmd @@ -35,7 +35,7 @@ dat <- qof_asthma_cohort_dist_gs_sf |> # Uploading a subset to save upload time -dat <- dat |> dplyr::filter(lsoa11cd %in% sample(levels(factor(dat$lsoa11cd)), 35)) +#dat <- dat |> dplyr::filter(lsoa11cd %in% sample(levels(factor(dat$lsoa11cd)), 35)) dat$row_id <- seq_len(nrow(dat)) # needs a unique row id ``` @@ -122,8 +122,8 @@ ds.summarystat(x = "D$has_asthma", y = "D$lsoa11cd", ```{r} ds.groupcor(x = "D$has_asthma", y = "D$distance_local_greenspace", - index = "D$lsoa11cd", datasources = conns, type = "combine", - names_region = "Liverpool") + index = "D$lsoa11cd", datasources = conns, type = "split", + names_region = names_region) ``` ```{r} From 67225f0e1b1c676d142bf88b90b6e32a649e928b Mon Sep 17 00:00:00 2001 From: Theertha Kariyathan Date: Mon, 12 Oct 2026 15:52:52 +0100 Subject: [PATCH 2/2] update opal function check --- vignettes/Opal_functioncheck.Rmd | 108 +++++++++++++++++++++++++++++-- 1 file changed, 102 insertions(+), 6 deletions(-) diff --git a/vignettes/Opal_functioncheck.Rmd b/vignettes/Opal_functioncheck.Rmd index 16b2ed4..d5051b7 100644 --- a/vignettes/Opal_functioncheck.Rmd +++ b/vignettes/Opal_functioncheck.Rmd @@ -3,6 +3,7 @@ title: "R Notebook" output: html_notebook --- +Build package locally. ```{r, eval=FALSE} # r console @@ -21,9 +22,106 @@ devtools::build() # go to administration --> R see if the profiles appear # go to administration --> Datashield activate profile and initialize ``` -This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. When you execute code within the notebook, the results appear beneath the code. -Try executing this chunk by clicking the *Run* button within the chunk or by placing your cursor inside it and pressing *Cmd+Shift+Enter*. +Create example data set. + +```{r, eval=FALSE} +# retrieve LSOA boundaries for Cheshire and Merseyside +# to reduce coverage, remove names from the `within_names` +cm_lsoa_sf <- boundr::bounds( + "lsoa", + within_level = "lad", + within_names = c( + "Cheshire East", + "Cheshire West and Chester", + "Halton", + "Knowsley", + "Liverpool", + "Sefton", + "St. Helens", + "Warrington", + "Wirral" + ), + lookup_year = 2011, # can be either 2011 or 2021 + opts = boundr::boundr_options(resolution = "BFC") +) + +liverpool_sf <- cm_lsoa_sf |> + dplyr::filter(lad22nm %in% "Liverpool") + +uprn_gs_sf <- "https://pldr.org/download/2k6r3/n12/UPRN_2_1_greenspace_distances_with_coords.csv" |> + readr::read_csv() |> + # convert to spatial object + sf::st_as_sf(coords = c("longitude", "latitude"), crs = 4326) |> + # add LSOA details + sf::st_join(cm_lsoa_sf) + + + +qof_asthma_tbl <- "https://pldr.org/download/e6nzv/ng1/QOF_4_03_Asthma_LSOA.csv" |> + readr::read_csv() |> + # filter out CM LSOAs + dplyr::filter(lsoa11 %in% cm_lsoa_sf$lsoa11cd) |> + # filter out last year + dplyr::filter(year == 2024) + + + +# set seed for reproducibility +set.seed(20250603) +qof_asthma_cohort_tbl <- qof_asthma_tbl |> + purrr::pmap(function(lsoa11, den, num, ...) { + # recalculate pop to be a multiple of 2 + pop <- ifelse(den %% 2 != 0, den + 1, den) + # extract UPRNs for this LSOA + uprn_tbl <- uprn_gs_sf |> + sf::st_drop_geometry() |> + dplyr::filter(lsoa11cd == lsoa11) |> + dplyr::mutate( + dist_gs_dec = dplyr::ntile(distance_local_greenspace, 10) + ) |> + dplyr::arrange(dist_gs_dec) + + # create cohort table + cohort_tbl <- tibble::tibble( + lsoa11cd = lsoa11, + sex = rep(c(0, 1), each = pop / 2), + has_asthma = FALSE, + uprn = sample(uprn_tbl$UPRN, size = pop, replace = TRUE) + ) + + # randomly assigned the outcome to `num` patients + idx <- sample(seq_len(pop), ceiling(num), replace = FALSE) + cohort_tbl$has_asthma[idx] <- TRUE + + # randomly re-assigned UPRNs with higher distances to patients with outcome + ## create vector of probabilities, so further distances are favoured + p <- c(0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.1, 0.2, 0.2, 0.2) + idx_uprn_tile <- sample(1:10, ceiling(num), replace = TRUE, prob = p) + ## subset UPRNs with match tile as per `idx_uprn_tile` + cohort_tbl$uprn[idx] <- purrr::map(idx_uprn_tile, function(x) { + # subset based on tile decile + aux <- dplyr::filter(uprn_tbl, dist_gs_dec == x) + # extract UPRN + aux$UPRN[sample(seq_len(nrow(aux)), 1)] + }) |> + purrr::list_c() + + # return final cohort table + return(cohort_tbl) + }) |> + purrr::list_c() + + +qof_asthma_cohort_dist_gs_sf <- qof_asthma_cohort_tbl |> + dplyr::left_join(uprn_gs_sf, by = c("uprn" = "UPRN", "lsoa11cd")) |> + dplyr::arrange(lsoa11cd, uprn) |> + sf::st_as_sf() + + +``` + +Subset the data for needed coloumns. ```{r, eval=FALSE} library(opalr) @@ -40,6 +138,7 @@ dat$row_id <- seq_len(nrow(dat)) # needs a unique row id ``` +Creating project and uploading table. ```{r} library(opalr) @@ -76,6 +175,7 @@ opal.table_save( ``` +Login ```{r} library(DSOpal) library(dsBase) @@ -143,8 +243,4 @@ fitdat <- ds.glmerSLMA(formula="D$has_asthma ~ D$sex + (1|D$lsoa11cd)", family="poisson",datasources=conns ) ``` -Add a new chunk by clicking the *Insert Chunk* button on the toolbar or by pressing *Cmd+Option+I*. - -When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the *Preview* button or press *Cmd+Shift+K* to preview the HTML file). -The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike *Knit*, *Preview* does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.