Current Progress - GLM Models
Anika Larson
12/6/2016
Loading libraries:
library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(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, unionlibrary(ggplot2)
library(TTR)
library(forecast)## Loading required package: timeDate## This is forecast 7.3NARMS_retailmeat_public <- read.csv("narms_public.csv")
d <- NARMS_retailmeat_publicSubsetting the data to contain only the genuses Salmonella and E coli, and renaming psec.
sec <- subset(d,(d$GENUS=="S"|d$GENUS=="EC"))Creating a year/month variable in sec:
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.98876404Subsetting psec.AXO to contain only proportions of isolates that were resistant to ceftriaxone each month (which groups Intermediate and Susceptible together by default):
sec.AXO <- subset(psec.AXO, AXO.SIR=="R")
head(sec.AXO)## Source: local data frame [6 x 4]
## Groups: ym [6]
##
## ym AXO.SIR n freq
## <S3: yearmon> <fctr> <int> <dbl>
## 1 Jan 2002 R 1 0.01063830
## 2 Feb 2002 R 6 0.06122449
## 3 Mar 2002 R 1 0.01123596
## 4 Apr 2002 R 1 0.01086957
## 5 May 2002 R 9 0.08571429
## 6 Jun 2002 R 2 0.02061856Plotting the prevalence of resistance over time:
sec.axo.rprev <- plot(sec.AXO$ym, sec.AXO$freq, main = "Prevalence of ceftriaxone resistance in isolates of Salmonella and E. coli from NARMS retail meat", xlab = "Time (months)", ylab = "Proportion of isolates resistant to ceftriaxone")
Note: if the prevalence of resistance over time is made into a time series object and smoothed, it doesn’t look like the change is at 2012 - it seems to plateau a few years earlier. Not sure what to do with this information.
#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")
Most recent attempt to create a generalized linear model, just evaluating the effect of the intervention with no covariates assuming that the trendlines before and after the intervention will have to intersect at the breakpoint:
Creating a date field, where each timepoint is treated as the 15th of the month (because isolates were received throughout) and a dummy variable for before/after the law restricting cephalosporins went into effect (assuming that it was in effect starting May 1, 2012), and creating a dateLaw object that is just the date of the law:
sec.AXO$by_date <- strptime(paste("15", sec.AXO$ym), "%d %b %Y")
sec.AXO$date <- as.Date(sec.AXO$by_date)
sec.AXO$dummy <- ifelse(sec.AXO$by_date<"2012-05-01",0,1)
dateLaw <- as.Date("2012-01-01")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))
summary(glm.sec.AXO)##
## Call:
## glm(formula = sec.AXO$freq ~ sec.AXO$date + pmax(sec.AXO$date -
## dateLaw, 0))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.061953 -0.022232 -0.007517 0.016976 0.117512
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.178e-01 3.300e-02 -6.601 5.95e-10 ***
## sec.AXO$date 2.130e-05 2.408e-06 8.846 1.78e-15 ***
## pmax(sec.AXO$date - dateLaw, 0) -7.108e-05 1.067e-05 -6.661 4.33e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.0009542413)
##
## Null deviance: 0.22500 on 159 degrees of freedom
## Residual deviance: 0.14982 on 157 degrees of freedom
## AIC: -653.7
##
## Number of Fisher Scoring iterations: 2Issues 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
Displaying the predictions of the GLM model against the actual prevalence data:
newdata <- data.frame(sec.AXO$date)
newdata$freq <- predict(glm.sec.AXO, newdata = newdata)
plot.glm.sec.AXO <- plot(sec.AXO$freq ~ sec.AXO$date, type = "l", main = "Frequency of Resistance to Cephalosporin Resistance in E. coli and Salmonella in Retail Meat Before and After 2012 Restriction", xlab = "Time", ylab = "Proportion of Isolates Resistant") + lines(newdata$freq ~ newdata$sec.AXO.date)
Merging the frequency data into the original dataset to include covariates:
sec$ym <- as.yearmon(paste(sec$Year, sec$Month, sep = "-"))
range(sec$ym)## [1] "Jan 2002" "Jul 2015"secm.AXO <- merge(sec,sec.AXO,by="ym")
merged.plot <- plot(secm.AXO$ym, secm.AXO$freq, main = "Prevalence of ceftriaxone resistance in isolates of Salmonella and E. coli from NARMS retail meat", xlab = "Time (months)", ylab = "Proportion of isolates resistant to ceftriaxone")
glm.secm.AXO <- glm(secm.AXO$freq ~ secm.AXO$date + pmax(secm.AXO$date - dateLaw, 0))
summary(glm.secm.AXO)##
## Call:
## glm(formula = secm.AXO$freq ~ secm.AXO$date + pmax(secm.AXO$date -
## dateLaw, 0))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.063884 -0.020273 -0.006654 0.018508 0.116165
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.309e-01 2.921e-03 -79.05 <2e-16
## secm.AXO$date 2.231e-05 2.132e-07 104.65 <2e-16
## pmax(secm.AXO$date - dateLaw, 0) -8.542e-05 1.044e-06 -81.79 <2e-16
##
## (Intercept) ***
## secm.AXO$date ***
## pmax(secm.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.0008875864)
##
## Null deviance: 27.248 on 19516 degrees of freedom
## Residual deviance: 17.320 on 19514 degrees of freedom
## AIC: -81754
##
## Number of Fisher Scoring iterations: 2newdatam <- data.frame(secm.AXO$date)
newdatam$freq <- predict(glm.secm.AXO, newdatam = newdatam)
plot.glm.secm.AXO <- plot(secm.AXO$freq ~ secm.AXO$date, type = "l", main = "Frequency of Resistance to Cephalosporin Resistance in E. coli and Salmonella in Retail Meat Before and After 2012 Restriction", xlab = "Time", ylab = "Proportion of Isolates Resistant") + lines(newdatam$freq ~ newdatam$secm.AXO.date)
However, this method produces a different set of estimations for the slopes at the line change, probably because it takes into account the different numbers of observations at each timepoint, which now each have the same value. It’s weighting those values unequally depending on the number of observations for the month. Is that appropriate?
Introducing serotype as a covariate:
glm.secm.AXO.sero <- glm(secm.AXO$freq ~ secm.AXO$date + pmax(secm.AXO$date - dateLaw, 0) + secm.AXO$SEROTYPE)
summary(glm.secm.AXO.sero)##
## Call:
## glm(formula = secm.AXO$freq ~ secm.AXO$date + pmax(secm.AXO$date -
## dateLaw, 0) + secm.AXO$SEROTYPE)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.080078 -0.021377 -0.006287 0.017516 0.121660
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -2.289e-01 2.947e-03 -77.686
## secm.AXO$date 2.212e-05 2.157e-07 102.553
## pmax(secm.AXO$date - dateLaw, 0) -8.465e-05 1.063e-06 -79.639
## secm.AXO$SEROTYPEAdelaide 1.414e-02 1.715e-02 0.824
## secm.AXO$SEROTYPEAgona 3.106e-04 3.251e-03 0.096
## secm.AXO$SEROTYPEAlachua 1.626e-02 9.907e-03 1.641
## secm.AXO$SEROTYPEAlbany -1.465e-04 4.592e-03 -0.032
## secm.AXO$SEROTYPEAlbert 6.348e-02 2.971e-02 2.137
## secm.AXO$SEROTYPEAnatum -6.108e-03 4.889e-03 -1.249
## secm.AXO$SEROTYPEBareilly -1.491e-02 1.213e-02 -1.229
## secm.AXO$SEROTYPEBerta 2.912e-03 4.251e-03 0.685
## secm.AXO$SEROTYPEBlockley -1.235e-02 1.329e-02 -0.930
## secm.AXO$SEROTYPEBraenderup 4.988e-03 6.068e-03 0.822
## secm.AXO$SEROTYPEBrandenburg 9.555e-03 5.027e-03 1.901
## secm.AXO$SEROTYPEBredeney 6.534e-04 7.674e-03 0.085
## secm.AXO$SEROTYPECerro 1.830e-03 1.329e-02 0.138
## secm.AXO$SEROTYPEChester -4.293e-03 1.329e-02 -0.323
## secm.AXO$SEROTYPEDerby 1.557e-03 4.168e-03 0.374
## secm.AXO$SEROTYPEDublin 1.562e-02 6.339e-03 2.463
## secm.AXO$SEROTYPEEnteritidis -1.666e-03 1.830e-03 -0.911
## secm.AXO$SEROTYPEGaminara 3.132e-03 2.101e-02 0.149
## secm.AXO$SEROTYPEGive -1.920e-02 2.971e-02 -0.646
## secm.AXO$SEROTYPEHaardt 2.803e-02 1.329e-02 2.109
## secm.AXO$SEROTYPEHadar -4.788e-03 1.642e-03 -2.916
## secm.AXO$SEROTYPEHeidelberg 2.011e-03 1.167e-03 1.724
## secm.AXO$SEROTYPEHindmarsh -1.580e-02 2.971e-02 -0.532
## secm.AXO$SEROTYPEHvittingfoss -2.110e-02 2.971e-02 -0.710
## secm.AXO$SEROTYPEI 4,12:d: - 1.994e-02 2.101e-02 0.949
## secm.AXO$SEROTYPEI 4,12:d:- 1.220e-03 5.832e-03 0.209
## secm.AXO$SEROTYPEI 4,12:i: - -1.363e-02 2.101e-02 -0.649
## secm.AXO$SEROTYPEI 4,12:non-motile 1.886e-02 2.971e-02 0.635
## secm.AXO$SEROTYPEI 4,12:r:- 1.817e-02 1.329e-02 1.367
## secm.AXO$SEROTYPEI 4,5,12:-:1,2 -9.225e-03 2.970e-02 -0.311
## secm.AXO$SEROTYPEI 4,5,12:d:- -1.957e-02 2.101e-02 -0.931
## secm.AXO$SEROTYPEI 4,[5],12:i:- 4.684e-03 3.534e-03 1.326
## secm.AXO$SEROTYPEI 4,5,12: i : - -1.930e-03 2.101e-02 -0.092
## secm.AXO$SEROTYPEI 4,5,12:r:- 5.232e-03 6.822e-03 0.767
## secm.AXO$SEROTYPEI 6,7:-1,5 -2.691e-02 2.971e-02 -0.906
## secm.AXO$SEROTYPEI 6,7:k:- -2.072e-03 1.715e-02 -0.121
## secm.AXO$SEROTYPEI 6,7:l,w:- 6.559e-03 2.971e-02 0.221
## secm.AXO$SEROTYPEI 8,20:-:z6 -3.692e-02 2.101e-02 -1.757
## secm.AXO$SEROTYPEIIIa 18:z4,z23:- 4.416e-03 2.762e-03 1.599
## secm.AXO$SEROTYPEIIIa 18:z4,z32:- -9.364e-03 8.245e-03 -1.136
## secm.AXO$SEROTYPEIIIa 35:z4,z23:- 1.544e-02 2.101e-02 0.735
## secm.AXO$SEROTYPEIIIa -:z4,z23:- -9.491e-03 2.970e-02 -0.320
## secm.AXO$SEROTYPEInfantis 1.558e-03 3.465e-03 0.450
## secm.AXO$SEROTYPEJaviana 3.132e-03 2.971e-02 0.105
## secm.AXO$SEROTYPEJohannesburg -1.249e-03 7.947e-03 -0.157
## secm.AXO$SEROTYPEKentucky 2.293e-03 1.305e-03 1.757
## secm.AXO$SEROTYPELitchfield 7.064e-03 9.905e-03 0.713
## secm.AXO$SEROTYPELiverpool 2.805e-02 2.101e-02 1.335
## secm.AXO$SEROTYPELivingstone -1.401e-02 1.329e-02 -1.054
## secm.AXO$SEROTYPELondon 2.048e-02 2.101e-02 0.975
## secm.AXO$SEROTYPEMbandaka -7.835e-03 3.811e-03 -2.056
## secm.AXO$SEROTYPEMeleagridis 5.198e-03 2.101e-02 0.247
## secm.AXO$SEROTYPEMinnesota -1.938e-02 1.123e-02 -1.726
## secm.AXO$SEROTYPEMontevideo -2.964e-05 3.231e-03 -0.009
## secm.AXO$SEROTYPEMuenchen -1.350e-03 4.348e-03 -0.310
## secm.AXO$SEROTYPEMuenster -3.283e-04 5.430e-03 -0.060
## secm.AXO$SEROTYPENewport 5.414e-04 4.435e-03 0.122
## secm.AXO$SEROTYPENorwich -3.205e-03 1.123e-02 -0.285
## secm.AXO$SEROTYPEOhio 6.647e-03 9.397e-03 0.707
## secm.AXO$SEROTYPEOranienburg -1.013e-02 1.485e-02 -0.682
## secm.AXO$SEROTYPEOrion 3.048e-02 2.101e-02 1.451
## secm.AXO$SEROTYPEOuakam -5.954e-03 2.101e-02 -0.283
## secm.AXO$SEROTYPEPanama -9.390e-03 1.485e-02 -0.632
## secm.AXO$SEROTYPEReading -1.528e-03 2.863e-03 -0.534
## secm.AXO$SEROTYPERough/ Nonmotile isolate 6.991e-02 2.971e-02 2.353
## secm.AXO$SEROTYPERough/ Nonmotile Isolate 3.006e-02 8.961e-03 3.355
## secm.AXO$SEROTYPESaintpaul 6.066e-03 1.516e-03 4.003
## secm.AXO$SEROTYPESandiego 1.372e-02 2.971e-02 0.462
## secm.AXO$SEROTYPESchwarzengrund -3.240e-04 2.723e-03 -0.119
## secm.AXO$SEROTYPESenftenberg -1.340e-04 3.043e-03 -0.044
## secm.AXO$SEROTYPESIIIa 40:z4,z23:- -2.542e-02 2.971e-02 -0.856
## secm.AXO$SEROTYPESinstorf 4.374e-02 2.970e-02 1.473
## secm.AXO$SEROTYPETennessee -2.174e-02 1.329e-02 -1.636
## secm.AXO$SEROTYPEThompson -1.927e-02 5.833e-03 -3.304
## secm.AXO$SEROTYPETyphimurium 9.419e-03 1.128e-03 8.346
## secm.AXO$SEROTYPETyphimurium var 5- 7.056e-03 4.799e-03 1.470
## secm.AXO$SEROTYPEUganda -1.257e-02 1.051e-02 -1.196
## secm.AXO$SEROTYPEUrbana -7.473e-03 2.971e-02 -0.252
## secm.AXO$SEROTYPEWorthington -4.815e-03 8.961e-03 -0.537
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## secm.AXO$date < 2e-16 ***
## pmax(secm.AXO$date - dateLaw, 0) < 2e-16 ***
## secm.AXO$SEROTYPEAdelaide 0.409757
## secm.AXO$SEROTYPEAgona 0.923887
## secm.AXO$SEROTYPEAlachua 0.100816
## secm.AXO$SEROTYPEAlbany 0.974553
## secm.AXO$SEROTYPEAlbert 0.032618 *
## secm.AXO$SEROTYPEAnatum 0.211601
## secm.AXO$SEROTYPEBareilly 0.219045
## secm.AXO$SEROTYPEBerta 0.493329
## secm.AXO$SEROTYPEBlockley 0.352499
## secm.AXO$SEROTYPEBraenderup 0.411067
## secm.AXO$SEROTYPEBrandenburg 0.057334 .
## secm.AXO$SEROTYPEBredeney 0.932153
## secm.AXO$SEROTYPECerro 0.890469
## secm.AXO$SEROTYPEChester 0.746615
## secm.AXO$SEROTYPEDerby 0.708775
## secm.AXO$SEROTYPEDublin 0.013770 *
## secm.AXO$SEROTYPEEnteritidis 0.362512
## secm.AXO$SEROTYPEGaminara 0.881469
## secm.AXO$SEROTYPEGive 0.518097
## secm.AXO$SEROTYPEHaardt 0.034932 *
## secm.AXO$SEROTYPEHadar 0.003548 **
## secm.AXO$SEROTYPEHeidelberg 0.084806 .
## secm.AXO$SEROTYPEHindmarsh 0.594839
## secm.AXO$SEROTYPEHvittingfoss 0.477568
## secm.AXO$SEROTYPEI 4,12:d: - 0.342742
## secm.AXO$SEROTYPEI 4,12:d:- 0.834284
## secm.AXO$SEROTYPEI 4,12:i: - 0.516590
## secm.AXO$SEROTYPEI 4,12:non-motile 0.525547
## secm.AXO$SEROTYPEI 4,12:r:- 0.171607
## secm.AXO$SEROTYPEI 4,5,12:-:1,2 0.756132
## secm.AXO$SEROTYPEI 4,5,12:d:- 0.351624
## secm.AXO$SEROTYPEI 4,[5],12:i:- 0.184981
## secm.AXO$SEROTYPEI 4,5,12: i : - 0.926805
## secm.AXO$SEROTYPEI 4,5,12:r:- 0.443113
## secm.AXO$SEROTYPEI 6,7:-1,5 0.364911
## secm.AXO$SEROTYPEI 6,7:k:- 0.903854
## secm.AXO$SEROTYPEI 6,7:l,w:- 0.825239
## secm.AXO$SEROTYPEI 8,20:-:z6 0.078858 .
## secm.AXO$SEROTYPEIIIa 18:z4,z23:- 0.109819
## secm.AXO$SEROTYPEIIIa 18:z4,z32:- 0.256076
## secm.AXO$SEROTYPEIIIa 35:z4,z23:- 0.462391
## secm.AXO$SEROTYPEIIIa -:z4,z23:- 0.749338
## secm.AXO$SEROTYPEInfantis 0.653032
## secm.AXO$SEROTYPEJaviana 0.916024
## secm.AXO$SEROTYPEJohannesburg 0.875110
## secm.AXO$SEROTYPEKentucky 0.078913 .
## secm.AXO$SEROTYPELitchfield 0.475735
## secm.AXO$SEROTYPELiverpool 0.181787
## secm.AXO$SEROTYPELivingstone 0.291999
## secm.AXO$SEROTYPELondon 0.329531
## secm.AXO$SEROTYPEMbandaka 0.039799 *
## secm.AXO$SEROTYPEMeleagridis 0.804537
## secm.AXO$SEROTYPEMinnesota 0.084356 .
## secm.AXO$SEROTYPEMontevideo 0.992680
## secm.AXO$SEROTYPEMuenchen 0.756256
## secm.AXO$SEROTYPEMuenster 0.951793
## secm.AXO$SEROTYPENewport 0.902828
## secm.AXO$SEROTYPENorwich 0.775359
## secm.AXO$SEROTYPEOhio 0.479355
## secm.AXO$SEROTYPEOranienburg 0.495183
## secm.AXO$SEROTYPEOrion 0.146746
## secm.AXO$SEROTYPEOuakam 0.776833
## secm.AXO$SEROTYPEPanama 0.527302
## secm.AXO$SEROTYPEReading 0.593401
## secm.AXO$SEROTYPERough/ Nonmotile isolate 0.018609 *
## secm.AXO$SEROTYPERough/ Nonmotile Isolate 0.000796 ***
## secm.AXO$SEROTYPESaintpaul 6.29e-05 ***
## secm.AXO$SEROTYPESandiego 0.644069
## secm.AXO$SEROTYPESchwarzengrund 0.905263
## secm.AXO$SEROTYPESenftenberg 0.964864
## secm.AXO$SEROTYPESIIIa 40:z4,z23:- 0.392166
## secm.AXO$SEROTYPESinstorf 0.140893
## secm.AXO$SEROTYPETennessee 0.101808
## secm.AXO$SEROTYPEThompson 0.000955 ***
## secm.AXO$SEROTYPETyphimurium < 2e-16 ***
## secm.AXO$SEROTYPETyphimurium var 5- 0.141518
## secm.AXO$SEROTYPEUganda 0.231571
## secm.AXO$SEROTYPEUrbana 0.801363
## secm.AXO$SEROTYPEWorthington 0.591079
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.000882295)
##
## Null deviance: 27.248 on 19516 degrees of freedom
## Residual deviance: 17.147 on 19435 degrees of freedom
## AIC: -81792
##
## Number of Fisher Scoring iterations: 2newdatam.sero <- data.frame(secm.AXO$date)
newdatam.sero$freq <- predict(glm.secm.AXO, newdatam.sero = newdatam.sero)
plot.glm.secm.AXO.sero <- plot(secm.AXO$freq ~ secm.AXO$date, type = "l", main = "Frequency of Resistance to Cephalosporin Resistance in E. coli and Salmonella in Retail Meat Before and After 2012 Restriction", xlab = "Time", ylab = "Proportion of Isolates Resistant") + lines(newdatam.sero$freq ~ newdatam.sero$secm.AXO.date)
Interesting. Did not affect significance of the intervention variable. Is there a way to quantify the total impact of serotype?
Okay trying meat type THEN DONE:
glm.secm.AXO.spec <- glm(secm.AXO$freq ~ secm.AXO$date + pmax(secm.AXO$date - dateLaw, 0) + secm.AXO$MEAT_TYPE)
summary(glm.secm.AXO.spec)##
## Call:
## glm(formula = secm.AXO$freq ~ secm.AXO$date + pmax(secm.AXO$date -
## dateLaw, 0) + secm.AXO$MEAT_TYPE)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.064699 -0.020972 -0.006149 0.017944 0.117021
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.295e-01 2.954e-03 -77.707 <2e-16
## secm.AXO$date 2.227e-05 2.136e-07 104.268 <2e-16
## pmax(secm.AXO$date - dateLaw, 0) -8.532e-05 1.045e-06 -81.625 <2e-16
## secm.AXO$MEAT_TYPEGB -1.564e-03 6.135e-04 -2.550 0.0108
## secm.AXO$MEAT_TYPEGT -1.013e-03 5.169e-04 -1.960 0.0500
## secm.AXO$MEAT_TYPEPC -1.693e-03 6.922e-04 -2.446 0.0144
##
## (Intercept) ***
## secm.AXO$date ***
## pmax(secm.AXO$date - dateLaw, 0) ***
## secm.AXO$MEAT_TYPEGB *
## secm.AXO$MEAT_TYPEGT *
## secm.AXO$MEAT_TYPEPC *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.0008872784)
##
## Null deviance: 27.248 on 19516 degrees of freedom
## Residual deviance: 17.312 on 19511 degrees of freedom
## AIC: -81758
##
## Number of Fisher Scoring iterations: 2newdatam.spec <- data.frame(secm.AXO$date)
newdatam.spec$freq <- predict(glm.secm.AXO.spec, newdatam.spec = newdatam.spec)
plot.glm.secm.AXO.spec <- plot(secm.AXO$freq ~ secm.AXO$date, type = "l", main = "Frequency of Resistance to Cephalosporin Resistance in E. coli and Salmonella in Retail Meat Before and After 2012 Restriction", xlab = "Time", ylab = "Proportion of Isolates Resistant") + lines(newdatam.spec$freq ~ newdatam.spec$secm.AXO.date, col = "blue")
Not clear where the issues with the trendlines are coming from in this model.
Repeating this process using ONLY Salmonella isolates:
sal <- subset(d,d$GENUS=="S")
sal$ym <- as.yearmon(paste(sal$Year, sal$Month, sep = "-"))
sal.AXO <- sal %>% group_by(ym, AXO.SIR) %>% summarise (n = n()) %>% mutate(freq = n / sum(n))
head(sal.AXO)## Source: local data frame [6 x 4]
## Groups: ym [5]
##
## ym AXO.SIR n freq
## <S3: yearmon> <fctr> <int> <dbl>
## 1 Jan 2002 S 18 1.0
## 2 Feb 2002 S 14 1.0
## 3 Mar 2002 S 8 1.0
## 4 Apr 2002 S 10 1.0
## 5 May 2002 R 2 0.2
## 6 May 2002 S 8 0.8sal.AXO <- subset(sal.AXO, AXO.SIR=="R")
head(sal.AXO)## Source: local data frame [6 x 4]
## Groups: ym [6]
##
## ym AXO.SIR n freq
## <S3: yearmon> <fctr> <int> <dbl>
## 1 May 2002 R 2 0.2000000
## 2 Jun 2002 R 1 0.1111111
## 3 Jul 2002 R 6 0.4000000
## 4 Sep 2002 R 2 0.2500000
## 5 Oct 2002 R 3 0.1500000
## 6 Nov 2002 R 1 0.0500000sal.axo.rprev <- plot(sal.AXO$ym, sal.AXO$freq, main = "Prevalence of ceftriaxone resistance in isolates of Salmonella from NARMS retail meat", xlab = "Time (months)", ylab = "Proportion of isolates resistant to ceftriaxone")
sal.AXO.ts <- ts(sal.AXO$freq, frequency = 12, start = c(2002, 1))
plot(sal.AXO.ts)
sal.AXO.ts.SMA12 <- SMA(sal.AXO.ts,n=12)
plot.ts(sal.AXO.ts.SMA12, main = "Prevalence of AXO resistance in retail meat - smoothed", ylab = "Proportion of isolates resistant")
sal.AXO$by_date <- strptime(paste("15", sal.AXO$ym), "%d %b %Y")
sal.AXO$date <- as.Date(sal.AXO$by_date)
sal.AXO$dummy <- ifelse(sal.AXO$by_date<"2012-05-01",0,1)
dateLaw <- as.Date("2012-01-01")
glm.sal.AXO <- glm(sal.AXO$freq ~ sal.AXO$date + pmax(sal.AXO$date - dateLaw, 0))
summary(glm.sal.AXO)##
## Call:
## glm(formula = sal.AXO$freq ~ sal.AXO$date + pmax(sal.AXO$date -
## dateLaw, 0))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.16405 -0.07290 -0.01692 0.05813 0.33433
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.658e-01 1.095e-01 -1.514 0.13223
## sal.AXO$date 2.463e-05 7.936e-06 3.103 0.00230 **
## pmax(sal.AXO$date - dateLaw, 0) -1.084e-04 3.378e-05 -3.211 0.00163 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.0089281)
##
## Null deviance: 1.4086 on 148 degrees of freedom
## Residual deviance: 1.3035 on 146 degrees of freedom
## AIC: -275.25
##
## Number of Fisher Scoring iterations: 2newdata.sal <- data.frame(sal.AXO$date)
newdata.sal$freq <- predict(glm.sal.AXO, newdata.sal = newdata.sal)
plot.glm.sal.AXO <- plot(sal.AXO$freq ~ sal.AXO$date, type = "l", main = "Frequency of Resistance to Cephalosporin Resistance in Salmonella in Retail Meat Before and After 2012 Restriction", xlab = "Time", ylab = "Proportion of Isolates Resistant") + lines(newdata.sal$freq ~ newdata.sal$sal.AXO.date)