Loading libraries:

library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(magrittr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(TTR)
library(forecast)
## Loading required package: timeDate
## This is forecast 7.3
library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2015). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2. http://CRAN.R-project.org/package=stargazer
library(stats)
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select

Loading data - Do not use the NARMS dataset currently available on their website. It doesn’t include the SIR column and requires a few more steps to use.

NARMS_retailmeat_public <- read.csv("narms_public.csv")
d <- NARMS_retailmeat_public

Subsetting the data to contain only the bacterial genuses Salmonella and E coli and renaming the subset sec.

sec <- subset(d,(d$GENUS=="S"|d$GENUS=="EC"))

Creating a year/month variable in sec from a combination of the variables Year and Month, which represent the year and month in which each sample was collected:

sec$ym <- as.yearmon(paste(sec$Year, sec$Month, sep = "-"))

Generating a table of proportion of isolates that were Susceptible, Resistant, or Intermediate to ceftriaxone (a cephalosporin drug, abbreviated AXO) by month. Creates the dataset sec.AXO (shorthand for preliminary salmonella e coli ceftriaxone) where ym is a yearmon variable and AXO.SIR is a categorical variable containing I, S and R with proportions indicating how many isolates that month were intermediate, susceptible, or resistant to ceftriaxone.

psec.AXO <- sec %>% group_by(ym, AXO.SIR) %>% summarise (n = n()) %>% mutate(freq = n / sum(n))
head(psec.AXO)
## Source: local data frame [6 x 4]
## Groups: ym [3]
## 
##              ym AXO.SIR     n       freq
##   <S3: yearmon>  <fctr> <int>      <dbl>
## 1      Jan 2002       R     1 0.01063830
## 2      Jan 2002       S    93 0.98936170
## 3      Feb 2002       R     6 0.06122449
## 4      Feb 2002       S    92 0.93877551
## 5      Mar 2002       R     1 0.01123596
## 6      Mar 2002       S    88 0.98876404

Subsetting psec.AXO to contain only proportions of isolates that were resistant to ceftriaxone each month (which groups Intermediate and Susceptible together by default). #This dataset - called psec.AXO for ceftriaxone - is used to create glm models for the graphs, because it generates a smooth set of predictions that can be plotted alongside montly prevalences on a graph. The equivalent for ampicillin is psec.AMP and the equivalent for chloramphenicol is psec.CHL.

psec.AXO <- subset(psec.AXO, AXO.SIR=="S")
head(psec.AXO)
## Source: local data frame [6 x 4]
## Groups: ym [6]
## 
##              ym AXO.SIR     n      freq
##   <S3: yearmon>  <fctr> <int>     <dbl>
## 1      Jan 2002       S    93 0.9893617
## 2      Feb 2002       S    92 0.9387755
## 3      Mar 2002       S    88 0.9887640
## 4      Apr 2002       S    91 0.9891304
## 5      May 2002       S    96 0.9142857
## 6      Jun 2002       S    95 0.9793814

Doing the same for sec.axo, which is a dataset that is similar except that it is also grouped by meat type. #This dataset - called sec.AXO, sec.AMP, or sec.CHL - is used to create the GLM models that generate the coefficients and p-values.

sec.AXO <- sec %>% group_by(ym, MEAT_TYPE, AXO.SIR) %>% summarise (n = n()) %>% mutate(freq = n / sum(n))
sec.AXO <- subset(sec.AXO, AXO.SIR=="S")
head(sec.AXO)
## Source: local data frame [6 x 5]
## Groups: ym, MEAT_TYPE [6]
## 
##              ym MEAT_TYPE AXO.SIR     n      freq
##   <S3: yearmon>    <fctr>  <fctr> <int>     <dbl>
## 1      Jan 2002        CB       S    26 0.9629630
## 2      Jan 2002        GB       S    20 1.0000000
## 3      Jan 2002        GT       S    35 1.0000000
## 4      Jan 2002        PC       S    12 1.0000000
## 5      Feb 2002        CB       S    31 0.9117647
## 6      Feb 2002        GB       S    23 1.0000000

Creating a date field in each:, where each timepoint is treated as the 15th of the month (isolates were collected throughout each month) and a dummy variable for before/after the law restricting cephalosporins went into effect (assuming that it was in effect starting April 5, 2012), and creating a dateLaw object that is just the date of the law:

psec.AXO$date <- strptime(paste("15", psec.AXO$ym), "%d %b %Y")
psec.AXO$date <- as.Date(psec.AXO$date)
psec.AXO$dummy <- ifelse(psec.AXO$date<"2012-04-05",0,1)
sec.AXO$date <- strptime(paste("15", sec.AXO$ym), "%d %b %Y")
sec.AXO$date <- as.Date(sec.AXO$date)
sec.AXO$dummy <- ifelse(sec.AXO$date<"2012-04-05",0,1)
dateLaw <- as.Date("2012-04-05")

Creating a GLM object that predicts the change in slope for a linear model after the date the law went into effect, with no covariates:

