+ "markdown": "---\ntitle: Functional programming examples\nformat: html\nexecute: \n echo: true\n---\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(ggplot2)\nlibrary(tidyr)\nlibrary(dplyr)\nlibrary(purrr)\n\ndat <- iris\n\n# functions\n\nremove_rownames <- function(x){\n rownames(x) <- NULL\n x\n}\n\nsum_lm <- function(x, conf.level = 0.95){\n xs <- data.frame(summary(x)$coefficients)\n xs$param <- rownames(xs)\n rownames(xs) <- NULL\n names(xs) <- c(\"b\", \"se\", \"t\", \"pval\", \"param\")\n cc <- data.frame(confint(x, level = conf.level))\n names(cc) <- c(\"ci.lb\", \"ci.ub\")\n rownames(cc) <- NULL\n cbind(\n xs[, c(\"param\", \"b\", \"se\", \"t\", \"pval\")],\n cc\n )\n}\n\n# fit a linear model for each Species\n\ndatl <- split(dat, dat$Species)\nfitl <- lapply(datl, function(d) lm(Sepal.Length ~ Petal.Width + Petal.Length, data = d))\nresl <- lapply(fitl, sum_lm)\n\nres_by_species <- do.call(rbind, resl)\nres_by_species <- remove_rownames(res_by_species)\nres_by_species$species <- rep(names(resl), sapply(resl, nrow))\n\nggplot(res_by_species, aes(x = param, \n y = b, \n ymin = ci.lb, \n ymax = ci.ub,\n color = species)) +\n geom_pointrange(position = position_dodge(width = 0.5))\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n\n```{.r .cell-code}\n# alternative version using other packages\n\nresl <- lapply(fitl, broom::tidy, conf.int = TRUE)\nres_by_species <- dplyr::bind_rows(resl, .id = \"species\")\n\nres_by_species\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n# A tibble: 9 × 8\n species term estimate std.error statistic p.value conf.low conf.high\n <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>\n1 setosa (Intercep… 4.25 0.411 10.3 1.13e-13 3.42 5.08 \n2 setosa Petal.Wid… 0.712 0.487 1.46 1.51e- 1 -0.268 1.69 \n3 setosa Petal.Len… 0.399 0.296 1.35 1.84e- 1 -0.196 0.994\n4 versicolor (Intercep… 2.38 0.449 5.30 3.04e- 6 1.48 3.28 \n5 versicolor Petal.Wid… -0.320 0.402 -0.795 4.30e- 1 -1.13 0.489\n6 versicolor Petal.Len… 0.934 0.169 5.52 1.44e- 6 0.594 1.27 \n7 virginica (Intercep… 1.05 0.514 2.05 4.63e- 2 0.0179 2.09 \n8 virginica Petal.Wid… 0.00706 0.179 0.0394 9.69e- 1 -0.354 0.368\n9 virginica Petal.Len… 0.995 0.0893 11.1 8.87e-15 0.815 1.17 \n```\n\n\n:::\n\n```{.r .cell-code}\n# bootstrapping\n\nboot <- function(data, B = 100){\n res <- vector(mode = \"list\", length = B)\n n <- nrow(data)\n for(i in 1:B){\n idx <- sample(x = 1:n, size = n, replace = TRUE)\n dataB <- data[idx, ]\n rownames(dataB) <- NULL\n res[[i]] <- dataB\n }\n return(res)\n}\n\nfit_lm <- function(data){\n lm(Sepal.Length ~ Petal.Width + Petal.Length, data = data)\n}\n\nbootl <- lapply(datl, boot, B = 100)\nfit_bootl <- lapply(bootl, function(x) lapply(x, fit_lm))\nres_bootl <- lapply(fit_bootl, function(x) lapply(x, sum_lm))\nres_bootl <- lapply(res_bootl, function(x) do.call(rbind, x))\nres_boot_by_species <- do.call(rbind, res_bootl)\nres_boot_by_species <- remove_rownames(res_boot_by_species)\nres_boot_by_species$species <- rep(names(res_bootl), sapply(res_bootl, nrow))\n\nggplot(res_boot_by_species, aes(x = param, y = b, fill = species)) +\n geom_boxplot()\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n\n```{.r .cell-code}\n# using nested tibbles\n\ndat |> \n group_by(Species) |> \n nest() |> \n mutate(boot = map(data, boot)) |> \n select(-data) |> \n unnest(boot) |> \n mutate(fit = map(boot, fit_lm)) |> \n mutate(res = map(fit, sum_lm)) |> \n select(-boot, -fit) |> \n unnest(res)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n# A tibble: 900 × 8\n# Groups: Species [3]\n Species param b se t pval ci.lb ci.ub\n <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>\n 1 setosa (Intercept) 4.27 0.369 11.6 2.20e-15 3.53 5.01 \n 2 setosa Petal.Width 1.57 0.596 2.64 1.12e- 2 0.376 2.77 \n 3 setosa Petal.Length 0.280 0.265 1.06 2.96e- 1 -0.253 0.813\n 4 setosa (Intercept) 3.62 0.460 7.87 4.02e-10 2.70 4.55 \n 5 setosa Petal.Width 0.457 0.575 0.796 4.30e- 1 -0.699 1.61 \n 6 setosa Petal.Length 0.847 0.353 2.40 2.05e- 2 0.136 1.56 \n 7 setosa (Intercept) 4.17 0.473 8.81 1.62e-11 3.22 5.12 \n 8 setosa Petal.Width 0.297 0.469 0.632 5.30e- 1 -0.648 1.24 \n 9 setosa Petal.Length 0.591 0.348 1.70 9.63e- 2 -0.110 1.29 \n10 setosa (Intercept) 4.22 0.431 9.78 6.55e-13 3.35 5.08 \n# ℹ 890 more rows\n```\n\n\n:::\n:::\n",
0 commit comments