R

Personal R tricks and code snippets.

echo "R_BROWSER=zen-browser" >> ~/.Renviron

library(readxl) df <- read_excel("supplement-v2-for_BNX.xlsx", sheet = "Case List", skip = 5, n_max = 26)

df <- distinct(df)

df %>% arrange( desc(value) ) %>% distinct()

it cannot be done with distinct, one has to group and summarize

df %>% group_by(id) %>% summarize( mean_value = mean(value) )

df <- bind_rows( df1, df2, df3, df4 )

library(lubridate)

rodne_cislo_to_dob <- function( rc ) {
	# Run df %>% mutate( dob = map_vec( rc, rodne_cislo_to_dob ) )
	if (is.na(rc)) {
		return(NA);
	}

	rc_year  <- as.numeric(str_sub(rc, 1, 2));
	rc_month <- as.numeric(str_sub(rc, 3, 4));
	rc_day   <- as.numeric(str_sub(rc, 5, 6));
	date <- Sys.Date();

	new_format <- FALSE;
	if (rc_month > 70) {
		rc_month <- rc_month - 70;
		new_format <- TRUE;
	} else if (rc_month > 50) {
		rc_month <- rc_month - 50;
	} else if (rc_month > 20) {
		rc_month <- rc_month - 20;
		new_format <- TRUE;
	}

	if (rc_day > 50 ) {
		rc_day <- rc_day - 50;
	}

	# 9mistne ci 10mistne rodne cislo?
	# 10mistna rodna cisla byla zavedena v roce 1954
	if (str_length(rc) == 9) {
		# pred rokem 1954
		add_year <- 1900;
	} else if (str_length(rc) == 10) {
		if (new_format) {
			add_year <- 2000;
		} else {
			# po roku 1954
			this_year <- as.numeric(date, 3, 4);
			if (rc_year < 54) {
				# pri 1900 by osobe bylo vice nez 100 let
				# spise bude add_year 2000
				add_year <- 2000;
			} else {
				add_year <- 1900;
			}
		}
	} else {
		warn("Neocekavana delka rodneho cisla!!!");
		return(NA);
	}

	rc_year <- rc_year + add_year;
	year(date)  <- rc_year;
	month(date) <- rc_month;
	day(date)   <- rc_day;
	return(date);
}

 df %>% mutate( dob = map_vec( rc, rodne_cislo_to_dob ) )

source("script.R")
savehistory("script.R")

library(lintr)
lint(filename="analyze.R")

colnames(df) <- make.names(str_trim(colnames(df)))

df <- df %>%
  mutate(
    new_column = str_glue("{column1} ({column2}–{column3}) [{column4}±{column5}]")
  )

render

umrti %>%
    ggplot(aes(y = fct_infreq(oddeleni), fill = oddeleni)) +
    geom_bar() +
    xlab("počet úmrtí") +
    ylab("oddělení") +
    ggtitle("Počet úmrtí 2021–2025") +
    theme(legend.position = "none")

library(forcats)
df %>%
  ggplot( aes( x = fct_rev(fct_reorder(angiolog, groin.to.open, median)), y = groin.to.open)) +
  geom_boxplot() +
  ggtitle("Groin to open podle angiologa") +
  ylab("groin to open [min]") +
  xlab("angiolog")

library(forcats)
df %>%
    ggplot(aes(y = fct_infreq(oddeleni), fill = typ)) +
    geom_bar() +
    xlab("typ propuštění") +
    ylab("oddělení") +
    ggtitle("Propuštění 2021–2025")

g + theme( axis.text.x = element_text( angle = 90, vjust = 0.5, hjust=1) )

library(ggpubr)

ggboxplot( df, x = "group", y = "age", color = "group", add = "jitter") +
    stat_compare_means()

To change the test: stat_compare_means( method = 't.test' ). Method can be t.test, wilcox.test, anova, kruskal.test.

To get pairwise comparisons:

