## —- fig.align=”center”, echo=FALSE, fig.width = 8, fig.height = 5, out.width = “1\\textwidth”—-
# When fixing rate (lambda) and changing shape (r) for Gamma Distribution,
# When the shape (r) increases, based on the formula,
# the mean increases (shift to the right),
# the variance increases
# the skewness decreases 偏度
# the excess kurtosis decreases 峰度
### F Distribution 画图，画df1相同，df2不同，及df1不同，df2相同的图
# df1 = q numerator 分子，df2 = n-K denominator 分母
# Plot 1: Fix df2 and changing df1
par(mfrow=c(1,2)) # par(mfrow=c(nrows, ncols))确定组图位置
curve(expr = df(x = x, df1 = 3, df2 = 5),
xlab = “”, ylab = “”, main = “”,
lwd = 2, col = 1, xlim = c(0, 4),
ylim = c(0, 1))
for (i in 1:2) {
curve(expr = df(x = x, df1 = 3,
df2=c(15,500)[i]),
lwd = 2, col = (2:3)[i], add = TRUE)
}
legend(x = “topright”, legend = c(“F(3,5)”, “F(3,15)”, “F(3,500)”),
lwd = 2, col = 1:3)
curve(expr = df(x = x, df1 = 1, df2 = 30),
xlab = “”, ylab = “”, main = “”,
lwd = 2, col = 1, xlim = c(0, 4),
ylim = c(0, 1))
for (i in 1:2) {
curve(expr = df(x = x, df1 = c(3,15)[i],
df2=30),
lwd = 2, col = (4:6)[i], add = TRUE)
}
legend(x = “topright”, legend = c(“F(1,30)”, “F(3,30)”, “F(15,30)”),
lwd = 2, col = 4:6)
## —- fig.align=”center”, echo=FALSE, fig.width = 8, fig.height = 5, out.width = “1\\textwidth”—-
## as the degrees of freedom increase, the shape of the t distribution comes closer to that of a standard normal bell curve.
## In case of small degrees of freedom values, we find the distribution to have heavier tails than a standard normal.
# plot the standard normal density
curve(dnorm(x),
xlim = c(-4, 4),
xlab = “”,
lty = 2,
ylab = “”,
main = “”)
# plot the t density for M=2
curve(dt(x, df = 2),
xlim = c(-4, 4),
col = 2,
add = T)
# plot the t density for M=4
curve(dt(x, df = 4),
xlim = c(-4, 4),
col = 3,
add = T)
# plot the t density for M=25 # 25 degrees of freedom
curve(dt(x, df = 25),
xlim = c(-4, 4),
col = 4,
add = T)
# add a legend
legend(“topright”, bty=”n”,
c(“N(0, 1)”, expression(t[2]), expression(t[2]), expression(t[25])),
col = 1:4,
lty = c(2, 1, 1, 1))
## —- fig.align=”center”, echo=FALSE, fig.width = 8, fig.height = 5, out.width = “1\\textwidth”—-
## critical value for F-test
par(mfrow=c(1,1)) # par(mfrow=c(nrows, ncols))确定组图位置
library(“scales”) # scales包用于可视化
curve(expr = df(x = x, df1 = 9, df2 = 120),
xlab = “”, ylab = “”, main = “”,
lwd = 2, col = 1, xlim = c(0, 4),
ylim = c(0, 1), xaxs=”i”,yaxs=”i”)
alpha <- 0.05 # confidence level
q <- qf(p = 1-alpha, df1 = 9, df2 = 120) # qf gives the quantile function, 即F表中对应值c_alpha
xx <- seq(0,q,len=25)
yy <- df(x = xx, df1 = 9, df2 = 120) # df gives the density, pf gives the distribution function
polygon(x = c(xx,rev(xx)), y = c(yy, rep(0,length(yy))), border = NA, col = alpha("green", 0.25))
# polygon根据顶点的x和y坐标绘制多边形
# rev不改变行/列，仅倒置向量中的元素顺序
q <- qf(p = 1-alpha, df1 = 9, df2 = 120)
xx <- seq(q,4,len=25)
yy <- df(x = xx, df1 = 9, df2 = 120)
polygon(x = c(xx,rev(xx)), y=c(yy, rep(0,length(yy))), border = NA, col = alpha("red", 0.25))
lines(x=c(0,q-0.02),y=c(0,0), col="darkgreen", lwd=10)
lines(x=c(q+0.02,4),y=c(0,0), col="red", lwd=10)
legend("topright", pch=c(22,NA, 22, NA), lty=c(NA,1,NA,1), lwd=c(NA,4,NA,4), cex=1,
col = c(alpha("green", 0.25),"darkgreen",alpha("red", 0.25),"red"),
legend=c(expression(1-alpha==~"95% of"~F['9,120']),"Non-Rejection Region",
expression(alpha==~"5% of"~F['9,120']),"Rejection Region"),
bty="n", pt.bg = c(alpha("green", 0.25),alpha("green", 0.25),alpha("red", 0.25),alpha("red", 0.25)))
curve(expr = df(x = x, df1 = 9, df2 = 120), lwd = 2, col = 1, add=TRUE, from = 0, to = 4)
lines(x=c(q,q), y=c(0,.6),lwd=2,lty=2)
text(x = q, y = .65, labels = expression(c[alpha]==1.9588))
## ---------------------------
## Reading the F-Table. 读F表
df1 <- 9 # numerator df 分子
df2 <- 120 # denominator df 分母
alpha <- 0.05 # significance level
## Critical value c_alpha (= (1-alpha) quantile):
c_alpha <- qf(p = 1-alpha, df1 = df1, df2 = df2)
c_alpha
## ---------------------------
## Changing the significance level from α = 0.05 to α = 0.01 makes the critical value cα larger
## and, therefore, the rejection region smaller (fewer Type I errors)
alpha <- 0.01
## Critical value c_\alpha = 1-\alpha-quantile:
c_alpha <- qf(p = 1-alpha, df1 = df1, df2 = df2)
c_alpha
## ---- fig.align="center", echo=FALSE, fig.width = 8, fig.height = 5, out.width = "1\\textwidth"----
## t distribution 画图, two-sided t-test, 双侧检验
library("scales")
curve(expr = dt(x = x, df = 12),
xlab = "", ylab = "", main = "",
lwd = 2, col = 1, xlim = c(-5, 5),
ylim = c(0, .6), xaxs="i",yaxs="i")
alpha <- 0.05/2
q <- qt(p = 1-alpha, df=12) # 读表，t值
xx1 <- seq(-5,-q,len=25)
yy1 <- dt(x = xx1, df = 12)
xx2 <- seq(q,5,len=25)
yy2 <- dt(x = xx2, df = 12)
xx3 <- seq(-q,q,len=25)
yy3 <- dt(x = xx3, df = 12)
polygon(x = c(xx1,rev(xx1)), y=c(yy1, rep(0,length(yy1))), border = NA, col = alpha("red", 0.25))
polygon(x = c(xx2,rev(xx2)), y=c(yy2, rep(0,length(yy2))), border = NA, col = alpha("red", 0.25))
polygon(x = c(xx3,rev(xx3)), y=c(yy3, rep(0,length(yy3))), border = NA, col = alpha("green", 0.25))
legend("topright", pch=c(22,22), lty=c(NA,NA), lwd=c(NA,NA), cex=1,
col = c(alpha("green", 0.25),alpha("red", 0.25)),
legend=c(expression("95% of"~t['12']),
expression("5% of"~t['12'])),
bty="n", pt.bg = c(alpha("green", 0.25),alpha("red", 0.25)))
curve(expr = dt(x = x, df= 12), lwd = 2, col = 1, add=TRUE)
lines(x=c(q,q), y=c(0,.35),lwd=2,lty=2)
lines(x=c(-q,-q), y=c(0,.35),lwd=2,lty=2)
text(x = q, y = .45, labels = expression(c[alpha/2]==2.18))
text(x = -q, y = .45, labels = expression(-c[alpha/2]==-2.18)) # 方括号内为脚标
## ---- fig.align="center", echo=FALSE, fig.width = 8, fig.height = 5, out.width = "1\\textwidth"----
## t distribution 画图, one-sided t-test, 单侧检验, 左侧
library("scales")
alpha <- 0.05
q <- qt(p = 1-alpha, df=12)
xx1 <- seq(-5,-q,len=25)
yy1 <- dt(x = xx1, df = 12)
xx2 <- seq(-q,5,len=25)
yy2 <- dt(x = xx2, df = 12)
##
xx3 <- seq(q,5,len=25)
yy3 <- dt(x = xx3, df = 12)
xx4 <- seq(-5,q,len=25)
yy4 <- dt(x = xx4, df = 12)
par(mfrow=c(1,2))
curve(expr = dt(x = x, df = 12),
xlab = "", ylab = "", main = "",
lwd = 2, col = 1, xlim = c(-5, 5),
ylim = c(0, .6), xaxs="i",yaxs="i")
polygon(x = c(xx1,rev(xx1)), y=c(yy1, rep(0,length(yy1))), border = NA, col = alpha("red", 0.25))
polygon(x = c(xx2,rev(xx2)), y=c(yy2, rep(0,length(yy2))), border = NA, col = alpha("green", 0.25))
lines(x=c(-q,-q), y=c(0,.35),lwd=2,lty=2)
text(x = -q, y = .45, labels = expression(-c[alpha]==-1.78))
legend("topright", pch=c(22,22), lty=c(NA,NA), lwd=c(NA,NA), cex=1,
col = c(alpha("green", 0.25),alpha("red", 0.25)),
legend=c(expression("95% of"~t['12']),
expression("5% of"~t['12'])),
bty="n", pt.bg = c(alpha("green", 0.25),alpha("red", 0.25)))
###########
## t distribution 画图, one-sided t-test, 单侧检验, 右侧
curve(expr = dt(x = x, df = 12),
xlab = "", ylab = "", main = "",
lwd = 2, col = 1, xlim = c(-5, 5),
ylim = c(0, .6), xaxs="i",yaxs="i")
polygon(x = c(xx3,rev(xx3)), y=c(yy3, rep(0,length(yy3))), border = NA, col = alpha("red", 0.25))
polygon(x = c(xx4,rev(xx4)), y=c(yy4, rep(0,length(yy4))), border = NA, col = alpha("green", 0.25))
lines(x=c(q,q), y=c(0,.35),lwd=2,lty=2)
text(x = q, y = .45, labels = expression(c[alpha]==1.78))
legend("topright", pch=c(22,22), lty=c(NA,NA), lwd=c(NA,NA), cex=1,
col = c(alpha("green", 0.25),alpha("red", 0.25)),
legend=c(expression("95% of"~t['12']),
expression("5% of"~t['12'])),
bty="n", pt.bg = c(alpha("green", 0.25),alpha("red", 0.25)))
par(mfrow=c(1,1))
## ---------------------------
## Reading the t-Table. 读t表
df <- 16 # degrees of freedom
alpha <- 0.05 # significance level
## One-sided critical value (= (1-alpha) quantile):
c_oneSided <- qt(p = 1-alpha, df = df)
c_oneSided
## Two-sided critical value (= (1-alpha/2) quantile):
c_twoSided <- qt(p = 1-alpha/2, df = df)
## lower critical value
-c_twoSided
## upper critical value
c_twoSided
## ---- fig.align="center", echo=FALSE, fig.width = 8, fig.height = 5, out.width = "1\\textwidth"----
library("scales") # transperent color
mean.alt <- 3
x <- seq(-4, 4, length=1000)
hx <- dt(x, df=12)
alpha <- 0.05
plot(x, hx, type="n", xlim=c(-6, 8), ylim=c(0, 0.65),
ylab = "", xlab = "", axes=T)
#axis(1)
xfit2 <- x + mean.alt
yfit2 <- dt(x, df=12)
## Print null hypothesis area
polygon(c(min(x), x, max(x)),
c(0, hx, 0),
col =alpha("grey", 0.5),
border=alpha("grey", 0.9))
## The alternative hypothesis area
## The red area: Type II error
lb <- min(xfit2)
ub <- round(qt(1-alpha, df=12),2)
i <- xfit2 >= lb & xfit2 <= ub
polygon(c(lb,xfit2[i],ub), c(0,yfit2[i],0), col=alpha("red", 0.25), border=alpha("red", 0.25))
## The green area: Power
i <- xfit2 >= ub
polygon(c(ub,xfit2[i],max(xfit2)), c(0,yfit2[i],0),
col=alpha(“green”, 0.25),
border=alpha(“green”, 0.25))
text(x=0+.25,y=.425, expression(t[12]), pos=2)
text(x=3-3/5,y=.425, expression(t[12]~”+”~3), pos=4)
legend(x=-6.5,y=.65, title=NULL, bty=”n”,
c(“Null Distribution”,”Type II Error”, expression(paste(“Power”))),
fill=c(alpha(“grey”, 0.5), alpha(“red”, 0.25), alpha(“green”, 0.5)), horiz=FALSE)
## —————————
## generate data from the following fully specified data generating process:
## Yi = beta_1 + beta_2*X_i2 + beta_3*X_i3 + eps_i
## beta = (beta_1, beta_2, beta_3)’ = (2,3,4)’
## X_i2 ～ U[2,10], X_i3 ～ U[12,22], eps_i ～ N(0,3^2), n = 10
## myDataGenerator()
## Function to generate artificial data
## If X=NULL: new X variables are generated,
## If the user gives X variables, the sampling of new Y variables is conditionally on the given X variables.
myDataGenerator <- function(n, beta, X=NULL, sigma=3){
if(is.null(X)){ # is.null returns TRUE if its argument's value is NULL and FALSE otherwise.
X <- cbind(rep(1,n),
runif(n,2,10),
runif(n,12,22))
}
eps <- rnorm(n, sd=sigma)
Y <- X %*% beta + eps
data <- data.frame("Y"=Y, "X_1"=X[,1], "X_2"=X[,2], "X_3"=X[,3])
return(data)
}
## Define a true beta vector
beta_true <- c(2,3,4)
## Check:
## Generate Y and X data
test_data <- myDataGenerator(n = 10, beta=beta_true)
## Generate new Y data conditionally on X
X_cond <- cbind(test_data$X_1, test_data$X_2, test_data$X_3)
test_data_new <- myDataGenerator(n = 10, beta = beta_true, X = X_cond)
## compare
round(head(test_data,3),2) # new Y, new X # round将x舍入为指定位的小数, head列出某个对象的开始部分
#> Y X_1 X_2 X_3
#> 1 93.24 1 2.42 20.19
#> 2 101.12 1 4.85 20.43
#> 3 83.82 1 7.20 15.00
round(head(test_data_new, 3), 2) # New Y, conditionally on X, head列出某个对象的开始部分
#> Y X_1 X_2 X_3
#> 1 83.74 1 2.42 20.19
#> 2 95.42 1 4.85 20.43
#> 3 81.67 1 7.20 15.00
## —- fig.align=”center”, echo=TRUE, fig.width = 8, fig.height = 5, out.width = “1\\textwidth”—-
## estimators beta_hat_k|X should be normal distributed conditionally on X
## 如果不conditionally on x, beta的分布和他的理论分布会差很多
## check the distribution using a Monte Carlo simulation for the case of beta_hat_2|X with a small sample size of n = 10.
set.seed(123)
n <- 10 # a small sample size
beta_true <- c(2,3,4) # true data vector
sigma <- 3 # true variance of the error term
## Let's generate a data set from our data generating process
mydata <- myDataGenerator(n = n, beta=beta_true)
X_cond <- cbind(mydata$X_1, mydata$X_2, mydata$X_3)
## True mean and variance of the true normal distribution of beta_hat_2|X:
# true mean
beta_true_2 <- beta_true[2]
# true variance
var_true_beta_2 <- sigma^2 * diag(solve(t(X_cond) %*% X_cond))[2]
## Let's generate 5000 realizations from beta_hat_2 conditionally on X and check whether their
## distribution is close to the true normal distribution
rep <- 5000 # MC replications
beta_hat_2 <- rep(NA, times=rep)
##
for(r in 1:rep){
MC_data <- myDataGenerator(n = n,
beta = beta_true,
X = X_cond) # 这里控制了X不变，conditionally on X
lm_obj <- lm(Y ~ X_2 + X_3, data = MC_data)
beta_hat_2[r] <- coef(lm_obj)[2]
}
## Compare
## True beta_2 versus average of beta_hat_2 estimates
beta_true_2
#> [1] 3
round(mean(beta_hat_2), 4)
#> [1] 3.0091
## True variance of beta_hat_2 versus empirical variance of beta_hat_2 estimates
round(var_true_beta_2, 4)
#> [1] 0.416
round(var(beta_hat_2), 4)
#> [1] 0.4235
## True normal distribution of beta_hat_2 versus empirical density of beta_hat_2 estimates
library(“scales”)
curve(expr = dnorm(x, mean = beta_true_2,
sd=sqrt(var_true_beta_2)),
xlab=””,ylab=””, col=gray(.2), lwd=3, lty=1,
xlim=range(beta_hat_2), ylim=c(0,1.1))
lines(density(beta_hat_2, bw = bw.SJ(beta_hat_2)),
col=alpha(“blue”,.5), lwd=3)
legend(“topleft”, lty=c(1,1), lwd=c(3,3),
col=c(gray(.2), alpha(“blue”,.5)), bty=”n”, legend=
c(expression(
“Theoretical Gaussian Density of”~hat(beta)[2]~’|’~X),
expression(
“Empirical Density Estimation based on MC realizations from”~
hat(beta)[2]~’|’~X)))
## —- fig.align=”center”, echo=TRUE, fig.width = 8, fig.height = 5, out.width = “1\\textwidth”—-
## However, what would happen if we would sample unconditionally on X?
## How does the distribution of beta_hat_2 would then look like?
set.seed(123)
## Let’s generate 5000 realizations from beta_hat_2
## WITHOUT conditioning on X
beta_hat_2_uncond <- rep(NA, times=rep)
##
for(r in 1:rep){
MC_data <- myDataGenerator(n = n,
beta = beta_true) # 这里没有X = X_cond，没有控制X，总是产生new Y和new X
lm_obj <- lm(Y ~ X_2 + X_3, data = MC_data)
beta_hat_2_uncond[r] <- coef(lm_obj)[2]
}
## Compare
## True beta_2 versus average of beta_hat_2 estimates
beta_true_2
#> [1] 3
round(mean(beta_hat_2_uncond), 4)
#> [1] 2.9973
## True variance of beta_hat_2 versus
## empirical variance of beta_hat_2 estimates
round(var_true_beta_2, 4)
#> [1] 0.416
round(var(beta_hat_2_uncond), 4)
#> [1] 0.2521
## Plot
curve(expr = dnorm(x, mean = beta_true_2,
sd=sqrt(var_true_beta_2)),
xlab=””,ylab=””, col=gray(.2), lwd=3, lty=1,
xlim=range(beta_hat_2_uncond), ylim=c(0,1.1))
lines(density(beta_hat_2_uncond, bw=bw.SJ(beta_hat_2_uncond)),
col=alpha(“blue”,.5), lwd=3)
legend(“topleft”, lty=c(1,1), lwd=c(3,3),
col=c(gray(.2), alpha(“blue”,.5)), bty=”n”, legend=
c(expression(
“Theoretical Gaussian Density of”~hat(beta)[2]~’|’~X),
expression(
“Empirical Density Estimation based on MC realizations from”~
hat(beta)[2])))
## Not good.
## Since we do not condition on X, the realizations of X affect the distribution
## of beta_hat and our theoretical Gaussian distribution result in Equation
## (4.1) does not apply anymore.
## —————————
## Testing Multiple Parameters
## H0两种形式: beta_2 = 3, beta_3 = 4 or H0: R*beta − r = 0
suppressMessages(library(“car”)) # for linearHyothesis()
# ?linearHypothesis
## Estimate the linear regression model parameters
lm_obj <- lm(Y ~ X_2 + X_3, data = MC_data)
## Option 1: H0是beta_2 = 3, beta_3 = 4的形式
car::linearHypothesis(model = lm_obj,
hypothesis.matrix = c("X_2=3", "X_3=4"))
#> Linear hypothesis test
#> Hypothesis:
#> X_2 = 3
#> X_3 = 4
#> Model 1: restricted model
#> Model 2: Y ~ X_2 + X_3
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 9 35.371
#> 2 7 27.357 2 8.0147 1.0254 0.4069
## Option 2: H0是R*beta − r = 0的形式
# R是beta前面系数，r是beta等于多少，q是受限条件的个数，k是beta的个数
R <- rbind(c(0,1,0),
c(0,0,1)) # c() create rows
car::linearHypothesis(model = lm_obj,
hypothesis.matrix = R, # q=2，k=3，所以R有两行三列
rhs = c(3,4)) # r, right-hand-side vector for hypothesis
#> Linear hypothesis test
#> Hypothesis:
#> X_2 = 3
#> X_3 = 4
#> Model 1: restricted model
#> Model 2: Y ~ X_2 + X_3
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 9 35.371
#> 2 7 27.357 2 8.0147 1.0254 0.4069
## F: F statistic, 1.0254, one-sided test, larger than 1 is large, 1.0254接近1，还行
## Pr(>F): p value, 大于等于F statistic的可能性，40.69%, 相当高了，不能拒绝H0
## —————————
## in repeated samples we should nevertheless observe alpha*100%
## type I errors (false rejections of H0)
## Let’s generate 5000 F-test decisions and check
## whether the empirical rate of type I errors is
## close to the theoretical significance level.
# 检查p值
rep <- 5000 # MC replications
F_test_pvalues <- rep(NA, times=rep)
##
for(r in 1:rep){
## generate new MC_data conditionally on X_cond
MC_data <- myDataGenerator(n = n,
beta = beta_true,
X = X_cond)
lm_obj <- lm(Y ~ X_2 + X_3, data = MC_data)
## save the p-value
p <- linearHypothesis(lm_obj,
c("X_2=3", "X_3=4"))$`Pr(>F)`[2] # select 2. element of p value vector
F_test_pvalues[r] <- p # save p value in the r th ran, so total 5000 p values
}
##
signif_level <- 0.05
rejections <- F_test_pvalues[F_test_pvalues < signif_level] # subselect
# 上面是all the cases of rejections, so much shorter than the rep
round(length(rejections)/rep, 3) # length(rejections)计算how many cases of rejections were there
# round保留三位小数
#> [1] 0.05
##
signif_level <- 0.01
rejections <- F_test_pvalues[F_test_pvalues < signif_level]
round(length(rejections)/rep, 3)
#> [1] 0.009
# type I error rate is approximately equal to the chosen significance level α
# 蒙特卡罗模拟次数越多，越接近选定的显著性水平，因为大数定律
# Script page 116
## —————————
## how well the F test detects certain violations of the null hypothesis.
## testing the following incorrect null hypothesis using Monte Carlo
## H0: beta_2 = 4 and beta_3 = 4
## HA: beta_2 ≠ 4 and/or beta_3 ≠ 4
set.seed(321)
rep <- 5000 # MC replications
F_test_pvalues <- rep(NA, times=rep)
##
for(r in 1:rep){
## generate new MC_data conditionally on X_cond
MC_data <- myDataGenerator(n = n,
beta = beta_true,
X = X_cond)
lm_obj <- lm(Y ~ X_2 + X_3, data = MC_data)
## save p-values of all rep-many tests
F_test_pvalues[r] <- linearHypothesis(lm_obj,
c("X_2=4", "X_3=4"))$`Pr(>F)`[2]
}
##
signif_level <- 0.05
rejections <- F_test_pvalues[F_test_pvalues < signif_level]
length(rejections)/rep
#> [1] 0.3924
# we can now reject the (false) null hypothesis in approximately 39% of
# all resamplings from the true data generating process.
## we can never use an insignificant test result (p-value≥ α)
## as a justification to accept the null hypothesis.
##
signif_level <- 0.10
rejections <- F_test_pvalues[F_test_pvalues < signif_level]
length(rejections)/rep
#> 0.5528
## —————————
car::linearHypothesis(lm_obj, c(“X_2=4”, “X_3=4”))
#> Linear hypothesis test
#> Hypothesis:
#> X_2 = 4
#> X_3 = 4
#> Model 1: restricted model
#> Model 2: Y ~ X_2 + X_3
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 9 141.245
#> 2 7 32.572 2 108.67 11.677 0.005889 **
#> —
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 看Pr(>F)，比significant level大，已经可以拒绝H0
## 但F-test不能判断H0的哪个部分不对，我们只知道至少有一个关于beta的假设不对
## —————————
## Computing confidence intervals
signif_level = 0.05
## 95% CI for beta_2
confint(lm_obj, parm = “X_2”, level = 1 – signif_level)
#> 2.5 % 97.5 %
#> X_2 1.370315 3.563536
## 95% CI for beta_3
confint(lm_obj, parm = “X_3″, level = 1 – signif_level)
#> 2.5 % 97.5 %
#> X_3 3.195389 4.695134
## —- fig.align=”center”, echo=TRUE, fig.width = 8, fig.height = 5, out.width = “1\\textwidth”—-
## using these two-sided confidence intervals to do hypothesis tests.
## check whether the confidence interval CI1−α contains the hypothetical value 4 or not
## In case of 4 ∈ CI1−α, we cannot reject the null hypothesis.
## In case of 4 ̸∈ CI1−α, we reject the null hypothesis.
## under the null hypothesis, CI1−α falsely rejects the null hypothesis in only α · 100% of resamplings.
## H0: beta_3 = 4 versus HA: beta_3 ≠ 4
## using these two-sided confidence intervals to do hypothesis tests
## Let’s generate 1000 CIs
set.seed(123)
signif_level <- 0.05
rep <- 5000 # MC replications
confint_m <- matrix(NA, nrow=2, ncol=rep)
##
for(r in 1:rep){
## generate new MC_data conditionally on X_cond
MC_data <- myDataGenerator(n = n,
beta = beta_true,
X = X_cond)
lm_obj <- lm(Y ~ X_2 + X_3, data = MC_data)
## save the p-value
CI <- confint(lm_obj, parm = "X_2", level = 0.95)
# confidence intervals, parm指示哪一个变量的置信区间
confint_m[,r] <- CI # save each confidence interval to each col r. r th col contains r th CI
}
##
inside_CI <- confint_m[1,] <= beta_true_2 & # 1. row, all col. # 判断5000个confint_m[1,] 是否小于等于true value，得到5000个true or false
beta_true_2 <= confint_m[2,] # 判断true value是否小于等于5000个confint_m[2,]，即模拟值是否cover true value
# inside is true, outside is false
## CI-lower, CI-upper, beta_true_2 inside?
head(cbind(t(confint_m), inside_CI)) # head列出某个对象的开始部分
#> inside_CI
#> [1,] 0.8555396 3.639738 1
#> [2,] 0.9143542 3.270731 1
#> [3,] 1.9336526 4.984167 1
#> [4,] 1.9985874 3.812695 1
#> [5,] 3.0108642 5.621791 0
#> [6,] 2.0967675 4.716398 1
round(length(inside_CI[inside_CI == FALSE])/rep, 2) # count how many 0=FALSE above
#> [1] 0.05
##
nCIs <- 100
plot(x=0,y=0,type="n",xlim=c(0,nCIs),ylim=range(confint_m[,1:nCIs]),
ylab="", xlab="Resamplings", main="Confidence Intervals")
for(r in 1:nCIs){
if(inside_CI[r]==TRUE){
lines(x=c(r,r), y=c(confint_m[1,r], confint_m[2,r]),
lwd=2, col=gray(.5,.5))
}else{
lines(x=c(r,r), y=c(confint_m[1,r], confint_m[2,r]),
lwd=2, col="darkred")
}
}
axis(4, at=beta_true_2, labels = expression(beta[2]))
abline(h=beta_true_2)