Skip to content
Snippets Groups Projects
Commit 224a77e8 authored by Aurelien Brionne's avatar Aurelien Brionne
Browse files

eggbiomechnaic upgrade

parent 06e1e098
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
DESCRIPTION
+ 2
− 2
View file @ 224a77e8
Package: eggbiomechanic
Version: 1.2.2
Version: 1.2.3
Title: Analysis of egg mechanical properties (eggshell, vitellines membranes).
Author: Aurelien Brionne [aut, cre]
Maintainer: Aurelien Brionne <aurelien.brionne@inra.fr>
......@@ -7,12 +7,12 @@ Depends: R (>= 3.6)
Imports: data.table,
DiagrammeR,
DT,
emmeans,
nparcomp,
nlme,
ggplot2,
gridExtra,
methods,
multcomp,
multcompView,
reshape,
zoo
......
NAMESPACE
+ 1
− 1
View file @ 224a77e8
......@@ -24,7 +24,6 @@ importFrom(data.table,data.table)
importFrom(data.table,fread)
importFrom(data.table,rbindlist)
importFrom(data.table,setorderv)
importFrom(emmeans,emmeans)
importFrom(ggplot2,aes)
importFrom(ggplot2,element_blank)
importFrom(ggplot2,facet_wrap)
......@@ -39,6 +38,7 @@ importFrom(ggplot2,theme)
importFrom(ggplot2,theme_bw)
importFrom(gridExtra,grid.arrange)
importFrom(methods,setClass)
importFrom(multcomp,glht)
importFrom(multcompView,multcompLetters)
importFrom(nlme,lme)
importFrom(zoo,rollmean)
NEWS
+ 4
− 0
View file @ 224a77e8
......@@ -4,3 +4,7 @@ CHANGES IN VERSION 1.2
o S3 to S4 class
o methods upgrade
o var_test print option
CHANGES IN VERSION 1.2
------------------------
o use glht instead emeans in Stats_lme
R/stat_diagram.R
+ 1
− 1
View file @ 224a77e8
......@@ -16,7 +16,7 @@ stat_diagram=function(){
B1(Fit an analysis of variance model</br><small>stats::aov</small>)
B2(Compute Tukey Honest Significant Differences</br><small>stats::TukeyHSD</small>)
C1(Fit an analysis of variance model</br>with linear mixed-effects </br><small>nlme::lme, stats::anova</small>)
C2(Estimated marginal means</br><small>emmeans::emmeans</small>)
C2(Compute Tukey Significant Differences</br><small>multcomp::glht</small>)
A-->B1
B1-->B2
A-->C1
......
R/stat_lm.R
+ 4
− 1
View file @ 224a77e8
......@@ -140,11 +140,14 @@ stat_lm=function(Data,f,print=TRUE){
x
)
# rename p value column
names(Data)[5]<-"p.value"
# remove NA
Data<-na.exclude(Data)
# extract pvalues
pvalues=unlist(Data[,"p adj",with=FALSE])
pvalues=unlist(Data[,"p.value",with=FALSE])
names(pvalues)<-Data$contrast
# pvalues to letters
......
R/stat_lme.R
+ 70
− 50
View file @ 224a77e8
......@@ -4,7 +4,7 @@
#' @importFrom ggplot2 ggplot aes geom_text geom_abline theme_bw theme labs element_blank
#' @importFrom multcompView multcompLetters
#' @importFrom nlme lme
#' @importFrom emmeans emmeans
#' @importFrom multcomp glht
#' @param Data a \code{\link[data.table]{data.table}}.
#' @param f The \code{\link[stats]{formula}}.
#' @param random The random effect from \code{\link[nlme]{lme}}.
......@@ -13,9 +13,9 @@
#' @details
#' This function is used in internal by the \code{\link{var_test}} function in order to perform Linear mixed-effects models statistics:
#' \enumerate{
#' \item analyse of variance using \code{\link[nlme]{lme}}.
#' \item analyse of variance using \code{\link[stats]{anova}} on \code{\link[nlme]{lme}} fitting model.
#' \item residues studie using qqplot.
#' \item post Hoc test with \code{\link[emmeans]{emmeans}}.
#' \item Tukey post Hoc test with \code{\link[multcomp]{glht}}.
#' }
#' @return a \code{\link[base]{list}} and text output.
#' @examples
......@@ -38,15 +38,25 @@ stat_lme=function(Data,f,random,method,print=TRUE){
fii<-formula(fi)
fi=sub("~","",fi)
randomi=paste("(",sub("^.+\\|","",random),")",sep="")
randomi=paste("(",random,")",sep="")
random<-formula(random)
# init list object
Results<-list()
Results<-list()
# lm
Lme=lme(f,data=Data,random=random,method=method)
# if print is TRUE
if(print==TRUE){
# add anova results
cat("\n\nused lme model:.\n\n")
# display model summary
print(summary(Lme))
}
# compute standardized residues
standardized_residues=residuals(Lme,method="pearson")
......@@ -102,77 +112,87 @@ stat_lme=function(Data,f,random,method,print=TRUE){
print(Results$anova_test)
}
# EMM tests
EMM<-emmeans(Lme,fii)
# contrast to do
type<-unlist(strsplit(fi,"[[:punct:]]"))
# EMM plot
EMMplot=plot(EMM,comparisons =TRUE)+
labs(title="EMM comparisons")+
theme_bw()+
theme(
panel.grid.major=element_blank(),
panel.grid.minor=element_blank()
# build list object
contr<-as.list(rep("Tukey",length(type)))
names(contr)<-type
# convert to mcp object used by glht
attr(contr,"interaction_average")<-FALSE
attr(contr,"covariate_average")<- FALSE
attr(contr,"class")<-"mcp"
# Tukey post hoc test
tukey<-glht(
Lme,
linfct =contr
)
# compute post hoc test summary
tukey<-summary(tukey)
# compute confidance intervals
tukey<-confint(tukey)
# qqplot in list
plots=list(
qqplot=qqplot,
EMMplot=EMMplot
qqplot=qqplot
)
# init list object
Results<-c(Results,plots=list(plots))
# contrast tests
contrasts_test<-pairs(EMM)
# contrast tests
# global contrast tests
contrasts_test<-data.table(
summary(contrasts_test)
contrast=gsub("[[:space:]]","",row.names(tukey$confint)),
data.table(tukey$confint),
p.value=tukey$test$pvalues
)
# rename columns 2
# rename secodn column
names(contrasts_test)[2]<-"diff"
# contrast tests adds
contrasts_test[,
`:=`(
contrast=gsub(" ","",contrast),
lwr=diff-SE,
upr=diff+SE
# format output
tukey=lapply(unique(gsub(":.+$","",contrasts_test$contrast)),function(x){
# contrast tests
Data<-contrasts_test[contrast%like%x]
# remove type
Data[,contrast:=gsub("^.+:","",contrast)]
# extract pvalues
pvalues=unlist(Data[,"p.value",with=FALSE])
names(pvalues)<-Data$contrast
# pvalues to letters
pvalues<-multcompLetters(pvalues)
# pvalues to letters in table
pvalues<-data.table(
group=factor(names(pvalues$Letters)),
text=pvalues$Letters
)
]
# extract pvalues
pvalues=unlist(contrasts_test[,"p.value",with=F])
names(pvalues)<-contrasts_test$contrast
# pvalues to letters
pvalues<-multcompLetters(pvalues)
# pvalues to letters in table
pvalues<-data.table(
group=factor(names(pvalues$Letters)),
text=pvalues$Letters
)
# return results
contrasts_test=list(
# results in a list
list(
stats=contrasts_test,
stats=Data,
letters=pvalues
)
)
names(contrasts_test)<-paste(fi,randomi,sep="")
})
names(tukey)<-paste(unique(gsub(":.+$","",contrasts_test$contrast)),randomi)
# add TukeyHSD results
Results<-c(Results,contrasts_test=list(contrasts_test))
Results<-c(Results,contrasts_test=list(tukey))
# print or not
if(print==TRUE){
# return results
cat("\n\nResults of EMM contrasts.\n\n")
cat("\n\nResults of Tukey contrasts.\n\n")
# return results
print(Results$contrasts_test)
......
R/var_test.R
+ 10
− 41
View file @ 224a77e8
......@@ -147,7 +147,7 @@ setMethod(
Summary<-lapply(Class,function(x){
# tested class
x<-gsub("\\(.+$","",x)
x<-gsub(" \\(.+$","",x)
# if interaction
if(length(grep(":|\\*|\\+",x))==1){
......@@ -204,18 +204,18 @@ setMethod(
var=x
# tested class
x<-gsub("\\(.+$","",x)
x<-gsub(" \\(.+$","",x)
# remove group
if("group"%in%names(Data)){Data[,group:=NULL]}
# if interaction
if(length(grep(":|\\*",x))==1){
if(length(grep(":|\\*|\\+",x))==1){
# add a group column for interaction
Data[,
"group":=paste(as.character(.SD),collapse=","),
.SDcols=strsplit(x,":|\\*")[[1]],
.SDcols=strsplit(x,":|\\*|\\+")[[1]],
by=1:nrow(Data)
]
......@@ -243,8 +243,8 @@ setMethod(
x=var,
y=x,
title="density plot",
fill=Class,
colour=Class
fill=var,
colour=var
)+
theme_bw()+
theme(
......@@ -278,7 +278,7 @@ setMethod(
setorderv(pvalues,"group")
# init anova tag
pval=Stats$anova_test$pval
pval=Stats$anova_test[group==x,pval]
# anova tag
text<-rep("ns",length(pval))
......@@ -288,7 +288,7 @@ setMethod(
text[pval<0.001]<-"p<0.001"
# add class
text=paste(paste(Stats$anova_test$group,text),collapse="\n")
text=paste(x,text,collapse="\n")
# box plot
bp<-ggplot(Data,aes(group,value,fill=group,colour=group)) +
......@@ -297,8 +297,8 @@ setMethod(
labs(
y=var,
title="box plot",
fill=Class,
colour=Class
fill=var,
colour=var
)+
theme_bw()+
theme(
......@@ -357,36 +357,6 @@ setMethod(
# remove special characters
var=gsub("/","-",var)
# display in 2 columns
if(length(Stats$plots)>1){
# print or not
if(print==TRUE){
# return for html
grid.arrange(
Stats$plots$qqplot,
Stats$plots$EMMplot,
dp,
bp,
dm,
ncol=1
)
}
# return in pdf by gene
if(!is.null(figs)){
pdf(paste(figs,"_",prefix,paste(variable,var,sep="_"),".pdf",sep=""),width=20)
print(Stats$plots$qqplot)
print(Stats$plots$EMMplot)
print(dp)
print(bp)
print(dm)
dev.off()
}
}else{
# print or not
if(print==TRUE){
......@@ -409,7 +379,6 @@ setMethod(
print(dm)
dev.off()
}
}
})
# return results in list
......
README.md
+ 1
− 1
View file @ 224a77e8
......@@ -24,7 +24,7 @@ remotes::install_gitlab(
devtools::build("eggbiomechanic")
# install package (from R console)
install.packages("eggbiomechanic_1.2.2.tar.gz", repos = NULL, type = "source")
install.packages("eggbiomechanic_1.2.3.tar.gz", repos = NULL, type = "source")
```
## Quick overview
......
man/stat_lme.Rd
+ 2
− 2
View file @ 224a77e8
......@@ -26,9 +26,9 @@ This function perform the statistical tests according a chain.
\details{
This function is used in internal by the \code{\link{var_test}} function in order to perform Linear mixed-effects models statistics:
\enumerate{
\item analyse of variance using \code{\link[nlme]{lme}}.
\item analyse of variance using \code{\link[stats]{anova}} on \code{\link[nlme]{lme}} fitting model.
\item residues studie using qqplot.
\item post Hoc test with \code{\link[emmeans]{emmeans}}.
\item Tukey post Hoc test with \code{\link[multcomp]{glht}}.
}
}
\examples{
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment