Как настроить данные для использования термина взаимодействия в ежедневных временных рядах

Question

У меня есть ежедневная смертность по всем причинам и стратифицирована различными категориями заболеваний, и мне хотелось выяснить, изменяются ли ассоциации между смертью и pm10 различными категории болезней. В данных примера смертность представляет собой общую ежедневную смерть, cvd те, кто умер из-за болезни сердца, а «другие» представляют собой все смертельные случаи, которые не связаны с сердечными заболеваниями. Чтобы моделировать ассоциации между pm10 и различными результатами, я использую следующие сценарии.

m1<-gam(death ~ pm10 + s(trend)+ s(temp), data=df1, na.action=na.omit, family=poisson) 
m2<-gam(cvd ~ pm10 + s(trend)+ s(temp), data=df1, na.action=na.omit, family=poisson) 
m3<-gam(others ~ pm10 + s(trend)+ s(temp), data=df1, na.action=na.omit, family=poisson) 

Для каждого дня смертельных исходов, из которых большинство из-за ССЗ. 1 января 1987 года было 130 смертей (65 случаев смерти от ССЗ и 65 из других причин). Моя цель - выяснить, есть ли разница в смерти в группе ССЗ и другие причины, связанные с воздействием ТЧ10. Вопрос исследования: имеет ли коэффициент смертности у пациентов с сердечно-сосудистыми заболеваниями и другими людьми, когда они подвергаются воздействию ТЧ10. В стратифицированном анализе я мог бы разделить данные на группу CVD и другую группу. Но в этой задаче мне интересно запустить модель с термином взаимодействия. Но я не мог понять, как это сделать. Я подумал о том, чтобы развернуть каждую строку дважды и создать фиктивную переменную для обеих групп (1 для других и 0 для CVD) и один столбец (newdeath), который содержит две строки для каждого дня, чтобы представлять смерть из-за других против CVD. При этом настройки (набор данных df2 показано ниже) Я хотел запустить следующий код:

minter<-gam(newdeath~ pm10*dummy + s(trend)+ s(temp), data=df2, na.action=na.omit, family=poisson) 

Однако я не уверен, действительно ли дает мне эту форму данные и модель для достижения того, что я хотел.

Следующий код будет производить выборки набора данных

library(mgcv) 
require(dlnm) 
df <- chicagoNMMAPS 
df <- chicagoNMMAPS 
df1 <- df[,c("date","dow","death","cvd","temp","pm10")] 
df1$trend<-seq(dim(df1)[1]) 
df1$others<-df1$death-df1$cvd # all other non-CVD deaths 

Я рассмотрел следующую дату установить, чтобы решить эту проблему, но не уверены, правильно ли это.

> dput(df2) 
structure(list(date = structure(c(6209, 6209, 6210, 6210, 6211, 
6211, 6212, 6212, 6213, 6213), class = "Date"), dow = structure(c(5L, 
5L, 6L, 6L, 7L, 7L, 1L, 1L, 2L, 2L), .Label = c("Sunday", "Monday", 
"Tuesday", "Wednesday", "Thursday", "Friday", "Saturday"), class = "factor"), 
    death = c(130L, 130L, 150L, 150L, 101L, 101L, 135L, 135L, 
    126L, 126L), cvd = c(65L, 65L, 73L, 73L, 43L, 43L, 72L, 72L, 
    64L, 64L), temp = c(-0.277777777777778, -0.277777777777778, 
    0.555555555555556, 0.555555555555556, 0.555555555555556, 
    0.555555555555556, -1.66666666666667, -1.66666666666667, 
    0, 0), pm10 = c(26.956073733, 26.956073733, NA, NA, 32.838694951, 
    32.838694951, 39.9560737332, 39.9560737332, NA, NA), trend = c(1L, 
    1L, 2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L), newdeath = c(65L, 65L, 
    77L, 73L, 58L, 43L, 63L, 72L, 62L, 64L), dummy = c(1, 0, 
    1, 0, 1, 0, 1, 0, 1, 0)), datalabel = "Written by R.    ", time.stamp = "24 Mar 2015 00:00", .Names = c("date", 
"dow", "death", "cvd", "temp", "pm10", "trend", "newdeath", "dummy" 
), formats = c("%dD_m_Y", "%9.0g", "%9.0g", "%9.0g", "%9.0g", 
"%9.0g", "%9.0g", "%9.0g", "%9.0g"), types = c(255L, 253L, 253L, 
253L, 255L, 255L, 253L, 253L, 254L), val.labels = c("", "dow", 
"", "", "", "", "", "", ""), var.labels = c("date", "dow", "death", 
"cvd", "temp", "pm10", "trend", "others", ""), row.names = c("1", 
"2", "3", "4", "5", "6", "7", "8", "9", "10"), version = 12L, label.table = structure(list(
    dow = structure(1:7, .Names = c("Sunday", "Monday", "Tuesday", 
    "Wednesday", "Thursday", "Friday", "Saturday"))), .Names = "dow"), class = "data.frame") 

Answer 1

Я согласен с @BondedDust ... кажется, что ваши данные уже слишком много агрегируются, чтобы ответить на ваш вопрос, да еще и с инструментами, которые вы собираетесь. То, что вы могли бы измерить, однако, соотношение между PM10 и доля смертей от болезней сердца:

df1$prop.cvd <- df1$cvd/df1$death 

А затем визуализировать

plot(prop.cvd ~ pm10, data = df1) 

И, возможно, использовать эту переменную в качестве переменной отклика в регрессионной модели , так вот - но тогда заданный вопрос отличается от вашего и не учитывает никакой задержки в результате pm10 и времени. Для этого вам понадобятся другие инструменты, но я не в состоянии помочь вам здесь. Возможно, задание вопроса о Cross-Validated может помочь вам двигаться дальше.

model <- glm(prop.cvd ~ pm10 + temp, data = df1) 
summary(model) 

Call: 
glm(formula = prop.cvd ~ pm10 + temp, data = df1) 

Deviance Residuals: 
    Min  1Q Median  3Q  Max 
-0.1632 -0.0351 -0.0014 0.0349 0.3323 

Coefficients: 
       Estimate Std. Error t value Pr(>|t|)  
(Intercept) 4.435e-01 1.529e-03 290.051 <2e-16 *** 
pm10   5.294e-05 4.161e-05 1.272 0.203  
temp  -6.436e-04 7.437e-05 -8.654 <2e-16 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for gaussian family taken to be 0.002732995) 

    Null deviance: 13.498 on 4862 degrees of freedom 
Residual deviance: 13.282 on 4860 degrees of freedom 
    (251 observations deleted due to missingness) 
AIC: -14898 

Number of Fisher Scoring iterations: 2