2017-05-21 14 views
0

Ich muss eine Simulationsstudie in R codieren. Ich habe X1, ..., X15 ~ N (0,1) erklärende Variablen und Y ~ N (2 + 2 * X1 + 0.8 * X2-1.2 * X15, 1) und ich muss n = 100 Werte simulieren und wiederhole das für iter = 100 mal. Dann muss ich für jedes erstellte lineare Modell das beste Submodell mit stepAIC finden. Ich schrieb den folgenden Code:Wie führe ich eine Simulationsstudie mit stepAIC durch

set.seed(1234) 
sim <- function (sd) { 
n <- 100 
p <- 15 
X <- matrix(rnorm(n * p), n, p) 
mu <- 2 + 2*X[,1] + 0.8*X[,2] - 1.2*X[,15] 
Y <- matrix(rnorm(100, mu,sd)) 
sim<-data.frame(Y,X) 
r<- lm(Y~., data = sim) 
library(MASS) 
r0<-lm(Y~1, data=sim) 
res<-stepAIC(r0,k=2,direction="forward", scope=list(lower=~1, upper=r)) 
return(res$coefficients) 
} 

sim(1) 
oo1<- lapply(1:100, sim) 

Wie ich einen unerfahrenen R-User bin, denke ich, dass ich etwas falsch mache. Der Zweck der Studie besteht darin, herauszufinden, ob innerhalb der 100 besten Teilmodelle (gemäß stepAIC) Modelle existieren, die den echten finden können (Y = 2 + 2 * X1 + 0.8 * X2-1.2 * X15 + e) . Falls ich die falschen Dinge mache, könnte ich Hilfe/Hinweise bekommen, um sie richtig zu implementieren?

Antwort

0

Hier ist eine Arbeitsversion des Codes:

library(MASS) 
set.seed(1234) 

sim <- function(sd, n, p) { 
    X <- matrix(rnorm(n * p), n, p) 
    mu <- 2 + 2*X[,1] + 0.8*X[,2] - 1.2*X[,p] 
    Y <- rnorm(n, mean=mu, sd=rep(sd,n)) 
    df <- data.frame(Y,X) 
    r <- lm(Y~., data=df) 
    r0 <- lm(Y~1, data=df) 
    res <- stepAIC(r0, k=2, direction="forward", 
      scope=list(lower=~1, upper=r), trace=F) 
    return(res$coefficients) 
} 

n <- 100 
p <- 15 
sim(1, n, p) 
oo1 <- lapply(1:100, sim, n, p) 
Verwandte Themen