forked from equipe22/BioQuality
-
Notifications
You must be signed in to change notification settings - Fork 0
/
fun_breakpoints.R
50 lines (48 loc) · 1.69 KB
/
fun_breakpoints.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
# This function loads into the memory the mobile quantiles ####
loadquant <- function(nombio,type="Median", directory=paste0(rept2,"MovingQuantiles/") ){
kk <- nombio
nomfichier <- sub(pattern = ":",replacement = ".",x = nombio,perl = F)
if(file.exists(paste0(directory,nomfichier,"/",nomfichier,".",type,".csv"))){
movingquantile <- read.csv2(paste0(directory,nomfichier,"/",nomfichier,".",type,".csv") )
return(movingquantile)
}else{
stop("Compute moving quantiles before breakpoint detection")
}
}
# Breakpoint detection function derived from the breakpoint package ####
detectbreakpoint <- function(inputData){
# Input : moving quantile
ansmean2 <- cpt.mean(inputData$value, method="PELT")
nbp <- ncpts(ansmean2)
place <- inputData$start[cpts(ansmean2)]
#
if(nbp>0){
toint <- c(1,cpts(ansmean2),nrow(inputData))
tocalc <- length(toint)-1
meanlist <- rep(NA,tocalc)
# If multiple breakpoints detected, compute the mean by interval
for(i in 1:tocalc){
meanlist[i] <- mean(inputData$value[toint[i]:toint[i+1]], na.rm = T)
}
# list of computed means
meanlist <- meanlist[which(!is.na(meanlist))]
proplist <- rep(NA, length(meanlist)-1)
if(length(meanlist)>1){
for(j in 1:(length(meanlist)-1) ){
proplist[j] <- round(100*(meanlist[j]-meanlist[j+1])/meanlist[j],1)
}
}else{
proplist <- 0
}
}else{
proplist <- NA
meanlist <- NA
maxchange <- NA
}
# Output :
# nbp : number of breakpoint
# place : positions of breakpoint
# meanlist : list of means
# maxchange : max variation of means between intervals
return(list(nbp=nbp, place=place,meanlist=meanlist,maxchange=max(proplist, na.rm=T)))
}