glm.sec.AXO <- glm(sec.AXO$freq ~ sec.AXO$date + pmax(sec.AXO$date - dateLaw, 0) + sec.AXO$MEAT_TYPE)
summary(glm.sec.AXO)
## 
## Call:
## glm(formula = sec.AXO$freq ~ sec.AXO$date + pmax(sec.AXO$date - 
##     dateLaw, 0) + sec.AXO$MEAT_TYPE)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.251966  -0.021846   0.007757   0.033788   0.110491  
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1.067e+00  2.565e-02  41.603  < 2e-16 ***
## sec.AXO$date                    -1.474e-05  1.852e-06  -7.959 8.01e-15 ***
## pmax(sec.AXO$date - dateLaw, 0)  6.670e-05  9.015e-06   7.399 4.36e-13 ***
## sec.AXO$MEAT_TYPEGB              1.157e-01  5.489e-03  21.083  < 2e-16 ***
## sec.AXO$MEAT_TYPEGT              7.225e-02  5.440e-03  13.281  < 2e-16 ***
## sec.AXO$MEAT_TYPEPC              1.158e-01  5.458e-03  21.223  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.002404379)
## 
##     Null deviance: 3.1244  on 641  degrees of freedom
## Residual deviance: 1.5292  on 636  degrees of freedom
## AIC: -2041.7
## 
## Number of Fisher Scoring iterations: 2
confint(glm.sec.AXO)
## Waiting for profiling to be done...
##                                         2.5 %        97.5 %
## (Intercept)                      1.016922e+00  1.117476e+00
## sec.AXO$date                    -1.837256e-05 -1.111152e-05
## pmax(sec.AXO$date - dateLaw, 0)  4.903053e-05  8.436859e-05
## sec.AXO$MEAT_TYPEGB              1.049637e-01  1.264800e-01
## sec.AXO$MEAT_TYPEGT              6.158668e-02  8.291145e-02
## sec.AXO$MEAT_TYPEPC              1.051460e-01  1.265422e-01

Limitations with this model: 1. it treats the effect of the intervention as something that increases linearly starting from the date the law was passed. This is probably fine for the first few years, but would not be as reasonable for a long-term evaluation. 2. it assumes linear relationships between the predictors and outcomes. See notes in later sections - smoothing the prevalence data over time generates a trend prior to the implementation of the law that does not look linear. Instead it seems to increase linearly up to a point and then plateau in 2011 (maybe when law was first discussed). 3. it assumes observations are independent

Making a graph to display on the poster:

glm.psec.AXO <- glm(psec.AXO$freq ~ psec.AXO$date + pmax(psec.AXO$date - dateLaw, 0))
summary(glm.psec.AXO)
## 
## Call:
## glm(formula = psec.AXO$freq ~ psec.AXO$date + pmax(psec.AXO$date - 
##     dateLaw, 0))
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.116198  -0.018643   0.007962   0.022407   0.063244  
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)
## (Intercept)                       1.216e+00  3.256e-02  37.341  < 2e-16
## psec.AXO$date                    -2.126e-05  2.368e-06  -8.978 7.19e-16
## pmax(psec.AXO$date - dateLaw, 0)  8.819e-05  1.085e-05   8.126 1.13e-13
##                                     
## (Intercept)                      ***
## psec.AXO$date                    ***
## pmax(psec.AXO$date - dateLaw, 0) ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.0009905636)
## 
##     Null deviance: 0.24554  on 162  degrees of freedom
## Residual deviance: 0.15849  on 160  degrees of freedom
## AIC: -659.96
## 
## Number of Fisher Scoring iterations: 2
confint(glm.psec.AXO)
## Waiting for profiling to be done...
##                                          2.5 %        97.5 %
## (Intercept)                       1.152105e+00  1.279749e+00
## psec.AXO$date                    -2.590572e-05 -1.662165e-05
## pmax(psec.AXO$date - dateLaw, 0)  6.691811e-05  1.094566e-04
newdata <- data.frame(psec.AXO$date)
newdata$freq <- predict(glm.psec.AXO, newdata = newdata)
plot.glm.psec.AXO <- plot(psec.AXO$freq ~ psec.AXO$date, type = "p", main = "Frequency of Susceptibility to Cephalosporin in E. coli and Salmonella in Retail Meat Before and After 2012 Restriction", xlab = "Time", ylab = "Proportion of Isolates susceptible") + lines(newdata$freq ~ newdata$psec.AXO.date, col = "blue") + abline(v=as.Date("2012-04-05"),col=3,lty=3)

Repeating with GLM predictions fit before and after the law:

Generating estimates:

before = sec.AXO$date<"2012-04-05"
sec.AXO.bef <- subset(sec.AXO, (before))
sec.AXO.aft <- subset(sec.AXO, (!before))
glm.sec.AXO.MT.bef <- glm(sec.AXO.bef$freq ~ sec.AXO.bef$date + sec.AXO.bef$MEAT_TYPE)
glm.sec.AXO.MT.aft <- glm(sec.AXO.aft$freq ~ sec.AXO.aft$date + sec.AXO.aft$MEAT_TYPE)
summary(glm.sec.AXO.MT.bef)
## 
## Call:
## glm(formula = sec.AXO.bef$freq ~ sec.AXO.bef$date + sec.AXO.bef$MEAT_TYPE)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.216733  -0.020483   0.007288   0.032727   0.114273  
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              1.067e+00  2.698e-02  39.569  < 2e-16 ***
## sec.AXO.bef$date        -1.530e-05  1.965e-06  -7.787 4.16e-14 ***
## sec.AXO.bef$MEAT_TYPEGB  1.253e-01  6.006e-03  20.862  < 2e-16 ***
## sec.AXO.bef$MEAT_TYPEGT  8.393e-02  6.006e-03  13.974  < 2e-16 ***
## sec.AXO.bef$MEAT_TYPEPC  1.228e-01  6.006e-03  20.438  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.00221835)
## 
##     Null deviance: 2.4767  on 491  degrees of freedom
## Residual deviance: 1.0803  on 487  degrees of freedom
## AIC: -1603.4
## 
## Number of Fisher Scoring iterations: 2
confint(glm.sec.AXO.MT.bef)
## Waiting for profiling to be done...
##                                 2.5 %        97.5 %
## (Intercept)              1.014509e+00  1.120250e+00
## sec.AXO.bef$date        -1.915113e-05 -1.144913e-05
## sec.AXO.bef$MEAT_TYPEGB  1.135233e-01  1.370660e-01
## sec.AXO.bef$MEAT_TYPEGT  7.215427e-02  9.569693e-02
## sec.AXO.bef$MEAT_TYPEPC  1.109790e-01  1.345217e-01
summary(glm.sec.AXO.MT.aft)
## 
## Call:
## glm(formula = sec.AXO.aft$freq ~ sec.AXO.aft$date + sec.AXO.aft$MEAT_TYPE)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.234795  -0.022435   0.009315   0.032775   0.092503  
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             1.880e-01  2.084e-01   0.902 0.368558    
## sec.AXO.aft$date        4.400e-05  1.298e-05   3.389 0.000903 ***
## sec.AXO.aft$MEAT_TYPEGB 8.439e-02  1.243e-02   6.790 2.68e-10 ***
## sec.AXO.aft$MEAT_TYPEGT 3.588e-02  1.194e-02   3.006 0.003124 ** 
## sec.AXO.aft$MEAT_TYPEPC 9.440e-02  1.211e-02   7.795 1.15e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.002813973)
## 
##     Null deviance: 0.64758  on 149  degrees of freedom
## Residual deviance: 0.40803  on 145  degrees of freedom
## AIC: -448.38
## 
## Number of Fisher Scoring iterations: 2
confint(glm.sec.AXO.MT.aft)
## Waiting for profiling to be done...
##                                 2.5 %       97.5 %
## (Intercept)             -2.204876e-01 5.964467e-01
## sec.AXO.aft$date         1.855593e-05 6.945041e-05
## sec.AXO.aft$MEAT_TYPEGB  6.002816e-02 1.087462e-01
## sec.AXO.aft$MEAT_TYPEGT  1.248445e-02 5.928475e-02
## sec.AXO.aft$MEAT_TYPEPC  7.066580e-02 1.181409e-01
newdata.AXO.MT.bef <- data.frame(sec.AXO.bef$date)
newdata.AXO.MT.bef$freq <- predict(glm.sec.AXO.MT.bef, newdata = newdata.AXO.MT.bef)
newdata.AXO.MT.aft <- data.frame(sec.AXO.aft$date)
newdata.AXO.MT.aft$freq <- predict(glm.sec.AXO.MT.aft, newdata = newdata.AXO.MT.aft)
before.p = psec.AXO$date<"2012-04-05"
psec.AXO.bef <- subset(psec.AXO, (before.p))
psec.AXO.aft <- subset(psec.AXO, (!before.p))
glm.psec.AXO.bef <- glm(psec.AXO.bef$freq ~ psec.AXO.bef$date)
glm.psec.AXO.aft <- glm(psec.AXO.aft$freq ~ psec.AXO.aft$date)
newdata.AXO.bef <- data.frame(psec.AXO.bef$date)
newdata.AXO.bef$freq <- predict(glm.psec.AXO.bef, newdata = newdata.AXO.bef)
newdata.AXO.aft <- data.frame(psec.AXO.aft$date)
newdata.AXO.aft$freq <- predict(glm.psec.AXO.aft, newdata = newdata.AXO.aft)

plot.glm.sec.AXO.ALL <- plot(psec.AXO$freq ~ psec.AXO$date, type = "p", main = "Monthly Prevalence of Ceftriaxone Susceptibility in E. coli and Salmonella isolates in U.S. Retail Meat Before and After 2012 Restriction", xlab = "Date Isolates Collected, 2002-2015", ylab = "Proportion of Isolates Susceptible by Month", pch=20, ylim=c(.5,1))  + lines(newdata.AXO.bef$freq ~ newdata.AXO.bef$psec.AXO.bef, col = "red") + lines(newdata.AXO.aft$freq ~ newdata.AXO.aft$psec.AXO.aft, col = "red") + lines(newdata$freq ~ newdata$psec.AXO.date, col = "blue") + abline(v=as.Date("2012-04-05"),col=3,lty=3)

Repeating for Ampicillin:

Generating a table of proportion of isolates that were Susceptible, Resistant, or Intermediate to ampicillin (abbreviated AMP) by month. Creates the dataset sec.AMP (shorthand for preliminary salmonella e coli ampicillin) where ym is a yearmon variable and AMP.SIR is a categorical variable containing I, S and R with proportions indicating how many isolates that month were intermediate, susceptible, or resistant to ampicillin.

psec.AMP <- sec %>% group_by(ym, AMP.SIR) %>% summarise (n = n()) %>% mutate(freq = n / sum(n))
head(psec.AMP)
## Source: local data frame [6 x 4]
## Groups: ym [2]
## 
##              ym AMP.SIR     n       freq
##   <S3: yearmon>  <fctr> <int>      <dbl>
## 1      Jan 2002       I     3 0.03191489
## 2      Jan 2002       R    17 0.18085106
## 3      Jan 2002       S    74 0.78723404
## 4      Feb 2002       I     2 0.02040816
## 5      Feb 2002       R    12 0.12244898
## 6      Feb 2002       S    84 0.85714286

Subsetting psec.AMP to contain only proportions of isolates that were resistant to ampicillin each month (which groups Intermediate and Susceptible together by default):

psec.AMP <- subset(psec.AMP, AMP.SIR=="S")
head(psec.AMP)
## Source: local data frame [6 x 4]
## Groups: ym [6]
## 
##              ym AMP.SIR     n      freq
##   <S3: yearmon>  <fctr> <int>     <dbl>
## 1      Jan 2002       S    74 0.7872340
## 2      Feb 2002       S    84 0.8571429
## 3      Mar 2002       S    74 0.8314607
## 4      Apr 2002       S    74 0.8043478
## 5      May 2002       S    81 0.7714286
## 6      Jun 2002       S    81 0.8350515

Doing the same for sec.axo, including meat type:

