Как использовать функцию репликации в R, чтобы повторить функцию

Question

У меня проблема при использовании replicate для повторения функции.

Я попытался использовать bootstrap для соответствия модели квадратичного c, используя концентрацию в качестве предиктора и Total_lignin в качестве ответа и собираясь сообщить оценку максимума с соответствующей стандартной ошибкой.

Моя идея состоит в том, чтобы создать функцию под названием bootFun, которая, по сути, делала все в пределах одной итерации a для l oop. bootFun принял только набор данных для предсказателя и ответ для использования (оба имени переменных в кавычках).

Однако SD равно 0, что неверно. Я не знаю, где не то место. Не могли бы вы помочь мне с этим?

# Load the libraries
library(dplyr)
library(tidyverse)

# Read the .csv and only use M.giganteus and S.ravennae.
dat <- read_csv('concentration.csv') %>% 
  filter(variety == 'M.giganteus' | variety == 'S.ravennae') %>% 
  arrange(variety)
# Check the data
head(dat)

# sample size
n <- nrow(dat)

# A function to do one iteration
bootFun <- function(dat, pred, resp){
  # Draw the sample size from the dataset
  sample <- sample_n(dat, n, replace = TRUE)

  # A quadratic model fit
  formula <- paste0('resp', '~', 'pred', '+', 'I(pred^2)')
  fit <- lm(formula, data = sample)

  # Derive the max of the value of concentration
  max <- -fit$coefficients[2]/(2*fit$coefficients[3])

  return(max)
}

max <- bootFun(dat = dat, pred = 'concentration', resp = 'Total_lignin' )

# Iterated times
N <- 5000

# Use 'replicate' function to do a loop
maxs <- replicate(N, max)

# An estimate of the max of predictor and corresponding SE
mean(maxs)
sd(maxs)

Answer 1

Базовый пакет boot, функция boot, может упростить задачу повторного вызова функции bootstrap. Первый аргумент должен быть набором данных, второй аргумент - аргументом индекса, который пользователь не устанавливает , и другие аргументы также могут быть переданы ему. В этом случае этими другими аргументами являются имена предиктора и ответа.

library(boot)

bootFun <- function(dat, indices, pred, resp){
  # Draw the sample size from the dataset
  dat.sample <- dat[indices, ]

  # A quadratic model fit
  formula <- paste0(resp, '~', pred, '+', 'I(', pred, '^2)')
  formula <- as.formula(formula)
  fit <- lm(formula, data = dat.sample)

  # Derive the max of the value of concentration
  max <- -fit$coefficients[2]/(2*fit$coefficients[3])

  return(max)
}

N <- 5000
set.seed(1234)  # Make the bootstrap results reproducible

results <- boot(dat, bootFun, R = N, pred = 'concentration', resp = 'Total_lignin')

results
#
#ORDINARY NONPARAMETRIC BOOTSTRAP
#
#
#Call:
#boot(data = dat, statistic = bootFun, R = N, pred = "concentration", 
#    resp = "Total_lignin")
#
#
#Bootstrap Statistics :
#      original        bias    std. error
#t1* -0.4629808 -0.0004433889  0.03014259
#


results$t0  # this is the statistic, not bootstrapped
#concentration 
#   -0.4629808

mean(results$t)  # bootstrap value
#[1] -0.4633233

Обратите внимание, что подобрать многочлен, функция poly намного проще, чем явно записывать полиномиальные члены один за другим.

formula <- paste0(resp, '~ poly(', pred, ',2, raw = TRUE)')

Проверить распределение статистики начальной загрузки c.

op <- par(mfrow = c(1, 2))
hist(results$t)
qqnorm(results$t)
qqline(results$t)
par(op)

Тестовые данные

set.seed(2020)  # Make the results reproducible
x <- cumsum(rnorm(100))
y <- x + x^2 + rnorm(100)
dat <- data.frame(concentration = x, Total_lignin = y)