Как указать все уровни для функции R {stats} predic () в нелинейной смешанной модели после подбора модели с помощью пакета medrc - PullRequest
Новая
       114

Как указать все уровни для функции R {stats} predic () в нелинейной смешанной модели после подбора модели с помощью пакета medrc

0 голосов
/ 18 июня 2020

У меня есть 3 испытания (испытание: e1, e2, e3), 2 продукта / испытание (продукты: A, B), 5 оценок / продукт (0,1,10,100,1000), всего 6 кривых (кривая : c1, ..., c6). После подбора нелинейной смешанной модели я хочу построить все кривые и кривые, полученные на основе модели, на одной диаграмме. Это ссылка (пакет medr c в github): https://doseresponse.github.io/medrc/articles/medrc.html

Это код для генерации нелинейной смешанной модели. 

#packages
library(drc)
library(medrc)
library(dplyr)
library(tidyr)

#my data
trial <- c("e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1",
           "e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1",
           "e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2",
           "e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2",
           "e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3",
           "e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3")

curve <- c("c1","c1","c1","c1","c1","c1","c1","c1","c1","c1","c1","c1","c1","c1","c1",
           "c2","c2","c2","c2","c2","c2","c2","c2","c2","c2","c2","c2","c2","c2","c2",
           "c3","c3","c3","c3","c3","c3","c3","c3","c3","c3","c3","c3","c3","c3","c3",
           "c4","c4","c4","c4","c4","c4","c4","c4","c4","c4","c4","c4","c4","c4","c4",
           "c5","c5","c5","c5","c5","c5","c5","c5","c5","c5","c5","c5","c5","c5","c5",
           "c6","c6","c6","c6","c6","c6","c6","c6","c6","c6","c6","c6","c6","c6","c6")

rates <- c(.1,.1,.1,1,1,1,10,10,10,100,100,100,1000,1000,1000,
           .1,.1,.1,1,1,1,10,10,10,100,100,100,1000,1000,1000,
           .1,.1,.1,1,1,1,10,10,10,100,100,100,1000,1000,1000,
           .1,.1,.1,1,1,1,10,10,10,100,100,100,1000,1000,1000,
           .1,.1,.1,1,1,1,10,10,10,100,100,100,1000,1000,1000,
           .1,.1,.1,1,1,1,10,10,10,100,100,100,1000,1000,1000)

product <- c("A","A","A","A","A","A","A","A","A","A","A","A","A","A","A",
             "B","B","B","B","B","B","B","B","B","B","B","B","B","B","B",
             "A","A","A","A","A","A","A","A","A","A","A","A","A","A","A",
             "B","B","B","B","B","B","B","B","B","B","B","B","B","B","B",
             "A","A","A","A","A","A","A","A","A","A","A","A","A","A","A",
             "B","B","B","B","B","B","B","B","B","B","B","B","B","B","B")

resp <- c(.295,.3232,.3015,.155,.1501,.1483,.0511,.036,.0445,.0021,.0022,.0035,.0015,.0025,.0009,         
      .312,.3373,.2994,.265,.2501,.2482,.1022,.103,.1142,.0220,.0198,.0159,.0036,.0099,.0100,
      .289,.3122,.3093,.141,.1612,.1398,.0722,.022,.0581,.0019,.0015,.0011,.0018,.0009,.0014,
      .325,.3451,.2952,.267,.2412,.2398,.1125,.109,.1019,.0554,.0547,.0118,.0029,.0075,.0078,
      .294,.2452,.2991,.121,.1925,.1485,.0871,.025,.0658,.0019,.0019,.0010,.0025,.0019,.0008,
      .285,.3412,.3069,.124,.1861,.1958,.1276,.132,.1985,.0325,.0201,.0225,.0031,.0089,.0094)


data.test <- data.frame(trial,curve,rates,product,resp) #my data frame

#my model
m1 <- medrm(resp ~ rates, 
            curveid=b + c + d + e ~ product, 
            data = data.test, 
            fct=LL.4(), 
            random = c + d ~ 1|trial,
            start=NULL)

