Chapter 3
Emily Runnion
2025-01-24
#input libraries
library(readr)
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(tidyr)
library(ordinal)## Warning: package 'ordinal' was built under R version 4.4.3##
## Attaching package: 'ordinal'## The following object is masked from 'package:dplyr':
##
## slicelibrary(RVAideMemoire)## Warning: package 'RVAideMemoire' was built under R version 4.4.3## *** Package RVAideMemoire v 0.9-83-7 ***library(ggplot2)
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectlibrary(emmeans)## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'library(cowplot)
library(Hmisc)##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:dplyr':
##
## src, summarize## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(reshape2)##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smithslibrary(stats)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodeinput data
pol_pes_22 <- read.csv("pol_pes_22.csv")
pol_pes_22$location <- as.factor(pol_pes_22$location)
pol_pes_22$pesticide <- as.factor(pol_pes_22$pesticide)
pol_pes_22$gh <- as.factor(pol_pes_22$gh)
pol_pes_22$week <- as.factor(pol_pes_22$week)
pol_pes_22$hives <- as.double(pol_pes_22$hives)
pol_pes_22$joules <- as.double(pol_pes_22$joules)
pol_pes_22$ordered_score <- as.ordered(pol_pes_22$score)
pol_pes_22$week_number <- as.numeric(pol_pes_22$week)
#scale the variables to see if this helps
pol_pes_22$hives_scaled <- scale(pol_pes_22$hives)
pol_pes_22$joules_scaled <- scale(pol_pes_22$joules)
pol_pes_22$nt_sc <- scale(pol_pes_22$night_temp)
pol_pes_22$dt_sc <- scale(pol_pes_22$day_temp)
pol_pes_22$hrsl_sc <- scale(pol_pes_22$hrs_arfc_light)look at all model possibilities
# everything model with regular variables not scaled
#m_all <- clmm(ordered_score ~ hives + week + night_temp + day_temp + hrs_arfc_light + joules + bay_applied_fungicide + bay_applied_insecticide + bay_applied_biological + (1|location), data = pol_pes_22, Hess = TRUE)
#drop1(m_all) this won't even run because model is so overfit
#summary(m_all)select coefficients
#get rid of day temp, artificial light, and joules
m_all1 <- clmm(ordered_score ~ hives + week + night_temp + bay_applied_fungicide + bay_applied_insecticide + bay_applied_biological + (1|location), data = pol_pes_22, Hess = TRUE)
drop1(m_all1, test = "Chisq")## Single term deletions
##
## Model:
## ordered_score ~ hives + week + night_temp + bay_applied_fungicide +
## bay_applied_insecticide + bay_applied_biological + (1 | location)
## Df AIC LRT Pr(>Chi)
## <none> 13108
## hives 1 13206 100.83 < 2.2e-16 ***
## week 12 15938 2854.39 < 2.2e-16 ***
## night_temp 1 13131 25.78 3.824e-07 ***
## bay_applied_fungicide 1 13202 96.20 < 2.2e-16 ***
## bay_applied_insecticide 1 13106 0.01 0.9156603
## bay_applied_biological 1 13120 15.00 0.0001073 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m_all2 <- update(m_all1, .~. -bay_applied_insecticide)
drop1(m_all2, test = "Chisq")## Single term deletions
##
## Model:
## ordered_score ~ hives + week + night_temp + bay_applied_fungicide +
## bay_applied_biological + (1 | location)
## Df AIC LRT Pr(>Chi)
## <none> 13106
## hives 1 13204 100.82 < 2.2e-16 ***
## week 12 15999 2917.69 < 2.2e-16 ***
## night_temp 1 13129 25.78 3.832e-07 ***
## bay_applied_fungicide 1 13201 97.03 < 2.2e-16 ***
## bay_applied_biological 1 13119 15.05 0.0001045 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# see what happens with week as a number not a factor
m_all_wn <- clmm(ordered_score ~ hives + week_number + night_temp + bay_applied_fungicide + bay_applied_insecticide + bay_applied_biological + (1|location), data = pol_pes_22, Hess = TRUE)
anova(m_all1, m_all2, test = "Chisq")## 'test' argument ignored in anova.clm## Likelihood ratio tests of cumulative link models:
##
## formula:
## m_all2 ordered_score ~ hives + week + night_temp + bay_applied_fungicide + bay_applied_biological + (1 | location)
## m_all1 ordered_score ~ hives + week + night_temp + bay_applied_fungicide + bay_applied_insecticide + bay_applied_biological + (1 | location)
## link: threshold:
## m_all2 logit flexible
## m_all1 logit flexible
##
## no.par AIC logLik LR.stat df Pr(>Chisq)
## m_all2 19 13106 -6533.8
## m_all1 20 13108 -6533.8 0.0112 1 0.9157anova(m_all1, m_all2, m_all_wn, test = "Chisq")## 'test' argument ignored in anova.clm## Likelihood ratio tests of cumulative link models:
##
## formula:
## m_all_wn ordered_score ~ hives + week_number + night_temp + bay_applied_fungicide + bay_applied_insecticide + bay_applied_biological + (1 | location)
## m_all2 ordered_score ~ hives + week + night_temp + bay_applied_fungicide + bay_applied_biological + (1 | location)
## m_all1 ordered_score ~ hives + week + night_temp + bay_applied_fungicide + bay_applied_insecticide + bay_applied_biological + (1 | location)
## link: threshold:
## m_all_wn logit flexible
## m_all2 logit flexible
## m_all1 logit flexible
##
## no.par AIC logLik LR.stat df Pr(>Chisq)
## m_all_wn 9 15055 -7518.7
## m_all2 19 13106 -6533.8 1969.9107 10 <2e-16 ***
## m_all1 20 13108 -6533.8 0.0112 1 0.9157
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#removing insecticide does lower the AIC but the model doesn't fit significantly better
#yes week should be a factor test for importance of random
# model with no random effect
m_no_rand <- clm(ordered_score ~ hives + week + night_temp + bay_applied_fungicide + bay_applied_insecticide + bay_applied_biological, data = pol_pes_22, Hess = TRUE)
summary(m_no_rand)## formula:
## ordered_score ~ hives + week + night_temp + bay_applied_fungicide + bay_applied_insecticide + bay_applied_biological
## data: pol_pes_22
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 9260 -6878.69 13795.38 7(0) 4.01e-09 3.1e+07
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## hives -0.011961 0.001701 -7.030 2.07e-12 ***
## week31 0.572906 0.226110 2.534 0.0113 *
## week32 2.531594 0.152783 16.570 < 2e-16 ***
## week33 4.164175 0.192853 21.593 < 2e-16 ***
## week34 4.456901 0.195171 22.836 < 2e-16 ***
## week35 5.879650 0.215466 27.288 < 2e-16 ***
## week36 5.304095 0.202036 26.253 < 2e-16 ***
## week37 4.946386 0.212621 23.264 < 2e-16 ***
## week38 5.164880 0.187790 27.503 < 2e-16 ***
## week39 3.788695 0.228828 16.557 < 2e-16 ***
## week40 4.407633 0.262933 16.763 < 2e-16 ***
## week41 4.577075 0.249399 18.352 < 2e-16 ***
## week42 4.845731 0.178754 27.108 < 2e-16 ***
## night_temp -0.018779 0.020185 -0.930 0.3522
## bay_applied_fungicideTRUE -1.537952 0.088427 -17.392 < 2e-16 ***
## bay_applied_insecticideTRUE -0.003791 0.079647 -0.048 0.9620
## bay_applied_biologicalTRUE 0.006576 0.093921 0.070 0.9442
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -1.0271 0.4862 -2.113
## 2|3 2.3648 0.4863 4.863# do we want random?
anova(m_all1, m_no_rand, test = "Chisq")## 'test' argument ignored in anova.clm## Likelihood ratio tests of cumulative link models:
##
## formula:
## m_no_rand ordered_score ~ hives + week + night_temp + bay_applied_fungicide + bay_applied_insecticide + bay_applied_biological
## m_all1 ordered_score ~ hives + week + night_temp + bay_applied_fungicide + bay_applied_insecticide + bay_applied_biological + (1 | location)
## link: threshold:
## m_no_rand logit flexible
## m_all1 logit flexible
##
## no.par AIC logLik LR.stat df Pr(>Chisq)
## m_no_rand 19 13795 -6878.7
## m_all1 20 13108 -6533.8 689.86 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#yes we want random plot interactions
There are no instances where fungicide + insecticide + biological were applied all at once - in fact, each type is ever only applied independently, so there should be no interactions between the types that we need to take into account (multivariate analysis possible?)
ggplot(pol_pes_22, aes(x = ordered_score, y = hives)) +
geom_boxplot(size = .75) +
geom_jitter(alpha = .5) +
facet_grid(bay_applied_fungicide ~ bay_applied_biological, margins = TRUE) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
ggplot(pol_pes_22, aes(x = ordered_score, y = hives)) +
geom_boxplot(size = .75) +
geom_jitter(alpha = .5) +
facet_grid(bay_applied_fungicide ~ bay_applied_insecticide, margins = TRUE) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
ggplot(pol_pes_22, aes(x = ordered_score, y = hives)) +
geom_boxplot(size = .75) +
geom_jitter(alpha = .5) +
facet_grid(bay_applied_insecticide ~ bay_applied_biological, margins = TRUE) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
estimate an ordered logistic regression model
m <- polr(ordered_score ~ hives + week + night_temp + bay_applied_fungicide + bay_applied_insecticide + bay_applied_biological, data = pol_pes_22, Hess = TRUE)
summary(m)## Call:
## polr(formula = ordered_score ~ hives + week + night_temp + bay_applied_fungicide +
## bay_applied_insecticide + bay_applied_biological, data = pol_pes_22,
## Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## hives -0.011960 0.001703 -7.02223
## week31 0.572879 0.226073 2.53404
## week32 2.531561 0.152747 16.57354
## week33 4.164131 0.192984 21.57754
## week34 4.456872 0.195299 22.82079
## week35 5.879605 0.215721 27.25562
## week36 5.304049 0.202206 26.23087
## week37 4.946302 0.212749 23.24951
## week38 5.164818 0.187947 27.48020
## week39 3.788582 0.228992 16.54459
## week40 4.407473 0.263093 16.75255
## week41 4.576946 0.249561 18.33998
## week42 4.845628 0.178888 27.08751
## night_temp -0.018785 0.020174 -0.93112
## bay_applied_fungicideTRUE -1.537936 0.088440 -17.38955
## bay_applied_insecticideTRUE -0.003770 0.079647 -0.04733
## bay_applied_biologicalTRUE 0.006557 0.093919 0.06981
##
## Intercepts:
## Value Std. Error t value
## 1|2 -1.0272 0.4860 -2.1136
## 2|3 2.3647 0.4859 4.8668
##
## Residual Deviance: 13757.38
## AIC: 13795.38(ctable <- coef(summary(m)))## Value Std. Error t value
## hives -0.011960343 0.001703211 -7.02223244
## week31 0.572878749 0.226073330 2.53403950
## week32 2.531560814 0.152747117 16.57354231
## week33 4.164131450 0.192984488 21.57754484
## week34 4.456871767 0.195298718 22.82079375
## week35 5.879604638 0.215720835 27.25561786
## week36 5.304048835 0.202206334 26.23087389
## week37 4.946301603 0.212748612 23.24951288
## week38 5.164817975 0.187946851 27.48020486
## week39 3.788582178 0.228992276 16.54458505
## week40 4.407472984 0.263092686 16.75254848
## week41 4.576945714 0.249561174 18.33997506
## week42 4.845628310 0.178887932 27.08750811
## night_temp -0.018784746 0.020174352 -0.93112016
## bay_applied_fungicideTRUE -1.537935868 0.088440211 -17.38955454
## bay_applied_insecticideTRUE -0.003769514 0.079646798 -0.04732787
## bay_applied_biologicalTRUE 0.006556873 0.093918672 0.06981437
## 1|2 -1.027191189 0.485985975 -2.11362311
## 2|3 2.364720482 0.485890893 4.86677260p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
## combined table
(ctable <- cbind(ctable, "p value" = p))## Value Std. Error t value p value
## hives -0.011960343 0.001703211 -7.02223244 2.183509e-12
## week31 0.572878749 0.226073330 2.53403950 1.127560e-02
## week32 2.531560814 0.152747117 16.57354231 1.082563e-61
## week33 4.164131450 0.192984488 21.57754484 2.919648e-103
## week34 4.456871767 0.195298718 22.82079375 2.850794e-115
## week35 5.879604638 0.215720835 27.25561786 1.425646e-163
## week36 5.304048835 0.202206334 26.23087389 1.181518e-151
## week37 4.946301603 0.212748612 23.24951288 1.438764e-119
## week38 5.164817975 0.187946851 27.48020486 3.027444e-166
## week39 3.788582178 0.228992276 16.54458505 1.751675e-61
## week40 4.407472984 0.263092686 16.75254848 5.425157e-63
## week41 4.576945714 0.249561174 18.33997506 3.969792e-75
## week42 4.845628310 0.178887932 27.08750811 1.381859e-161
## night_temp -0.018784746 0.020174352 -0.93112016 3.517914e-01
## bay_applied_fungicideTRUE -1.537935868 0.088440211 -17.38955454 9.900070e-68
## bay_applied_insecticideTRUE -0.003769514 0.079646798 -0.04732787 9.622519e-01
## bay_applied_biologicalTRUE 0.006556873 0.093918672 0.06981437 9.443414e-01
## 1|2 -1.027191189 0.485985975 -2.11362311 3.454747e-02
## 2|3 2.364720482 0.485890893 4.86677260 1.134354e-06# gets odds ratios
exp(coef(m))## hives week31
## 0.9881109 1.7733648
## week32 week33
## 12.5731151 64.3367784
## week34 week35
## 86.2173787 357.6678056
## week36 week37
## 201.1495849 140.6538072
## week38 week39
## 175.0056003 44.1936971
## week40 week41
## 82.0618297 97.2170117
## week42 night_temp
## 127.1831674 0.9813906
## bay_applied_fungicideTRUE bay_applied_insecticideTRUE
## 0.2148241 0.9962376
## bay_applied_biologicalTRUE
## 1.0065784#confidence interval
(ci <- confint(m))## Waiting for profiling to be done...## 2.5 % 97.5 %
## hives -0.01530556 -0.008634376
## week31 0.13127979 1.017815346
## week32 2.23550845 2.834688451
## week33 3.78888259 4.545030769
## week34 4.07744359 4.842659605
## week35 5.46104446 6.305811488
## week36 4.91129825 5.703420548
## week37 4.53276315 5.366391242
## week38 4.79963684 5.535951585
## week39 3.34164735 4.239190977
## week40 3.89341127 4.924845637
## week41 4.08982272 5.067979691
## week42 4.49877836 5.199672102
## night_temp -0.05838937 0.020802942
## bay_applied_fungicideTRUE -1.71192030 -1.365239566
## bay_applied_insecticideTRUE -0.15980096 0.152471581
## bay_applied_biologicalTRUE -0.17726125 0.190897437#odds ratios and confidence intervals
exp(cbind(OR = coef(m), ci))## OR 2.5 % 97.5 %
## hives 0.9881109 0.9848110 0.9914028
## week31 1.7733648 1.1402868 2.7671429
## week32 12.5731151 9.3512353 17.0250953
## week33 64.3367784 44.2069752 94.1633245
## week34 86.2173787 58.9944629 126.8061579
## week35 357.6678056 235.3431033 547.7458973
## week36 201.1495849 135.8156220 299.8914415
## week37 140.6538072 93.0152215 214.0888771
## week38 175.0056003 121.4662979 253.6490414
## week39 44.1936971 28.2656517 69.3517220
## week40 82.0618297 49.0780196 137.6680891
## week41 97.2170117 59.7293019 158.8530707
## week42 127.1831674 89.9072299 181.2128128
## night_temp 0.9813906 0.9432826 1.0210208
## bay_applied_fungicideTRUE 0.2148241 0.1805188 0.2553195
## bay_applied_insecticideTRUE 0.9962376 0.8523134 1.1647094
## bay_applied_biologicalTRUE 1.0065784 0.8375609 1.2103353Every 1 unit increase in hives the odds of seeing a higher score is multiplied 0.98 times
When fungicide was applied, the odds of seeing scors greater than 1 are 0.21 times higher than when there wasn’t fungicide When fungiicde was not applied the odds of seeing a score = 1 is 0.21 times higher than when there was fungicide
Proportional odds assumption
sf <- function(y) {
c('Y>=1' = qlogis(mean(y >= 1)),
'Y>=2' = qlogis(mean(y >= 2)),
'Y>=3' = qlogis(mean(y >= 3)))
}
(s <- with(pol_pes_22, summary(as.numeric(ordered_score) ~ hives + week + night_temp + bay_applied_fungicide + bay_applied_insecticide + bay_applied_biological, fun=sf)))## as.numeric(ordered_score) N= 9260
##
## +-----------------------+-----------+----+----+------------+------------+
## | | | N|Y>=1| Y>=2| Y>=3|
## +-----------------------+-----------+----+----+------------+------------+
## | hives| [ 70,150)|2795| Inf| 0.204646844|-1.278010401|
## | | [150,165)|2207| Inf| 3.612314886| 0.298527305|
## | | [165,190)|2398| Inf| 2.364693479|-0.185685989|
## | | [190,220]|1860| Inf| 2.580101162|-0.645137961|
## +-----------------------+-----------+----+----+------------+------------+
## | week| 30| 816| Inf|-1.968509981|-6.008813185|
## | | 31| 396| Inf|-0.451985124|-4.359269648|
## | | 32| 593| Inf|-0.003372684|-2.739046301|
## | | 33| 594| Inf| 4.768988271|-0.395654238|
## | | 34| 594| Inf| 4.768988271|-0.094346015|
## | | 35| 594| Inf| Inf| 1.262529731|
## | | 36| 594| Inf| 6.385194399| 0.708337346|
## | | 37| 383| Inf| 5.249652195| 0.257284595|
## | | 38|1116| Inf| 4.365824660| 0.474665642|
## | | 39|1126| Inf| 1.569474656|-1.003743513|
## | | 40|1056| Inf| 2.337453919|-0.850907973|
## | | 41| 696| Inf| 2.824122372|-0.443931389|
## | | 42| 702| Inf| 2.754570217| 0.188591170|
## +-----------------------+-----------+----+----+------------+------------+
## | night_temp|[ 6.9,14.6)|2986| Inf| 2.099184020|-0.779315106|
## | |[14.6,18.1)|2686| Inf| 2.022095774|-0.007446051|
## | |[18.1,19.1)|2004| Inf| 0.486414754|-0.907922482|
## | |[19.1,22.7]|1584| Inf| 1.688795240|-0.103627990|
## +-----------------------+-----------+----+----+------------+------------+
## | bay_applied_fungicide| No|7059| Inf| 2.654182513|-0.043355706|
## | | Yes|2201| Inf|-0.155697688|-2.608315235|
## +-----------------------+-----------+----+----+------------+------------+
## |bay_applied_insecticide| No|8224| Inf| 1.392995574|-0.567793543|
## | | Yes|1036| Inf| 4.290459441| 0.435318071|
## +-----------------------+-----------+----+----+------------+------------+
## | bay_applied_biological| No|8306| Inf| 1.546059331|-0.487647686|
## | | Yes| 954| Inf| 1.384984603|-0.138586163|
## +-----------------------+-----------+----+----+------------+------------+
## | Overall| |9260| Inf| 1.528687884|-0.450672589|
## +-----------------------+-----------+----+----+------------+------------+s[, 4] <- s[, 4] - s[, 3]
s[, 3] <- s[, 3] - s[, 3]
s # print## as.numeric(ordered_score) N= 9260
##
## +-----------------------+-----------+----+----+----+---------+
## | | | N|Y>=1|Y>=2| Y>=3|
## +-----------------------+-----------+----+----+----+---------+
## | hives| [ 70,150)|2795| Inf| 0|-1.482657|
## | | [150,165)|2207| Inf| 0|-3.313788|
## | | [165,190)|2398| Inf| 0|-2.550379|
## | | [190,220]|1860| Inf| 0|-3.225239|
## +-----------------------+-----------+----+----+----+---------+
## | week| 30| 816| Inf| 0|-4.040303|
## | | 31| 396| Inf| 0|-3.907285|
## | | 32| 593| Inf| 0|-2.735674|
## | | 33| 594| Inf| 0|-5.164643|
## | | 34| 594| Inf| 0|-4.863334|
## | | 35| 594| Inf| | -Inf|
## | | 36| 594| Inf| 0|-5.676857|
## | | 37| 383| Inf| 0|-4.992368|
## | | 38|1116| Inf| 0|-3.891159|
## | | 39|1126| Inf| 0|-2.573218|
## | | 40|1056| Inf| 0|-3.188362|
## | | 41| 696| Inf| 0|-3.268054|
## | | 42| 702| Inf| 0|-2.565979|
## +-----------------------+-----------+----+----+----+---------+
## | night_temp|[ 6.9,14.6)|2986| Inf| 0|-2.878499|
## | |[14.6,18.1)|2686| Inf| 0|-2.029542|
## | |[18.1,19.1)|2004| Inf| 0|-1.394337|
## | |[19.1,22.7]|1584| Inf| 0|-1.792423|
## +-----------------------+-----------+----+----+----+---------+
## | bay_applied_fungicide| No|7059| Inf| 0|-2.697538|
## | | Yes|2201| Inf| 0|-2.452618|
## +-----------------------+-----------+----+----+----+---------+
## |bay_applied_insecticide| No|8224| Inf| 0|-1.960789|
## | | Yes|1036| Inf| 0|-3.855141|
## +-----------------------+-----------+----+----+----+---------+
## | bay_applied_biological| No|8306| Inf| 0|-2.033707|
## | | Yes| 954| Inf| 0|-1.523571|
## +-----------------------+-----------+----+----+----+---------+
## | Overall| |9260| Inf| 0|-1.979360|
## +-----------------------+-----------+----+----+----+---------+plot(s, which = 1:3, pch = 1:3, xlab = 'logit', main = ' ',
xlim = range(as.matrix(s[, 3:4]), finite = TRUE))
# fungicide might actually be our best bet at non-nominal data #obtain predicted probabilites
newdat <- expand.grid(
hives = seq(min(pol_pes_22$hives), max(pol_pes_22$hives), length.out = 5),
week = factor(1:13, levels = levels(pol_pes_22$week)), # ensure same factor levels
night_temp = seq(min(pol_pes_22$night_temp), max(pol_pes_22$night_temp), length.out = 5),
bay_applied_fungicide = c(TRUE, FALSE),
bay_applied_insecticide = c(TRUE, FALSE),
bay_applied_biological = c(TRUE, FALSE)
)
newdat <- cbind(newdat, predict(m, newdat, type = "probs"))
##oh my god i think it worked
head(newdat)## hives week night_temp bay_applied_fungicide bay_applied_insecticide
## 1 70.0 <NA> 6.9 TRUE TRUE
## 2 107.5 <NA> 6.9 TRUE TRUE
## 3 145.0 <NA> 6.9 TRUE TRUE
## 4 182.5 <NA> 6.9 TRUE TRUE
## 5 220.0 <NA> 6.9 TRUE TRUE
## 6 70.0 <NA> 6.9 TRUE TRUE
## bay_applied_biological 1 2 3
## 1 TRUE NA NA NA
## 2 TRUE NA NA NA
## 3 TRUE NA NA NA
## 4 TRUE NA NA NA
## 5 TRUE NA NA NA
## 6 TRUE NA NA NAlnewdat <- melt(newdat, id.vars = c("hives", "week", "night_temp", "bay_applied_fungicide", "bay_applied_insecticide", "bay_applied_biological"),
variable.name = "Level", value.name="Probability")
## view first few rows
head(lnewdat)## hives week night_temp bay_applied_fungicide bay_applied_insecticide
## 1 70.0 <NA> 6.9 TRUE TRUE
## 2 107.5 <NA> 6.9 TRUE TRUE
## 3 145.0 <NA> 6.9 TRUE TRUE
## 4 182.5 <NA> 6.9 TRUE TRUE
## 5 220.0 <NA> 6.9 TRUE TRUE
## 6 70.0 <NA> 6.9 TRUE TRUE
## bay_applied_biological Level Probability
## 1 TRUE 1 NA
## 2 TRUE 1 NA
## 3 TRUE 1 NA
## 4 TRUE 1 NA
## 5 TRUE 1 NA
## 6 TRUE 1 NAlnewdat$week <- as.numeric(lnewdat$week)
ggplot(lnewdat, aes(x = week, y = Probability, colour = Level)) +
geom_line() + facet_wrap(vars(bay_applied_fungicide))## Warning: Removed 7800 rows containing missing values or values outside the scale range
## (`geom_line()`).
ggplot(lnewdat, aes(x = bay_applied_fungicide, y = Probability, color = Level, group = Level)) +
geom_boxplot()## Warning: Removed 7800 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
run multivariate anlysis mixed model with nominal effects
m1 <- clmm2(location = ordered_score ~ bay_applied_fungicide,
nominal = ~ week,
random = location,
data = pol_pes_22,
Hess = T,
nAGQ = 10, link = "probit", method = "nlminb")
summary(m1)## Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
## quadrature approximation with 10 quadrature points
##
## Call:
## clmm2(location = ordered_score ~ bay_applied_fungicide, nominal = ~week,
## random = location, data = pol_pes_22, Hess = T, link = "probit",
## nAGQ = 10, method = "nlminb")
##
## Random effects:
## Var Std.Dev
## location 0.1537035 0.3920504
##
## Location coefficients:
## Estimate Std. Error z value Pr(>|z|)
## bay_applied_fungicideTRUE -0.5048 0.0512 -9.8593 < 2.22e-16
##
## No scale coefficients
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2.(Intercept) 0.8254 0.0897 9.1990
## 2|3.(Intercept) 2.5668 0.2460 10.4358
## 1|2.week31 -0.4673 0.1022 -4.5730
## 2|3.week31 -0.2843 0.2976 -0.9553
## 1|2.week32 -1.3318 0.0811 -16.4204
## 2|3.week32 -1.4335 0.2513 -5.7046
## 1|2.week33 -3.2892 0.1833 -17.9397
## 2|3.week33 -2.3136 0.2484 -9.3127
## 1|2.week34 -3.3171 0.1871 -17.7264
## 2|3.week34 -2.5064 0.2483 -10.0930
## 1|2.week35 -6.2164 18.5366 -0.3354
## 2|3.week35 -3.3816 0.2500 -13.5244
## 1|2.week36 -4.0358 0.3494 -11.5502
## 2|3.week36 -3.0376 0.2491 -12.1947
## 1|2.week37 -3.4868 0.2619 -13.3159
## 2|3.week37 -2.7697 0.2521 -10.9870
## 1|2.week38 -3.4347 0.1423 -24.1411
## 2|3.week38 -2.8951 0.2455 -11.7935
## 1|2.week39 -2.1268 0.0887 -23.9722
## 2|3.week39 -2.0416 0.2447 -8.3416
## 1|2.week40 -2.6489 0.0969 -27.3318
## 2|3.week40 -2.1418 0.2447 -8.7514
## 1|2.week41 -2.7412 0.1167 -23.4852
## 2|3.week41 -2.3154 0.2473 -9.3643
## 1|2.week42 -2.8087 0.1143 -24.5767
## 2|3.week42 -2.7984 0.2463 -11.3637
##
## log-likelihood: -6513.224
## AIC: 13082.45
## Condition number of Hessian: 1200367.94m2 <- clmm2(location = ordered_score ~ bay_applied_fungicide,
nominal = ~ week + hives,
random = location,
data = pol_pes_22,
Hess = T,
nAGQ = 10, link = "probit", method = "nlminb")## Warning in rho$updateU(rho): Non finite negative log-likelihood
## at iteration 96summary(m2)## Warning in sqrt(diag(vc)): NaNs produced## Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
## quadrature approximation with 10 quadrature points
##
## Call:
## clmm2(location = ordered_score ~ bay_applied_fungicide, nominal = ~week +
## hives, random = location, data = pol_pes_22, Hess = T, link = "probit",
## nAGQ = 10, method = "nlminb")
##
## Random effects:
## Var Std.Dev
## location 0.1591919 0.3989886
##
## Location coefficients:
## Estimate Std. Error z value Pr(>|z|)
## bay_applied_fungicideTRUE -0.5067 0.0513 -9.8774 < 2.22e-16
##
## No scale coefficients
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2.(Intercept) 0.7249 0.1864 3.8891
## 2|3.(Intercept) 1.5396 0.2686 5.7326
## 1|2.week31 -0.5178 0.1251 -4.1399
## 2|3.week31 -0.7761 0.3083 -2.5174
## 1|2.week32 -1.4115 0.1080 -13.0719
## 2|3.week32 -2.0995 0.2681 -7.8306
## 1|2.week33 -3.3963 0.2040 -16.6443
## 2|3.week33 -3.0428 0.2665 -11.4166
## 1|2.week34 -3.4162 0.2084 -16.3950
## 2|3.week34 -3.2343 0.2664 -12.1395
## 1|2.week35 -7.3571 NaN NaN
## 2|3.week35 -4.2018 0.2706 -15.5271
## 1|2.week36 -4.1405 0.3617 -11.4483
## 2|3.week36 -3.8605 0.2698 -14.3080
## 1|2.week37 -3.6058 0.2870 -12.5634
## 2|3.week37 -3.6774 0.2753 -13.3600
## 1|2.week38 -3.5328 0.1838 -19.2186
## 2|3.week38 -3.7857 0.2686 -14.0920
## 1|2.week39 -2.2329 0.1782 -12.5314
## 2|3.week39 -3.1677 0.2762 -11.4667
## 1|2.week40 -2.7760 0.2131 -13.0252
## 2|3.week40 -3.5073 0.2863 -12.2504
## 1|2.week41 -2.8614 0.2344 -12.2097
## 2|3.week41 -3.6058 0.2858 -12.6179
## 1|2.week42 -2.8562 0.1531 -18.6578
## 2|3.week42 -3.5095 0.2646 -13.2620
## 1|2.hives 0.0013 0.0017 0.7558
## 2|3.hives 0.0118 0.0011 10.5855
##
## log-likelihood: -6454.514
## AIC: 12969.03
## Condition number of Hessian: 903604111751.62m3 <- clmm2(location = ordered_score ~ bay_applied_fungicide,
nominal = ~ week,
random = location,
data = pol_pes_22,
Hess = T,
nAGQ = 10)
summary(m3)## Warning in sqrt(diag(vc)): NaNs produced## Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
## quadrature approximation with 10 quadrature points
##
## Call:
## clmm2(location = ordered_score ~ bay_applied_fungicide, nominal = ~week,
## random = location, data = pol_pes_22, Hess = T, nAGQ = 10)
##
## Random effects:
## Var Std.Dev
## location 0.4537316 0.673596
##
## Location coefficients:
## Estimate Std. Error z value Pr(>|z|)
## bay_applied_fungicideTRUE -0.8901 0.0894 -9.9519 < 2.22e-16
##
## No scale coefficients
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2.(Intercept) 1.3710 0.1567 8.7483
## 2|3.(Intercept) 5.4584 0.1223 44.6204
## 1|2.week31 -0.7193 0.1746 -4.1195
## 2|3.week31 -0.9269 0.4546 -2.0387
## 1|2.week32 -2.2427 0.1410 -15.9034
## 2|3.week32 -3.4249 0.2355 -14.5430
## 1|2.week33 -6.2958 0.4729 -13.3145
## 2|3.week33 -5.0208 0.1347 -37.2836
## 1|2.week34 -6.2990 0.4728 -13.3226
## 2|3.week34 -5.3480 0.1336 -40.0181
## 1|2.week35 -31.4440 NaN NaN
## 2|3.week35 -6.8378 0.1448 -47.2220
## 1|2.week36 -8.0249 1.0120 -7.9296
## 2|3.week36 -6.2301 0.1031 -60.4222
## 1|2.week37 -6.7804 0.7246 -9.3570
## 2|3.week37 -5.7835 0.1513 -38.2298
## 1|2.week38 -6.2038 0.3031 -20.4689
## 2|3.week38 -5.9940 0.1225 -48.9272
## 1|2.week39 -3.5737 0.1557 -22.9455
## 2|3.week39 -4.5938 0.1280 -35.9005
## 1|2.week40 -4.4775 0.1729 -25.9020
## 2|3.week40 -4.7602 0.1080 -44.0868
## 1|2.week41 -4.7253 0.2181 -21.6611
## 2|3.week41 -5.0213 0.1317 -38.1335
## 1|2.week42 -4.7776 0.2104 -22.7093
## 2|3.week42 -5.8284 0.1349 -43.1954
##
## log-likelihood: -6530.368
## AIC: 13116.74
## Condition number of Hessian: 1.028937e+12AIC(m1, m2, m3)## df AIC
## m1 28 13082.45
## m2 30 12969.03
## m3 28 13116.74#Hess = T, nAGQ = 10, link = "probit", method = "nlminb", threshold = "equidistant"
# 'arg' should be one of “logistic”, “probit”, “cloglog”, “loglog”, “cauchit”, “Aranda-Ordaz”, “log-gamma”
# 'arg' should be one of “ucminf”, “nlminb”, “model.frame”
# arg' should be one of “flexible”, “symmetric”, “equidistant”#I’m trying to figure out how to do the formula manually here but it’s not going well
# Intercept (week 1)
eta <- -0.8901
theta <- c(1.37, 5.45)
# f_j values: theta - eta
f1 <- theta[1] - eta
f2 <- theta[2] - eta
# category probabilities
P1 <- 1 / (1 + exp(-f1))
P2 <- 1 / (1 + exp(-f2)) - P1
P3 <- 1 - 1 / (1 + exp(-f2))
round(c(score_1 = P1, score_2 = P2, score_3 = P3), 4)## score_1 score_2 score_3
## 0.9055 0.0927 0.0018# week 2
eta <- -0.8901
theta <- c(-0.7193, -0.9269)
# f_j values: theta - eta
f1 <- theta[1] - eta
f2 <- theta[2] - eta
#category probabilities
P1 <- 1 / (1 + exp(-f1))
P2 <- 1 / (1 + exp(-f2)) - P1
P3 <- 1 - 1 / (1 + exp(-f2))
round(c(score_1 = P1, score_2 = P2, score_3 = P3), 4)## score_1 score_2 score_3
## 0.5426 -0.0518 0.5092