quizAnwers

peter — Mar 16, 2013, 12:30 AM

5/25

[1] 0.2

5/55

[1] 0.09091

####

set.seed(3343)
pValues = rep(NA,100)
for(i in 1:100){
  z = rnorm(20)
  x = rnorm(20)
  y = rnorm(20,mean=0.5*x)
  pValues[i] = summary(lm(y ~ x))$coef[2,4]
}
# Controls FWER
sum(p.adjust(pValues,method="bonferroni") < 0.1)

[1] 7

# Controls FDR
sum(p.adjust(pValues,method="BH") < 0.1)

[1] 61

####

# generate the b_i

####

b0 = 1
b1 = 2

# case 1
set.seed(17)
beta1 = rep(NA,1000)
for (j in 1:1000){
  e <- rnorm(100)
  x <- rnorm(100)
  x
  threshold = quantile(x,.80)
  for (i in 1:20) {
    if (x[i] > threshold) x[i] = NA
  }
  x
  y = b0 + b1*x + e

  beta1[j] = lm(y ~ x)$coefficients[2]
}
mean(beta1)

[1] 1.996

# case 2
beta1 = rep(NA,1000)
for (j in 1:1000){
  e <- rnorm(100)
  x <- rnorm(100)
  y = b0 + b1*x + e
  threshold = quantile(y,.80)
  for (i in 1:20) {
    if (y[i] > threshold) y[i] = NA
  }

  beta1[j] = lm(y ~ x)$coefficients[2]
}
mean(beta1)

[1] 1.978

####

library(MASS)

b0 = 1
b1 = 2

# case 1
set.seed(17)
beta1 = rep(NA,1000)
for (j in 1:1000){
e <- rnorm(100)
x <- rnorm(100)
x
threshold = quantile(x,.80)
for (i in 1:20) {
  if (x[i] > threshold) x[i] = NA
}
x
y = b0 + b1*x + e

beta1[j] = rlm(y ~ x,na.action=na.omit)$coefficients[2]
}
mean(beta1)

[1] 1.995

# case 2
beta1 = rep(NA,1000)
for (j in 1:1000){
e <- rnorm(100)
x <- rnorm(100)
y = b0 + b1*x + e
threshold = quantile(y,.80)
for (i in 1:20) {
  if (y[i] > threshold) y[i] = NA
}

beta1[j] = rlm(y ~ x)$coefficients[2]
}
mean(beta1)

[1] 1.977