r 置信区间T2 Bonferroni Hotelling Bootstrap

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了r 置信区间T2 Bonferroni Hotelling Bootstrap相关的知识,希望对你有一定的参考价值。

####
#### Perspiration from 20 healthy females was analyzed. Three components, X_1 = sweat 
#### rate, X_2 = sodium content, and X_3 = potassium content, were measured.
####

## Load data.
X <- read.table("T5-1.DAT", header = FALSE)
colnames(X) <- c("Sweat", "Sodium", "Potassium")
attach(X)

n <- nrow(X)
p <- 3

## Summary statistics.
x_bar <- colMeans(X)
S <- var(X)

##
## Check normality.
##

pdf("figures/sweat_QQ.pdf", width = 7, height = 3.5)
par(mfrow = c(1, 3))
qqnorm(Sweat, main = "Sweat"); qqline(Sweat)
qqnorm(Sodium, main = "Sodium"); qqline(Sodium)
qqnorm(Potassium, main = "Potassium"); qqline(Potassium)
dev.off()

pdf("figures/sweat_pairs.pdf")
pairs(X)
dev.off()

##
## 95% confidence ellipse.
##

## The ellipse is centered at the sample mean vector, with axes equal to the eigenvectors 
## of S.
ee <- eigen(S)
lambda <- ee$values
ee <- ee$vectors

## The scaled F percentile that defines the desired squared distance.
scaled_F <- ((p * (n - 1)) / (n - p)) * qf(0.95, p, n - p)

## The half-lengths of the ellipse.
sqrt((lambda / n) * scaled_F)

##
## 95% T^2 simultaneous confidence intervals for the mean components.
##

x_bar[1] + c(-1, 1) * sqrt(scaled_F) * sqrt(S[1, 1] / n)
x_bar[2] + c(-1, 1) * sqrt(scaled_F) * sqrt(S[2, 2] / n)
x_bar[3] + c(-1, 1) * sqrt(scaled_F) * sqrt(S[3, 3] / n)

##
## 95% Bonferroni simultaneous confidence intervals for the mean components.
##

x_bar[1] + c(-1, 1) * qt(1 - 0.05 / (2 * p), n - 1) * sqrt(S[1, 1] / n)
x_bar[2] + c(-1, 1) * qt(1 - 0.05 / (2 * p), n - 1) * sqrt(S[2, 2] / n)
x_bar[3] + c(-1, 1) * qt(1 - 0.05 / (2 * p), n - 1) * sqrt(S[3, 3] / n)

##
## Hotelling's T^2 test of H_0: mu' = [3.5, 54.5, 8.8]. This is equivalent to checking 
## whether mu_0 is inside the 95% confidence ellipse from above. It is, so, as expected, 
## we fail to reject H_0.
##

mu_0 <- c(4.0, 45.0, 10.0)
T2 <- n * t(x_bar - mu_0) %*% solve(S) %*% (x_bar - mu_0)
p_value <- 1 - pf(((n - p) / (p * (n - 1))) * T2, p, n - p)
or 
library(ICSNP)
HotellingsT2(dta2,mu=c(4,45,10))

##
## Bootstrap test, using Lambda = (|S| / |S_0|)^{n/2} as the test statistic. Under the 
## assumption of multivariate normality, this equals the likelihood ratio statistic. If 
## n were "big", and assuming we were happy assuming multivariate normality, we could 
## compute a p-value via the large-sample chi-square result for likelihood ratio tests. 
## However, we can *always* use the bootstrap, even if n is small and normality does not 
## hold. 
##

set.seed(101)

## Our observed value of Lambda.
S_0_f <- function(X, mu_0) {
  matrix(rowSums(apply(X, 1, function(x) { (x - mu_0) %*% t(x - mu_0) })), nrow = p) / 
    (n - 1)
}

Lambda <- (det(S) / det(S_0_f(X, mu_0))) ^ (n / 2)

B <- 500
Lambda_b <- rep(NA, B)
X_0 <- scale(X, center = TRUE, scale = FALSE) + matrix(rep(mu_0, each = n), nrow = n)
for(b in 1:B) {
  cat(".")

  ## Create bootstrap sample.
  X_b <- X_0[sample(1:n, replace = TRUE), ]

  ## Compute Lambda.
  S_b <- var(X_b)
  S_0_b <- S_0_f(X_b, mu_0)
  Lambda_b[b] <- (det(S_b) / det(S_0_b)) ^ (n / 2)
}

## We again fail to reject H_0.
p_value_boot <- mean(Lambda_b <= Lambda)

以上是关于r 置信区间T2 Bonferroni Hotelling Bootstrap的主要内容,如果未能解决你的问题,请参考以下文章

r语言怎么计算回归模型的置信区间

R与python中的ACF置信区间:为啥它们不同?

置信椭圆与R画法

你如何计算 Python 中 Pearson's r 的置信区间?

R语言编写自定义函数计算R方使用自助法Bootstrapping估计多元回归模型的R方的置信区间可视化获得的boot对象估计单个统计量的置信区间分别使用分位数法和BCa法

R中线性回归模型后的置信区间[关闭]