Чтобы построить график: 

#plotting
pdata <- data.test%>%
  group_by(curve, product) %>%
  expand(rates=exp(seq(-3, 10, length=50)))
#pdata$resp_ind <- predict(m1, newdata=pdata)
pdata$resp <- predict(m1, newdata=pdata, level=0)

ggplot(data.test, aes(x=log(rates), y=resp, 
                      colour=product, group=curve, shape=product)) +
  geom_point() +
  geom_line(data=pdata) +
  #geom_line(data=pdata, aes(y=resp_ind), linetype=2) +
  theme_bw() +
  scale_x_continuous("DOSE", 
                     breaks=log(c(.1, 1, 10, 100, 1000)), 
                     labels=c(.1, 1, 10, 100, 1000))

Обратите внимание, что две строки кода закомментированы. При извлечении данных прогноза для каждой кривой я не могу указать уровни, т.е. кривые, которые дали случайные компоненты. Что мне не хватает? 

pdata$resp_ind <- predict(m1, newdata=pdata)

приводит к ошибке: 

Error in predict.nlme(object$fit, newdata = newdata, level = level) : 
cannot evaluate groups for desired levels on 'newdata'

Поэтому я не могу использовать эту строку кода для построения каждой линии кривой 

geom_line(data=pdata, aes(y=resp_ind), linetype=2) +

Это похожие вопросы, но я получаю средний тренд с кодом: 

pdata$resp <- predict(m1, newdata=pdata, level=0)

Я хотел указать уровни, чтобы получить все кривые. R : lme, невозможно оценить группы для желаемых уровней в 'newdata' https://stats.stackexchange.com/questions/58031/prediction-on-mixed-effect-models-what-to-do-with-random-effects

1 Ответ

0 голосов
/ 18 июня 2020

Я могу определить проблему и поделюсь тем, что обнаружил.

Чтобы код графика работал так, как в вопросе, строка случайного фактора в модели должна иметь вместо этого curve из trial

    #my model
    m1 <- medrm(resp ~ rates, 
            curveid=b + c + d + e ~ product, 
            data = data.test, 
            fct=LL.4(), 
            random = c + d ~ 1|curve,
            start=NULL)


#plotting
pdata <- data.test%>%
  group_by(curve, product) %>%
  expand(rates=exp(seq(-3, 10, length=50)))
pdata$resp_ind <- predict(m1, newdata=pdata)
pdata$resp <- predict(m1, newdata=pdata, level=0)

ggplot(data.test, aes(x=log(rates), y=resp,
                      colour=product, group=curve, shape=product)) +
  geom_point() +
  geom_line(data=pdata) +
  geom_line(data=pdata, aes(y=resp_ind), linetype=2) +
  theme_bw() +
  scale_x_continuous("DOSE",
                     breaks=log(c(.1, 1, 10, 100, 1000)),
                     labels=c(.1, 1, 10, 100, 1000))

Другие модели с другими случайными параметрами с trial должны иметь trial в group_by для построения: 

#my model
m2 <- medrm(resp ~ rates,
            curveid=b + c + d + e ~ product,
            data = data.test,
            fct=LL.4(),
            random = c + d ~ 1|trial/curve,
            start=NULL)

#plotting
pdata <- data.test%>%
  group_by(trial, curve, product) %>%
  expand(rates=exp(seq(-3, 10, length=50)))
pdata$resp_ind <- predict(m2, newdata=pdata)
pdata$resp <- predict(m2, newdata=pdata, level=0)

ggplot(data.test, aes(x=log(rates), y=resp,
                      colour=product, group=curve, shape=product)) +
  geom_point() +
  geom_line(data=pdata) +
  geom_line(data=pdata, aes(y=resp_ind), linetype=2) +
  theme_bw() +
  scale_x_continuous("DOSE",
                     breaks=log(c(.1, 1, 10, 100, 1000)),
                     labels=c(.1, 1, 10, 100, 1000))

Правильная модель для использование зависит от каждого случая, и это другой предмет.

...