Смешанная модель на каждом уровне факторов

Question

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

Вот пример. «x» и «y» - переменные, а «fac 2» - фактор, который должен моделироваться как случайный эффект. Что я хочу сделать, так это применить модель к каждому уровню «fac1».

Вот как выглядит фрейм данных:

> head(df)
   x      y fac1 fac2
1 14 1.0328    A    1
2 18 1.0205    A    1
3 22 1.9262    A    1
4 26 2.3026    A    1
5 30 2.5159    A    1
6 34 2.6633    A    1

Вот модель

lp<- function(x, a, b, c){
  ifelse(x < c, a + b * x, a + b * c)
}

library(nlme)
f1<- nlme(y ~ lp(x, a, b, c),
           random = a+b+c~1|fac2,
           fixed = a+b+c,
           data = df,
           start = c(a=1, b=0.05, c=40,40))

Я хочу получить что-то подобное, но для каждого уровня «fac 1» (в данном случае это «A» и «B»):

> fixef(f1)
         a          b          c 
-0.4653545  0.1025822 30.3707837 
> ranef(f1)
            a             b          c
1 -0.02172505  3.580473e-13  0.4892547
2 -0.16928799  1.110698e-12  1.5177167
3 -0.20562406  1.196632e-12  1.6351403
4 -0.34252338  1.742883e-12  2.3815664
5  0.29974892 -1.608822e-12 -2.1983787
6  0.17999633 -8.927181e-13 -1.2198569
7  0.17491778 -1.083529e-12 -1.4805913
8  0.08449744 -8.231909e-13 -1.1248512

В случае, если вам нужны данные примера:

df<-structure(list(x = c(14L, 18L, 22L, 26L, 30L, 34L, 38L, 40L, 
42L, 14L, 18L, 22L, 26L, 30L, 34L, 38L, 40L, 42L, 14L, 18L, 22L, 
26L, 30L, 34L, 38L, 42L, 44L, 14L, 18L, 22L, 26L, 30L, 34L, 38L, 
40L, 42L, 10L, 14L, 18L, 22L, 26L, 30L, 34L, 36L, 38L, 42L, 10L, 
14L, 18L, 22L, 26L, 30L, 34L, 38L, 42L, 10L, 14L, 18L, 22L, 26L, 
30L, 34L, 37L, 38L, 42L, 10L, 14L, 18L, 22L, 26L, 30L, 34L, 36L, 
38L, 42L), y = c(1.0328, 1.0205, 1.9262, 2.3026, 2.5159, 2.6633, 
2.7435, 2.855, 2.6624, 0.881, 1.0738, 1.6, 2.1519, 2.3339, 2.4908, 
2.7169, 2.7106, 2.7731, 0.6859, 1.1838, 1.6867, 1.9957, 2.3212, 
2.5685, 2.6384, 2.557, 2.7263, 0.6374, 1.062, 1.5332, 1.8987, 
2.0687, 2.5393, 2.5393, 2.5241, 2.4515, 0.7777, 1.3255, 1.7045, 
2.2217, 2.4302, 2.7138, 2.7788, 2.733, 2.7156, 2.741, 0.6405, 
1.1825, 1.5864, 2.0159, 2.3993, 2.7801, 2.7167, 2.8142, 2.6103, 
0.6804, 1.0521, 1.7468, 1.8914, 2.4697, 2.773, 2.6471, 2.4977, 
2.76, 2.595, 0.7479, 1.0025, 1.4848, 1.8616, 2.3183, 2.6273, 
2.5209, 2.6643, 2.4964, 2.4766), fac1 = structure(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, 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), .Label = c("A", "B"), class = "factor"), 
    fac2 = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 
    5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
    7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L), .Label = c("1", "2", "3", "4", "5", 
    "6", "7", "8"), class = "factor")), row.names = c(NA, -75L
), class = "data.frame")