sec.AMP <- sec %>% group_by(ym, MEAT_TYPE, AMP.SIR) %>% summarise (n = n()) %>% mutate(freq = n / sum(n))
sec.AMP <- subset(sec.AMP, AMP.SIR=="S")
head(sec.AMP)
## Source: local data frame [6 x 5]
## Groups: ym, MEAT_TYPE [6]
## 
##              ym MEAT_TYPE AMP.SIR     n      freq
##   <S3: yearmon>    <fctr>  <fctr> <int>     <dbl>
## 1      Jan 2002        CB       S    20 0.7407407
## 2      Jan 2002        GB       S    16 0.8000000
## 3      Jan 2002        GT       S    29 0.8285714
## 4      Jan 2002        PC       S     9 0.7500000
## 5      Feb 2002        CB       S    30 0.8823529
## 6      Feb 2002        GB       S    21 0.9130435

Creating a date field in each:, where each timepoint is treated as the 15th of the month (because isolates were collected throughout) and a dummy variable for before/after the law restricting cephalosporins went into effect (assuming that it was in effect starting April 5, 2012), and creating a dateLaw object that is just the date of the law:

psec.AMP$date <- strptime(paste("15", psec.AMP$ym), "%d %b %Y")
psec.AMP$date <- as.Date(psec.AMP$date)
psec.AMP$dummy <- ifelse(psec.AMP$date<"2012-04-05",0,1)
sec.AMP$date <- strptime(paste("15", sec.AMP$ym), "%d %b %Y")
sec.AMP$date <- as.Date(sec.AMP$date)
sec.AMP$dummy <- ifelse(sec.AMP$date<"2012-04-05",0,1)
dateLaw <- as.Date("2012-04-05")

Creating a GLM object that predicts the change in slope for a linear model after the date the law went into effect, with no covariates:

glm.sec.AMP <- glm(sec.AMP$freq ~ sec.AMP$date + pmax(sec.AMP$date - dateLaw, 0) + sec.AMP$MEAT_TYPE)
summary(glm.sec.AMP)
## 
## Call:
## glm(formula = sec.AMP$freq ~ sec.AMP$date + pmax(sec.AMP$date - 
##     dateLaw, 0) + sec.AMP$MEAT_TYPE)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.35340  -0.06264   0.00293   0.06900   0.36109  
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1.070e+00  5.345e-02  20.023  < 2e-16 ***
## sec.AMP$date                    -2.269e-05  3.860e-06  -5.879 6.68e-09 ***
## pmax(sec.AMP$date - dateLaw, 0)  1.298e-04  1.888e-05   6.873 1.50e-11 ***
## sec.AMP$MEAT_TYPEGB              1.739e-01  1.144e-02  15.212  < 2e-16 ***
## sec.AMP$MEAT_TYPEGT             -2.059e-01  1.135e-02 -18.140  < 2e-16 ***
## sec.AMP$MEAT_TYPEPC              8.187e-02  1.137e-02   7.200 1.71e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.0104353)
## 
##     Null deviance: 19.6115  on 640  degrees of freedom
## Residual deviance:  6.6264  on 635  degrees of freedom
## AIC: -1097.6
## 
## Number of Fisher Scoring iterations: 2
confint(glm.sec.AMP)
## Waiting for profiling to be done...
##                                         2.5 %        97.5 %
## (Intercept)                      9.654525e-01  1.174966e+00
## sec.AMP$date                    -3.025617e-05 -1.512583e-05
## pmax(sec.AMP$date - dateLaw, 0)  9.276802e-05  1.667792e-04
## sec.AMP$MEAT_TYPEGB              1.515354e-01  1.963609e-01
## sec.AMP$MEAT_TYPEGT             -2.281743e-01 -1.836758e-01
## sec.AMP$MEAT_TYPEPC              5.958050e-02  1.041553e-01

Issues with this model: 1. it treats the effect of the intervention as something that increases linearly starting from the date the law was passed. This is probably fine for the first few years, but would not be as reasonable for a long-term evaluation. 2. it assumes linear relationships between the predictors and outcomes. 3. it assumes observations are independent

Making a graph to display on the poster:

glm.psec.AMP <- glm(psec.AMP$freq ~ psec.AMP$date + pmax(psec.AMP$date - dateLaw, 0))
summary(glm.psec.AMP)
## 
## Call:
## glm(formula = psec.AMP$freq ~ psec.AMP$date + pmax(psec.AMP$date - 
##     dateLaw, 0))
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.302482  -0.044700  -0.000805   0.039872   0.200831  
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)
## (Intercept)                       1.228e+00  6.599e-02  18.616  < 2e-16
## psec.AMP$date                    -3.615e-05  4.800e-06  -7.532 3.45e-12
## pmax(psec.AMP$date - dateLaw, 0)  1.466e-04  2.199e-05   6.665 4.05e-10
##                                     
## (Intercept)                      ***
## psec.AMP$date                    ***
## pmax(psec.AMP$date - dateLaw, 0) ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.004067785)
## 
##     Null deviance: 0.89877  on 162  degrees of freedom
## Residual deviance: 0.65085  on 160  degrees of freedom
## AIC: -429.71
## 
## Number of Fisher Scoring iterations: 2
confint(glm.psec.AMP)
## Waiting for profiling to be done...
##                                          2.5 %        97.5 %
## (Intercept)                       1.099080e+00  1.357747e+00
## psec.AMP$date                    -4.555859e-05 -2.674481e-05
## pmax(psec.AMP$date - dateLaw, 0)  1.034771e-04  1.896796e-04
newdata <- data.frame(psec.AMP$date)
newdata$freq <- predict(glm.psec.AMP, newdata = newdata)
plot.glm.psec.AMP <- plot(psec.AMP$freq ~ psec.AMP$date, type = "p", main = "Frequency of Susceptibility to Cephalosporin in E. coli and Salmonella in Retail Meat Before and After 2012 Restriction", xlab = "Time", ylab = "Proportion of Isolates susceptible") + lines(newdata$freq ~ newdata$psec.AMP.date, col = "blue") + abline(v=as.Date("2012-04-05"),col=3,lty=3)

Repeating with added before / after predictions:

Generating estimates:

