Simulating Minirhyzotron data
The data is simulated as a mixture between a Uniform background distribution and N gaussian distributions where N follows a uniform distribution and each gaussian is parameterized by mu and sigma. Mu also follows a uniform distribution and sigma follows a half-normal distribution. Pixels are assigned to the background or the gaussian mixture component according to theta. Pixels assigned to the gaussian mixture are randomly assigned to each of the N gaussian distributions. See the rRhyzoDist function below. We also define functions that will generate single value traits from the MV frequencies.
Code
rRhyzoDist <- function(n, theta = 0.3, u1_max = 20, u2_max = 5500, sd = 200, abs_max = 5500) {
#* split n_pixels based on theta into background and gaussians
n_unif_pixels <- ceiling(n * theta)
n_gauss_pixels <- floor(n * (1 - theta))
#* background is uniform
background <- runif(n_unif_pixels, 0, u2_max)
#* simulate a number of gaussians randomly between 1 and u1_max
n_gaussians <- runif(1, 1, u1_max)
#* each gaussian has a mean that is uniform between 1 and u2_max
mu_is <- lapply(seq_len(n_gaussians), function(i) {
runif(1, 1, u2_max)
})
#* each gaussian has a sigma that is half-normal based on sd
sd_is <- lapply(seq_len(n_gaussians), function(i) {
extraDistr::rhnorm(1, sd)
})
#* assign pixels randomly to gaussians
index <- sample(seq_len(n_gaussians), size = n_gauss_pixels, replace = TRUE)
px_is <- lapply(seq_len(n_gaussians), function(i) {
sum(index == i)
})
#* draws n_pixels time from each gaussian
d <- unlist(lapply(seq_len(n_gaussians), function(i) {
rnorm(px_is[[i]], mu_is[[i]], sd_is[[i]])
}))
#* combine data
d <- c(d, background)
#* make sure no gaussians return data out of bounds
d[d < 0] <- runif(sum(d < 0), 0, abs_max)
d[d >= abs_max] <- runif(sum(d >= abs_max), 0, abs_max)
return(d)
}
lastNonZeroBin <- function(d) {
max(d[d$value > 0, "label"])
}
tubeAngleToDepth <- function(x, theta) {
sin(theta) * x
}
mv_mean <- function(d) {
weighted.mean(d$label, d$value)
}
mv_median <- function(d) {
median(rep(d$label, d$value))
}
mv_std <- function(d) {
sd(rep(d$label, d$value))
}
sv_from_mv <- function(df, theta) { # note this should also return mean/median/std
metaCols <- colnames(df)[-which(grepl("value|label|trait", colnames(df)))]
out <- aggregate(as.formula(paste0("value ~ ",
paste(metaCols, collapse = "+"))),
data = df, sum, na.rm = TRUE)
colnames(out)[ncol(out)] <- "area"
out$max_pixel <- unlist(lapply(split(df, interaction(df[, metaCols])), lastNonZeroBin))
out$height <- tubeAngleToDepth(out$max_pixel, 0.35)
out$mean_x_frequencies <- unlist(lapply(split(df, interaction(df[, metaCols])), mv_mean))
out$median_x_frequencies <- unlist(lapply(split(df, interaction(df[, metaCols])), mv_median))
out$std_x_frequencies <- unlist(lapply(split(df, interaction(df[, metaCols])), mv_std))
return(out)
}The simulated data looks realistic based on limited test data available at time of writing.
Code
set.seed(123)
ex <- do.call(rbind, lapply(1:20, function(rep) {
n_total_pixels <- runif(1, 100, 3000)
x <- rRhyzoDist(n = n_total_pixels)
h <- hist(x, plot = FALSE, breaks = seq(0, 5500, 20))
breaks <- h$breaks[-1]
counts <- h$counts
data.frame(
rep = as.character(rep),
value = counts, label = breaks,
trait = "x_frequencies"
)
}))
pcv.joyplot(ex, "x_frequencies", group = c("rep"))
For the rest of this vignette we will use simulated data. The first simulated data df assumes that roots can appear and disappear from the minirhyzotron images over time. The second simulated data assumes that once a root is seen along the minirhyzotron that it will stay visible.
Code
n_times <- 5
parameters <- data.frame(
time = c(1:n_times),
n_min = rep(0, n_times),
n_max = seq(2000, 5500, length.out = n_times),
theta = rep(0.3, n_times),
u1_max = seq(10, 20, length.out = n_times),
u2_max = seq(2000, 5500, length.out = n_times),
u2_max_noise = seq(400, 200, length.out = n_times),
sd = seq(150, 250, length.out = n_times)
)
set.seed(123)
df <- do.call(rbind, lapply(seq_len(nrow(parameters)), function(time) {
pars <- parameters[parameters$time == time, ]
do.call(rbind, lapply(1:10, function(rep) {
n_total_pixels <- runif(1, pars$n_min, pars$n_max)
u2_max_iter <- ceiling(rnorm(1, pars$u2_max, pars$u2_max_noise))
x <- rRhyzoDist(
n = n_total_pixels, theta = pars$theta,
u1_max = pars$u1_max,
u2_max = u2_max_iter, sd = pars$sd
)
h <- hist(x, plot = FALSE, breaks = seq(0, 5500, 20))
breaks <- h$breaks[-1]
counts <- h$counts
data.frame(
rep = as.character(rep),
time = as.character(time),
value = counts, label = breaks,
trait = "x_frequencies"
)
}))
}))
df$rep <- factor(df$rep, levels = seq_along(unique(df$rep)), ordered = TRUE)
sv <- sv_from_mv(df)Code
n_times <- 5
parameters <- data.frame(
time = c(1:n_times),
n_min_new_px = rep(50, n_times),
n_max_new_px = seq(1000, 2500, length.out = n_times),
theta = rep(0.3, n_times),
u1_max = seq(10, 20, length.out = n_times),
mean_added_depth = seq(log(200), log(1200), length.out = n_times),
added_depth_noise = rep(0.1, n_times),
sd = seq(150, 250, length.out = n_times)
)
set.seed(567)
df2 <- do.call(rbind, lapply(1:10, function(rep) {
dList <- list()
add_area <- 2000
previous_max_depth <- 2000
d <- numeric(0)
for (time in seq_len(nrow(parameters))) {
pars <- parameters[parameters$time == time, ]
n_total_pixels <- add_area + runif(1, pars$n_min_new_px, pars$n_max_new_px)
max_depth <- ceiling(previous_max_depth + rlnorm(1, pars$mean_added_depth, pars$added_depth_noise))
x <- rRhyzoDist(
n = n_total_pixels, theta = pars$theta,
u1_max = pars$u1_max,
u2_max = max_depth, sd = pars$sd
)
d <- append(d, x)
h <- hist(d, plot = FALSE, breaks = seq(0, 5500, 20))
breaks <- h$breaks[-1]
counts <- h$counts
dList[[time]] <- data.frame(
rep = as.character(rep),
time = as.character(time),
value = counts, label = breaks,
trait = "x_frequencies"
)
add_area <- 0
previous_max_depth <- max_depth
}
do.call(rbind, dList)
}))
df2$rep <- factor(df2$rep, levels = seq_along(unique(df2$rep)), ordered = TRUE)
sv2 <- sv_from_mv(df2)
