Как написать цикл for для расчета остатка модели из набора данных, содержащего NA?

Question

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

#dataframe
loc <- c("Loc1", "Loc2", "Loc3", "Loc4")
day <- as.numeric(c(1, 14, 20, 31, 37, 59))
empty <- expand.grid(loc,  day) 
empty <- empty %>% arrange(Var1,Var2)

response <- as.numeric(c(4398,NA, 6000.00,9234,11680,12395
                         ,2000,4273,8000,NA,NA,12762
                         ,2300,4000.00,5161,8682,12000.00,13388
                         ,NA,6225,6547,9441,7999,8688))
resp.data <- cbind(empty, response)
names(resp.data) <- c("loc", "day", "response")  

Вот что я сделал.

# run loop to calculate residuals from a linear fit
residuals <- as.data.frame(matrix(nrow = 6, ncol = 4)) 
for (i in seq_along(unique(resp.data$loc))) { 
        data_loc <- resp.data %>% filter(loc == unique(resp.data$loc)[i]) 
        model_loc <- lm(data = data_loc, 
                                  response ~ day) 
        temp <- c(resid(model_loc)) 
        if (length(temp)<6){
                temp <- c(rep('na',6-length(temp)), temp)  
        }
        residuals[i] <- temp
}

Проблема, с которой я столкнулся, заключается в том, что данные наблюдений имеют некоторые NA, и поэтому я не смогу рассчитать остаток для этого конкретного наблюдения. Я пришел с решением, но если не работал, потому что NA остатков не совпадают с NA наблюдаемых данных. Это мой результат.

# getting the final dataset with the residuals 
residuals <- residuals %>% rename_at(vars(names(residuals)), ~ unique(resp.data$loc)) %>%
        gather(key = "loc", value = "res")

resp.data$res <- residuals$res 

    loc day response               res
1  Loc1   1     4398                na
2  Loc1  14       NA  35.7766491917869
3  Loc1  20     6000 -1271.46657929227
4  Loc1  31     9234  278.234709480122
5  Loc1  37    11680  1805.52632153779
6  Loc1  59    12395 -848.071100917431
7  Loc2   1     2000                na
8  Loc2  14     4273                na
9  Loc2  20     8000 -672.182985553773
10 Loc2  31       NA -760.310593900481
11 Loc2  37       NA  1876.93820224719
12 Loc2  59    12762 -444.444622792938
13 Loc3   1     2300  274.745821042281
14 Loc3  14     4000 -806.877089478858
15 Loc3  20     5161 -929.703048180924
16 Loc3  31     8682  237.616027531956
17 Loc3  37    12000  2271.79006882989
18 Loc3  59    13388 -1047.57177974435
19 Loc4   1       NA                na
20 Loc4  14     6225 -561.709846254499
21 Loc4  20     6547 -567.168138698069
22 Loc4  31     9441  1726.49165848872
23 Loc4  37     7999 -42.9666339548574
24 Loc4  59     8688 -554.647039581289

Может ли кто-нибудь дать мне несколько советов здесь?

Заранее большое спасибо.

Answer 1

1) Для каждого подмножества выполните регрессию, используя na.action = na.exclude, вычислите их остатки, добавьте их к подмножеству и соберите все вместе.

library(dplyr)
resp.data %>%
  group_by(loc) %>%
  do(mutate(., resid = resid(lm(response ~ day, ., na.action = na.exclude)))) %>%
  ungroup

2) или без dplyr:

do.call("rbind", by(resp.data, resp.data$loc, function(x) {
   cbind(x, resid = resid(lm(response ~ day, x, na.action = na.exclude)))
}))

3) Другой подход заключается в вычислении остатков и их добавлении. Это работает здесь, но может быть немного более хрупким, поскольку предполагает, что вычисленный остаточный вектор находится в том же порядке, что и кадр входных данных.

reg.list <- by(resp.data, resp.data$loc, lm, formula = response ~ day,
  na.action = na.exclude)
transform(resp.data, resid = c(sapply(reg.list, resid)))