0

I want to perform a bootstrap to get a better estimation of the beta estimates, simulated with the simex method.

The basic code is below. Value1 and Value2 should be repeated.

library(simex)
library(data.table) 

value1 <- rnorm(24, mean = 1, sd = 0.3) 
value2 <- rnorm(24,mean = 1.41, sd = 0.5)
group1 <- as.data.frame(value1)
group1$Treat <- 1
setnames(group1, "value1", "Value") 
group2 <-  as.data.frame(value2)
group2$Treat <- 2
setnames(group2, "value2", "Value") 
simulated_data <- rbind(group1, group2)

naive_model <- lm(Treat~Value, data = simulated_data, x = T, y = T)
simex_simulated_data <- simex(naive_model, SIMEXvariable="Value", measurement.error = 0.12, B=100, fitting.method="quadratic")

EDIT:

This is how far I get, not far. It seems thaat the function doesn`t like x=T, however its needed for Simex. 

> getRegrnaive <- function(dat, idx) {   bsFit <- lm(Value~Treat,
> subset=idx, data=dat)   coef(bsFit) }
> 
> getRegr <- function(dat, idx) {   bsFit <- lm(Value~Treat, subset=idx,
> data=dat, x = T, y = T)   simex_boot <- simex(bsFit,
> SIMEXvariable="Value", measurement.error = 0.12, B=100,
> fitting.method="quadratic")   coef(simex_boot) }
> 
> #------------------------ 
>nR <- 999 
>bsRegrsim <- boot(simulated_data, statistic=getRegr, R=nR) 
> #---------------------------------------- 
>bsRegrnaiv <- boot(simulated_data, statistic=getRegrnaive, R=nR)

or something like this?

getRegr <- function(simulated_data, indices) { d <- simulated_data[indices,] bsFit <- lm(Value~Treat, data=d, x = T, y = T) simex_boot <- simex(bsFit, SIMEXvariable="Value", measurement.error = 0.12, B=100, fitting.method="quadratic")
coef(simex_boot) } summary(simex_boot)

however it is not doing what I want... the problem (inducing variability of the outcome) is the initial pupulation which I want to replicate....

so, I did this: 

f <- function (n=8) {
  sderrormanual <- 0.12
  sderrormacro <- 0.02
  sdmeasuredtreatment <- 0.3
  sdmeasuredcontrol <- 0.5
  #n <-8
  errormanual <- rnorm (n*2, mean = 0, sd = sderrormanual) 
  errormacro <- rnorm (n*2, mean = 0, sd = sderrormacro) 
  value1 <- rnorm(n, mean = 1, sd = 0.18) #normalverteilung anhand der system. Analyse von Felix, ausgehend von einem errechneten, "wahren" SD
  value2 <- rnorm(n, mean = 14.1, sd = 0.38)#normalverteilung anhand der system. Analyse von Felix, ausgehend von einem errechneten, "wahren" SD
  group1 <- as.data.frame(value1)
  group1$Treat <- 1
  setnames(group1, "value1", "Value") 
  group2 <-  as.data.frame(value2)
  group2$Treat <- 2
  setnames(group2, "value2", "Value") 
  simulated_data <- rbind(group1, group2)
  simulated_data$errormanual <- errormanual 
  simulated_data$errormacro <- errormacro 
  simulated_data$Valueerrormanual <- simulated_data$Value + simulated_data$errormanual
  simulated_data$Valueerrormacro <- simulated_data$Value + simulated_data$errormacro
  d <- simulated_data
  bsFit <- lm(Valueerrormanual~Treat, data=d)
}

But than

out <- replicate(10, f(), simplify = "array")

results in some crazy estimates, which have nothing todo with the original result from

bsFit <- lm(Valueerrormanual~Treat, data=d)

I don`t understand

Err... it was a typo

rbru
  • 11
  • 4

1 Answers1

0
getRegr <- function(simulated_data, indices) {
  d <- simulated_data[indices,]
  bsFit <- lm(Value~Treat, data=d, x = T, y = T)
  simex_boot <- simex(bsFit, SIMEXvariable="Value", measurement.error = 0.12, B=100, fitting.method="quadratic")
  coef(simex_boot)
}
summary(simex_boot)
rbru
  • 11
  • 4