comparisons <- list( c("group1", "group2"), c("group2", "group3"), c("group1", "group3" ))
ggboxplot( df, x = "group", y = "age", color = "group", add = "jitter" ) +
    stat_compare_means( comparisons = comparisons, label.y = c(58, 60, 62) ) +
    stat_compare_means( method = 'kruskal.test', label.y = 66 )

First stat_compare_means does pairwise comparison, second stat_compare_means does global test.

df %>% 
    ggplot( aes( x = as.factor(rok) ) ) + 
    geom_bar() +
    stat_count(geom = "text", size = 4, aes(label = ..count..), position=position_stack(vjust=1.15)) +
    facet_grid(sluzba~.) +
    ggtitle("Počet rekanalizací (IVT+MT) v jednotlivých letech dle pracovní doby") +
    ylab("počet") +
    xlab("rok")

Error in dyn.load(file, DLLpath = DLLpath, ...) : unable to load shared object '/home/pavel/R/x86_64-pc-linux-gnu-library/4.4/stringi/libs/stringi.so'

to solve: install.packages("stringi")

getOption('timeout')
options(timeout=120)

unable to load shared object '/home/pavel/R/x86_64-pc-linux-gnu-library/4.5/xml2/libs/xml2.so'

sudo pacman -S libxml2-legacy

df <- df %>%
    mutate(
        tier = ifelse(
            `STIM TYPE` == "SHAM",
            "SHAM",
            substr( `Study ID`, start = 1, stop = 1)
        )
    ) %>%
    drop_na(`Study ID`)

Descriptive statistics by group:

library(skimr)
df %>% group_by(group) %>% skim()

 df %>% 
    summarize(
              "DTG_mean" = mean(door.to.groin, na.rm = TRUE),
              "DTG_sd" = sd(door.to.groin, na.rm = TRUE),
              "DTG_median" = median(door.to.groin, na.rm = TRUE),
              "DTG_Q1" = quantile(door.to.groin, p = 0.25, na.rm = TRUE),
              "DTG_Q3" = quantile(door.to.groin, p = 0.75, na.rm = TRUE),

              "GTO_mean" = mean(groin.to.open, na.rm = TRUE),
              "GTO_sd" = sd(groin.to.open, na.rm = TRUE),
              "GTO_median" = median(groin.to.open, na.rm = TRUE),
              "GTO_Q1" = quantile(groin.to.open, p = 0.25, na.rm = TRUE),
              "GTO_Q3" = quantile(groin.to.open, p = 0.75, na.rm = TRUE),

              "DTN_mean" = mean(door.to.needle, na.rm = TRUE),
              "DTN_sd" = sd(door.to.needle, na.rm = TRUE),
              "DTN_median" = median(door.to.needle, na.rm = TRUE),
              "DTN_Q1" = quantile(door.to.needle, p = 0.25, na.rm = TRUE),
              "DTN_Q3" = quantile(door.to.needle, p = 0.75, na.rm = TRUE),
    ) %>%
    mutate(
        door_to_groin = str_glue("{DTG_median} ({DTG_Q1}–{DTG_Q3}) [{DTG_mean}±{DTG_sd}]"),
        groin_to_open = str_glue("{GTO_median} ({GTO_Q1}–{GTO_Q3}) [{GTO_mean}±{GTO_sd}]"),
        door_to_needle = str_glue("{DTN_median} ({DTN_Q1}–{DTN_Q3}) [{DTN_mean}±{DTN_sd}]"),
    ) %>% 
    select(door_to_groin, groin_to_open, door_to_needle)