before = sec.AMP$date<"2012-04-05"
sec.AMP.bef <- subset(sec.AMP, (before))
sec.AMP.aft <- subset(sec.AMP, (!before))
glm.sec.AMP.MT.bef <- glm(sec.AMP.bef$freq ~ sec.AMP.bef$date + sec.AMP.bef$MEAT_TYPE)
glm.sec.AMP.MT.aft <- glm(sec.AMP.aft$freq ~ sec.AMP.aft$date + sec.AMP.aft$MEAT_TYPE)
summary(glm.sec.AMP.MT.bef)
## 
## Call:
## glm(formula = sec.AMP.bef$freq ~ sec.AMP.bef$date + sec.AMP.bef$MEAT_TYPE)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.34444  -0.06247   0.00835   0.06761   0.25431  
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              1.056e+00  5.759e-02  18.338  < 2e-16 ***
## sec.AMP.bef$date        -2.225e-05  4.195e-06  -5.304 1.72e-07 ***
## sec.AMP.bef$MEAT_TYPEGB  1.835e-01  1.282e-02  14.307  < 2e-16 ***
## sec.AMP.bef$MEAT_TYPEGT -1.813e-01  1.282e-02 -14.136  < 2e-16 ***
## sec.AMP.bef$MEAT_TYPEPC  8.118e-02  1.282e-02   6.331 5.53e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.01011254)
## 
##     Null deviance: 13.9878  on 491  degrees of freedom
## Residual deviance:  4.9248  on 487  degrees of freedom
## AIC: -857.03
## 
## Number of Fisher Scoring iterations: 2
confint(glm.sec.AMP.MT.bef)
## Waiting for profiling to be done...
##                                 2.5 %        97.5 %
## (Intercept)              9.432613e-01  1.169027e+00
## sec.AMP.bef$date        -3.047302e-05 -1.402859e-05
## sec.AMP.bef$MEAT_TYPEGB  1.583264e-01  2.085920e-01
## sec.AMP.bef$MEAT_TYPEGT -2.064037e-01 -1.561381e-01
## sec.AMP.bef$MEAT_TYPEPC  5.605172e-02  1.063173e-01
summary(glm.sec.AMP.MT.aft)
## 
## Call:
## glm(formula = sec.AMP.aft$freq ~ sec.AMP.aft$date + sec.AMP.aft$MEAT_TYPE)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.37495  -0.06097  -0.00284   0.05803   0.41269  
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             -9.623e-01  4.010e-01  -2.400 0.017689 *  
## sec.AMP.aft$date         1.105e-04  2.498e-05   4.422 1.92e-05 ***
## sec.AMP.aft$MEAT_TYPEGB  1.442e-01  2.378e-02   6.065 1.10e-08 ***
## sec.AMP.aft$MEAT_TYPEGT -2.843e-01  2.300e-02 -12.361  < 2e-16 ***
## sec.AMP.aft$MEAT_TYPEPC  8.630e-02  2.317e-02   3.724 0.000281 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.01030168)
## 
##     Null deviance: 5.6076  on 148  degrees of freedom
## Residual deviance: 1.4834  on 144  degrees of freedom
## AIC: -251.98
## 
## Number of Fisher Scoring iterations: 2
confint(glm.sec.AMP.MT.aft)
## Waiting for profiling to be done...
##                                 2.5 %        97.5 %
## (Intercept)             -1.748201e+00 -0.1763136837
## sec.AMP.aft$date         6.150447e-05  0.0001594329
## sec.AMP.aft$MEAT_TYPEGB  9.762732e-02  0.1908464912
## sec.AMP.aft$MEAT_TYPEGT -3.294091e-01 -0.2392443894
## sec.AMP.aft$MEAT_TYPEPC  4.087743e-02  0.1317148656
newdata.AMP.MT.bef <- data.frame(sec.AMP.bef$date)
newdata.AMP.MT.bef$freq <- predict(glm.sec.AMP.MT.bef, newdata = newdata.AMP.MT.bef)
newdata.AMP.MT.aft <- data.frame(sec.AMP.aft$date)
newdata.AMP.MT.aft$freq <- predict(glm.sec.AMP.MT.aft, newdata = newdata.AMP.MT.aft)
before.p = psec.AMP$date<"2012-04-05"
psec.AMP.bef <- subset(psec.AMP, (before.p))
psec.AMP.aft <- subset(psec.AMP, (!before.p))
glm.psec.AMP.bef <- glm(psec.AMP.bef$freq ~ psec.AMP.bef$date)
glm.psec.AMP.aft <- glm(psec.AMP.aft$freq ~ psec.AMP.aft$date)
newdata.AMP.bef <- data.frame(psec.AMP.bef$date)
newdata.AMP.bef$freq <- predict(glm.psec.AMP.bef, newdata = newdata.AMP.bef)
newdata.AMP.aft <- data.frame(psec.AMP.aft$date)
newdata.AMP.aft$freq <- predict(glm.psec.AMP.aft, newdata = newdata.AMP.aft)

plot.glm.sec.AMP.ALL <- plot(psec.AMP$freq ~ psec.AMP$date, type = "p", main = "Frequency of Susceptibility to Ampicillin in E. coli and Salmonella in U.S. Retail Meat", xlab = "Date Isolates Collected, 2002-2015", ylab = "Proportion of Susceptible Isolates by Month", pch=20, ylim=c(.5,1))  + lines(newdata.AMP.bef$freq ~ newdata.AMP.bef$psec.AMP.bef, col = "red") + lines(newdata.AMP.aft$freq ~ newdata.AMP.aft$psec.AMP.aft, col = "red") + lines(newdata$freq ~ newdata$psec.AMP.date, col = "blue") + abline(v=as.Date("2012-04-05"),col=3,lty=3)

Repeating for Chloramphenicol:

Generating a table of proportion of isolates that were Susceptible, Resistant, or Intermediate to chloramphenicol (abbreviated CHL) by month. Creates the dataset sec.CHL (shorthand for preliminary salmonella e coli chloramphenicol) where ym is a yearmon variable and CHL.SIR is a categorical variable containing I, S and R with proportions indicating how many isolates that month were intermediate, susceptible, or resistant to chloramphenicol.

psec.CHL <- sec %>% group_by(ym, CHL.SIR) %>% summarise (n = n()) %>% mutate(freq = n / sum(n))
head(psec.CHL)
## Source: local data frame [6 x 4]
## Groups: ym [2]
## 
##              ym CHL.SIR     n       freq
##   <S3: yearmon>  <fctr> <int>      <dbl>
## 1      Jan 2002       I     2 0.02127660
## 2      Jan 2002       R     2 0.02127660
## 3      Jan 2002       S    90 0.95744681
## 4      Feb 2002       I     2 0.02040816
## 5      Feb 2002       R     2 0.02040816
## 6      Feb 2002       S    94 0.95918367

Subsetting psec.CHL to contain only proportions of isolates that were resistant to chloramphenicol each month (which groups Intermediate and Susceptible together by default):

psec.CHL <- subset(psec.CHL, CHL.SIR=="S")
head(psec.CHL)
## Source: local data frame [6 x 4]
## Groups: ym [6]
## 
##              ym CHL.SIR     n      freq
##   <S3: yearmon>  <fctr> <int>     <dbl>
## 1      Jan 2002       S    90 0.9574468
## 2      Feb 2002       S    94 0.9591837
## 3      Mar 2002       S    87 0.9775281
## 4      Apr 2002       S    89 0.9673913
## 5      May 2002       S   104 0.9904762
## 6      Jun 2002       S    96 0.9896907

Doing the same for sec.axo, including meat type:

sec.CHL <- sec %>% group_by(ym, MEAT_TYPE, CHL.SIR) %>% summarise (n = n()) %>% mutate(freq = n / sum(n))
sec.CHL <- subset(sec.CHL, CHL.SIR=="S")
head(sec.CHL)
## Source: local data frame [6 x 5]
## Groups: ym, MEAT_TYPE [6]
## 
##              ym MEAT_TYPE CHL.SIR     n      freq
##   <S3: yearmon>    <fctr>  <fctr> <int>     <dbl>
## 1      Jan 2002        CB       S    26 0.9629630
## 2      Jan 2002        GB       S    19 0.9500000
## 3      Jan 2002        GT       S    34 0.9714286
## 4      Jan 2002        PC       S    11 0.9166667
## 5      Feb 2002        CB       S    33 0.9705882
## 6      Feb 2002        GB       S    23 1.0000000

Creating a date field in each:, where each timepoint is treated as the 15th of the month (because isolates were collected throughout) and a dummy variable for before/after the law restricting cephalosporins went into effect (assuming that it was in effect starting April 5, 2012), and creating a dateLaw object that is just the date of the law:

psec.CHL$date <- strptime(paste("15", psec.CHL$ym), "%d %b %Y")
psec.CHL$date <- as.Date(psec.CHL$date)
psec.CHL$dummy <- ifelse(psec.CHL$date<"2012-04-05",0,1)
sec.CHL$date <- strptime(paste("15", sec.CHL$ym), "%d %b %Y")
sec.CHL$date <- as.Date(sec.CHL$date)
sec.CHL$dummy <- ifelse(sec.CHL$date<"2012-04-05",0,1)
dateLaw <- as.Date("2012-04-05")

Creating a GLM object that predicts the change in slope for a linear model after the date the law went into effect, with no covariates:

glm.sec.CHL <- glm(sec.CHL$freq ~ sec.CHL$date + pmax(sec.CHL$date - dateLaw, 0) + sec.CHL$MEAT_TYPE)
summary(glm.sec.CHL)
## 
## Call:
## glm(formula = sec.CHL$freq ~ sec.CHL$date + pmax(sec.CHL$date - 
##     dateLaw, 0) + sec.CHL$MEAT_TYPE)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.31186  -0.01994   0.00677   0.02888   0.06988  
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      9.506e-01  2.574e-02  36.928  < 2e-16 ***
## sec.CHL$date                     1.852e-06  1.859e-06   0.996 0.319606    
## pmax(sec.CHL$date - dateLaw, 0)  1.256e-05  9.047e-06   1.388 0.165650    
## sec.CHL$MEAT_TYPEGB             -1.826e-02  5.508e-03  -3.314 0.000971 ***
## sec.CHL$MEAT_TYPEGT             -1.903e-02  5.459e-03  -3.486 0.000524 ***
## sec.CHL$MEAT_TYPEPC             -4.227e-02  5.478e-03  -7.718 4.61e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.002421326)
## 
##     Null deviance: 1.7082  on 641  degrees of freedom
## Residual deviance: 1.5400  on 636  degrees of freedom
## AIC: -2037.2
## 
## Number of Fisher Scoring iterations: 2
confint(glm.sec.CHL)
## Waiting for profiling to be done...
##                                         2.5 %        97.5 %
## (Intercept)                      9.001607e-01  1.001069e+00
## sec.CHL$date                    -1.791782e-06  5.494802e-06
## pmax(sec.CHL$date - dateLaw, 0) -5.175195e-06  3.028718e-05
## sec.CHL$MEAT_TYPEGB             -2.905254e-02 -7.460573e-03
## sec.CHL$MEAT_TYPEGT             -2.973100e-02 -8.331205e-03
## sec.CHL$MEAT_TYPEPC             -5.300879e-02 -3.153732e-02

Issues with this model: 1. it treats the effect of the intervention as something that increases linearly starting from the date the law was passed. This is probably fine for the first few years, but would not be as reasonable for a long-term evaluation. 2. it assumes linear relationships between the predictors and outcomes. 3. it assumes observations are independent

Making a graph to display on the poster:

