[SOLVED] CS ## – fig.align=center, echo=FALSE, fig.width = 8, fig.height = 5, out.width = 1\textwidth-

## – 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 df1df2df1df2
# 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, Fc_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)) # polygonxy# 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. Fdf1 <- 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) # txx1<- 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. tdf<- 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 # roundx, 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) # Xconditionally 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)
## Lets 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_condXnew Ynew 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: H0beta_2 = 3, beta_3 = 4car::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: H0R*beta r = 0
# Rbetarbetaqkbeta
R <- rbind(c(0,1,0), c(0,0,1)) # c() create rowscar::linearHypothesis(model = lm_obj, hypothesis.matrix = R, # q=2k=3Rrhs = 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.02541
## Pr(>F): p value, F statistic40.69%, H0

##
## in repeated samples we should nevertheless observe alpha*100%
## type I errors (false rejections of H0)
## Lets 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 levelH0
## F-testH0beta

##
## 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
## Lets 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, parmconfint_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. # 5000confint_m[1,] true value5000true or falsebeta_true_2 <= confint_m[2,] # true value5000confint_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)