polarolouis/anova-phylogenetique-projet-msv
1
0
Fork 0
mirror of https://github.com/Polarolouis/anova-phylogenetique-projet-msv.git synced 2026-06-17 18:25:25 +02:00
Code Issues Projects Releases Packages Wiki Activity

🐛Fixing forgotten F test for pvalues, missing measurement error

Browse source
This commit is contained in:
Louis Lacoste 2024-01-24 11:03:37 +01:00
parent 39df26a712
commit 73448ad446
1 changed files with 9 additions and 7 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

14
R/utils.R
View file

@ -89,7 +89,6 @@ ddf_species <- function(fitlm, phylo) {
#' This code computes the satterthwaite approximation to obtain degrees of #' This code computes the satterthwaite approximation to obtain degrees of
#' freedom for a given tree #' freedom for a given tree
ddf_satterthwaite_sum <- function(fit_phylolm, phylo, REML = FALSE) { ddf_satterthwaite_sum <- function(fit_phylolm, phylo, REML = FALSE) {
if (is.null(fit_phylolm$sigma2_error) || fit_phylolm$sigma2_error == 0) { if (is.null(fit_phylolm$sigma2_error) || fit_phylolm$sigma2_error == 0) {
#  In this case we return the classical df #  In this case we return the classical df
return(ddf_species(fit_phylolm, phylo)) return(ddf_species(fit_phylolm, phylo))
@ -243,7 +242,7 @@ compute_trait_values <- function(
infere_anova_phyloanova <- function(y, groups, tree, stoch_process = "BM") { infere_anova_phyloanova <- function(y, groups, tree, stoch_process = "BM") {
#  The fits #  The fits
fit_anova <- lm(y ~ groups) fit_anova <- lm(y ~ groups)
fit_phylolm <- phylolm(y ~ groups, phy = tree, model = stoch_process) fit_phylolm <- phylolm(y ~ groups, phy = tree, model = stoch_process, measurement_error = TRUE)
return(list(anova = fit_anova, phyloanova = fit_phylolm)) return(list(anova = fit_anova, phyloanova = fit_phylolm))
} }
@ -257,7 +256,8 @@ infere_anova_phyloanova <- function(y, groups, tree, stoch_process = "BM") {
pvalues_from_fits <- function( pvalues_from_fits <- function(
fit_anova, fit_anova,
fit_phylolm, tree, fit_phylolm, tree,
tested_method = c("vanilla", "satterthwaite", "lrt")) { tested_method = c("vanilla", "satterthwaite", "lrt"),
REML = FALSE) {
#  For sanity test #  For sanity test
match.arg(tested_method) match.arg(tested_method)
@ -281,19 +281,21 @@ pvalues_from_fits <- function(
switch(tested_method, switch(tested_method,
"vanilla" = { "vanilla" = {
pvalue_phylolm <- compute_F_statistic( F_stat <- compute_F_statistic(
r_squared = fit_phylolm$r.squared, r_squared = fit_phylolm$r.squared,
df1 = df1, df1 = df1,
df2 = df2 df2 = df2
) )
pvalue_phylolm <- 1 - pf(F_stat, df1, df2)
}, },
"satterthwaite" = { "satterthwaite" = {
df2 <- ddf_satterthwaite_sum(fit_phylolm = fit_phylolm, phylo = tree, REML = REML) df2 <- ddf_satterthwaite_sum(fit_phylolm = fit_phylolm, phylo = tree, REML = REML)$ddf
pvalue_phylolm <- compute_F_statistic( F_stat <- compute_F_statistic(
r_squared = fit_phylolm$r.squared, r_squared = fit_phylolm$r.squared,
df1 = df1, df1 = df1,
df2 = df2 df2 = df2
) )
pvalue_phylolm <- 1 - pf(F_stat, df1, df2)
}, },
"lrt" = { "lrt" = {
h0_phylolm <- phylolm(fit_phylolm$y ~ 1, h0_phylolm <- phylolm(fit_phylolm$y ~ 1,