fill in NA values with predicted values from linear regression for grouped data in R-CodePudding

I want to fill in NAs in the prop_below5 variable for those depths (depth_round) where I don't have measurements using a linear interpolation and keep the actual measured value when available. I want to fit a linear regression to each Cruise/Cast data group.

When I subset the data selecting the first Cruise and Cast and apply lm() and predict() I get correct values:

OS2017 <- data[data$Cruise=="OS1701" & data$Cast==1,]
OS2017$prop_below5_interp <- ifelse(is.na(OS2017$prop_below5), predict(lm(OS2017$prop_below5 ~ OS2017$depth_round), OS2017), OS2017$prop_below5)
head(OS2017)

  Cruise Cast depth_round prop_below5 prop_below5_interp
1 OS1701    1           0          NA         0.04478375
2 OS1701    1           1  0.04142012         0.04142012
3 OS1701    1           2          NA         0.04351924
4 OS1701    1           3          NA         0.04288698
5 OS1701    1           4          NA         0.04225472
6 OS1701    1           5          NA         0.04162247

When I try to implement this within dplyr to do the same thing for each Cruise/Cast group I don't get the same answer...

data <- data %>% group_by(Cruise, Cast) %>% mutate(prop_below5_interp=ifelse(is.na(prop_below5), predict(lm(prop_below5 ~ depth_round,.), .), prop_below5))
data <- as.data.frame(data)
head(data)

  Cruise Cast depth_round prop_below5 prop_below5_interp
1 OS1701    1           0          NA         0.29926186
2 OS1701    1           1  0.04142012         0.04142012
3 OS1701    1           2          NA         0.29258267
4 OS1701    1           3          NA         0.28924308
5 OS1701    1           4          NA         0.28590349
6 OS1701    1           5          NA         0.28256390
>

How do I do this properly within dplyr?

Here is a snippet of data:

structure(list(Cruise = c("OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", "OS1701", 
"OS1701", "OS1701", "OS1701", "OS1701", "OS1701"), Cast = c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 4L, 
4L, 4L, 4L), depth_round = c(0L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 
8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 
21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 
34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 
47L, 48L, 49L, 50L, 51L, 52L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 
8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 
21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 
34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 0L, 1L, 2L, 3L), 
    prop_below5 = c(NA, 0.0414201183431953, NA, NA, NA, NA, NA, 
    NA, NA, NA, 0.0436507936507936, NA, NA, NA, NA, NA, NA, NA, 
    NA, NA, 0.0296803652968037, NA, NA, NA, NA, NA, NA, NA, NA, 
    NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
    NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
    NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
    NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
    NA, NA, NA, NA, NA, NA, NA, NA, 0.4, NA, NA)), row.names = c(53L, 
33L, 48L, 36L, 51L, 38L, 23L, 12L, 4L, 43L, 31L, 41L, 46L, 49L, 
45L, 35L, 9L, 50L, 25L, 26L, 40L, 14L, 42L, 52L, 44L, 30L, 19L, 
16L, 5L, 1L, 10L, 7L, 18L, 17L, 8L, 32L, 20L, 34L, 11L, 21L, 
2L, 24L, 15L, 6L, 3L, 47L, 22L, 39L, 37L, 13L, 27L, 28L, 29L, 
84L, 58L, 85L, 86L, 93L, 89L, 75L, 77L, 91L, 66L, 67L, 81L, 82L, 
54L, 92L, 71L, 95L, 61L, 63L, 73L, 87L, 74L, 62L, 64L, 65L, 56L, 
60L, 59L, 70L, 72L, 96L, 88L, 76L, 90L, 79L, 78L, 80L, 94L, 57L, 
55L, 68L, 69L, 83L, 98L, 99L, 126L, 134L), class = "data.frame")
>

CodePudding user response：

If your question is "why doesn't dplyr work" I can't help you there. But what you ask is easily achievable using data.table.

library(data.table)
setDT(df)
#   count of non-NA values by group
df[, N:=sum(!is.na(prop_below5)), by=.(Cruise, Cast)]
#   need at least 2 non-NA values or lm will fail
df[N > 1
   , interp:=predict(lm(prop_below5~depth_round, data=.SD[!is.na(prop_below5)]), newdata=.SD)
   , by=.(Cruise, Cast)]

The code above identifies groups with insufficient non-NA values to use lm(...), then runs the interpolation/extrapolation on the groups where it is possible. In the data you provided, all the groups other than the first one have either 0 or 1 non-NA value.