## —- 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 levelq <- qf(p = 1-alpha, df1 = 9, df2 = 120) # qf gives the quantile function, 即F表中对应值c_alphaxx<- seq(0,q,len=25)yy<- df(x = xx, df1 = 9, df2 = 120) # df gives the density, pf gives the distribution functionpolygon(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/2q<- 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.05q<- 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 valuec_twoSided## —- fig.align=”center”, echo=FALSE, fig.width = 8, fig.height = 5, out.width = “1\textwidth”—-library(“scales”) # transperent colormean.alt <- 3x<- seq(-4, 4, length=1000)hx <- dt(x, df=12)alpha <- 0.05plot(x, hx, type=”n”, xlim=c(-6, 8), ylim=c(0, 0.65),ylab = “”, xlab = “”, axes=T)#axis(1)xfit2 <- x + mean.altyfit2 <- dt(x, df=12)## Print null hypothesis areapolygon(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 errorlb <- min(xfit2)ub <- round(qt(1-alpha, df=12),2)i <- xfit2 >= lb & xfit2 <= ubpolygon(c(lb,xfit2[i],ub), c(0,yfit2[i],0), col=alpha(“red”, 0.25), border=alpha(“red”, 0.25))## The green area: Poweri <- 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 + epsdata <- data.frame(“Y”=Y, “X_1″=X[,1], “X_2″=X[,2], “X_3″=X[,3])return(data)}## Define a true beta vectorbeta_true <- c(2,3,4)## Check:## Generate Y and X datatest_data <- myDataGenerator(n = 10, beta=beta_true) ## Generate new Y data conditionally on XX_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)## compareround(head(test_data,3),2) # new Y, new X # round将x舍入为指定位的小数, head列出某个对象的开始部分#>Y X_1X_2 X_3
#> 193.24 1 2.42 20.19
#> 2 101.12 1 4.85 20.43
#> 383.82 1 7.20 15.00
round(head(test_data_new, 3), 2) # New Y, conditionally on X, head列出某个对象的开始部分
#>Y X_1X_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 sizebeta_true <- c(2,3,4) # true data vectorsigma <- 3# true variance of the error term## Let’s generate a data set from our data generating processmydata<- 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 meanbeta_true_2 <- beta_true[2] # true variancevar_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 distributionrep<- 5000 # MC replicationsbeta_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 Xlm_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 estimatesbeta_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 Xlm_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 estimatesbeta_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.DfRSS Df Sum of SqF Pr(>F)
#> 19 35.371
#> 27 27.35728.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 rowscar::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.DfRSS Df Sum of SqF Pr(>F)
#> 19 35.371
#> 27 27.35728.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 replicationsF_test_pvalues<- rep(NA, times=rep)##for(r in 1:rep){## generate new MC_data conditionally on X_condMC_data <- myDataGenerator(n= n,beta = beta_true,X= X_cond)lm_obj<- lm(Y ~ X_2 + X_3, data = MC_data)## save the p-valuep <- 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.05rejections <- F_test_pvalues[F_test_pvalues < signif_level] # subselect # 上面是all the cases of rejections, so much shorter than the repround(length(rejections)/rep, 3) # length(rejections)计算how many cases of rejections were there# round保留三位小数#> [1] 0.05
##
signif_level <-0.01rejections <- 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 replicationsF_test_pvalues<- rep(NA, times=rep)##for(r in 1:rep){## generate new MC_data conditionally on X_condMC_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 testsF_test_pvalues[r] <- linearHypothesis(lm_obj,c(“X_2=4”, “X_3=4”))$`Pr(>F)`[2]
}
##
signif_level <-0.05rejections <- 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.10rejections <- 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 SqF Pr(>F)
#> 19 141.245
#> 2732.5722108.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.05rep<- 5000 # MC replicationsconfint_m<- matrix(NA, nrow=2, ncol=rep)##for(r in 1:rep){## generate new MC_data conditionally on X_condMC_data <- myDataGenerator(n= n,beta = beta_true,X= X_cond)lm_obj<- lm(Y ~ X_2 + X_3, data = MC_data)## save the p-valueCI <- 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 falsebeta_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 <- 100plot(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)
Reviews
There are no reviews yet.