1. Observations potentially influential
tn<-600
trials.data<-trials.data%>%filter(to_next>=tn)
dim(trials.data)
## [1] 6258 22
plot(trials.data$visit.length, xlab = "Observation", ylab = "Visit length (min)",
main = "Visit length in Distant replacement")
which(trials.data$visit.length>100)
## [1] 5758 6111
a<-trials.data[which(trials.data$visit.length>100),]%>%select(.,Date:Consumed,visit.length,to_next, trial.day, Loc_trial:follower_trial)
kable(a,caption = "Observations potentially influential")%>%
kable_styling(bootstrap_options = c("striped"),full_width = F, position = "left",font_size = 14)
Observations potentially influential
|
Date
|
Entry
|
Exit
|
Consumed
|
visit.length
|
to_next
|
trial.day
|
Loc_trial
|
Eartag_trial
|
follower_trial
|
|
2019-11-04
|
08:42:13
|
10:42:13
|
0.107
|
120.0000
|
79357
|
39
|
Loc_1_t_7
|
1454_t_7
|
1453_t_7
|
|
2019-11-04
|
08:35:32
|
10:38:30
|
0.366
|
122.9667
|
79592
|
39
|
Loc_2_t_7
|
1467_t_7
|
1455_t_7
|
2. Estimating variance component and residuals check with and without observations potentially influential
trials.data<-trials.data[-c(5758,6111),]
# M1 fit
m1<-lmer(visit.length ~ Loc_trial -1 + hour_entry+ wt_median+
(1|Eartag_trial),data=trials.data)
a0<-round(as.matrix(as.data.frame(VarCorr(m1))[,4],1,2),3)
colnames(a0)<-"Estimate"
rownames(a0)<-c("Eartag M1","Residual M1")
resm1wo<-residuals(m1)%>%as.data.frame()
colnames(resm1wo)<-"residual_M1600wo"
# M2 fit
m2<-lmer(visit.length ~ Loc_trial -1 + hour_entry+ wt_median+
(1|Eartag_trial)+(1|follower_trial),data=trials.data)
## boundary (singular) fit: see ?isSingular
b0<-round(as.matrix(as.data.frame(VarCorr(m2))[,4],1,3),3)
colnames(b0)<-"Estimate"
rownames(b0)<-c("Eartag M2","Follower M2","Residual M2")
resm2wo<-residuals(m2)%>%as.data.frame()
colnames(resm2wo)<-"residual_M2600wo"
resid.wo<-cbind(resm1wo,resm2wo)
2.1. Residual check: Q-Q plot and summary
#M1 all observations
resid.all %>% ggplot(aes(sample =residual_M0600)) +stat_qq_band(bandType = "pointwise", fill = "#8DA0CB", alpha = 0.4) + stat_qq_line(colour = "#8DA0CB") +stat_qq_point() + ggtitle("Normal Q-Q plot Model 1 all observations") +
xlab("Normal quantiles") +ylab("Visit lenght measurements quantiles") + theme_light()
# M2 all observations
resid.all %>% ggplot(aes(sample =residual_M1600)) +stat_qq_band(bandType = "pointwise", fill = "#8DA0CB", alpha = 0.4) + stat_qq_line(colour = "#8DA0CB") +stat_qq_point() + ggtitle("Normal Q-Q plot Model 2 all observations") +
xlab("Normal quantiles") +ylab("Visit lenght measurements quantiles") + theme_light()
# M1 without observations
resid.wo %>% ggplot(aes(sample =residual_M1600wo)) +stat_qq_band(bandType = "pointwise", fill = "#1DA0CB", alpha = 0.4) + stat_qq_line(colour = "#8DA0CB") +stat_qq_point() + ggtitle("Normal Q-Q plot Model 1 without observations") +
xlab("Normal quantiles") +ylab("Visit lenght measurements quantiles") + theme_light()
# M2 without observations
resid.wo %>% ggplot(aes(sample =residual_M1600wo)) +stat_qq_band(bandType = "pointwise", fill = "#1DA0CB", alpha = 0.4) + stat_qq_line(colour = "#8DA0CB") +stat_qq_point() + ggtitle("Normal Q-Q plot Model 2 without observations") +
xlab("Normal quantiles") +ylab("Visit lenght measurements quantiles") + theme_light()
resisum<-round(rbind(summary(resid.all$residual_M0600),
summary(resid.all$residual_M1600),
summary(resid.wo$residual_M1600wo),
summary(resid.wo$residual_M2600wo)), 4)
nx<-as.matrix(c("M1","M2","M1","M2"))
colnames(nx)<-"Model"
nwo<-as.matrix(c(rep("With Obs",2), rep("Without Obs", 2)))
colnames(nwo)<-"Subset"
resisum<-cbind(nwo,nx,resisum)
kable(resisum,caption = "Summary Residuals Distant replacement")%>%kable_styling(bootstrap_options = c("striped"),full_width = F, position = "left",font_size = 14)%>%
column_spec(1,italic = T, bold = T)%>%
collapse_rows(columns = 1)
Summary Residuals Distant replacement
|
Subset
|
Model
|
Min.
|
1st Qu.
|
Median
|
Mean
|
3rd Qu.
|
Max.
|
|
With Obs
|
M1
|
-19.708
|
-4.4465
|
-0.7556
|
0
|
3.5369
|
112.4495
|
|
M2
|
-19.708
|
-4.4465
|
-0.7556
|
0
|
3.5369
|
112.4495
|
|
Without Obs
|
M1
|
-19.9284
|
-4.3751
|
-0.7454
|
0
|
3.5282
|
47.4163
|
|
M2
|
-19.9284
|
-4.3751
|
-0.7454
|
0
|
3.5282
|
47.4163
|
3. Estimated variance components
kable(m,caption = "Estimated of variance components Distant replacement")%>%kable_styling(bootstrap_options = c("striped"),full_width = F, position = "left",font_size = 14)%>%add_header_above(c(" "=2,"Without points "= 1, "With points"=1))%>%
column_spec(1,italic = T, bold = T)%>%
collapse_rows(columns = 1)
Estimated of variance components Distant replacement
|
|
Without points
|
With points
|
|
Model
|
Variance
|
Estimate
|
Estimate
|
|
M1
|
Eartag
|
11.186
|
10.64
|
|
Residual
|
43.482
|
47.767
|
|
M2
|
Eartag
|
11.186
|
10.64
|
|
Follower
|
0
|
0
|
|
Residual
|
43.482
|
47.767
|
4. Model comparison
kable(modc,caption = "Model comparison Distant replacement")%>%kable_styling(bootstrap_options = c("striped"),full_width = F, position = "left",font_size = 14)%>%
column_spec(1,italic = T, bold = T)%>%
column_spec(8:10, color = "white")%>%
row_spec(c(2,4), color = "black")%>%
collapse_rows(columns = 1)
Model comparison Distant replacement
|
Subset
|
Model
|
Df
|
AIC
|
BIC
|
logLik
|
deviance
|
Chisq
|
Chi Df
|
Pr(>Chisq)
|
|
Without points
|
M1
|
40
|
41719.20
|
41988.86
|
-20819.60
|
41639.20
|
NA
|
NA
|
NA
|
|
M2
|
41
|
41721.20
|
41997.60
|
-20819.60
|
41639.20
|
0
|
1
|
1
|
|
With points
|
M1
|
40
|
42303.49
|
42573.15
|
-21111.74
|
42223.49
|
NA
|
NA
|
NA
|
|
M2
|
41
|
42305.49
|
42581.89
|
-21111.74
|
42223.49
|
0
|
1
|
1
|