高校数学の質問スレ（医者・東大卒専用） Part438 (882ﾚｽ)
上下前次1-新
抽出解除 必死ﾁｪｯｶｰ(本家) (べ) 自ID ﾚｽ栞 あぼーん
826: 2025/05/05(月) 19:07:36.02 ID:LDY1RQtT(1/4)調 AAS
# This R function, riskratio.boot, calculates the risk ratio and its Highest Density Interval (HDI)
# using a bootstrap method. It takes the number of events and the total number of observations
# for two groups as input.
riskratio.boot <- function(r1, r2, n1, n2, nboot = 5000, conf.level = 0.95, verbose = FALSE){
# Combine the number of events and total observations for both groups.
r <- c(r1, r2)
n <- c(n1, n2)
# Create a 2x2 matrix representing the contingency table (number of events and non-events).
m <- cbind(r, n - r) ; m
# Perform bootstrap resampling to estimate the risk ratio distribution.
rr <- replicate(nboot, {
# Simulate the number of events in the first group by sampling with replacement
# from a population with the observed event rate.
sample(c(rep(1, r1), rep(0, n1 - r1)), n1, replace = TRUE) |> sum() -> R1
# Simulate the number of events in the second group similarly.
sample(c(rep(1, r2), rep(0, n2 - r2)), n2, replace = TRUE) |> sum() -> R2
# Calculate the risk ratio from the bootstrapped event counts.
(R1 / n1) / (R2 / n2)
})
# If verbose is TRUE, plot the posterior distribution of the risk ratio.
if(verbose){
source('plotpost.R') # Assuming 'plotpost.R' is a script for plotting posterior distributions.
plotpost(rr, compVal = 1, showCurve = 1, lwd = 4) # Plot the distribution with a comparison value of 1.
}
# Calculate the mean of the bootstrapped risk ratios.
b_mean <- mean(rr)
# Calculate the Highest Density Interval (HDI) of the bootstrapped risk ratios.
b_ci <- HDInterval::hdi(rr, credMass = conf.level)
# Return a list containing the bootstrap mean and the HDI.
list(b_mean = b_mean, b_ci = b_ci)
}
# Example usage of the riskratio.boot function.
riskratio.boot(244, 282, 2345, 2333)
827: 2025/05/05(月) 19:08:42.51 ID:LDY1RQtT(2/4)調 AAS
Description:
The `riskratio.boot` function in R estimates the risk ratio between two groups and provides its Highest Density Interval (HDI) using a bootstrap resampling approach. It takes the counts of events and the total number of observations for each of the two groups as input.
Usage:
riskratio.boot(r1, r2, n1, n2, nboot = 5000, conf.level = 0.95, verbose = FALSE)
Arguments:
* `r1`: The number of events in the first group.
* `r2`: The number of events in the second group.
* `n1`: The total number of observations in the first group.
* `n2`: The total number of observations in the second group.
* `nboot`: The number of bootstrap replicates to perform (default is 5000). A larger number of replicates generally provides more stable estimates.
* `conf.level`: The credibility level for the Highest Density Interval (HDI) (default is 0.95, corresponding to a 95% HDI).
* `verbose`: A logical value indicating whether to display a plot of the bootstrap distribution of the risk ratio (default is `FALSE`). If `TRUE`, it assumes a script named `plotpost.R` is available in the working directory for plotting.
828: 2025/05/05(月) 19:09:02.55 ID:LDY1RQtT(3/4)調 AAS
Details:
The function works by simulating the event outcomes in each group through bootstrap resampling. For each group, it draws `n1` (or `n2`) samples with replacement from a hypothetical population that has the observed proportion of events (`r1/n1` or `r2/n2`). The number of events in each resampled set (`R1` and `R2`) is then used to calculate a bootstrapped risk ratio `(R1/n1) / (R2/n2)`. This process is repeated `nboot` times to generate a distribution of risk ratios. The function then calculates the mean of this distribution and its Highest Density Interval (HDI), which represents the most credible range for the true risk ratio given the data and the bootstrap procedure.
If `verbose` is set to `TRUE`, the function attempts to plot the distribution of the bootstrapped risk ratios using a script named `plotpost.R`. This requires that the `plotpost.R` script exists in the current working directory and is capable of handling the vector of bootstrapped risk ratios.
Value:
The function returns a list with the following components:
* `b_mean`: The mean of the bootstrapped risk ratios.
* `b_ci`: The Highest Density Interval (HDI) of the bootstrapped risk ratios, as a vector with the lower and upper bounds.
Note:
The `verbose = TRUE` option depends on an external script `plotpost.R`. If you intend to use this option, ensure that the script is available and correctly implemented for plotting posterior-like distributions. The HDI is calculated using the `hdi` function from the `HDInterval` package, so this package must be installed (`install.packages("HDInterval")`) and loaded (`library(HDInterval)`) if you intend to use the default behavior.
829: 2025/05/05(月) 19:23:25.28 ID:LDY1RQtT(4/4)調 AAS
> riskratio.boot(244,282,2345,2333,nboot=10000)
$b_mean
[1] 0.8641319
$b_ci
lower upper
0.7312212 1.0106121
attr(,"credMass")
[1] 0.95
上下前次1-新書関写板覧索設栞歴
ｽﾚ情報 赤ﾚｽ抽出 画像ﾚｽ抽出 歴の未読ｽﾚ
ぬこの手 ぬこTOP 0.044s