Inter-rater and intra-rater reliability
Thomas Weissensteiner
2024-12-09
Note: This document would benefit from a few additions, improvements, and possibly corrections. I ran out of time but will post an updated version after the webinar at Rpubs/thomas-weissensteiner
# - r libraries and data - #
library(readxl) # version 1.4.3 # read_excel
library(dplyr) # version 1.1.4 # general grammar
library(purrr) # version 1.0.2 # list_rbind
library(tidyr) # version 1.3.1 # pivot_longer, pivot_wider
library(ggplot2) # version 3.4.0
library(ggside) # version 0.3.1 # geom_ysidedensity (density side plots)
library(ggiraph) # version 0.8.10 # girafe, geom_..._interactive (interactive plots)
library(patchwork) # version 1.3.0 # plot_spacer, plot_layout, plot_annotation (assembly of interactive plots)
library(eirasagree) # version 5.1 # script for visual and statistical assessments of accuracy, precision and agreement
library(kableExtra) # rendering tables in HTML format
setwd("C:/Users/ThomasWeissensteiner/OneDrive/Documents/portfolio/WWchallenges/WWchallenge_56")
# Data were downloaded manually into the working directory from https://github.com/VIS-SIG/Wonderful-Wednesdays/tree/master/data/2024/2024-11-13
exampleData <- read_excel("goniometer.xlsx")
1. Distribution of measured values
1.1. Raw data
The data is based on a paper by Eliasziw et al., 1994
1. One therapist uses two goniometers (i.e., electro
vs. manual) to measure knee joint angles. Twenty-nine subjects (n=29)
were measured three consecutive times (m=3) on each goniometer (t=2) at
one joint position (i.e., full passive extension).
# - Chunk_1 - #
# - Requires R-libraries and object exampleData (r libraries and data) - #
## Data structure
exampleData %>% .[1:10,]## # A tibble: 10 × 4
## id time rater y
## <dbl> <chr> <chr> <dbl>
## 1 1 1 manual -2
## 2 1 1 electro 2
## 3 1 2 manual 0
## 4 1 2 electro 1
## 5 1 3 manual 1
## 6 1 3 electro 1
## 7 2 1 manual 16
## 8 2 1 electro 12
## 9 2 2 manual 16
## 10 2 2 electro 14## Plot of measurements in individual subjects vs replicate number
p1 <-
exampleData %>%
mutate(
time = as.numeric(time),
rater = factor(rater),
id = factor(id)
) %>%
ggplot(
aes(
x = time,
y = y, rater,
color = id,
data_id = id)
) +
geom_line(color = "lightgrey") +
geom_point_interactive(
aes(
tooltip = id)
) +
facet_wrap(.~rater) +
scale_x_continuous(
breaks = 1:3
) +
xlab("Replicate") +
ylab("Measurement") +
labs(
title = "Consecutive replicate measurements in individual patients",
caption = "Hover over a point to highlight an individual patient"
) +
theme_light() +
theme(
plot.title = element_text(hjust = 0.5),
strip.background = element_blank(),
strip.text.x = element_text_interactive(
colour="black",
size = 10
),
legend.position = "none",
plot.caption = element_text(hjust=0.5)
)
girafe(
ggobj = p1,
options = list(
opts_hover(css = ''), # CSS code of selected line
opts_hover_inv(css = "opacity:0.05;"), # CSS code of all other lines
opts_sizing(rescale = FALSE)
)
)
* Manual measurements are spread over a greater range and tend
to be somewhat higher than the values obtained with the electro
method
* Raw data show no obvious trends across repeated measurements for
either method
Note: The R packages ggiraph and patchwork work well together for generating interactive multiplots. However, there is presently no error-free method for making interactive line plots with ggiraph. As a substitute, I used geom_point_interactive to add colour to the lines
1. 2. Means and their variance across consecutive replicate measurements
# - Chunk_3c - #
# - Requires R-libraries (r libraries and data) and object exampleData_int (chunk_1) - #
# Standard error function
std <- function(x) sd(x)/sqrt(length(x))
# Calculate intra-subject replicate means
replicates <- exampleData %>% .$time %>% unique
exampleData_replMean <-
exampleData %>%
pivot_wider(
names_from = "time",
values_from = "y") %>%
rowwise() %>%
mutate(
replicate_Mean = mean(
c_across(
all_of(replicates)
)
),
replicate_SD = sd(
c_across(
all_of(replicates)
)
),
replicate_SE = std(
c_across(
all_of(replicates)
)
)
) %>%
pivot_longer(
cols = replicates,
names_to = "time",
values_to = "y"
)
## Table with summary statistics of the original sample data
# Generate summary table
exampleData_replMean_summary <-
exampleData_replMean %>%
{ rbind(
group_by (., rater, time) %>% # for each method and replicate, calculate mean and SE across subjects
summarise(
Mean = mean(y),
SD = sd(y),
SE = std(y)
),
group_by (., rater) %>% # for each method, calculate means across subjects for mean and SE across replicates
summarise(
Mean = mean(replicate_Mean),
SD = sd(replicate_Mean),
SE = std(replicate_Mean)
)
)
} %>%
mutate(
errorbar_pos = rep(- 0.06, n()),
) %>% # defines vertical position of errorbars and their mean in the facets
replace(is.na(.), "replicate mean")
exampleData_replMean_summary %>%
select(!errorbar_pos) %>%
kable(., "html") %>%
kable_styling(
bootstrap_options = c(
"striped", "hover", "responsive", full_width = FALSE
)
) | rater | time | Mean | SD | SE |
|---|---|---|---|---|
| electro | 1 | 0.0000000 | 7.086204 | 1.3158750 |
| electro | 2 | 0.0344828 | 7.480846 | 1.3891582 |
| electro | 3 | 0.1034483 | 7.148151 | 1.3273783 |
| manual | 1 | 1.5862069 | 7.218613 | 1.3404628 |
| manual | 2 | 1.4827586 | 7.471951 | 1.3875064 |
| manual | 3 | 1.2413793 | 7.452972 | 1.3839822 |
| electro | replicate mean | 0.0459770 | 7.109729 | 0.7622430 |
| manual | replicate mean | 1.4367816 | 7.263108 | 0.7786869 |
## Plot distribution, mean and SE of adjusted measurements
# Change over consecutive replicates
p2 <-
exampleData_replMean %>%
bind_rows(
exampleData_replMean %>%
filter(time == "1") %>% # select a unique row for each subject's replicate mean (avoid generating one row for each replicate)
mutate(
.keep = "none",
id = id,
rater = rater,
y = replicate_Mean)
) %>%
select(id, time, rater, y) %>% # selected columns must have unique rows across all column pairs
replace(is.na(.), "replicate mean") %>%
ggplot(
aes(x = y)
) +
geom_density( # show smoothed distribution of measurement values
aes(y*0.7,
fill = rater
),
alpha = 0.3,
color = "grey",
position = position_nudge(y = 0.02)
) +
geom_point_interactive( # show data points of all subjects
aes(
y = -0.02, # defines vertical position of the data point rows in the facets
color = rater,
data_id = id,
tooltip = id
),
size = 3,
alpha = 0.5
) +
geom_errorbarh ( # show SE bar
data = exampleData_replMean_summary,
inherit.aes = FALSE,
aes(
xmin = Mean - SE,
xmax = Mean + SE,
y = errorbar_pos - 0.025,
color = rater
),
height = 0.025 # height of bar caps
) +
geom_errorbarh ( # show SD bar
data = exampleData_replMean_summary,
inherit.aes = FALSE,
aes(
xmin = Mean - SD,
xmax = Mean + SD,
y = errorbar_pos,
color = rater
),
height = 0.04 # height of bar caps
) +
geom_point( # add means to SE bars
data = exampleData_replMean_summary,
aes(
x = Mean,
y = errorbar_pos,
color = rater
),
size = 2
) +
facet_wrap(
.~ time,
dir = "v",
strip.position = "left",
nrow = 4
) +
labs(
title = "Inter-method comparison of mean, SD and SE \nof consecutive replicate measurements",
y = "Replicate",
x = "Measurement",
tag = "A"
) +
scale_fill_manual( values = c("blue", "red")) +
scale_color_manual(values = c("blue", "red")) +
guides(
fill = guide_legend(title="Method"),
color = guide_legend(title="Method")
) +
theme_light() +
theme(
plot.title = element_text(
size = 12,
hjust = 0.5),
axis.text.y = element_blank(),
strip.background = element_rect(
fill = "white",
color = "darkgrey"
),
axis.ticks.y = element_blank(),
strip.text.y = element_text(
colour="black",
size = 12
)
)
# Variance vs. value (replicate mean)
p3 <-
exampleData_replMean %>%
filter(time == "1") %>% # select a unique row for each subject's replicate mean (avoid generating one row for each replicate)
ggplot(
aes(
x = replicate_Mean,
y = replicate_SD,
colour = rater,
data_id = id,
tooltip = id
)
) +
geom_point_interactive( # show data points of for replicate means of all subjects
size = 3,
alpha = 0.5
) +
labs(
title = "Intra-subject replicate standard deviation vs. mean",
caption = "Hover over points to highlight individual patients",
x = "replicate mean",
y = "SD",
tag = "B"
) +
scale_color_manual(values = c("blue", "red")) +
theme_light() +
theme(
plot.title = element_text(
size = 12,
hjust = 0.5),
plot.caption = element_text(
size = 9,
hjust = 0.5),
legend.position = "none"
)
girafe(
ggobj =
p2 +
plot_spacer() +
p3 +
plot_layout(
nrow = 3,
heights = unit(c(18, 0.2, 5), "cm"
)
) +
plot_annotation(
title = "Distributions and means of original measurements\n",
theme = theme(
plot.title = element_text(
size = 14,
hjust = 0.5
)
)
),
height_svg = 12,
width_svg = 6,
options = list(
opts_hover(css = "stroke:black; stroke-width:2px;"), # CSS code of selected line
opts_hover_inv(css = "opacity:0.1;"), # CSS code of all other lines
opts_sizing(rescale = FALSE) # allows to zoom in and out of saved html image
)
)
* Means of manual measurements are higher at all time points (A)
* The distribution of measurement values appears to be similar and
approximately normal for both methods at all time points (B) * Variance
across replicates has a similar dependency on replicate means for both
methods, being greatest around mean ~ 2.5
The standard deviations and errors for the replicate means in this
figure are not entirely correct: Overall variance should account for
both within-sample and between-sample variation, but shown here is only
the latter which is an underestimate. Because replicates within a sample
are also correlated, more appropriate methods for their calculations
would be * Mixed-effects model or variance components analysis (e.g. R
package lme4) * Nested ANOVA (Analysis of Variance)
1. 3. Distribution of differences between measurements taken at various timepoints
# - Chunk_4d - #
# - Requires R-libraries (r libraries and data) and object exampleData_int (chunk_1) - #
exampleData_methodDiff <-
exampleData_replMean %>%
select(id, time, rater, y) %>% # selected columns must have unique rows across all column pairs
pivot_wider(
names_from = rater,
values_from = y
) %>%
mutate(
time = factor(time),
difference = manual - electro
)
replicate_summary <-
exampleData_replMean %>%
select(id, time, rater, y) %>%
pivot_wider(
names_from = rater,
values_from = y
) %>%
split(.$id) %>%
lapply(
function(x) summarise(x,
replMean_manual = mean(manual),
replSD_manual= sd(manual),
replMean_electro = mean(electro),
replSD_electro= sd(electro),
)
) %>%
list_rbind(names_to = "id") %>%
mutate(
id = as.double(id),
replMean_diff = replMean_manual - replMean_electro,
replMean_diff_SD = sqrt(replSD_manual^2 + replSD_electro^2),
replMean_diff_CI =
qt(0.975, df = length(replicates)+ 1 # 3 replicates = 4 degrees of freedom: t-distribution = 2.766
) *
replMean_diff_SD/sqrt(length(replicates)
)
)
# Histogram plot of measurement differences
p4 <-
exampleData_methodDiff %>%
select(id, time, difference) %>%
bind_rows(
replicate_summary %>%
mutate(time = NA) %>%
select(id, time, replMean_diff) %>%
mutate(
.keep = "none",
id = id,
time = time,
difference = replMean_diff
)
) %>%
select(id, time, difference) %>%
mutate(
time = as.character(time),
) %>%
replace(is.na(.), "replicate mean") %>%
ggplot(
aes(difference)
) +
geom_histogram(
bins = 12,
position = "stack",
color = "darkgrey"
) +
facet_wrap(
.~ time,
dir = "v",
strip.position = "left",
nrow = 4
) +
labs(
title = "Histogram plots",
x = "\nmanual - electro",
y = "count"
) +
theme_light() +
theme(
plot.title = element_text(
size = 12,
hjust = 0.5),
axis.text.y = element_blank(),
strip.background = element_rect(
fill = "white",
color = "darkgrey"
),
strip.placement = "outside",
strip.text.y = element_text(
colour="black",
size = 12
),
axis.ticks = element_blank(),
axis.text.x = element_text(vjust = -2)
)
# qq-plot of measurement differences (unscaled)
p5 <-
exampleData_methodDiff %>%
select(id, time, difference) %>%
bind_rows(
replicate_summary %>%
mutate(time = NA) %>%
select(id, time, replMean_diff) %>%
mutate(
.keep = "none",
id = id,
time = time,
difference = replMean_diff
)
) %>%
select(id, time, difference) %>%
mutate(
time = as.character(time),
) %>%
replace(is.na(.), "replicate mean") %>%
ggplot(
aes(sample = difference)
) +
geom_qq(
dparams = list(
mean = mean(exampleData_methodDiff$difference),
sd = sd(exampleData_methodDiff$difference)
),
size = 3,
color = "black",
alpha = 0.3,
) +
geom_abline(
linetype = "dashed",
slope = 1,
intercept = 0
) +
facet_wrap(
.~ time,
dir = "v",
strip.position = "left",
nrow = 4
) +
labs(
title = "qq-plots",
x = "\nmanual - electro (theoretical data)",
y = "manual - electro"
) +
scale_x_continuous(
breaks = seq(-2.5, 7.5, 0.5)
) +
theme_light() +
theme(
plot.title = element_text(
size = 12,
hjust = 0.5),
strip.background = element_blank(),
strip.text.x = element_blank(),
axis.ticks = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_text(vjust = 6),
axis.title.x = element_text(vjust = 14)
)
p4 + p5 +
plot_layout(
ncol = 2,
heights = c(12, 12)
) +
plot_annotation(
title = "Distribution of inter-method differences \nfor consecutive replicate measurements and their means",
theme = theme(
plot.title = element_text(
size = 14,
hjust = 0.5
)
)
)
- Differences between values obtained by the two methods at the three
time points appear to have normal distributions
2. Method agreement
2. 1. Bland-Altman plots for replicate means and individual time points
# - Chunk_5 - #
# - Requires R-libraries (r libraries and data) and object exampleData_int (chunk_1) - #
## Bland-Altman plot with differences for each individual replicate (original measurement values)
exampleData_BA <-
exampleData_replMean %>%
select(id, time, rater, y) %>%
pivot_wider(
names_from = "rater",
values_from = "y"
) %>%
mutate(
method_mean = (manual + electro)/2,
method_diff = manual - electro
)
summaryData_BA <-
exampleData_BA %>%
group_by(time) %>%
summarise(
Mean = mean(method_diff),
SD = sd(method_diff),
LoA_up = Mean + 1.96* SD,
LoA_low = Mean - 1.96* SD
)
p6 <-
exampleData_BA %>%
ggplot(
aes(
x= method_mean,
y= method_diff,
colour = time,
data_id = id,
tooltip = id
)
) +
geom_point_interactive(
size = 3,
alpha = 0.5
) +
geom_ysidedensity(
aes(fill = time),
alpha = 0.1
) +
geom_hline(
data = summaryData_BA,
aes(
yintercept = summaryData_BA$Mean,
color = time
),
show.legend = FALSE
) +
geom_hline(
data = summaryData_BA,
aes(
yintercept = summaryData_BA$LoA_up,
color = time
),
alpha = 0.7,
show.legend = FALSE
) +
geom_hline(
data = summaryData_BA,
aes(
yintercept = summaryData_BA$LoA_low,
color = time
),
alpha = 0.7,
show.legend = FALSE
) +
scale_color_manual(
values = c("#F8766D", "#00BA38", "#619CFF", "darkgrey")
) +
scale_y_continuous(
limits = c(
-5, # min(exampleData_BA$method_diff),
7 # max(exampleData_BA$method_diff)
),
breaks = seq(-7, 7, 1)
) +
labs(
tag = "B",
x = "(manual+electro)/2",
y = "manual-electro"
) +
guides(
theme(
legend.text = element_text(size = 13),
color = guide_legend(title="Mean and limits of agreement for individual replicate datasets", position = "top"),
fill = guide_legend(title="Mean and limits of agreement for individual replicate datasets", position = "top")
),
color = FALSE
) +
theme_light() +
theme_ggside_void()
## Bland-Altman plot with differences between inter-replicate means (original measurement values)
# Calculate inter-replicate mean and SD for each sample
exampleData_BA_means <-
exampleData_BA %>%
select(id, time, method_mean, method_diff) %>%
split(.$id) %>%
lapply(
function(x) summarise(x,
Mean_methodMean = mean(method_mean),
SD_methodMean = sd(method_mean),
Mean_methodDiff = mean(method_diff),
SD_methodDiff = sd(method_diff)
)
) %>%
list_rbind(names_to = "id")
# Calculate summary stats of inter-replicate means and LoAs
SD_means <-
function(x) {
sqrt(sum(x^2)/(length(x)-1))
}
summaryData_BA_means <-
exampleData_BA_means %>%
summarise(
Mean = mean(Mean_methodDiff),
SD = SD_means(Mean_methodDiff)
) %>%
mutate(
LoA_up = Mean + 1.96*SD,
LoA_low = Mean - 1.96*SD
)
# Bland-Altman plot
p7 <-
exampleData_BA_means %>%
ggplot(
aes(
x= Mean_methodMean,
y= Mean_methodDiff,
data_id = id,
tooltip = id
)
) +
geom_point_interactive(
color = "black",
alpha = 0.3,
size = 3
) +
geom_ysidedensity(
inherit.aes = FALSE,
aes(y = Mean_methodDiff)
) +
geom_hline(
data = summaryData_BA_means,
aes(yintercept = summaryData_BA_means$Mean)
) +
geom_hline(
data = summaryData_BA_means,
aes(yintercept = summaryData_BA_means$LoA_up),
) +
geom_hline(
data = summaryData_BA_means,
aes(yintercept = summaryData_BA_means$LoA_low)
) +
scale_y_continuous(
limits = c(
-5, # min(exampleData_BA$method_diff), # same y-scale limits as p6
7 # max(exampleData_BA$method_diff)
),
breaks = seq(-7, 7, 1)
) +
labs(
tag = "A",
title = "Mean and limits of agreement for averaged replicate measurements ",
x = "(manual+electro)/2",
y = "manual-electro",
caption =
summaryData_BA_means %>%
round(digits = 2) %>%
{ paste("Mean = ", .[1],
"\n upper LoA = ", .[3],
"\n lower LoA = ", .[4]
)
}
) +
theme_light() +
theme(
plot.title = element_text(
size = 11,
hjust = 0.3
),
legend.position = "none"
) +
theme_ggside_void()
girafe(
ggobj = p7 + plot_spacer() + p6 +
plot_layout(
nrow = 3,
heights = c(4, 0.2, 4)
) +
plot_annotation(
title =
"Bland-Altman plots with replicate measurements\n",
caption = "Hover over points to highlight individual patients",
theme = theme(
plot.title = element_text(
face = "bold",
size = 12,
hjust = 0.5,
vjust = 3
),
plot.caption = element_text(
size = 9,
hjust = 0.5
)
)
),
# height_svg = 9, # svg = 9 for saving to individual image to html, omit for knitting
options = list(
opts_hover(css = ''), # CSS code of selected line
opts_hover_inv(css = "opacity:0.05;"), # CSS code of all other lines
opts_sizing(rescale = FALSE)
)
)
2. 2. Probability of Agreement plot
This method 2 can take replicate measurements as
input.
## Probability agreement (NT Stevens, 2019)
# Script customised for analysis with own dataset
# Choose acceptable difference for which the probability of agreement should be evalualated
delta <- 3
# Reformat the input dataframe (check variable names and formats)
rr <-
exampleData %>%
rename(
Subject = id,
Repeat = time,
MS = rater,
Measurements = y
) %>%
mutate(
MS =
case_when(
MS == "manual" ~ 1,
MS == "electro" ~ 2
),
across(
.cols = everything(),
as.double
)
)
# Calculate number of subjects
n_subj <- rr$Subject %>% max
############################################################################
## Relevant functions required for the Probability of Agreement analysis ##
############################################################################
## Author: Nathaniel T. Stevens
## Assistant Professor, University of Waterloo
## Date: March 2019
## Analysis script for the probabilioty of agreement methodology
## https://github.com/vandainacio/Comparison-of-Six-Agreement-Methods-for-Continuous-Repeated-Measures-Data/blob/master/PAgreement_Analysis_Script.R
## Note that the dataframe must have four columns with headers 'Subject', 'MS', 'Repeat' and 'Measurements',
## and N = sum{m_ij} rows. The 'Subject' column contains subject codes and the entries of 'MS' must be 1's and
## 2's indicating which measurement system a given row corresponds to. The entries of 'Repeat' indicate which
## repeated measurement a given row corrsponds to, and the entries of the 'Measurements' column correspond to
## the repeated measurements by the given system on each of the n subjects. Note that the observations must be
## ordered by subject with repeated measurements appearing consecutively with MS1 measurements appearing first
## and MS2 second.
#############################
## Log-Likelihood Function ##
#############################
log.likelihood <- function(param, n, r1vec, r2vec, data){
# Calculate the log-likelihood contribution for each subject and then sum them all
mu <- param[1]
alpha <- param[2]
beta <- param[3]
sp <- param[4]
s1 <- param[5]
s2 <- param[6]
l <- rep(0, n)
for(i in 1:n){
r1 <- r1vec[i]
r2 <- r2vec[i]
Yi1 <- subset(data, subset = (data$Subject==i & data$MS==1))$Measurements
Yi2 <- subset(data, subset = (data$Subject==i & data$MS==2))$Measurements
A <- sum((Yi1-rep(mu, r1))^2)
B <- sum(Yi1-rep(mu, r1))^2
C <- sum((Yi2-rep(alpha + beta*mu, r2))^2)
D <- sum(Yi2-rep(alpha + beta*mu, r2))^2
E <- sum(Yi1-rep(mu, r1)) * sum(Yi2-rep(alpha + beta*mu, r2))
# Log-likelihood contribution for subject i
l[i] <- -(r1 / 2 + r2 / 2) * log(2 * pi) - r1 * log(s1) - r2 * log(s2) - log1p(sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) / 2 - A / s1 ^ 2 / 2 + sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * B / s1 ^ 4 / 2 - C / s2 ^ 2 / 2 + sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * D / s2 ^ 4 / 2 + sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * E / s1 ^ 2 / s2 ^ 2
}
# The log-likelihood function:
l <- sum(l)
return(-l) # Return negative because optim() is a minimizer
}
##############################
## Log-Likelihood Gradients ##
##############################
gradient <- function(param, n, r1vec, r2vec, data){
# Calculate the gradients for each parameter / subject combination and sum over all subjects
# for the gradients for each parameter.
mu <- param[1]
alpha <- param[2]
beta <- param[3]
sp <- param[4]
s1 <- param[5]
s2 <- param[6]
g <- matrix(0, nrow = 6, ncol = n)
for(i in 1:n){
r1 <- r1vec[i]
r2 <- r2vec[i]
Yi1 <- subset(data, subset = (data$Subject==i & data$MS==1))$Measurements
Yi2 <- subset(data, subset = (data$Subject==i & data$MS==2))$Measurements
A <- sum((Yi1-rep(mu, r1))^2)
C <- sum((Yi2-rep(alpha + beta*mu, r2))^2)
F <- sum(Yi1-rep(mu, r1))
G <- sum(Yi2-rep(alpha + beta*mu, r2))
B <- F^2
D <- G^2
E <- F*G
# Derivative with respect to mu:
g[1,i] <- F / s1 ^ 2 - r1 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * F / s1 ^ 4 + beta * G / s2 ^ 2 - r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 3 * G / s2 ^ 4 - r1 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * G / s1 ^ 2 / s2 ^ 2 - r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * F / s1 ^ 2 / s2 ^ 2
# Derivative with respect to alpha:
g[2,i] <- G / s2 ^ 2 - r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * G / s2 ^ 4 - r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * F / s1 ^ 2 / s2 ^ 2
# Derivative with respect to beta:
g[3,i] <- -beta * r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) / s2 ^ 2 - beta * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 * B / s1 ^ 4 + mu * G / s2 ^ 2 + (-2 * G * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * mu * r2 - 2 * D * beta ^ 3 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + 2 * D * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta) / s2 ^ 4 / 2 + (-F * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * mu * r2 - 2 * E * beta ^ 2 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + E * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2))) / s1 ^ 2 / s2 ^ 2
# Derivative with respect to sp:
g[4,i] <- -sp * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) + sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * B / s1 ^ 4 - sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * B / s1 ^ 4 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) + sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * D / s2 ^ 4 - sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * D / s2 ^ 4 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) + 2 * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * E / s1 ^ 2 / s2 ^ 2 - 2 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * E / s1 ^ 2 / s2 ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)
# Derivative with respect to s1:
g[5,i] <- -r1 / s1 + sp ^ 2 * r1 / s1 ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) + A / s1 ^ 3 + sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * B / s1 ^ 7 * r1 - 2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * B / s1 ^ 5 + sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * D / s2 ^ 4 * r1 / s1 ^ 3 + 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * E / s1 ^ 5 / s2 ^ 2 * r1 - 2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * E / s1 ^ 3 / s2 ^ 2
# Derivative with respect to s2:
g[6,i] <- -r2 / s2 + sp ^ 2 * r2 * beta ^ 2 / s2 ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) + sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * B / s1 ^ 4 * r2 * beta ^ 2 / s2 ^ 3 + C / s2 ^ 3 + sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 4 * D / s2 ^ 7 * r2 - 2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * D / s2 ^ 5 + 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 * E / s1 ^ 2 / s2 ^ 5 * r2 - 2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * E / s1 ^ 2 / s2 ^ 3
}
# The gradients:
g <- apply(X = g, MARGIN = 1, FUN = sum)
return(-g) # Return negative because optim() is a minimizer
}
###############################
## Fisher Information Matrix ##
###############################
fisher.info <- function(param, n, r1vec, r2vec){
# Calculate the Fisher Information Matrix associated with the 6 parameters
# Expected values of the various sums of squares
mu <- param[1]
alpha <- param[2]
beta <- param[3]
sp <- param[4]
s1 <- param[5]
s2 <- param[6]
dldmu2 <- rep(0, n); dldmualpha <- rep(0, n); dldmubeta <- rep(0, n); dldmusp <- rep(0, n); dldmus1 <- rep(0, n); dldmus2 <- rep(0, n); dldalpha2 <- rep(0, n); dldalphabeta <- rep(0, n); dldalphasp <- rep(0, n); dldalphas1 <- rep(0, n); dldalphas2 <- rep(0, n); dldbeta2 <- rep(0, n); dldbetasp <- rep(0, n); dldbetas1 <- rep(0, n); dldbetas2 <- rep(0, n); dldsp2 <- rep(0, n); dldsps1 <- rep(0, n); dldsps2 <- rep(0, n); dlds12 <- rep(0, n); dlds1s2 <- rep(0, n); dlds22 <- rep(0, n);
for(i in 1:n){
r1 <- r1vec[i]
r2 <- r2vec[i]
A <- r1*(sp^2+s1^2)
B <- r1*(r1*sp^2+s1^2)
C <- r2*((beta^2)*(sp^2)+s2^2)
D <- r2*(r2*(beta^2)*(sp^2)+s2^2)
E <- (r1*r2)*beta*sp^2
F <- 0
G <- 0
# Second partial derivatives for each subject
dldmu2[i] <- -r1 / s1 ^ 2 + r1 ^ 2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) / s1 ^ 4 - r2 * beta ^ 2 / s2 ^ 2 + r2 ^ 2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 4 / s2 ^ 4 + 2 * r1 * r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 / s1 ^ 2 / s2 ^ 2
dldmualpha[i] <- -r2 * beta / s2 ^ 2 + r2 ^ 2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 3 / s2 ^ 4 + r1 * r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 2
dldmubeta[i] <- 2 * r1 * beta * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 * F / s1 ^ 4 + (-beta * mu * r2 + G) / s2 ^ 2 - (-sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 3 * mu * r2 ^ 2 - 2 * G * beta ^ 4 * r2 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + 3 * G * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * r2) / s2 ^ 4 - (-sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * mu * r1 * r2 - 2 * G * beta ^ 2 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 * r1 + G * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * r1) / s1 ^ 2 / s2 ^ 2 - (-2 * F * beta ^ 3 * r2 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + 2 * F * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * r2) / s1 ^ 2 / s2 ^ 2
dldmusp[i] <- -2 * r1 * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * F / s1 ^ 4 + 2 * r1 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * F / s1 ^ 4 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) - 2 * r2 * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 3 * G / s2 ^ 4 + 2 * r2 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 * G / s2 ^ 4 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) - 2 * r1 * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * G / s1 ^ 2 / s2 ^ 2 + 2 * r1 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * G / s1 ^ 2 / s2 ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) - 2 * r2 * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * F / s1 ^ 2 / s2 ^ 2 + 2 * r2 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * F / s1 ^ 2 / s2 ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)
dldmus1[i] <- -2 * F / s1 ^ 3 - 2 * r1 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * F / s1 ^ 7 + 4 * r1 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * F / s1 ^ 5 - 2 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 * G / s2 ^ 4 * r1 / s1 ^ 3 - 2 * r1 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * G / s1 ^ 5 / s2 ^ 2 + 2 * r1 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * G / s1 ^ 3 / s2 ^ 2 - 2 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * F / s1 ^ 5 / s2 ^ 2 * r1 + 2 * r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * F / s1 ^ 3 / s2 ^ 2
dldmus2[i] <- -2 * r1 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * F / s1 ^ 4 * r2 * beta ^ 2 / s2 ^ 3 - 2 * beta * G / s2 ^ 3 - 2 * r2 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 5 * G / s2 ^ 7 + 4 * r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 3 * G / s2 ^ 5 - 2 * r1 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 * G / s1 ^ 2 / s2 ^ 5 * r2 + 2 * r1 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * G / s1 ^ 2 / s2 ^ 3 - 2 * r2 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 4 * F / s1 ^ 2 / s2 ^ 5 + 2 * r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * F / s1 ^ 2 / s2 ^ 3
dldalpha2[i] <- -r2 / s2 ^ 2 + r2 ^ 2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 / s2 ^ 4
dldalphabeta[i] <- -r2 * mu / s2 ^ 2 - (-sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * mu * r2 ^ 2 - 2 * G * beta ^ 3 * r2 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + 2 * G * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * r2) / s2 ^ 4 - (-2 * F * beta ^ 2 * r2 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + F * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * r2) / s1 ^ 2 / s2 ^ 2
dldalphasp[i] <- -2 * r2 * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * G / s2 ^ 4 + 2 * r2 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * G / s2 ^ 4 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) - 2 * r2 * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * F / s1 ^ 2 / s2 ^ 2 + 2 * r2 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * F / s1 ^ 2 / s2 ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)
dldalphas1[i] <- -2 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * G / s2 ^ 4 * r1 / s1 ^ 3 - 2 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * F / s1 ^ 5 / s2 ^ 2 * r1 + 2 * r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * F / s1 ^ 3 / s2 ^ 2
dldalphas2[i] <- -2 * G / s2 ^ 3 - 2 * r2 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 4 * G / s2 ^ 7 + 4 * r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * G / s2 ^ 5 - 2 * r2 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 * F / s1 ^ 2 / s2 ^ 5 + 2 * r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * F / s1 ^ 2 / s2 ^ 3
dldbeta2[i] <- -(-2 * beta ^ 2 * r2 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * r2) / s2 ^ 2 - r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * (-4 * beta ^ 2 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2))) / s2 ^ 2 * B / s1 ^ 4 - r2 * mu ^ 2 / s2 ^ 2 + (8 * G * beta ^ 3 * r2 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 * mu - 2 * D * beta ^ 2 * r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * (-4 * beta ^ 2 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2))) / s2 ^ 2 - 8 * D * beta ^ 2 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 - 8 * G * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * mu * r2 + 2 * D * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) + 2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * mu ^ 2 * r2 ^ 2) / s2 ^ 4 / 2 + (4 * F * beta ^ 2 * r2 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 * mu - 2 * E * beta * r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * (-4 * beta ^ 2 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2))) / s2 ^ 2 - 4 * E * beta * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 - 2 * r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * mu * F) / s1 ^ 2 / s2 ^ 2
dldbetasp[i] <- -2 * beta * r2 * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) / s2 ^ 2 + 2 * beta * r2 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) - 4 * beta * r2 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 * B / s1 ^ 4 + 4 * beta * r2 * sp ^ 5 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 / s2 ^ 2 * B / s1 ^ 4 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) + (-4 * G * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * mu * r2 + 4 * G * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * mu * r2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) - 8 * D * beta ^ 3 * r2 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + 8 * D * beta ^ 3 * r2 * sp ^ 5 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 / s2 ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) + 4 * D * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta - 4 * D * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) / s2 ^ 4 / 2 + (-2 * F * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * mu * r2 + 2 * F * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * mu * r2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) - 8 * E * beta ^ 2 * r2 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + 8 * E * beta ^ 2 * r2 * sp ^ 5 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 / s2 ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) + 2 * E * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) - 2 * E * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) / s1 ^ 2 / s2 ^ 2
dldbetas1[i] <- -2 * beta * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 * r1 / s1 ^ 3 - 4 * beta * r2 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 / s2 ^ 2 * B / s1 ^ 7 * r1 + 4 * beta * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 * B / s1 ^ 5 + (-4 * G * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * mu * r2 * r1 / s1 ^ 3 - 8 * D * beta ^ 3 * r2 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 / s2 ^ 2 * r1 / s1 ^ 3 + 4 * D * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * r1 / s1 ^ 3) / s2 ^ 4 / 2 + (-2 * F * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * mu * r2 * r1 / s1 ^ 3 - 8 * E * beta ^ 2 * r2 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 / s2 ^ 2 * r1 / s1 ^ 3 + 2 * E * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * r1 / s1 ^ 3) / s1 ^ 2 / s2 ^ 2 - 2 * (-F * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * mu * r2 - 2 * E * beta ^ 2 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + E * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2))) / s1 ^ 3 / s2 ^ 2
dldbetas2[i] <- -2 * beta ^ 3 * r2 ^ 2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 5 + 2 * beta * r2 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) / s2 ^ 3 - 4 * beta ^ 3 * r2 ^ 2 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 / s2 ^ 5 * B / s1 ^ 4 + 2 * beta * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 3 * B / s1 ^ 4 - 2 * mu * G / s2 ^ 3 + (-4 * G * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 4 * mu * r2 ^ 2 / s2 ^ 3 - 8 * D * beta ^ 5 * r2 ^ 2 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 / s2 ^ 5 + 8 * D * beta ^ 3 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 3) / s2 ^ 4 / 2 - 2 * (-2 * G * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * mu * r2 - 2 * D * beta ^ 3 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + 2 * D * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta) / s2 ^ 5 + (-2 * F * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 * mu * r2 ^ 2 / s2 ^ 3 - 8 * E * beta ^ 4 * r2 ^ 2 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 / s2 ^ 5 + 6 * E * beta ^ 2 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 3) / s1 ^ 2 / s2 ^ 2 - 2 * (-F * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * mu * r2 - 2 * E * beta ^ 2 * r2 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + E * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2))) / s1 ^ 2 / s2 ^ 3
dldsp2[i] <- -(r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) + 2 * sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 + 1 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * B / s1 ^ 4 - 5 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * B / s1 ^ 4 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) + 4 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * B / s1 ^ 4 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) ^ 2 + 1 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * D / s2 ^ 4 - 5 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * D / s2 ^ 4 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) + 4 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 2 * D / s2 ^ 4 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) ^ 2 + 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * E / s1 ^ 2 / s2 ^ 2 - 10 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * E / s1 ^ 2 / s2 ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) + 8 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * beta * E / s1 ^ 2 / s2 ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) ^ 2
dldsps1[i] <- 2 * sp * r1 / s1 ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) - 2 * sp ^ 3 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * r1 / s1 ^ 3 + 4 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * B / s1 ^ 7 * r1 - 4 * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * B / s1 ^ 5 - 4 * sp ^ 5 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * B / s1 ^ 7 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) * r1 + 4 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * B / s1 ^ 5 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) + 4 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * D / s2 ^ 4 * r1 / s1 ^ 3 - 4 * sp ^ 5 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 2 * D / s2 ^ 4 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) * r1 / s1 ^ 3 + 8 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * E / s1 ^ 5 / s2 ^ 2 * r1 - 4 * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * E / s1 ^ 3 / s2 ^ 2 - 8 * sp ^ 5 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * beta * E / s1 ^ 5 / s2 ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) * r1 + 4 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * E / s1 ^ 3 / s2 ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)
dldsps2[i] <- 2 * sp * r2 * beta ^ 2 / s2 ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) - 2 * sp ^ 3 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * r2 * beta ^ 2 / s2 ^ 3 + 4 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * B / s1 ^ 4 * r2 * beta ^ 2 / s2 ^ 3 - 4 * sp ^ 5 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * B / s1 ^ 4 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) * r2 * beta ^ 2 / s2 ^ 3 + 4 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 4 * D / s2 ^ 7 * r2 - 4 * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * D / s2 ^ 5 - 4 * sp ^ 5 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 4 * D / s2 ^ 7 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) * r2 + 4 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * D / s2 ^ 5 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) + 8 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 * E / s1 ^ 2 / s2 ^ 5 * r2 - 4 * sp / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * E / s1 ^ 2 / s2 ^ 3 - 8 * sp ^ 5 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 3 * E / s1 ^ 2 / s2 ^ 5 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2) * r2 + 4 * sp ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * E / s1 ^ 2 / s2 ^ 3 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)
dlds12[i] <- r1 / s1 ^ 2 - 3 * sp ^ 2 * r1 / s1 ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) + 2 * sp ^ 4 * r1 ^ 2 / s1 ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 - 3 * A / s1 ^ 4 + 4 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * B / s1 ^ 10 * r1 ^ 2 - 11 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * B / s1 ^ 8 * r1 + 10 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * B / s1 ^ 6 + 4 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 2 * D / s2 ^ 4 * r1 ^ 2 / s1 ^ 6 - 3 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * D / s2 ^ 4 * r1 / s1 ^ 4 + 8 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * beta * E / s1 ^ 8 / s2 ^ 2 * r1 ^ 2 - 14 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * E / s1 ^ 6 / s2 ^ 2 * r1 + 6 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * E / s1 ^ 4 / s2 ^ 2
dlds1s2[i] <- 2 * sp ^ 4 * r1 / s1 ^ 3 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * r2 * beta ^ 2 / s2 ^ 3 + 4 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * B / s1 ^ 7 * r1 * r2 * beta ^ 2 / s2 ^ 3 - 4 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * B / s1 ^ 5 * r2 * beta ^ 2 / s2 ^ 3 + 4 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 4 * D / s2 ^ 7 * r1 / s1 ^ 3 * r2 - 4 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * D / s2 ^ 5 * r1 / s1 ^ 3 + 8 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 3 * E / s1 ^ 5 / s2 ^ 5 * r1 * r2 - 4 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta * E / s1 ^ 5 / s2 ^ 3 * r1 - 4 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 * E / s1 ^ 3 / s2 ^ 5 * r2 + 4 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * E / s1 ^ 3 / s2 ^ 3
dlds22[i] <- r2 / s2 ^ 2 - 3 * sp ^ 2 * r2 * beta ^ 2 / s2 ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) + 2 * sp ^ 4 * r2 ^ 2 * beta ^ 4 / s2 ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 + 4 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * B / s1 ^ 4 * r2 ^ 2 * beta ^ 4 / s2 ^ 6 - 3 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * B / s1 ^ 4 * r2 * beta ^ 2 / s2 ^ 4 - 3 * C / s2 ^ 4 + 4 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 6 * D / s2 ^ 10 * r2 ^ 2 - 11 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 4 * D / s2 ^ 8 * r2 + 10 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta ^ 2 * D / s2 ^ 6 + 8 * sp ^ 6 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 5 * E / s1 ^ 2 / s2 ^ 8 * r2 ^ 2 - 14 * sp ^ 4 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 * E / s1 ^ 2 / s2 ^ 6 * r2 + 6 * sp ^ 2 / (1 + sp ^ 2 * (r1 / s1 ^ 2 + r2 * beta ^ 2 / s2 ^ 2)) * beta * E / s1 ^ 2 / s2 ^ 4
}
# Second partial derivatives
dldmu2 <- sum(dldmu2)
dldmualpha <- sum(dldmualpha)
dldmubeta <- sum(dldmubeta)
dldmusp <- sum(dldmusp)
dldmus1 <- sum(dldmus1)
dldmus2 <- sum(dldmus2)
dldalpha2 <- sum(dldalpha2)
dldalphabeta <- sum(dldalphabeta)
dldalphasp <- sum(dldalphasp)
dldalphas1 <- sum(dldalphas1)
dldalphas2 <- sum(dldalphas2)
dldbeta2 <- sum(dldbeta2)
dldbetasp <- sum(dldbetasp)
dldbetas1 <- sum(dldbetas1)
dldbetas2 <- sum(dldbetas2)
dldsp2 <- sum(dldsp2)
dldsps1 <- sum(dldsps1)
dldsps2 <- sum(dldsps2)
dlds12 <- sum(dlds12)
dlds1s2 <- sum(dlds1s2)
dlds22 <-sum(dlds22)
# Create the matrix
J <- matrix(0, nrow = 6, ncol = 6)
J[1,] <- c(dldmu2, dldmualpha, dldmubeta, dldmusp, dldmus1, dldmus2)
J[2,] <- c(dldmualpha, dldalpha2, dldalphabeta, dldalphasp, dldalphas1, dldalphas2)
J[3,] <- c(dldmubeta, dldalphabeta, dldbeta2, dldbetasp, dldbetas1, dldbetas2)
J[4,] <- c(dldmusp, dldalphasp, dldbetasp, dldsp2, dldsps1, dldsps2)
J[5,] <- c(dldmus1, dldalphas1, dldbetas1, dldsps1, dlds12, dlds1s2)
J[6,] <- c(dldmus2, dldalphas2, dldbetas2, dldsps2, dlds1s2, dlds22)
return(-J) #change the sign since we want the negative of the expected value instead of the expected value
}
############################
## P(agreement) Gradients ##
############################
delta.method <-
function(delta, p, alpha, beta, s1, s2){
# Function that calculates the gradient vector for theta(p)
D <- matrix(0, nrow = 6)
# Derivative with respect to mu
D[1,] <- 0
# Derivative with respect to alpha
D[2,] <- -pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha - delta - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 2) * sqrt(2) * (s1 ^ 2 + s2 ^ 2) ^ (-0.1e1 / 0.2e1) / 2 + pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha + delta - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 2) * sqrt(2) * (s1 ^ 2 + s2 ^ 2) ^ (-0.1e1 / 0.2e1) / 2
# Derivative with respect to beta
D[3,] <- -pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha - delta - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 2) * sqrt(2) * p * (s1 ^ 2 + s2 ^ 2) ^ (-0.1e1 / 0.2e1) / 2 + pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha + delta - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 2) * sqrt(2) * p * (s1 ^ 2 + s2 ^ 2) ^ (-0.1e1 / 0.2e1) / 2
# Derivative with respect to sp
D[4,] <- 0
# Derivative with respect to s1
D[5,] <- pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha - delta - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 2) * sqrt(2) * (p * beta + alpha - delta - p) * (s1 ^ 2 + s2 ^ 2) ^ (-0.3e1 / 0.2e1) * s1 / 2 - pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha + delta - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 2) * sqrt(2) * (p * beta + alpha + delta - p) * (s1 ^ 2 + s2 ^ 2) ^ (-0.3e1 / 0.2e1) * s1 / 2
# Derivative with respect to s2
D[6,] <- pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha - delta - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 2) * sqrt(2) * (p * beta + alpha - delta - p) * (s1 ^ 2 + s2 ^ 2) ^ (-0.3e1 / 0.2e1) * s2 / 2 - pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha + delta - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 2) * sqrt(2) * (p * beta + alpha + delta - p) * (s1 ^ 2 + s2 ^ 2) ^ (-0.3e1 / 0.2e1) * s2 / 2
return(D)
}
######################
## Modified QQ-Plot ##
######################
mod.qqplot <- function(data, n){
## This function creates a modified QQ-plot as described in Stevens et al. (2015) but for unbalanced repeated measurements
## Inputs:
# data = a dataframe with four columns (Subject, MS, Repeat, Measurements) and N = sum{r_ij} rows. The entries of Measurements
# correspond to the repeated measurements by the given system on each of the n subjects. Note that the observations
# are ordered by subject with repeated measurements appearing consecutively with MS1 measurements appearing first and MS2 second
# n = the number of subjects used in the study
## Sample call: mod.qqplot(data = rr, n = 21)
mm1 <- rep(0, n) # vector of subject-means (MS1)
mm2 <- rep(0, n) # vector of subject-means (MS2)
for(i in 1:n){
mm1[i] <- mean(subset(data, (subset = data$Subject==i & data$MS==1))$Measurements)
mm2[i] <- mean(subset(data, (subset = data$Subject==i & data$MS==2))$Measurements)
}
mu1 <- mean(mm1)
mu2 <- mean(mm2)
sd1 <- sd(mm1)
sd2 <- sd(mm2)
rand.norm1 <- matrix(0, nrow = n, ncol = 50)
rand.norm2 <- matrix(0, nrow = n, ncol = 50)
for (i in 1:50){
rand.norm1[,i] <- rnorm(n, mu1, sd1)
rand.norm2[,i] <- rnorm(n, mu2, sd2)
}
par(mfrow=c(1,2))
# MS1
blackpts1 <- qqnorm(mm1, ylim = range(mm1, rand.norm1), main = "QQ-Plot of MS1 Part Averages vs Standard Normal", pch = 19, cex = 0.75)
for (i in 1:50){
graypts1 <- qqnorm(rand.norm1[,i], plot.it = FALSE)
points(graypts1, col = "grey87", cex = 0.75)
}
points(blackpts1, col = "black", pch = 19, cex = 0.75)
qqline(mm1, col = "red")
# MS2
blackpts2 <- qqnorm(mm2, ylim = range(mm2, rand.norm2), main = "QQ-Plot of MS2 Part Averages vs Standard Normal", pch = 19, cex = 0.75)
for (i in 1:50){
graypts2 <- qqnorm(rand.norm2[,i], plot.it = FALSE)
points(graypts2, col = "grey87", cex = 0.75)
}
points(blackpts2, col = "black", pch = 19, cex = 0.75)
qqline(mm2, col = "red")
}
########################
## Repeatabiltiy Plot ##
########################
repeat.plot <- function(data, n){
## This function creates a repeatability plot as described in Stevens et al. (2015) but for unbalanced repeated measurements
## Inputs:
# data = a dataframe with four columns (Subject, MS, Repeat, Measurements) and N = sum{r_ij} rows. The entries of Measurements
# correspond to the repeated measurements by the given system on each of the n subjects. Note that the observations
# are ordered by subject with repeated measurements appearing consecutively with MS1 measurements appearing first and MS2 second
# n = the number of subjects used in the study
## Sample call: repeat.plot(data = rr, n = 21)
mm1 <- rep(0, n) # vector of subject-means (MS1)
mm2 <- rep(0, n) # vector of subject-means (MS2)
r1vec <- rep(0, n) # vector of repeated measurement numbers (MS1)
r2vec <- rep(0, n) # vector of repeated measurement numbers (MS2)
for(i in 1:n){
mm1[i] <- mean(subset(data, subset = (data$Subject==i & data$MS==1))$Measurements)
mm2[i] <- mean(subset(data, subset = (data$Subject==i & data$MS==2))$Measurements)
r1vec[i] <- max(subset(data, subset = (data$Subject==i & data$MS==1))$Repeat)
r2vec[i] <- max(subset(data, subset = (data$Subject==i & data$MS==2))$Repeat)
}
mm1 <- rep(mm1, r1vec) # vector of MS1 subject-means repeated r_i1 times
res1 <- data[data$MS==1,]$Measurements - mm1
mm2 <- rep(mm2, r2vec) # vector of MS1 subject-means repeated r_i2 times
res2 <- data[data$MS==2,]$Measurements - mm2
par(mfrow=c(1,2))
plot(mm1, res1, pch = 19, cex = 0.5, main = "MS1 Residual Measurements vs Subject-Average", ylab = "Observed Residuals", xlab = "Subject Averages", xlim = c(min(mm1,mm2),max(mm1,mm2)), ylim = c(min(res1,res2),max(res1,res2)))
abline(h = 0, col = "red")
plot(mm2, res2, pch = 19, cex = 0.5, main = "MS2 Residual Measurements vs Subject-Average", ylab = "Observed Residuals", xlab = "Subject Averages", xlim = c(min(mm1,mm2),max(mm1,mm2)), ylim = c(min(res1,res2),max(res1,res2)))
abline(h = 0, col = "red")
}
####################################
## Repeated Measures Scatter Plot ##
####################################
rep.scatter.plot <- function(data, n){
Y1 <- subset(data, subset = data$MS==1)$Measurements
Y2 <- subset(data, subset = data$MS==2)$Measurements
sub_means1 <- rep(0, n)
sub_means2 <- rep(0, n)
for(i in 1:n){
sub_i1 <- subset(data, subset = (data$Subject==i & data$MS==1))
sub_means1[i] <- mean(sub_i1$Measurements)
sub_i2 <- subset(data, subset = (data$Subject==i & data$MS==2))
sub_means2[i] <- mean(sub_i2$Measurements)
}
par(mfrow = c(1,1))
plot(sub_means1, sub_means2, ylim = c(min(sub_means1, sub_means2), max(sub_means1, sub_means2)), xlim = c(min(sub_means1, sub_means2), max(sub_means1, sub_means2)), xlab = "MS1 Subject Averages", ylab = "MS2 Subject Averages", main = "Repeated Measures Scatterplot", pch = 16)
abline(a = 0, b = 1, col = "red", lty = 2)
}
###########################
## Ordinary Scatter Plot ## (this only works when number of measurements per subject is the same across devices)
###########################
scatter.plot <- function(data, n){
Y1 <- subset(data, subset = data$MS==1)$Measurements
Y2 <- subset(data, subset = data$MS==2)$Measurements
par(mfrow = c(1,1))
plot(Y1, Y2, ylim = c(min(Y1, Y2), max(Y1, Y2)), xlim = c(min(Y1, Y2), max(Y1, Y2)), xlab = "MS1 Readings", ylab = "MS2 Readings", main = "Scatterplot", pch = 16)
abline(a = 0, b = 1, col = "red", lty = 2)
}
###########################
## P(agreement) Function ##
###########################
prob.agree <- function(n, delta, data){
## This function conducts the probability of agreement analysis for the
## comparison of two homoscedastic measurement systems with possibly
## unbalanced measurements.
## Sample call: prob.agree(n = 21, delta = 5, data = rr)
## Inputs:
# * n = the number of parts/subjects being measured by each system
# * delta = the the maximum allowable difference between measurement by two systems
# * data = a dataframe with four columns (Subject, MS, Repeat, Measurements) and N = sum{r_ij} rows. The entries of Measurements
# correspond to the repeated measurements by the given system on each of the n subjects. Note that the observations
# are ordered by subject with repeated measurements appearing consecutively with MS1 measurements appearing first and MS2 second
# Make sure the correct inputs are inputted
if(is.null(n) == TRUE){
print("You must enter the number of subjects, n.")
}else if(is.null(delta) == TRUE){
print("You must define the clinically acceptable difference by inputing delta.")
}else if(is.null(data) == TRUE){
print("You must enter data for both measurement systems.")
}else{ # Function inputs are correct, therefore we can proceed
# Get initial guesses for parameters (needed for likelihood maximization)
Y1 <- subset(data, subset = data$MS==1)$Measurements
Y2 <- subset(data, subset = data$MS==2)$Measurements
sub_means1 <- rep(0, n)
sub_means2 <- rep(0, n)
sub_vars1 <- rep(0, n)
sub_vars2 <- rep(0, n)
r1vec <- rep(0, n)
r2vec <- rep(0, n)
for(i in 1:n){
sub_i1 <- subset(data, subset = (data$Subject==i & data$MS==1))
sub_means1[i] <- mean(sub_i1$Measurements)
sub_vars1[i] <- var(sub_i1$Measurements)
r1vec[i] <- max(sub_i1$Repeat)
sub_i2 <- subset(data, subset = (data$Subject==i & data$MS==2))
sub_means2[i] <- mean(sub_i2$Measurements)
sub_vars2[i] <- var(sub_i2$Measurements)
r2vec[i] <- max(sub_i2$Repeat)
}
muinit <- mean(c(Y1,Y2)) # overall mean for both MS observations
alphainit <- mean(Y2) - mean(Y1)
betainit <- 1
s1init <- sqrt(mean(sub_vars1)) # repeatability for MS1 (pooled across parts)
s2init <- sqrt(mean(sub_vars2)) # repeatability for MS2 (pooled across parts)
spinit <- mean(c(sqrt(var(Y1) - s1init^2), sqrt(var(Y2) - s2init^2)))
# Saving initial parameters in matrix
init.param <- t(as.matrix(c(muinit,alphainit,betainit,spinit,s1init,s2init)))
# Maximize the Likelihood
options(warn = -1)
maximize <- optim(init.param, fn = log.likelihood, gr = gradient, n, r1vec, r2vec, data, method = "BFGS", hessian = T)
# Get MLEs and standard errors
param.hat <- as.matrix(t(maximize$par)) # save parameters estimates in matrix
muhat <- param.hat[1]
alphahat <- param.hat[2]
betahat <- param.hat[3]
sphat <- param.hat[4]
s1hat <- param.hat[5]
s2hat <- param.hat[6]
#Get Inverse of Fisher Information to get Assymptotic Variances
I <- fisher.info(param.hat, n, r1vec, r2vec)
inverseI <- solve(I)
# Square root of diagonal of Inverse Fisher information matrix: find Standard Errors
standard.errors <- as.matrix(sqrt(diag(inverseI)))
# standard.errors <- as.matrix(sqrt(diag(solve(maximize$hessian))))
# Construct estimate of P(agreement) and get its standard error for CIs
p <- seq(muhat-3*sphat, muhat+3*sphat, 0.1)
thetahat <- matrix(0, nrow = length(p))
SEtheta <- matrix(0, nrow = length(p))
for(i in 1:length(p)){
k1hat <- (-delta-alphahat-p[i]*(betahat-1))/(sqrt((s1hat^2)+(s2hat^2)))
k2hat <- (delta-alphahat-p[i]*(betahat-1))/(sqrt((s1hat^2)+(s2hat^2)))
thetahat[i] <- pnorm(k2hat)-pnorm(k1hat) # estimated theta (evaluated at MLEs)
D <- delta.method(delta = delta,p = p[i], alpha = alphahat, beta = betahat, s1 = s1hat,s2 = s2hat)
SEtheta[i] <- sqrt(t(D) %*% inverseI %*% D) # Standard error of Probability of Agreement
}
LCL <- thetahat-1.96*SEtheta
LCL[LCL < 0] = 0
UCL <- thetahat+1.96*SEtheta
LCL[LCL > 1] = 1
# Probability of Agreement plot
par(mfrow=c(1,1))
plot(
p,
thetahat,
main = bquote(
plain("Probability of Agreement Plot with ") ~ delta == .(delta)),
ylim = c(0,1),
col = "red",
type = "l",
ylab = "θ(s)", # probability of agreement of measurements (thetahat)
xlab = "s" # MLE of true value of measurement
)
points(
p,
thetahat - 1.96 * SEtheta,
type = "l",
col = "blue",
lty = 2
)
points(
p,
thetahat + 1.96 * SEtheta,
type = "l",
col = "blue",
lty = 2
)
labels <-
c("θ(s)","95% CI")
legend(
"bottomright",
inset = 0.02,
labels,
col = c("red", "blue"),
lty = c(7,2)
)
}
# Order of return: (mu, alpha, beta, sp, s1, s2)
return(
list(
Estimates = param.hat,
St.Errors = standard.errors,
PAgree = data.frame(
True.Val = p,
PA = thetahat,
SE = SEtheta)
)
)
}
## Construct relevant plots:
## Scatter plot (ignoring repeated measurements)
scatter.plot(data = rr, n = n_subj )## Scatter plot (accounting for repeated measurements)
rep.scatter.plot(data = rr, n = n_subj)## Repeatability plot
repeat.plot(data = rr, n = n_subj )## Modified QQ-plot
mod.qqplot(data = rr, n = n_subj )## Probability of agreement analysis:
p <- prob.agree(n = n_subj , delta = delta, data = rr)## Inspect parameter estimates and standard errors. (Order of return: mu, alpha, beta, sigma_s, sigma_1, sigma_2).
p$Estimates
p$St.Errors
## Inspect probability of agreement estimates.
p$PAgree## Probability agreement (NT Stevens, 2019)
# customized agreement plot
p$PAgree %>%
ggplot(
aes(
x = True.Val,
y = PA
)
) +
geom_ribbon(
aes(
ymin = PA - 1.96 * SE,
ymax = PA + 1.96 * SE
),
color = "lightgrey",
alpha = 0.2
) +
geom_line() +
scale_y_continuous(
limits = c(0, 1.05) # adjust upper limit so that all of 95%CI band is included
) +
xlab("True value (ML estimate)") +
ylab("Probability of agreement") +
labs(
title = paste ("Probability of agreement \nbetween manual and electro measurements"),
subtitle = paste ("Maximum allowed inter-rater difference:", paste(delta)),
caption = "Code: Parker RA et al., BMC Med Res Methodol 2020; 20:154"
) +
theme_light() +
theme(
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5)
)
2.3. Decomposition of accuracy, precision, and agreement, with statistical tests for equivalence
Both method appear to have good equivalence according to the
classical Bland-Altman plot and Probability of Agreement analysis. These
methods emphasize clinical significance but lack a clear null hypothesis
on method equivalence. Instead they are relying on visual inspection to
draw what what might be regarded as subjective conclusions.
Recently, a novel strategy was proposed based on the Bland and Altman method but enabling statistical decision-making. Causes of non-equivalence are identified in three steps 3:
# - Chunk_6 - #
# - Requires R-libraries (r libraries and data) and object exampleData_int (chunk_1)
## Assumption: manual is the new/putative method, electro is the reference
# Optional adjustment of new method for bias, using parameters obtained from Deming regression with the original data (see "regression" plot obtained with original data)
# C <- 1.3982
# m <- 0.9908
# default
C <- 0
m <- 1
## Run eirasagree
exampleData_replMean %>%
filter(time == "1") %>%
select(
id,
rater,
replicate_Mean # use "y" instead of "replicate_Mean" if raw replicated measurements should be used
) %>%
pivot_wider(
names_from = "rater",
values_from = "replicate_Mean"
) %>%
mutate(electro = (C + electro * m) ) %>%
AllStructuralTests(
reference.cols = {grep("electro", colnames(.) )},
newmethod.cols = {grep("manual", colnames(.) )},
alpha = 0.05,
out.format = "txt",
B = 5e3
)- Accuracy This test resembles the Bland-Altman plot in that differences between measurements obtained from the same subjects (yi − xi) are charted against the centered values of the reference measurement (xi − mean(x)). A simple linear regression is then applied. Under the null hypothesis of method equivalence the y-intercept, representing the mean difference between measurements, is 0. The null hypothesis is rejected if bootstrapping shows 0 outside the 95% confidence interval for the y-intercept.
- Precision The slope of the Bland-Altman plot is used to detect unequal precision, as regression of 𝑦 −𝑥 against 𝑥 +𝑦 will not be 0 when the true ratio of variability of measurement errors λ ≠ 1. If a horizontal line cannot be fitted into the 95% confidence band defined by the functional regression, the null hypothesis of equal variability of measurement errors is rejected.
- Agreement This test applies the Deming regression to verify if two measurement techniques measure the same values in the same subjects […] depends on λ, estimated in the previous step, to compute the true X and Y values before the computation of the regression estimates. When the value of lambda is not assumed to be 1, it affects the band width. Analytically, the null hypothesis is rejected if α ≠ 0 and β ≠ 1. Since these two parameters are estimated together, the Bonferroni correction is applied to control the probability of a type I error and preserve test power (effective significance level is 2.5%). Graphically, two alternative statistical approaches were implemented for the bisector line agreement: the assessments of the 95% prediction ellipse and of the 95% confidence band of regression, both done by bootstrapping. In the first one, the null hypothesis is to be rejected if (β, α) is not inside the ellipse, in the second one, if the bisector line cannot fit inside the band. These methods test intercept and slope together and provide stronger statistical power than an independent assessment of these two factors.
eirasagree also calculates slope and intercept of the Deming regression line:
Decision by 95% confidence band: avg{manual - electro} = 1.3763 +
0.0184 {(electro + manual)/2}
This information can be used to adjust method results, and test
whether method equivalence can be achieved by a simple data
transformation.
3. Adjustment of electro measurements for bias improves the agreement of the means
3. 1. Replicate means of adjusted measurements and their variance
# - Chunk_7 - #
# - Requires R-libraries (r libraries and data) and object exampleData_int (chunk_1) - #
C <- 1.3982
m <- 0.9908
# default
# C <- 0
# m <- 1
## Adjust data for bias of the new method
exampleData_replMean_adj <-
exampleData_replMean %>%
select(id, time, rater, y) %>% # selected columns must have unique rows across all column pairs
pivot_wider( # adjust electro values for constant and proportional bias
names_from = "rater",
values_from = "y"
) %>%
mutate(electro = C + electro * m ) %>%
pivot_longer(
cols = c(manual, electro),
names_to = "rater",
values_to = "y"
) %>%
pivot_wider( # calculate replicate means and SE
names_from = "time",
values_from = "y") %>%
rowwise() %>%
mutate(
replicate_Mean = mean(
c_across(
all_of(replicates)
)
),
replicate_SE = std(
c_across(
all_of(replicates)
)
)
) %>%
pivot_longer(
cols = replicates,
names_to = "time",
values_to = "y"
)
# Standard error function
std <- function(x) sd(x)/sqrt(length(x))
# Calculate intra-subject replicate means
replicates <- exampleData_replMean_adj %>% .$time %>% unique
exampleData_replMean_adj <-
exampleData_replMean_adj %>%
pivot_wider(
names_from = "time",
values_from = "y") %>%
rowwise() %>%
mutate(
replicate_Mean = mean(
c_across(
all_of(replicates)
)
),
replicate_SD = sd(
c_across(
all_of(replicates)
)
),
replicate_SE = std(
c_across(
all_of(replicates)
)
)
) %>%
pivot_longer(
cols = replicates,
names_to = "time",
values_to = "y"
)
# Table with summary statistics of the original sample data
exampleData_replMean_summary_adj <-
exampleData_replMean_adj %>%
{ rbind(
group_by (., rater, time) %>% # for each method and replicate, calculate mean and SE across subjects
summarise(
Mean = mean(y),
SD = sd(y),
SE = std(y)
),
group_by (., rater) %>% # for each method, calculate means across subjects for mean and SE across replicates
summarise(
Mean = mean(replicate_Mean),
SD = sd(replicate_Mean),
SE = std(replicate_Mean)
)
)
} %>%
mutate(
errorbar_pos = rep(- 0.06, n()),
) %>% # defines vertical position of errorbars and their mean in the facets
replace(is.na(.), "replicate mean")
exampleData_replMean_summary_adj %>%
select(!errorbar_pos) %>%
kable(., "html") %>%
kable_styling(
bootstrap_options = c(
"striped", "hover", "responsive", full_width = FALSE
)
) | rater | time | Mean | SD | SE |
|---|---|---|---|---|
| electro | 1 | 1.398200 | 7.021011 | 1.3037690 |
| electro | 2 | 1.432365 | 7.412022 | 1.3763779 |
| electro | 3 | 1.500697 | 7.082388 | 1.3151664 |
| manual | 1 | 1.586207 | 7.218613 | 1.3404628 |
| manual | 2 | 1.482759 | 7.471951 | 1.3875064 |
| manual | 3 | 1.241379 | 7.452972 | 1.3839822 |
| electro | replicate mean | 1.443754 | 7.044320 | 0.7552303 |
| manual | replicate mean | 1.436782 | 7.263108 | 0.7786869 |
## Plot distribution, mean and SE of adjusted measurements
# Change over consecutive replicates
p2 <-
exampleData_replMean_adj %>%
bind_rows(
exampleData_replMean_adj %>%
filter(time == "1") %>% # select a unique row for each subject's replicate mean (avoid generating one row for each replicate)
mutate(
.keep = "none",
id = id,
rater = rater,
y = replicate_Mean)
) %>%
select(id, time, rater, y) %>% # selected columns must have unique rows across all column pairs
replace(is.na(.), "replicate mean") %>%
ggplot(
aes(x = y)
) +
geom_density( # show smoothed distribution of measurement values
aes(y*0.7,
fill = rater
),
alpha = 0.3,
color = "grey",
position = position_nudge(y = 0.02)
) +
geom_point_interactive( # show data points of all subjects
aes(
y = -0.02, # defines vertical position of the data point rows in the facets
color = rater,
data_id = id,
tooltip = id
),
size = 3,
alpha = 0.5
) +
geom_errorbarh ( # show SE bar
data = exampleData_replMean_summary_adj,
inherit.aes = FALSE,
aes(
xmin = Mean - SE,
xmax = Mean + SE,
y = errorbar_pos - 0.025,
color = rater
),
height = 0.025 # height of bar caps
) +
geom_errorbarh ( # show SD bar
data = exampleData_replMean_summary_adj,
inherit.aes = FALSE,
aes(
xmin = Mean - SD,
xmax = Mean + SD,
y = errorbar_pos,
color = rater
),
height = 0.04 # height of bar caps
) +
geom_point( # add means to SE bars
data = exampleData_replMean_summary_adj,
aes(
x = Mean,
y = errorbar_pos,
color = rater
),
size = 2
) +
facet_wrap(
.~ time,
dir = "v",
strip.position = "left",
nrow = 4
) +
labs(
title = "Inter-method comparison of mean, SD and SE \nof consecutive replicate measurements",
y = "Replicate",
x = "Measurement",
tag = "A"
) +
scale_fill_manual( values = c("blue", "red")) +
scale_color_manual(values = c("blue", "red")) +
guides(
fill = guide_legend(title="Method"),
color = guide_legend(title="Method")
) +
theme_light() +
theme(
plot.title = element_text(
size = 12,
hjust = 0.5),
axis.text.y = element_blank(),
strip.background = element_rect(
fill = "white",
color = "darkgrey"
),
axis.ticks.y = element_blank(),
strip.text.y = element_text(
colour="black",
size = 12
)
)
# Variance vs. value (replicate mean)
p3 <-
exampleData_replMean_adj %>%
filter(time == "1") %>% # select a unique row for each subject's replicate mean (avoid generating one row for each replicate)
ggplot(
aes(
x = replicate_Mean,
y = replicate_SD,
colour = rater,
data_id = id,
tooltip = id
)
) +
geom_point_interactive( # show data points of for replicate means of all subjects
size = 3,
alpha = 0.5
) +
labs(
title = "Intra-subject replicate standard deviation vs. mean",
caption = "Hover over points to highlight individual patients",
x = "replicate mean",
y = "SD",
tag = "B"
) +
scale_color_manual(values = c("blue", "red")) +
theme_light() +
theme(
plot.title = element_text(
size = 12,
hjust = 0.5),
plot.caption = element_text(
size = 9,
hjust = 0.5),
legend.position = "none"
)
girafe(
ggobj =
p2 +
plot_spacer() +
p3 +
plot_layout(
nrow = 3,
heights = unit(c(18, 0.2, 5), "cm"
)
) +
plot_annotation(
title = "Distributions and means of measurements, \nelectro method adjusted for bias",
theme = theme(
plot.title = element_text(
size = 14,
hjust = 0.5
)
)
),
height_svg = 12,
width_svg = 6,
options = list(
opts_hover(css = "stroke:black; stroke-width:2px;"), # CSS code of selected line
opts_hover_inv(css = "opacity:0.1;"), # CSS code of all other lines
opts_sizing(rescale = FALSE) # allows to zoom in and out of saved html image
)
)
* Sample and sample/replicate means become almost identical when
measurements with the electro method are adjusted for bias Note: The
electro was meant to be the reference method in eirasagree. In a later
version, I intend to change this and the following figure to showing
adjusted manual values
3. 2. Bland-Altman plot with adjusted measurements
# - Chunk_8 - #
# - Requires R-libraries (r libraries and data) and object exampleData_int (chunk_1) - #
## Bland-Altman plot with differences for each individual replicate (original measurement values)
exampleData_BA_adj <-
exampleData_replMean_adj %>%
select(id, time, rater, y) %>%
pivot_wider(
names_from = "rater",
values_from = "y"
) %>%
mutate(
method_mean = (manual + electro)/2,
method_diff = manual - electro
)
summaryData_BA_adj <-
exampleData_BA_adj %>%
group_by(time) %>%
summarise(
Mean = mean(method_diff),
SD = sd(method_diff),
LoA_up = Mean + 1.96* SD,
LoA_low = Mean - 1.96* SD
)
p9 <-
exampleData_BA_adj %>%
ggplot(
aes(
x= method_mean,
y= method_diff,
colour = time,
data_id = id,
tooltip = id
)
) +
geom_point_interactive(
size = 3,
alpha = 0.5
) +
geom_ysidedensity(
aes(fill = time),
alpha = 0.1
) +
geom_hline(
data = summaryData_BA_adj,
aes(
yintercept = summaryData_BA_adj$Mean,
color = time
),
show.legend = FALSE
) +
geom_hline(
data = summaryData_BA_adj,
aes(
yintercept = summaryData_BA_adj$LoA_up,
color = time
),
alpha = 0.7,
show.legend = FALSE
) +
geom_hline(
data = summaryData_BA_adj,
aes(
yintercept = summaryData_BA_adj$LoA_low,
color = time
),
alpha = 0.7,
show.legend = FALSE
) +
scale_color_manual(
values = c("#F8766D", "#00BA38", "#619CFF", "darkgrey")
) +
scale_y_continuous(
limits = c(
-5, # min(exampleData_BA_adj$method_diff),
7 # max(exampleData_BA_adj$method_diff)
),
breaks = seq(-7, 7, 1)
) +
labs(
tag = "B",
x = "(manual+electro)/2",
y = "manual-electro"
) +
guides(
theme(
legend.text = element_text(size = 13),
color = guide_legend(title="Mean and limits of agreement for individual replicate datasets", position = "top"),
fill = guide_legend(title="Mean and limits of agreement for individual replicate datasets", position = "top")
),
color = FALSE
) +
theme_light() +
theme_ggside_void()
## Bland-Altman plot with differences between inter-replicate means (original measurement values)
# Calculate inter-replicate mean and SD for each sample
exampleData_BA_means_adj <-
exampleData_BA_adj %>%
select(id, time, method_mean, method_diff) %>%
split(.$id) %>%
lapply(
function(x) summarise(x,
Mean_methodMean = mean(method_mean),
SD_methodMean = sd(method_mean),
Mean_methodDiff = mean(method_diff),
SD_methodDiff = sd(method_diff)
)
) %>%
list_rbind(names_to = "id")
# Calculate summary stats of inter-replicate means and LoAs
SD_means <-
function(x) {
sqrt(sum(x^2)/(length(x)-1))
}
summaryData_BA_means_adj <-
exampleData_BA_means_adj %>%
summarise(
Mean = mean(Mean_methodDiff),
SD = SD_means(Mean_methodDiff)
) %>%
mutate(
LoA_up = Mean + 1.96*SD,
LoA_low = Mean - 1.96*SD
)
# Bland-Altman plot
p10 <-
exampleData_BA_means_adj %>%
ggplot(
aes(
x= Mean_methodMean,
y= Mean_methodDiff,
data_id = id,
tooltip = id
)
) +
geom_point_interactive(
color = "black",
alpha = 0.3,
size = 3
) +
geom_ysidedensity(
inherit.aes = FALSE,
aes(y = Mean_methodDiff)
) +
geom_hline(
data = summaryData_BA_means_adj,
aes(yintercept = summaryData_BA_means_adj$Mean)
) +
geom_hline(
data = summaryData_BA_means_adj,
aes(yintercept = summaryData_BA_means_adj$LoA_up),
) +
geom_hline(
data = summaryData_BA_means_adj,
aes(yintercept = summaryData_BA_means_adj$LoA_low)
) +
scale_y_continuous(
limits = c(
-5, # min(exampleData_BA_adj$method_diff), # same y-scale limits as p6
7 # max(exampleData_BA_adj$method_diff)
),
breaks = seq(-7, 7, 1)
) +
labs(
tag = "A",
title = "Mean and limits of agreement for averaged replicate measurements ",
x = "(manual+electro)/2",
y = "manual-electro",
caption =
summaryData_BA_means_adj %>%
round(digits = 2) %>%
{ paste("Mean = ", .[1],
"\n upper LoA = ", .[3],
"\n lower LoA = ", .[4]
)
}
) +
theme_light() +
theme(
plot.title = element_text(
size = 11,
hjust = 0.3
),
legend.position = "none"
) +
theme_ggside_void()
girafe(
ggobj = p10 + plot_spacer() + p9 +
plot_layout(
nrow = 3,
heights = c(4, 0.2, 4)
) +
plot_annotation(
title =
"Bland-Altman plots with replicate measurements, \nelectro method adjusted for bias",
caption = "Hover over points to highlight individual patients",
theme = theme(
plot.title = element_text(
face = "bold",
size = 12,
hjust = 0.5,
vjust = 3
),
plot.caption = element_text(
size = 9,
hjust = 0.5
)
)
),
# height_svg = 9, # svg = 9 for saving to individual image to html, omit for knitting
options = list(
opts_hover(css = ''), # CSS code of selected line
opts_hover_inv(css = "opacity:0.05;"), # CSS code of all other lines
opts_sizing(rescale = FALSE)
)
)
4. References
Statistical methodology for the concurrent assessment of interrater and intrarater reliability: using goniometric measurements as an example. Eliasziw M, Young SL, Woodbury MG, Fryday-Field K. Phys Ther. 1994;74(8):777-88 https://doi.org/10.1093/ptj/74.8.777
Assessing agreement between two measurement systems: An alternative to the limits of agreement approach. Stevens NT, Steiner SH, MacKay RJ. Stat Methods Med Res. 2017;26(6):2487-2504 https://doi.org/10.1177/0962280215601133
Is the Bland-Altman plot method useful without inferences for accuracy, precision, and agreement? Silveira PSP, Vieira JE, Siqueira JO. Rev Saude Publica. 2024;58:01 https://doi.org/10.11606/s1518-8787.2024058005430