Question

R Code

Directions: All work has to be your own, you may not work in groups. Show all work. Submit your solutions in a pdf document on Moodle. Include your R code (which must be commented and properly indented) in the pdf file. Name this pdf file ‘your last name’-HW5.pdf. Also submit one text file with your R code, which must be commented and properly indented. You may only use ‘runif’ to generate random numbers; other random number generating functions have to be your own.

You will need to generate from the Gamma distribution in this homework. To do that, you may use the fact that if Z1, . . . , Zm are i.i.d. exponential with mean µ, then the sum of Zi (from i=1 to m) ~ Gamma(m, µ).

1. Let X be the yet-to-be measured diastolic blood pressure of a randomly selected person that takes a particular vitamin. Let Y be the yet-to-be measured diastolic blood pressure for a randomly selected person that does not take this particular vitamin. Suppose that X ~ Gamma(η = 79, ν = 1) and Y ~ Gamma(η = 80, ν = 1).

(a) Use R to graph the densities for X and Y in the same plot. Use dots for the density of Y.

(b) It can be shown that E(X) = 79 and var(X) = 79. Produce a simulation-based 99.5% approximate confidence interval for E((X− 79)2 ) = var(X). Use 500,000 replications. Is 79 in this interval?

(c) It can be shown that E(Y ) = 80 and var(Y ) = 80. Produce a simulation-based 99.5% approximate confidence interval for E((Y − 80)2 ) = var(Y ). Use 500,000 replications. Is 80 in this interval?

(d) Researchers plan to design an experiment to study the association between taking this vitamin and diastolic blood pressure. They plan to use a two-independent-samples model. Let the yet-tobe observed diastolic blood pressures for individuals that take this vitamin be X1, . . . , X_n1 , which are independent copies of X. Let the yet-to-be observed diastolic blood pressures for individuals that do not take this vitamin be Y1, . . . , Y_n2 , which are independent copies of Y. Assume that X1, . . . , X_n1, Y1, . . . , Y_n2 are independent. Use simulation to find the value of n1 = n2 such that the power of the two-independent-samples t-test of H0 : E(X) = E(Y ), HA : E(Y ) does not equal E(X), is roughly 90 % when the 5 % significance level is used.

2. The researchers were unhappy with the sample sizes required in problem 1 (d). They plan to find a set of n identical twins and have one of the two twins take the vitamin and have the other twin not take the vitamin. Let Xi be the yet-to-be observed diastolic blood pressure for the person taking the vitamin in the ith pair of twins. Let Yi be the yet-to-be observed diastolic blood pressure for the person not taking the vitamin in the ith pair of twins. Suppose that Xi = Ti + Ui and Yi = Ti + Vi , where T1, . . . , Tn are i.i.d Gamma(60, 1), U1, . . . , Un are i.i.d. Gamma(21, 1), and V1, . . . , Vn are i.i.d. Gamma(22, 1).

(a) Compare the distribution of X1 to the Gamma(81, 1) using a percentiles plot.

(b) Compare the distribution of Y2 to the Gamma(82, 1) using a percentiles plot.

(c) Construct a simulation based approximate 97 % confidence for the correlation between X1 and Y1.

(d) Use simulation to find the value of n such that the power of the paired t-test of H0 : E(X) = E(Y ), HA : E(Y ) does not equal E(X), is roughly 90 % when the 5 % significance level is used.

3. Suppose that Xi = µ + T + Ui , i = 1, . . . , n, where T ~ N(0, ρ) and U1, . . . , Un are i.i.d. N(0, 1 − ρ). Note that U is the same for all i, so the Xis are correlated. Use 10,000 replications and report the score approximate 98% confidence interval for the coverage probability of Xbar ± (t_0.99,n−1) * (S/√ n) for each (n, ρ) ∈ {20, 40, 80} × {0.3, 0.6, 0.8}. What happens to the coverage probability as ρ increases? What happens to the coverage probability as n increases? Set µ = 11 in the simulations.

Answer 1

1) x_seq = seq(from=20, t0 = 110 length_out = 1e3)

X_den.vals = dgamma(x=x.seq, shape = 79, scale = 1)

plot(x.seq, X_den.vals, type = "1", xlab = "x", ylab = "Density")

Y.den.vals = dgamma(x=x.seq, shape = 80, scale=1)

lines(x.seq,Y.den.vals,lty = 3)

2) set_seed(3701); reps = 5e5

x_list = rgamma(n = reps, shape=79,scale =1)

## CI for var(X) = E[(X-79)*2]

t.test((x.list-79)^2,conf.level=0.995)

$conf.int[1:2]

[1]78.67027 79.57529

the true value of 79 is in this interval

c) set_seed(3701); reps = 5e5

x_list = rgamma(n = reps, shape=80,scale =1)

## CI for var(X) = E[(X-79)*2]

t.test((x.list-80)^2,conf.level=0.995)

$conf.int[1:2]

[1]78.66698 80.58333

the true value of 80 is in this interval

d) Researchers plan to design an experiment to study the association between taking this vitamin and diastolic blood pressure. They plan to use a two-independent-samples model. Let the yet-to-be observed diastolic blood pressures for individuals thattake this vitamin be X1,...,Xn1, which are independent copies of X. Let the yet-to-be observed diastolic blood pressures for individuals that do not take this vitamin beY1,...,Yn2, which are independent copies of Y. Assume that X1n1,Y1n2are independent. Use simulation to find the value of n1=n2 such that the power of thetwo-independent-samples t-test of

H0:E(X)-E(Y) = 0,

Ha:E(X)-E(Y)6= 0,

is roughly 85% when a 1% significance level is used.