Я работаю над моделью SEIR для данных по кори, которая имеет следующие переменные и образцы наблюдений.
PN NAME FN HN AGE SEX PRO ERU CL DEAD IFTO SI C PR CA NI GE TD TM HNX HNY
1 Mueller 41 61 7.00 2 11.21 11.25 1 0.00 45 10 0 4 4 3 1 0 0.0 57 40
2 Mueller 41 61 6.00 2 11.23 11.27 1 0.00 45 12 0 4 4 3 1 3 40.3 57 40
3 Mueller 41 61 4.00 2 11.28 12.02 0 0.00 172 9 0 4 4 3 2 1 40.5 57 40
4 Seibold 61 62 13.00 1 11.27 11.28 2 0.00 180 10 0 1 1 1 1 3 40.7 66 41
Ниже приведен словарь данных.
# PN patient no.
# NAME patient name
# FN family index
# HN house no.
# AGE age in years
# SEX 0 = unknown
# 1 = male
# 2 = female
# PRO date of prodromes (month.day)
# ERU date of rash (month.day)
# CL class: 0 = preschool
# 1 = 1st class
# 2 = 2nd class
# DEAD 0.00 = survived
# month.day of death
# IFTO no. of patient who is the putative source
# of infection according to H. Oesterle (0 = unknown)
# SI serialinterval = no. of days between dates of prodromes
# of infection source and infected person
# C complications: 0 = no complicatons
# 1 = bronchopneumonia
# 2 = severe bronchitis
# 3 = lobar pneumonia
# 4 = pseudocroup
# 5 = cerebral edema
# PR duration of prodromes in days
# CA number of cases in family
# NI number of initial cases
# GE generation number of the case
# TD day of max. fever (days after rush)
# TM max. fever (degree celsius)
# HNX x coordinate of house (in 2.5m units)
# HNY y coordinate of house (in 2.5m units)
Здесь классы столбцов представляет класс (0,1 и 2).
Я хочу рассчитать максимальное логарифмическое правдоподобие с помощью функции на основе классов (0,1,2), на основе классов, которые можно назвать «классная» модель.
Во-первых, люди ведут себя так же, как и раньше, в отношении периода воздействия и периода заражения. Эта часть модели - это то, что я сделал.
Новая часть - это трансмиссия. В исходной модели
S (t + 1) ~ Binomial (S (t), q ^ I (t)).
В модели классных комнат я должен учитывать два параметра , q и q c где каждый человек находится в классе 1, классе 2 или классе 0 (= дошкольное учреждение).
Я хочу рассмотреть S_k (t) и I_k (t) для количества восприимчивые и инфекционные в классе k.
Первое, что мне нужно, это написать функцию для оценки этих количеств. Обратите внимание, что S (t) = S_0 (t) + S_1 (t) + S_2 (t), et c.
Идея состоит в том, что q c - это вероятность того, что восприимчивый в классе 1 (или 2) позволяет избежать заражения одним возбудителем того же класса. Такой восприимчивый человек также может быть заражен инфекциями двух других классов, но с вероятностью q.
Это означает, что (при условии обычной независимости)
S_1 (t + 1) ~ Binomial ( S_1 (t), (q c) ^ I_1 (t) q ^ {I_0 (t) + I_2 (t)})
Поскольку восприимчивый в классе 1 остается восприимчивым тогда и только тогда, когда они избегают заражение от (i) возбудителей I_1 (t) класса 1 и (ii) других возбудителей инфекции в популяции, из которых I_0 (t) + I_2 (t).
Имеется аналогичное выражение для S_2 (t + 1).
S_0 (t + 1) немного отличается, поскольку дети класса 0 не ходят в школу. Итак,
S_0 (t + 1) ~ Binomial (S_0 (t), q ^ I (t))
, где I (t) = I_0 (t) + I_1 (t) + I_2 (t).
Итак, я изо всех сил пытаюсь написать функцию для вероятности (или логарифмической вероятности), которая теперь является функцией двух параметров (q и q c). Я могу использовать optim в R, чтобы найти максимум.
Вот моя предыдущая функция правдоподобия журнала без учета классов.
ress=rep(0,50)
likelihood<-function(q)
{
for(i in 1:49){
ress[i+1]=(Current_dayStats$St[i]-Current_dayStats$St[i+1])*log(1-q^(Current_dayStats$It[i]))+(Current_dayStats$St[i+1])*log(q^(Current_dayStats$It[i]))
}
# ress[is.infinite(ress)]<-0
# ress[is.nan(ress)]<-0
# ress
return (sum(ress))
}
q=seq(0.1,0.9,0.1)
optim(par<-q,fn=likelihood)
q=seq(0.01,0.99,0.01)
which.max(sapply(q,likelihood))
MLE<-q[which.max(sapply(q,likelihood))]
MLE
plot(q,sapply(q, likelihood),type = "l")
Current_dayStats - это следующий фрейм данных. Также 50 - это количество дней. Общая численность населения - 180.
Пожалуйста, поделитесь своими любезными предложениями, как я могу двигаться дальше.
It Rt Et St
1 2 1 3 180
2 2 1 5 179
3 2 1 6 178
4 2 1 7 178
5 2 1 7 176
6 2 1 7 176