У меня проблема с циклом R, я хочу сгенерировать конкретный результат с прогнозами и доверительными интервалами, применяя линейную модель к разным выборкам

#r #dataframe #loops

#r #фрейм данных #циклы

Вопрос:

У меня есть совокупность с N наблюдениями, подобными этому

 *Y  X   ID 
……  ….. 1 
……  …   2 
……  ……. 3 
……  ….. . 
……. ……..    .*

Я сгенерировал этот код для взятия разных выборок и применил к ним линейную модель:

 N=1000
X=rnorm(N,2,1)
Y=8*X rnorm(N,0,1)
POP=cbind(X,Y)
POPULATION=as.data.frame(POP)
POPULATION$ID=seq.int(nrow(population))
J=10
n=100
PREDICTIONS=matrix(,nrow = n,ncol=J) 
for (i in 1:J) {
 SAMPLE=POPULATION[sample(nrow(POPULATION),size = n,replace = F),] 
  Y1=SAMPLE$Y 
  X1=SAMPLE$X 
  LM=lm(Y1~X1) 
  PREDICTIONS[,i]=as.array(predict(LM,SAMPLE)) 
}

Я хочу объединить прогнозы и доверительные интервалы с кадром данных совокупности. То есть я хочу что-то вроде этого :

 ID  Estimate1   LW  UP  Estimate2   LW  UP  …   ….  ….  
1   NA  NA  NA  8.25    4.3 5.7 NA  NA  NA  
2   3.5 1.2 4.2 NA  NA  NA  NA  NA  NA 
3   NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  ... ... . .
4   7.8 4.2 10.5    7.14    6.2 8.1 NA  NA  NA ....... 
5   .   .   .   .   .   .   .   .   . 
.   .   .   .   .   .   .   .   .   . 
.   .   .   .   .   .   .   .   .   .*

Как я могу настроить цикл, чтобы получить что-то подобное?

Ответ №1:

Вот как вы могли бы это сделать.

 set.seed(2) # with a seed your example is reproducible!

N <- 1000
X <- rnorm(N,2,1)
Y <- 8*X   rnorm(N,0,1)
POPULATION <- data.frame(X = X, Y = Y, ID = seq_len(N))
J <- 10
n <- 100

for (i in 1:J) {

 rows <- sample(nrow(POPULATION), size = n, replace = FALSE)

 SAMPLE <- POPULATION[rows,]
 LM <- lm(Y~X, SAMPLE)
 PR <- predict(LM, SAMPLE, interval = "confidence")
 cols <- paste(colnames(PR), i, sep = "_")
 POPULATION[rows,cols] <- asplit(PR,2)

}

head(POPULATION)[1:9]

#>           X         Y ID fit_1 lwr_1 upr_1    fit_2    lwr_2    upr_2
#> 1 1.1030855  9.290884  1    NA    NA    NA       NA       NA       NA
#> 2 2.1848492 18.433460  2    NA    NA    NA       NA       NA       NA
#> 3 3.5878453 27.755556  3    NA    NA    NA       NA       NA       NA
#> 4 0.8696243  6.995558  4    NA    NA    NA       NA       NA       NA
#> 5 1.9197482 14.527104  5    NA    NA    NA 15.57564 15.38493 15.76634
#> 6 2.1324203 17.616534  6    NA    NA    NA       NA       NA       NA

Однако, подобным образом, вы получите много недостающих данных POPULATION .

Вы уверены, что не хотите применять predict ко всем данным?

Вот так:

 for (i in 1:J) {
 
 rows <- sample(nrow(POPULATION), size = n, replace = FALSE)
 
 SAMPLE <- POPULATION[rows,] 
 LM <- lm(Y~X, SAMPLE) 
 PR <- predict(LM, POPULATION, interval = "confidence")
 cols <- paste(colnames(PR), i, sep = "_")
 POPULATION[,cols] <- asplit(PR,2)
 
}

head(POPULATION)[1:9]
#>          X         Y ID     fit_1     lwr_1     upr_1     fit_2     lwr_2     upr_2
#> 1 1.1030855  9.290884  1  8.782858  8.498869  9.066846  9.018652  8.741132  9.296172
#> 2 2.1848492 18.433460  2 17.395832 17.181911 17.609754 17.704131 17.516264 17.891998
#> 3 3.5878453 27.755556  3 28.566451 28.186621 28.946281 28.968784 28.610305 29.327263
#> 4 0.8696243  6.995558  4  6.924046  6.606492  7.241600  7.144193  6.829553  7.458833
#> 5 1.9197482 14.527104  5 15.285106 15.072070 15.498142 15.575636 15.384929 15.766343
#> 6 2.1324203 17.616534  6 16.978395 16.765726 17.191063 17.283179 17.096002 17.470357