Commit 751714e5 authored by de Freitas Costa, Eduardo's avatar de Freitas Costa, Eduardo
Browse files

Changes after meeting

parent ab216080
......@@ -133,7 +133,7 @@ expo_fa1<-c(0.002319,0,0.002319,0.002319,0.002319,0.002319,0.002319,0.002319,0.0
#Farm exposure is proportional to the number of farms 500 for chicken and 1000 for veal. see (https://www.statista.com/statistics/647188/total-number-of-veal-calves-in-the-netherlands/)
consu_ch<-34.5 #Chicken consumption per capta in grams/portion https://www.wateetnederland.nl/resultaten/voedingsmiddelen/consumptie/vlees (35% Dutch production)
consu_beef<-26.46 #Beef consumption per capta in grams/portion https://www.wateetnederland.nl/resultaten/voedingsmiddelen/consumptie/vlees (35% Dutch production)
consu_beef<-43 #Beef consumption per capta in grams/portion https://www.wateetnederland.nl/resultaten/voedingsmiddelen/consumptie/vlees (35% Dutch production)
consu_veg<-600 #Vegetables consumption per capta in grams/portion https://www.wateetnederland.nl/resultaten/voedingsmiddelen/consumptie/groenten
non_veg<-0.955 #Frequency of non-vegetarians in NL #RIVM 2017
prob_consu<-c(1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1) #Weekly consumption of 1 portion of chicken in NL #RIVM 2017
......@@ -199,7 +199,7 @@ heavy_prop<-c(0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.1,0.01,0.01,0.01,0.01,0.01,0.
0.01,0.001,0.1,0.01,0.01) #proportion of heavy ATM flocks No data, scenario created based on Luiken et al. (2019)
heavy_prop1<-c(0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.1,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,
0.01,0.001,0.1,0.01,0.01) #proportion of heavy ATM veal farms No data, scenario created based on Luiken et al. (2019)
0.01,0.001,0.1,0.01,0.01) #proportion of heavy ATM veal farms No data, scenario created based on guess
......@@ -245,11 +245,10 @@ conc_beef1<-5.68 # concentration o ESBL in feces log10(CFU/gram)
sig_conc_beef1<-0.67 # SD concentration o ESBL in feces log10(CFU/gram)
car_beef_me<-1.190*1000
car_beef_sd<-0.02*1000
car_beef_me<-380*1000 #Caracass weight in grams (Anita)
car_beef_sd<-15*1000 #Caracass sd in grams (guess)
red_carc<-2.51 #reduction in log10cfu/cm2
area_carc<-42000 #Area of a carcass in cm2
#bacterial growth at consumer phase parameters
......
......@@ -873,7 +873,7 @@ trns2$Trns1<-ifelse(trns2$Trns>=2,2,trns2$Trns)
summary(glmer(output~factor(treatment)+breed+factor(Sexe)+factor(Trns1)+(1|DF),
family = binomial(link="logit"),data=trns2,
family = binomial(link="logit"),data=trns1,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)) )->beta1.1)
#Shrink
......@@ -1201,6 +1201,11 @@ lsmeans(beta24.1, pairwise ~ factor(`24.2`), adjust="tukey", data=wk24.2)$lsmean
# ESBL ~ Cummulative regression (DF)
```{r}
#Creating function
back<-function(x){1/(1+exp(-x))}
#ESBL WK=trans~Atb treatments DF (Group)
data0<-subset(veal1[,c(3,4,9,10,50,51)],time %in%c(1))
......@@ -1226,22 +1231,25 @@ kable(table(data0$output,data0$days))
summary(glmer(output~factor(treatment)+breed+factor(Sexe)+days+(1|DairyFarm),
family = binomial(link="logit"),data=data0,control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)) )->model1.1)
kable(Anova(model1.1,type="III"))
model1.1
#Shrink
c<-15/16*pi/sqrt(3)
sr<-1/sqrt((sum((model1@theta)^2)/c^2)+1)
sr<-1/sqrt((sum((model1.1@theta)^2)/c^2)+1)
c(
exp(coef(summary(model1))[6:7,1]*sr-1.96*coef(summary(model1))[6:7,2]*sr),
exp(coef(summary(model1))[6:7,1]*sr),
exp(coef(summary(model1))[6:7,1]*sr+1.96*coef(summary(model1))[6:7,2]*sr)
exp(coef(summary(model1.1))[6,1]*sr-1.96*coef(summary(model1.1))[6,2]*sr),
exp(coef(summary(model1.1))[6,1]*sr),
exp(coef(summary(model1.1))[6,1]*sr+1.96*coef(summary(model1.1))[6,2]*sr)
)
lsmeans(model1.1, pairwise ~ days, adjust="tukey", data=data0)$lsmeans
```
......@@ -1277,10 +1285,6 @@ kable(table(data1$output,data1$days))
summary(glmer(output~factor(treatment)+breed+factor(Sexe)+days+(1|DairyFarm)+(1|VealFarm),
family = binomial(link="logit"),data=data1,control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)) )->model2.1)
kable(Anova(model2.1,type="III"))
#Shrink
c<-15/16*pi/sqrt(3)
......@@ -1288,12 +1292,14 @@ c<-15/16*pi/sqrt(3)
sr<-1/sqrt((sum((model2.1@theta)^2)/c^2)+1)
c(
exp(coef(summary(model2.1))[6:7,1]*sr-1.96*coef(summary(model2.1))[6:7,2]*sr),
exp(coef(summary(model2.1))[6:7,1]*sr),
exp(coef(summary(model2.1))[6:7,1]*sr+1.96*coef(summary(model2.1))[6:7,2]*sr)
exp(coef(summary(model2.1))[6,1]*sr-1.96*coef(summary(model2.1))[6,2]*sr),
exp(coef(summary(model2.1))[6,1]*sr),
exp(coef(summary(model2.1))[6,1]*sr+1.96*coef(summary(model2.1))[6,2]*sr)
)
lsmeans(model2.1, pairwise ~ days, adjust="tukey", data=data1)$lsmeans
```
......@@ -1338,11 +1344,12 @@ c<-15/16*pi/sqrt(3)
sr<-1/sqrt((sum((model3.1@theta)^2)/c^2)+1)
c(
exp(coef(summary(model3.1))[6:7,1]*sr-1.96*coef(summary(model3.1))[6:7,2]*sr),
exp(coef(summary(model3.1))[6:7,1]*sr),
exp(coef(summary(model3.1))[6:7,1]*sr+1.96*coef(summary(model3.1))[6:7,2]*sr)
exp(coef(summary(model3.1))[6,1]*sr-1.96*coef(summary(model3.1))[6,2]*sr),
exp(coef(summary(model3.1))[6,1]*sr),
exp(coef(summary(model3.1))[6,1]*sr+1.96*coef(summary(model3.1))[6,2]*sr)
)
lsmeans(model3.1, pairwise ~ days, adjust="tukey", data=data2)$lsmeans
......@@ -1386,12 +1393,12 @@ c<-15/16*pi/sqrt(3)
sr<-1/sqrt((sum((model4.1@theta)^2)/c^2)+1)
c(
exp(coef(summary(model4.1))[6:8,1]*sr-1.96*coef(summary(model4.1))[6:8,2]*sr),
exp(coef(summary(model4.1))[6:8,1]*sr),
exp(coef(summary(model4.1))[6:8,1]*sr+1.96*coef(summary(model4.1))[6:8,2]*sr)
exp(coef(summary(model4.1))[6,1]*sr-1.96*coef(summary(model4.1))[6,2]*sr),
exp(coef(summary(model4.1))[6,1]*sr),
exp(coef(summary(model4.1))[6,1]*sr+1.96*coef(summary(model4.1))[6,2]*sr)
)
lsmeans(model4.1, pairwise ~ days, adjust="tukey", data=data3)$lsmeans
```
......@@ -1432,12 +1439,12 @@ c<-15/16*pi/sqrt(3)
sr<-1/sqrt((sum((model5.1@theta)^2)/c^2)+1)
c(
exp(coef(summary(model5.1))[6:7,1]*sr-1.96*coef(summary(model5.1))[6:7,2]*sr),
exp(coef(summary(model5.1))[6:7,1]*sr),
exp(coef(summary(model5.1))[6:7,1]*sr+1.96*coef(summary(model5.1))[6:7,2]*sr)
exp(coef(summary(model5.1))[6,1]*sr-1.96*coef(summary(model5.1))[6,2]*sr),
exp(coef(summary(model5.1))[6,1]*sr),
exp(coef(summary(model5.1))[6,1]*sr+1.96*coef(summary(model5.1))[6,2]*sr)
)
lsmeans(model5.1, pairwise ~ days, adjust="tukey", data=data4)$lsmeans
```
......@@ -1472,24 +1479,26 @@ kable(table(data5$output,data5$days))
summary(glmer(output~factor(treatment)+breed+factor(Sexe)+days+(1|DairyFarm)+(1|VealFarm),
family = binomial(link="logit"),data=data5,control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)) )->model6.1)
kable(Anova(model6.1,type="III"))
table(data5$breed)
kable(Anova(model6.1,type="III"))
.#Shrink
#Shrink
c<-15/16*pi/sqrt(3)
sr<-1/sqrt((sum((model6.1@theta)^2)/c^2)+1)
c(
exp(coef(summary(model6.1))[6:7,1]*sr-1.96*coef(summary(model6.1))[6:7,2]*sr),
exp(coef(summary(model6.1))[6:7,1]*sr),
exp(coef(summary(model6.1))[6:7,1]*sr+1.96*coef(summary(model6.1))[6:7,2]*sr)
exp(coef(summary(model6.1))[6,1]*sr-1.96*coef(summary(model6.1))[6,2]*sr),
exp(coef(summary(model6.1))[6,1]*sr),
exp(coef(summary(model6.1))[6,1]*sr+1.96*coef(summary(model6.1))[6,2]*sr)
)
lsmeans(model6.1, pairwise ~ days, adjust="tukey", data=data5)$lsmeans
```
......
......@@ -163,7 +163,7 @@ set.seed(13)
#End step in slaughterhouse: CHILLING
Next_beef5[1]<-max(Next_beef1[1]-(10^2.51*500),0)
Next_beef5[1]<-max(Next_beef1[1]-(10^red_carc*area_carc),0)
################################
......@@ -586,7 +586,7 @@ set.seed(13)
#End step in slaughterhouse: CHILLING
Next5[t]<-max(Next1[t]-(10^2.51*500),0)
Next5[t]<-max(Next1[t]-(10^red_carc*area_carc),0)
#Beef slaughterhouse#
......
......@@ -76,7 +76,7 @@ knitr::opts_chunk$set(echo = TRUE) # code chunks show
```
```{r, include=FALSE}
#Importing data form veal calve studies
#Importing data form veal calves study
veal<-read_excel(here("Data","Raw","veal_calves2.xlsx"),sheet = "Data")
......@@ -249,7 +249,7 @@ veal4$DF1<-factor(veal4$DF1,labels = c("low","high"))
# ATB effect
## Gamma regression for aquiring
## Regression for aquiring
```{r, echo=FALSE}
......@@ -287,13 +287,13 @@ mean((veal_teste$right[veal_teste$atb_gr=="high"]-veal_teste$left[veal_teste$atb
b1 <- flexsurvreg(Surv(left, right, type="interval2") ~ atb_ind,
b1 <- flexsurvreg(Surv(left, right, type="interval2") ~ atb_ind, #Gamma regression
data = veal_teste,dist="gamma")
b1
fit_aft_exp <- ic_par(Surv(left, right, type="interval2") ~ atb_ind,data = veal_teste, model = "aft", dist = "exponential")
fit_aft_exp <- ic_par(Surv(left, right, type="interval2") ~ atb_ind,data = veal_teste, model = "aft", dist = "exponential") #Exponential
summary(fit_aft_exp)
......@@ -328,13 +328,13 @@ ggsurvplot(b1, data = veal4)+
# Group treatment
b2 <- flexsurvreg(Surv(left, right, type="interval2") ~ atb_gr, dist="gamma", data = veal4)
b2 <- flexsurvreg(Surv(left, right, type="interval2") ~ atb_gr, dist="gamma", data = veal4) # Gamma regression
b2
fit_aft_exp2 <- ic_par(Surv(left, right, type="interval2") ~ atb_gr,data = veal_teste, model = "aft", dist = "exponential")
fit_aft_exp2 <- ic_par(Surv(left, right, type="interval2") ~ atb_gr,data = veal_teste, model = "aft", dist = "exponential") #Exponential regression
summary(fit_aft_exp2)
......@@ -366,7 +366,7 @@ ggsurvplot(b2, data = veal4)+
```
## Gamma regression for fade
## Regression for fade
```{r, echo=FALSE}
......@@ -406,11 +406,11 @@ mean((veal_teste2$right[veal_teste2$atb_gr=="high"]-veal_teste2$left[veal_teste2
# Individual treatment
b3 <- flexsurvreg(Surv(left, right, type="interval2") ~ atb_ind, dist="gamma", data = veal6)
b3 <- flexsurvreg(Surv(left, right, type="interval2") ~ atb_ind, dist="gamma", data = veal6) #Gamma regression
b3
fit_aft_exp3 <- ic_par(Surv(left, right, type="interval2") ~ atb_ind,data = veal_teste2, model = "aft", dist = "exponential")
fit_aft_exp3 <- ic_par(Surv(left, right, type="interval2") ~ atb_ind,data = veal_teste2, model = "aft", dist = "exponential") #Exponential
summary(fit_aft_exp3)
......@@ -444,12 +444,12 @@ ggsurvplot(b2, data = veal6)+
# Group treatment
b4 <- flexsurvreg(Surv(left, right, type="interval2") ~ atb_gr, dist="gamma", data = veal6)
b4 <- flexsurvreg(Surv(left, right, type="interval2") ~ atb_gr, dist="gamma", data = veal6) # Gamma regression
b4
fit_aft_exp4 <- ic_par(Surv(left, right, type="interval2") ~ atb_gr,data = veal_teste2, model = "aft", dist = "exponential")
fit_aft_exp4 <- ic_par(Surv(left, right, type="interval2") ~ atb_gr,data = veal_teste2, model = "aft", dist = "exponential") #Exponential
summary(fit_aft_exp4)
......@@ -483,8 +483,7 @@ ggsurvplot(b4, data = veal6)+
# DFs effect
# DFs Effect
## Gamma regression for aquiring
```{r, echo=FALSE}
......@@ -553,82 +552,10 @@ ggsurvplot(b4, data = veal6)+
```
# DF separated
## Gamma regression for acquire
```{r}
b5 <- flexsurvreg(Surv(left, right, type="interval2") ~ factor(DF), dist="gamma", data = veal4)
b5
#Expected time to lose the carriership
mg5<-function(x){
rate.1<-( ((b5$res[2]))*exp(b3$res[3]*x))
mg3<-b3$res[1]*(1/rate.1)
mg3
}
averages3<-list(
#high ATM
high=mg3(1),
#low ATM
low=mg3(0),
p_value=pchisq((b3$res[3]^2/b3$cov[3,3]),1,lower.tail = F)
)
averages3
ggsurvplot(b5, data = veal4)+
ylab("Non-cariership probability")+
xlab("Time (weeks)")
```
## Gamma regression for fade
```{r}
b6 <- flexsurvreg(Surv(left, right, type="interval2") ~ factor(DF), dist="gamma", data = veal6)
b6
#Expected time to lose the carriership
mg6<-function(x){
rate.1<-( ((b4$res[2]))*exp(b4$res[3]*x))
mg4<-b4$res[1]*(1/rate.1)
mg4
}
averages4<-list(
#high ATM
high=mg4(1),
#low ATM
low=mg4(0),
p_value=pchisq((b4$res[3]^2/b4$cov[3,3]),1,lower.tail = F)
)
averages4
ggsurvplot(b6, data = veal6)+
ylab("Cariership probability")+
xlab("Time (weeks)")
```
# VF separated
# VFs Effect
## Gamma regression for acquire
```{r}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment