2017-05-21 5 views
0

Je dois coder une étude de simulation en R. J'ai X1, ..., X15 ~ N (0,1) variables explicatives et Y ~ N (2 + 2 * X1 + 0.8 * X2-1.2 * X15, 1) et j'ai besoin de simuler n = 100 valeurs et répéter que pour iter = 100 fois. Ensuite, pour chaque modèle linéaire créé, je dois trouver le meilleur sous-modèle, en utilisant stepAIC. J'ai écrit le code suivant:comment mener une étude de simulation en utilisant stepAIC

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) 

Comme je suis un inexpérimenté R-utilisateur, je pense que je fais quelque chose de mal. Le but de l'étude est de trouver si dans les 100 meilleurs sous-modèles (selon l'étape AIC), il y a des modèles qui peuvent trouver le réel (Y = 2 + 2 * X1 + 0.8 * X2-1.2 * X15 + e) . Au cas où je ferais les mauvaises choses, pourrais-je obtenir de l'aide/des conseils pour l'implémenter correctement?

Répondre

0

est ici une version de travail de votre code:

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)