glm.psec.CHL <- glm(psec.CHL$freq ~ psec.CHL$date + pmax(psec.CHL$date - dateLaw, 0))
summary(glm.psec.CHL)
## 
## Call:
## glm(formula = psec.CHL$freq ~ psec.CHL$date + pmax(psec.CHL$date - 
##     dateLaw, 0))
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.086964  -0.010504   0.001577   0.016902   0.036381  
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      9.376e-01  2.326e-02  40.314   <2e-16 ***
## psec.CHL$date                    1.704e-06  1.692e-06   1.007    0.315    
## pmax(psec.CHL$date - dateLaw, 0) 1.135e-05  7.751e-06   1.464    0.145    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.0005053225)
## 
##     Null deviance: 0.085811  on 162  degrees of freedom
## Residual deviance: 0.080852  on 160  degrees of freedom
## AIC: -769.68
## 
## Number of Fisher Scoring iterations: 2
confint(glm.psec.CHL)
## Waiting for profiling to be done...
##                                          2.5 %       97.5 %
## (Intercept)                       8.920339e-01 9.832027e-01
## psec.CHL$date                    -1.611892e-06 5.019146e-06
## pmax(psec.CHL$date - dateLaw, 0) -3.843645e-06 2.653898e-05
newdata <- data.frame(psec.CHL$date)
newdata$freq <- predict(glm.psec.CHL, newdata = newdata)
plot.glm.psec.CHL <- plot(psec.CHL$freq ~ psec.CHL$date, type = "p", main = "Frequency of Susceptibility to Cephalosporin in E. coli and Salmonella in Retail Meat Before and After 2012 Restriction", xlab = "Time", ylab = "Proportion of Isolates susceptible") + lines(newdata$freq ~ newdata$psec.CHL.date, col = "blue") + abline(v=as.Date("2012-04-05"),col=3,lty=3)

Repeating with added split:

Generating estimates:

before = sec.CHL$date<"2012-04-05"
sec.CHL.bef <- subset(sec.CHL, (before))
sec.CHL.aft <- subset(sec.CHL, (!before))
glm.sec.CHL.MT.bef <- glm(sec.CHL.bef$freq ~ sec.CHL.bef$date + sec.CHL.bef$MEAT_TYPE)
glm.sec.CHL.MT.aft <- glm(sec.CHL.aft$freq ~ sec.CHL.aft$date + sec.CHL.aft$MEAT_TYPE)
summary(glm.sec.CHL.MT.bef)
## 
## Call:
## glm(formula = sec.CHL.bef$freq ~ sec.CHL.bef$date + sec.CHL.bef$MEAT_TYPE)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.306740  -0.019206   0.009062   0.033268   0.076076  
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              9.163e-01  2.954e-02  31.022  < 2e-16 ***
## sec.CHL.bef$date         4.375e-06  2.152e-06   2.033   0.0426 *  
## sec.CHL.bef$MEAT_TYPEGB -1.533e-02  6.577e-03  -2.331   0.0202 *  
## sec.CHL.bef$MEAT_TYPEGT -1.385e-02  6.577e-03  -2.106   0.0357 *  
## sec.CHL.bef$MEAT_TYPEPC -4.388e-02  6.577e-03  -6.671  6.9e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.002660061)
## 
##     Null deviance: 1.4316  on 491  degrees of freedom
## Residual deviance: 1.2954  on 487  degrees of freedom
## AIC: -1514.1
## 
## Number of Fisher Scoring iterations: 2
confint(glm.sec.CHL.MT.bef)
## Waiting for profiling to be done...
##                                 2.5 %        97.5 %
## (Intercept)              8.584508e-01  9.742418e-01
## sec.CHL.bef$date         1.579291e-07  8.591946e-06
## sec.CHL.bef$MEAT_TYPEGB -2.821974e-02 -2.439540e-03
## sec.CHL.bef$MEAT_TYPEGT -2.673975e-02 -9.595479e-04
## sec.CHL.bef$MEAT_TYPEPC -5.676614e-02 -3.098594e-02
summary(glm.sec.CHL.MT.aft)
## 
## Call:
## glm(formula = sec.CHL.aft$freq ~ sec.CHL.aft$date + sec.CHL.aft$MEAT_TYPE)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.142515  -0.020195   0.004744   0.026727   0.066846  
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              3.406e-01  1.495e-01   2.278  0.02418 *  
## sec.CHL.aft$date         4.037e-05  9.315e-06   4.333 2.73e-05 ***
## sec.CHL.aft$MEAT_TYPEGB -2.649e-02  8.917e-03  -2.971  0.00348 ** 
## sec.CHL.aft$MEAT_TYPEGT -3.497e-02  8.566e-03  -4.083 7.32e-05 ***
## sec.CHL.aft$MEAT_TYPEPC -3.587e-02  8.689e-03  -4.128 6.14e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.001448548)
## 
##     Null deviance: 0.27328  on 149  degrees of freedom
## Residual deviance: 0.21004  on 145  degrees of freedom
## AIC: -547.98
## 
## Number of Fisher Scoring iterations: 2
confint(glm.sec.CHL.MT.aft)
## Waiting for profiling to be done...
##                                 2.5 %        97.5 %
## (Intercept)              4.758455e-02  6.337141e-01
## sec.CHL.aft$date         2.210832e-05  5.862381e-05
## sec.CHL.aft$MEAT_TYPEGB -4.396908e-02 -9.015140e-03
## sec.CHL.aft$MEAT_TYPEGT -5.176283e-02 -1.818481e-02
## sec.CHL.aft$MEAT_TYPEPC -5.290509e-02 -1.884290e-02
newdata.CHL.MT.bef <- data.frame(sec.CHL.bef$date)
newdata.CHL.MT.bef$freq <- predict(glm.sec.CHL.MT.bef, newdata = newdata.CHL.MT.bef)
newdata.CHL.MT.aft <- data.frame(sec.CHL.aft$date)
newdata.CHL.MT.aft$freq <- predict(glm.sec.CHL.MT.aft, newdata = newdata.CHL.MT.aft)
before.p = psec.CHL$date<"2012-04-05"
psec.CHL.bef <- subset(psec.CHL, (before.p))
psec.CHL.aft <- subset(psec.CHL, (!before.p))
glm.psec.CHL.bef <- glm(psec.CHL.bef$freq ~ psec.CHL.bef$date)
glm.psec.CHL.aft <- glm(psec.CHL.aft$freq ~ psec.CHL.aft$date)
newdata.CHL.bef <- data.frame(psec.CHL.bef$date)
newdata.CHL.bef$freq <- predict(glm.psec.CHL.bef, newdata = newdata.CHL.bef)
newdata.CHL.aft <- data.frame(psec.CHL.aft$date)
newdata.CHL.aft$freq <- predict(glm.psec.CHL.aft, newdata = newdata.CHL.aft)

