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())

Get test data

data("resCathData")
d <- resCathData
# reduce number of pots 
d <- d %>% 
  filter(cod_cuve %in% sample(unique(d$cod_cuve), 20))
head(d)

Log model for cathode resistance for all data

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

View cathode resistance and model for each pot

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

Comparing cathode resistance between groups

lm model

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)

lmer model

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")

Comparing average Cathode resistance over life

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))

Stan model

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


OlivierGranacher/modelCathRes documentation built on June 23, 2022, 3:03 p.m.