{{R kruskal.test(dNIHSS90 ~ tier, data = df) %>% print()

#! one-sided p-values!!! library(dunn.test) dunn.test(df\(dNIHSS90, df\)tier) dunn.test(df\(dNIHSS90, df\)tier, kw=TRUE, method="bonferroni")

#two sided p-values!!! library(FSA) dunnTest( dNIHSS90 ~ tier, df) dunnTest( dNIHSS90 ~ tier, df, method="bonferroni") }}

test <- df %>% group_by(group) %>% select( group, sex ) %>% table() %>% fisher.test()

library(broom)
tidy(test)$p.value

test <- aov(age ~ group, df)

library(broom)
p <- tidy(test)$p.value

`test <- kruskal.test( age ~ group, df )

library(broom) p <- tidy(test)$p.value`

fisher_p_values <- tibble(
    test = c( "humoral_month3", "humoral_month12", "cellular_month3", "cellular_month12"),
    p_value = c(
        response %>%
        filter( month == 3 ) %>%
        select( group, humoral_response ) %>%
        table() %>%
        fisher.test() %>%
        tidy() %>%
        select(p.value) %>%
        pull(),

        response %>%
        filter( month == 12 ) %>%
        select( group, humoral_response ) %>%
        table() %>%
        fisher.test() %>%
        tidy() %>%
        select(p.value)
        %>% pull(),

        response %>%
        filter( month == 3 ) %>%
        select( group, cellular_response ) %>%
        table() %>%
        fisher.test() %>%
        tidy() %>%
        select(p.value) %>%
        pull(),

        response %>%
        filter( month == 12 ) %>%
        select( group, cellular_response ) %>%
        table() %>%
        fisher.test() %>%
        tidy() %>%
        select(p.value) %>%
        pull()
    )) %>% mutate(
        Bonferroni = p.adjust( p_value, method = "bonferroni"),
        Benjamini_Hochberg = p.adjust( p_value, method = "BH" ),
    )

library(tidyverse)
library(readxl)
library(lubridate)
library(ggplot2)
library(lme4)
library(brms)
library(car)
library(skimr)
library(broom)
library(ggpubr)

print(version)
packageVersion("tidyverse")
packageVersion("readxl")
packageVersion("lubridate")
packageVersion("ggplot2")
packageVersion("lme4")
packageVersion("brms")
packageVersion("car")
packageVersion("skimr")
packageVersion("broom")
packageVersion("ggpubr")

df_long <- df %>%
    pivot_longer(
	cols = starts_with( "humoral_" ),
	names_to = "month",
	names_prefix = "humoral_M",
	values_to = "humoral_response"
    ) %>%
    mutate( month = as.numeric(month) )

df_joined <- df1 %>%
    inner_join(
        df2,
        by = join_by( subject_id == subject_id, month == month, group == group, sex == sex, age == age )
    ) %>% select(
           subject_id,
           group,
           sex,
           age,
           month,
           cellular_response,
           humoral_response,
   )

library(lme4)
model <- lmer( outcome ~ group + month + age + sex + (1 | subject_id), data = df )
Anova(model)

Subject_id as a random effect:

library(lme4)
model <- glmer( binary_outcome ~ group + month + age + sex + (1 | subject_id), data = df, family = binomial )
Anova(model)

library(brms)
model <- brm( outcome ~ group + sex + age + month + (1|subject_id), data = response )
plot(model)
posterior_interval(model)

library(brms)
model <- brm( binary_outcome ~ group + sex + age + month + (1|subject_id), family = bernoulli, data = response )
plot(model)
posterior_interval(model)

pedigree.csv:

id,fid,mid,sex,label
John Donn *1615,0,0,1,*28.6.1615 +25.12.1680
Anne Bowl *1630,0,0,2,*14.7.1626 +12.8.1682
Philip Donn *1652,John Donn *1615,Anne Bowl *1630,1,*3.2.1652

pedigree.R:

library(pedtools)

# Read data
x <- readPed("pedigree.csv", colSep = ",", header = TRUE, comment.char="#", colSkip = 5)
raw <- read.csv("pedigree.csv", sep=",", header = TRUE, comment.char="#")

# Plot data PDF
cairo_pdf("pedigree.pdf", width = 14, height = 10)
plot(x, proband = "", cex = 0.6, symbolsize = 0.8, margins = c(2, 3, 2, 3))
dev.off()

# Plot data PNG
png("pedigree.png", width = 4000, height = 1600, res = 200)
plot(x, proband = "", cex = 1.0, symbolsize = 0.8, margins = c(2, 3, 2, 3))
dev.off()