plot.glm.sec.CHL.ALL <- plot(psec.CHL$freq ~ psec.CHL$date, type = "p", main = "Monthly Prevalence of Susceptibility to Chloramphenicol in E. coli and Salmonella in Retail Meat Before and After 2012 Restriction", xlab = "Date Isolates Collected, 2002-2015", ylab = "Proportion of Isolates susceptible", pch=20, ylim=c(.5,1))  + lines(newdata.CHL.bef$freq ~ newdata.CHL.bef$psec.CHL.bef, col = "red") + lines(newdata.CHL.aft$freq ~ newdata.CHL.aft$psec.CHL.aft, col = "red") + lines(newdata$freq ~ newdata$psec.CHL.date, col = "blue") + abline(v=as.Date("2012-04-05"),col=3,lty=3)

Below are some additional attempts to visualize ceftriaxone over time as a time series object:

If the prevalence of resistance over time is made into a time series object and smoothed to average over a year, it doesn’t look like the change is at 2012 - it seems to plateau a few years earlier. This could be explained by when the law was first publicized - Dr. Call suggested this.

#creating time series object:
sec.AXO.ts <- ts(sec.AXO$freq, frequency = 12, start = c(2002, 1))
plot(sec.AXO.ts)

#smoothing, I know 12 is a lot but I was thinking it creates yearly moving averages:
sec.AXO.ts.SMA12 <- SMA(sec.AXO.ts,n=12)
plot.ts(sec.AXO.ts.SMA12, main = "Prevalence of AXO resistance in retail meat - smoothed", ylab = "Proportion of isolates resistant")

trying out ARIMA model. I don’t know if this is helpful for evaluating whether there was a change after the 2012 intervention but I am trying to use this to learn more about building an appropriate model for the data.

plot.ts(sec.AXO.ts)

sec.AXO.ts.diff.1 <- diff(sec.AXO.ts, differences=1)
plot.ts(sec.AXO.ts.diff.1)

acf(sec.AXO.ts.diff.1, lag.max=20)

I don’t know what’s going on with the x-axis (shouldn’t that be 1-20?) but if I am interpreting this right it looks like there is a significant correlation between residuals at lag 1, which I would think would make sense.

pacf(sec.AXO.ts.diff.1, lag.max=20) 

auto.arima(sec.AXO.ts.diff.1)
## Series: sec.AXO.ts.diff.1 
## ARIMA(1,0,1)(0,0,2)[12] with zero mean     
## 
## Coefficients:
##           ar1      ma1    sma1    sma2
##       -0.2279  -0.9766  0.2840  0.2664
## s.e.   0.0411   0.0078  0.0446  0.0360
## 
## sigma^2 estimated as 0.003303:  log likelihood=920.82
## AIC=-1831.65   AICc=-1831.55   BIC=-1809.33

It looks like the partial autocorrelation (which I do not fully understand how it is determined / how it is distinct from autocorrelation) is significant until lag 4, although it gets lower and lower from lag 1. Does that mean I can use an ARMA(0,1)? The auto.arima function says that the best fit is ARIMA(0,0,1)

sec.AXO.arima.diff.1 <- arima(sec.AXO.ts.diff.1, order=c(0,0,1))
sec.AXO.arima.forecasts.diff.1 <- forecast.Arima(sec.AXO.arima.diff.1, order=c(0,0,1))
## Warning in forecast.Arima(sec.AXO.arima.diff.1, order = c(0, 0, 1)): The
## non-existent order arguments will be ignored.
acf(sec.AXO.arima.forecasts.diff.1$residuals, lag.max=20)

sec.AXO.arima <- arima(sec.AXO.ts, order=c(0,0,1))
sec.AXO.arima.forecasts <- forecast.Arima(sec.AXO.arima, order=c(0,0,1))
## Warning in forecast.Arima(sec.AXO.arima, order = c(0, 0, 1)): The non-
## existent order arguments will be ignored.
acf(sec.AXO.arima.forecasts$residuals, lag.max=20)

I dont’ understand how to introduce the 2012 law into a model that uses a time series object, or how to introduce meat type as a covariate. It might also work to somehow reproduce this example from MacLeod and Yu’s paper:

“A brief example of a medical intervention analysis carried out using arima() will now be discussed. In a medical time series of monthly average creatinine clearances, a step intervention analysis model with a multiplicative seasonal ARIMA(0, 1, 1) (1, 0, 0)12 error term was fit. The intervention effect was found to be significant at 1%. To illustrate this finding, Figure 8 compares the forecasts before and after the intervention date. The forecasts are from a model fit to the pre-intervention series. The plot visually confirms the decrease in creatinine clearances after the intervention.” at http://www.stats.uwo.ca/faculty/aim/tsar/tsar.pdf

Medical intervention

Medical intervention

However, it would be helpful to have some kind of measure of the magnitude and uncertainty of the change, or the magnitude and uncertainty of the difference between the predicted and actual values, similar to the p-value for the coefficient of the predictor of the variable representing the law. I don’t know what that would be for a time series model.