knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(modelCathRes) library(dplyr, quietly = T, warn.conflicts = F) library(ggplot2) library(rstan) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores())
data("resCathData") d <- resCathData # reduce number of pots d <- d %>% filter(cod_cuve %in% sample(unique(d$cod_cuve), 20)) head(d)
Formula: $a + b.log(x)$
funLog <- findLogFunction(d, agebsq, rucv) # log fit curve curve <- ggplot2::stat_function(xlim = c(min(d$agebsq), max(d$agebsq)), fun = funLog, color = "darkred", size = 1.2) ggplot2::ggplot(d, ggplot2::aes(x = agebsq, y = rucv)) + ggplot2::geom_point(shape = 21, fill = NA, size = .9, alpha = .5) + ggplot2::coord_cartesian(ylim = c(.4, 1.0)) + curve
Same model $a + b.log(age)$ for all pots.
plotCathRes(d, cod_cuve, agebsq, rucv, ncol = 4) + labs(title = "Evolution of Cathode Resistance per pot") + coord_cartesian(ylim = c(.5, 1.0)) + curve
Using lm with log(age) function for age dependence and 0 intercept. Estimate in table is initial value of cathode resistance.
modelCathoResLm(d, group, agebsq, rucv)
Each group has own log coefficient and intercept.
# lmer model from lme4with formula rucv ~ 1 + (1 + log(agebsq) | group)) (modLmer <- modelCathoResLme(d, group, agebsq, rucv)) # Coefficients for each group (coefs <- coef(modLmer)$group) # Adding predictions to data d <- modelr::add_predictions(d, model = modLmer) # Graph of cathode resistance with age for each group modelCathRes::plotCathRes(d, group, agebsq, rucv, ncol = 4) + coord_cartesian(ylim = c(.5, 1.0)) + geom_line(aes(x = agebsq, y = pred), color = "darkorange")
Using a + b.log(age) model to estimate the mean resistance over total life
# Comparison of measured mean values and predicted mean values # calculateCathodeResMean(d, group, agebsq, rucv) d %>% group_by(group) %>% summarise(meanResValues = mean(rucv, na.rm = T)) %>% inner_join(calculateCathodeResMean(d, group, agebsq, rucv) , "group" ) %>% mutate(across(.cols = where(is.numeric), .fns = signif, dig = 3))
Temps d'exécution beaucoup trop long pour le modèle le plus simple...
# Test data n <- 1500 # nbs of obs per group dTest <- tibble(group = c(rep("A", n), rep("B", n)), rucv = c(rnorm(n, 0.45, .08), rnorm(n, 0.55, .12)) ) #stanData <- tidybayes::compose_data(select(d, c(group, rucv))) stanDataTest <- tidybayes::compose_data(dTest) names(stanDataTest) <- c("G", "N_G", "R", "N") samples <- rstan::sampling(object = rstan::stan_model("resGroupSimple.stan"), data = stanDataTest, chains = 2, iter = 1000) samples
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.