编辑代码

library(MASS)
L <- 10^6
Beta<-cbind(rep(1, 13),seq(13:1)^-1,seq(13:1)^-1)
beta<-Beta[,1]
p <- length(beta)
Rho<-c(0.3,0.6,0.9)
RR <- seq(0.1, 0.9, by = 0.1)
e<-rep(NA,L)
N<-c(200,400)
####数据构造####
#for (i in 1:3) { 
  rho<-Rho[1]
  #set.seed(123)  
  cov_matrix <- (1-rho)* diag(12) + rho * matrix(1, 12, 12)
  M <- mvrnorm(n = L, mu = rep(1, 12), Sigma = cov_matrix)
  X <- cbind(1, M)
  
  MSE.SAIC<-MSE.SBIC<-MSE.NAIC<-MSE.NBIC<-matrix(NA,2,9)
  for (r in 1:9) {
    Rsqure<-RR[r]
    sigma <- sqrt((1 / Rsqure-1) *t(beta[2:13])%*%cov_matrix%*%beta[2:13])
    Sigma<-c(sigma)*rowSums(X^2)/(2*p-1)
    #set.seed(101)  
    #for (t in 1:L) {
    #  e[t] <- rnorm(1, mean = 0, sd = Sigma[t])
    #}
    #Y<- X %*% beta+e
    mu<-X %*% beta
    XF <- X
    #set.seed(456)  
    for (s in 2:4) {
      indices <- X[, s] < 1
      col_indices <- 3 * (s - 1) + c(2, 3, 4)
      XF[indices, col_indices] <- NA
    }
    
    mse.SAIC<-mse.SBIC<-mse.NAIC<-mse.NBIC<-rep(NA,100)
    #mse2<-rep(NA,100)
    for (v in 1:2) {
      n<-N[v]
     for (u in 1:100) {
     
      #set.seed(789)  
      M1 <- mvrnorm(n = n, mu = rep(1, 12), Sigma = cov_matrix)
      X1 <- cbind(1, M1)
      
      #sigma1 <- sqrt((1 / Rsqure-1) * var(X1 %*% beta))
      Sigma1<-c(sigma)*rowSums(X1^2)/(2*p-1)
      e1<-rep(NA,n)
      #set.seed(202)  
      for (t in 1:n) {
        e1[t] <- rnorm(1, mean = 0, sd = Sigma1[t])
      }
      Y1<- X1 %*% beta+e1
      
      X1F <- X1
      for (s in 2:4) {
        indices <- X1[, s] < 1
        col_indices <- 3 * (s - 1) + c(2, 3, 4)
        X1F[indices, col_indices] <- NA
      }
      
      CMCV <- list(
        c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13),
        c(2, 3, 4, 8, 9, 10, 11, 12, 13),
        c(2, 3, 4, 5, 6, 7, 11, 12, 13),
        c(2, 3, 4, 5, 6, 7, 8, 9, 10),
        c(2, 3, 4, 11, 12, 13),
        c(2, 3, 4, 8, 9, 10),
        c(2, 3, 4, 5, 6, 7),
        c(2, 3, 4)
      )
      
      Model<-list("M1","M2","M3","M4","M5","M6","M7","M8")
      nw<-rep(NA,8)
      AIC<-BIC<-rep(NA,8)
####备选模型与权重选择####      
      
      for (j in 1:8) {
        complete_rows <- complete.cases(X1F[,CMCV[[j]]])
        nw[j]<-sum(complete_rows)
        X1CM<-X1F[complete_rows,CMCV[[j]]]
        Y1CM<-Y1[complete_rows,]
        Model[[j]]<-glm(Y1CM~X1CM)
        AIC[j]<-Model[[j]]$aic
        BIC[j]<-BIC(Model[[j]])
      }
      
      
      
      w_SAIC<-exp(-AIC/2)/sum(exp(-AIC/2))
      w_SBIC<-exp(-BIC/2)/sum(exp(-BIC/2))
      w_NAIC<-exp(-AIC/(2*nw))/sum(exp(-AIC/(2*nw)))
      w_NBIC<-exp(-BIC/(2*nw))/sum(exp(-BIC/(2*nw)))
####MSE####      
      XFCC<-XF[complete.cases(XF),]
      lCC<-length(XFCC[,1])
      Pred<-matrix(NA,lCC,8)
      for (j in 1:8) {
        Pred[,j]<- XFCC[,c(1,CMCV[[j]])]%*%coef(Model[[j]])
      }



      
      #Pred_SAIC<-cbind(Pred1,Pred2,Pred3,Pred4,Pred5,Pred6,Pred7,Pred8)%*%w_SAIC
      Pred_SAIC<-Pred%*%w_SAIC
      Pred_SBIC<-Pred%*%w_SBIC
      Pred_NAIC<-Pred%*%w_NAIC
      Pred_NBIC<-Pred%*%w_NBIC
      
      mse.SAIC[u]<-mean((Pred_SAIC-mu[complete.cases(XF)])^2)
      mse.SBIC[u]<-mean((Pred_SAIC-mu[complete.cases(XF)])^2)
      mse.NAIC[u]<-mean((Pred_NAIC-mu[complete.cases(XF)])^2)
      mse.NBIC[u]<-mean((Pred_NBIC-mu[complete.cases(XF)])^2)
      
      }
      
      
      
   
      MSE.SAIC[v,r]<-mean(mse.SAIC)
      MSE.SBIC[v,r]<-mean(mse.SBIC)
      MSE.NAIC[v,r]<-mean(mse.NAIC)
      MSE.NBIC[v,r]<-mean(mse.NBIC)
      
    
    }
    
    
    
    #mean(mse.SAIC)
    #median(mse.SAIC)
    #mean(mse.SBIC)
    #median(mse.SBIC)
    #mean(mse.NAIC)
    #median(mse.NAIC)
    #mean(mse.NBIC)
    #median(mse.NBIC)
    
    
    
  }
#}
  
  

  MSE.SAIC
  MSE.SBIC
  MSE.NAIC
  MSE.NBIC