Code-sources de fonctions personnelles de la librairie KefiR

bootreg <- function(reg,plot=T,analyse=T,conf.level=0.95,pval=0.05,iter=1000) {

  numind <- nrow(reg$model)

  indices <- c(1:numind ) # num if individus

  enregistrement <- 0

  erreur <- c() ; predictions <- c() ; verity <- c()

  for (i in c(1:iter)) {

    indices_training <- sample(indices,size=numind ,replace=T)

    indices_test <- setdiff(indices ,indices_training)

    if (length(indices_test) > 0) {

      training  <- reg$model[indices_training,]

      test   <-    reg$model[indices_test,]

      formula <- eval(reg$call[[2]])

      oldw <- getOption("warn")

      options(warn = -1)

      reg1 <- lm(formula=formula(formula),data=training,subset=indices_training)

      reg1 <- update(reg1,subset=indices_training)

      options(warn = oldw)

      if (enregistrement == 0) {

        pval_mdl <- pf(summary(reg1)$fstatistic[1],summary(reg1)$fstatistic[2],summary(reg1)$fstatistic[3],lower.tail=FALSE)

        names(pval_mdl) <- "p-value of model"

        if (length(c(pval_mdl,summary(reg1)[[4]][,4])) == (length(reg$coefficients)+1)) {

          enregistrement <- enregistrement + 1

          p_values <- c(pval_mdl,summary(reg1)[[4]][,4])  # p_values correspond aux valeurs Pr récupérés pour chaque variable et chaque cycle du bootstrap (ici 1000)

          # Pour extraction des p-values sur glm() : mettre cela : coef(summary(model2))[,4]

          coeff <- reg1$coefficients

          oldw <- getOption("warn")

          options(warn = -1)

          predictions <- c(predictions,predict(reg1,reg$model)[indices_test])

          options(warn = oldw)

          verity <- c(verity,test[,1])

        }

      } else {

        pval_mdl <- pf(summary(reg1)$fstatistic[1],summary(reg1)$fstatistic[2],summary(reg1)$fstatistic[3],lower.tail=FALSE)

        if (length(c(pval_mdl,summary(reg1)[[4]][,4])) == (length(reg$coefficients)+1)) {

          p_values <- rbind(p_values,c(pval_mdl,summary(reg1)[[4]][,4]))

          coeff <- rbind(coeff,reg1$coefficients) # Coeff correspond aux coefficients récupérés pour chaque variable et chaque cycle du bootstrap (ici 1000)

          oldw <- getOption("warn")

          options(warn = -1)

          predictions <- c(predictions,predict(reg1,reg$model)[indices_test])

          options(warn = oldw)

          verity <- c(verity,test[,1])

        }

      }

      # Test sur les predictions

      #predictions =c() ; verity = c()

    }}

  coeff <- na.omit(coeff) ; p_values <- na.omit(p_values)

  predverity <- data.frame(predictions ,verity); predverity <- na.omit(predverity)

  predictions  <- predverity[,1] ; verity  <- predverity[,2]

  confiance <- function(x,conf.level=0.99) { # seuil

    temp = sort(x) ; valeur_seuil = round(length(x)*conf.level)

    temp <- temp[valeur_seuil]

    return(temp)

  }

  mode <- function(x) {

    densite <- density(x)

    mode <- densite$x[which(densite$y==max(densite$y))]

    return(mode)

  }

  apply(coeff,2,median) -> coeff_median

  if (plot==T) {

    boxplot_Pr <- function(x,main="") {

      my_min <- min(c(apply(x,2,quantile)[2,]),0.0009)

      if (my_min <=0) {my_min -> 1e-20}

      boxplot(x,log="y",ylim=c(my_min,1),main=main)

      abline(h=0.05,col="red",lwd=2)

      abline(h=0.01,col="orange",lwd=2)

      abline(h=0.001,col="green",lwd=2)

    }

    hist_decrypt <- function(x,breaks=20,main="",sub="",conf.level=0.95) {

      densite <- density(x,na.rm=T) # créer une liste des densités

      hist(x,freq=F,col="#AAFFAA",ylim=c(0,max(densite$y)),breaks=breaks,main=main,cex=0.5,sub=sub) # il faut que freq=F

      lines(densite, col = "red",lwd=3) # Superposer la ligne

      abline(v=confiance(x,conf.level),col="black",lwd=5)

      abline(v=mean(x),col="red",lwd=3)

      abline(v=median(x),col="orange",lwd=2)

      abline(v=mode(x),col="green",lwd=1)

      legend("topright",col=c("black","red","orange","green"),

             c("p-value max in the confidence interval","mean","median","mode"),

             lwd=c(5,3,2,1))

    }

    layout(matrix(1:4,2,2))

    # p-value of the model

    #return(p_values)

    hist_decrypt(p_values[,1],sub="If the repartition is uniform:\ninsignificant model",breaks=20,

                 main="Distribution of the p-values of the model",conf.level=conf.level)

    # p-value of the model & coefficients

    boxplot_Pr(p_values,main="Distribution of the p-values of\nthe model and its coefficients")

    apply(coeff,2,function(x) {(x-median(x))/median(x)*100}) -> percent_coeff

    boxplot(percent_coeff,main="Fluctuation of coefficients (in %)",ylog=NA)

    abs(predictions-verity)-> temp

    by(temp,verity,confiance,conf.level=conf.level) -> CONF

    by(temp,verity,mean) -> MOY

    x <- as.numeric(names(MOY))

    plot(as.numeric(names(MOY)),CONF,type="l",lwd=3,col="red",

         xlim=c(min(x),max(x)),ylim=c(0,max(CONF)),

         xlab="Experimental values",ylab="Predictions",

         main="Average prediction error and\nmaximum error in the confidence interval")

    points(as.numeric(names(MOY)),MOY,type="l",lwd=2,col="black")

  }

  apply(p_values,2,median) -> p.values_median

  apply(p_values,2,confiance,conf.level=conf.level) -> p.values_max

  coeff_model <- c(NA,reg$coefficients)

  coeff_median <- c(NA,coeff_median)

  coeff_IC <- c(NA,apply(coeff,2,int.ech,conf.level=conf.level))

  synth <- rbind(p.values_median,p.values_max,coeff_model,coeff_median,coeff_IC)

  synth <- t(synth)

  rownames(synth)[1] <- "Model"

  if (analyse==T){return(synth)

  } else {

    if (max(synth[,2])>pval) {return(FALSE)

    } else {return(TRUE)}

  }

}

# Si j'ai un échantillon de tel taille, quel effet je serais capable de voir

check_ech <- function(n, mean,sd,conf.level=0.99,conf=0.95, iter=200,cut=100){

  # Version 02 - 17 novembre 2020

 # conf : confiance souhaitée sur les tests de Student

 # conf.level : fréquence de fois où l'on souhaite avoir la p-value souhaitée définie paar conf

 # 

 seq <- (10:-10)

 seq <- 10^(seq)

 #

 # conf

 conftest <- 1

 i <- 1

 res_min <- res_max <- 0

 #

 while ((conftest > conf) & (i < length(seq))) {

 #

  moy_sup = mean+mean*seq[i]

  moy_inf = mean-mean*seq[i]

  pval_sup <- c()

  pval_inf <- c()

  for (j in 1:iter) {

   pop_ref <- rnorm(n, mean = mean, sd = sd)

   pop_test_sup <- rnorm(n, mean = moy_sup, sd = sd)

   pop_test_inf <- rnorm(n, mean = moy_inf, sd = sd)

      pval_sup <- c(pval_sup, t.test(pop_ref,pop_test_sup)$p.value)

   pval_inf <- c(pval_inf, t.test(pop_ref,pop_test_inf)$p.value)

    }

  conftest_sup <- length(pval_sup[pval_sup<(1-conf.level)])/iter

  conftest_inf <- length(pval_inf[pval_inf<(1-conf.level)])/iter

  conftest <- min(conftest_sup,conftest_inf)

  if ((i > 1)&(conftest < conf)) {

   #cat("Limite de détection sur l'intervalle : ",seq[i-1]," à ",seq[i],"\n")

   res_min <- seq[i-1]

   res_max <- seq[i]

  }

  i <- i+1

 }

 if (conftest > 1) {cat("Il n'a pas été possible de trouver une variation significative avec un échantillon de taille : ",n,"\n")

 } else {

 if (res_min != res_max) {

  pas <- abs(mean*res_min-mean*res_max)/cut

  #cat("Le pas sera donc de : ",pas,"entre ",mean+mean*res_min," et ",mean+res_max*mean,"\n")

  i <- cut; conftest <- 1

  while ((i > -1) & (conftest > conf) ) {

   pvalsup <- c() ; pvalinf <- c()

   moysup = (mean+i*pas)

   moyinf = (mean-i*pas)

   #cat("Moyenne testée de ",moysup," et ",moyinf,"sur un pas de ",pas*i,"\n")

   for (j in 1:iter) {

    pop_ref <- rnorm(n, mean = mean, sd = sd)

    pop_test_sup <- rnorm(n, mean = moysup, sd = sd)

    pop_test_inf <- rnorm(n, mean = moyinf, sd = sd)

       pvalsup <- c(pvalsup, t.test(pop_ref,pop_test_sup)$p.value)

       pvalinf <- c(pvalinf, t.test(pop_ref,pop_test_inf)$p.value)

     }

   conftest_sup <- length(pvalsup[pvalsup<(1-conf.level)])/iter

   conftest_inf <- length(pvalinf[pvalinf<(1-conf.level)])/iter

   conftest <- min(conftest_sup,conftest_inf)

   if (conftest < conf) {   

    cat("Cet échantillon de ",n," individus est discriminant avec une p-value < à ",1-conf,"\n\tdans ",conf.level*100,"% des cas\n\tpour des différences de moyennes de +/- :\n")

    return((i-1)*pas)

   }

   #cat("Conftest : ",conftest, "pour pas[i]",i,"\n")

   i <- i-1

  }

  }}

}

cor.ech <- function(x,y,iter=100) {

  nb_return <- c()

 for (i in 1:40 ) {

 nb_final <- c()

 for (i in iter) {

 pvalue = 1

 nb = 0

 while (pvalue > 0.05) {

  nb <- nb+1

  indices <- 1:length(x)

  indices_temp <- sample(indices ,nb,replace=T)

  x1 <- c(x,x[indices_temp])

  y1 <- c(y,y[indices_temp])

  pval <- c() ; j=0

  for (i in 1:iter) {

   indices <- 1:length(x1)

   indices_temp <- sample(indices ,length(indices), replace=T)

   if (length(unique(x1[indices_temp])) > 2) {

    j=j+1

    x_temp <- x1[indices_temp]

    y_temp  <- y1[indices_temp]

    #print(x_temp)

    pval[j] <- cor.test(x_temp,y_temp)$p.value

   }

  }

  pvalue <- quantile(pval,probs=0.95)

  if (nb > 1000) {pvalue<-0}

 }

 nb_final <- c(nb_final,nb)

 }

 max(nb_return,(quantile(nb_final,probs=0.95,names=F)+length(x)))-> nb_return

 }

 cat("Pour avoir une corrélation significative,\nil aurait au moins fallu un nombre de valeurs de : \n")

 return(nb_return)

}

corrigraph <- function(data,pval=0.01,exclude=0.3, ampli=4) {

  # Fonction réalisée par Antoine Massé

 # Version 01

 # Janvier 2021

 # Merci de dire où vous l'avez trouvé !

 # pval : seuil de significativité à 0.95%

 # exclude : seuil d'exclusion des petites corrélations : valeur entre 0 et 1

 # ampli : Contraste de taille des vertices (ronds)

 # En projet : tenir compte des NA (et les permettre), permettre de mettre les Y et les X de chaque côté comme sur un réseau neuralnet.

 # Matrice de corrélation

 cor(data) -> cor_matrice

 # Matrice de p-values de ces corrélations

 require("psych")

 corr.test(data)$p -> pval_matrice # Matrice des p-values

 # Correction de la matrice de corrélation par les pvalues

 # Ne conserver que les corrélations significatives

 ifelse(pval_matrice<pval,1,0) -> pval_matrice 

 pval_matrice* cor_matrice -> mymat

 require(igraph)

 net <- graph_from_adjacency_matrix(mymat, weighted=T,mode="lower")

 net <- simplify(net, remove.multiple = T, remove.loops = TRUE) # élaguer les liens redondants

 net <- delete.edges(net, E(net)[ abs(weight) < exclude ])

 E(net)$colour <- ifelse(E(net)$weight<0,"red","blue")

 E(net)$weight <- abs(E(net)$weight)

 # Clusterisation

 clp <- cluster_optimal(net)

 class(clp)

 l <- layout_with_fr(net)

 plot(clp, net, layout = l,vertex.size=sqrt(betweenness(net))*ampli+10,vertex.color="yellow",

 edge.arrow.size =0,arrow.mode=0,edge.width=(abs(E(net)$weight)*2)^3, edge.color =E(net)$colour)

}

# Si j'ai un effet de tant, quel taille d'échantillon me permettra de le voir

# Version 02 - 17/11/2020

identify_ech <- function(mean1,mean2,sd1,sd2,iter=500,conf.level = 0.99,conf=0.95,nmin=10,nmax=1000) {

  # conf.level : confiance attendue sur le test de Student lorsque je compare 2 échantillons

  # conf : Pourcentage de fois où j'espère avoir une p-value qui réponde à mon test si il y a bien un effet lorsque je compare 2 échantillons

 conftest <- 0 ;  intervalle <- 1 ; test <- 0 ; nprec <- 0 ;

 nmint = nmin ; nmaxt = nmax

 while ((intervalle >= 1) & (nprec < nmax) & (test==0)) {

  #cat("boucle1\n")

  n <- nmint

  intervalle <- round((nmaxt-nmint)/9,0)

  if (intervalle < 1) {intervalle <-1}

  if(intervalle == 1) {test = test+1}

  cat("Pas testé : ",intervalle,"\n")

  while ((conftest < conf) & (n <= nmax)) {

   #cat("boucle2\n")

   pval <- c()

   for (i in 1:iter) {

    pop_ref <- rnorm(n, mean = mean1, sd = sd1)

    pop_test <- rnorm(n, mean = mean2, sd = sd2)

    pval <- c(pval,t.test(pop_ref,pop_test)$p.value)

   }

   #print(pval)

   conftest <- length(pval[pval<(1-conf.level)])/iter

   #cat("conftest",conftest,"\n")

   nprec_prec = n-intervalle

   nprec <- n

   n <- n+intervalle 

  }

  cat("Boucle suivante : de ",nprec_prec," à ",nprec,"\n")

  nmint = nprec_prec

  nmaxt = nprec

  confdef <- conftest

  conftest <- 0

 }

 if (nprec == nmax) {

  cat("Warning : La taille d'échantillon requise nécessite de dépasser la population max indiquée.\n")

 } else {cat("Pour détecter un effet dans les conditions expérimentales,\nla population doit avoir une taille minimale de : ",nprec,"\n")}

 cat("Confiance atteinte : ",confdef,"\n")

}

# Fonction à copier-coller dans R 

int.ech = function(vector,conf.level=0.95,na.rm=T) { # VERSION 2016 

     if (length(vector)==0) { cat("Erreur ! Le vecteur ",substitute(vector),"est vide.\n")} 

      else { s = var(vector,na.rm=na.rm) 

      n = length(vector)-sum(is.na(vector)) 

      ddl = n - 1 ; proba = (1-conf.level)*100 ; proba = (100-proba/2)/100 

      t_student = qt(proba, ddl) # quantile 

       intervalle = t_student * sqrt(s/n) 

      moyenne = mean(vector,na.rm=na.rm) ; return(intervalle) }} 

i

int.pop =function(vector,conf.level=0.95 ,sigma=c() ) { 

    # vector : un échantillon 

    # sigma : écart-type de la population 

    if (length(sigma)==0) {sigma=sd(vector)} 

    n = length(vector) 

    proba =(100-(((1-conf.level)*100)/2))/100 

    z = qnorm(proba) ; # Récupération des quantiles de la loi normale ==> Par exemple : si on veut englober 95% des valeurs : estimer la variabilité de la population en prenant le risque de négliger 5% de celle-ci, spot 2,5% en tête de distribution 

    intervalle = z*sigma/sqrt(n) 

    moyenne = mean(vector) 

    return(intervalle) } 

int.prop = function(proportion,n,conf.level) { 

    # Il faut donner une proportion entre 0 et 1 et une taille d'échantillon 

    # proportion peut aussi correspondre à une liste de proportions 

    nombre_proportions = length(proportion) 

     P = proportion 

    proba = (1-((1-conf.level)/2))  

    qnorm(proba) 

    IC = c(qnorm(proba)*sqrt(P*(1-P)/n)) 

    return(IC) } 

# Il suffit de copier coller cette commande meanbp telle quelle pour calculer ensuite des moyennes mobiles par itération

meanbp = function(vector, it=1000,na.rm=T) {

  #Methode cf. Mosteller et Tukey, 1977 : 

 #permet de pondérer les valeurs en fonctions de la distance qui les sépare de la médiane

 # Version 2018

if (na.rm==T) {vector <- vector[!is.na(vector)]}

    if (length(vector) >2) {

        EIQ = abs(quantile(vector,probs=0.75)-quantile(vector,probs=0.25))

  if (EIQ ==0) {moyenne = mean(vector)}

 else {

         Z = (vector - median(vector))/(3*EIQ)

        for (i in c(1:length(vector))) {

           if (abs(Z[i])>1) {Z[i]=1}        }

        w = (1-Z^2)^2

        moyenne = sum(vector*w)/sum(w)

        #Methode complément A.M. : permet de pondérer les valeurs en fonctions de la distance qui les sépare de la moyenne calculée précédement (1000 cycles maximum)

        for (j in c(1:it)) {

            moyenne_old = moyenne

                d = abs(vector-moyenne)

                EIQ = (quantile(d,probs=0.75))+(quantile(d,probs=0.25))

                  Z = (vector - moyenne)/(3*EIQ)

                for (i in c(1:length(vector))) {

                   if (abs(Z[i])>1) {Z[i]=1}

                }

                w = (1-Z^2)^2

                moyenne = sum(vector*w)/sum(w)

            variation = sd(vector)/(moyenne_old-moyenne)

            if (moyenne ==moyenne_old) { break } # si moyenne ne bouge plus

            if (variation > 1000) { break} # si le signal sur bruit est élevé

        }}   }

    else {moyenne = mean(vector)}

    return(moyenne)}

# Exemple

x <- rnorm(100,20,5) # simuler un échantillon d'une population dont la moyenne est de 20

meanbp(x,100) # 100 car on ne refera au maximum que 100 fois le calcul

# si on compare avec mean(x) ou median(x) : on voit que meanbp rend un résultat plus proche de la véritable moyenne qui est de 20

require("plotly")

require("htmlwidgets")

# VERSION 02 janvier 2021

# Mettre titre d'une colonne en variable

parco <- function(data,Y=c(),X=c(),save=F,file="file.html") {

  colnames(data)<-gsub("-",".",colnames(data))

  # data : une data.frame

  # Y : un numéro de colonne ou un titre

  # X : un vecteur contenant ou des numéros de colonnes ou des titres

  if (is.numeric(Y) == T) {

      mon_titre = paste("Variable :",colnames(data)[Y])

      j <- colnames(data)[Y]

  }else {

  Y<-gsub("-",".",Y)

      mon_titre = paste("Variable :",Y)

      j <- Y

  }

  if (length(X)==0) {X <- 1:ncol(data)}

  if (is.numeric(X) == T) {

   dims <- list() ; k<-0

      for (i in X) {

     k<-k+1

        dims[[k]]<- list(range = range(data[,i],na.rm=T),

       label = colnames(data)[i],values = formula(paste0("~",colnames(data)[i])))

      }

  }else {

   dims <- list() ; k<-0

      for (i in X) {

     k<-k+1

         dims[[k]]<- list(range = range(data[[i]],na.rm=T),

       label = i,values = formula(paste0("~",i)))

      }

  } 

  # Cette version a été réalisée par Antoine Massé

  # Site  :Aide à l'utilisation de R

  # Merci d'en indiquer la source si vous la recopiée

  # parallele coordinate plot sur variable

  fig <- data %>% plot_ly(type = 'parcoords',

    line = list(color = ~get(j),

    colorscale = 'Jet',

    showscale = TRUE,

    reversescale = TRUE,

    cmin = min( data[[j]]),

    cmax = max( data[[j]])),

    dimensions = dims) %>% layout(title = mon_titre)

  if (save==F) {print(fig)}

  if (save==T) {saveWidget(fig, file, selfcontained = F, libdir = "lib")}

}

install.packages("plotly")

library("plotly")

install.packages("plot3D");install.packages("rgl");install.packages("plot3Drgl")

library("plot3D") ; library("rgl");library("plot3Drgl")

pde <- function(x,y,z,xlim=c(),ylim=c(),zlim=c(),dim=45,xlab="",ylab="",

zlab="",main="3D plot",col=c("blue","cyan","yellow","red","#990000"),alpha=0.1,pch=16,

col.pt="blue",cex.pt=6,mode=1) {

  # limites des axes des graphiques

 if (length(xlim) == 0) {xlim <- range(x)}

 if (length(ylim) == 0) {ylim <- range(y)}

 if (length(zlim) == 0) {zlim <- range(z)}

 # dim=nbre de points de la nappe

 # Débogage

 if(cor(x,y)==1) {print("Les vecteurs ne doivent pas converger afin de se répartir sur une surface 2D.\n") }

 if((length(x)!=length(y)) | (length(x)!=length(z)) | (length(unique(paste(x,y)))<6) | (length(unique(x)) < 3) | (length(unique(y)) < 3)){

  cat("Ca ne peut pas marcher pas.\n")

 if((length(x)!=length(y)) | (length(x)!=length(z))) {print("Les tailles des vecteurs ne sont pas identiques.\n") }

 #if (length(unique(x))<6) {print("Il n'y a pas le minimum de 6 valeurs uniques en x.") }

 #if (length(unique(y))<6) {print("Il n'y a pas le minimum de 6 valeurs uniques en y.") }

 if (length(paste(x,y))<6) {print("Il n'y a pas le minimum de 6 combinaisons x et y différentes.") }

  if ((length(unique(x))) < 3) {print("Attention x varie trop peu !")}

  if ((length(unique(y))) < 3) {print("Attention y varie trop peu !")}

 }else{

 # Régression linéaire sur les résultats hédoniques en fonction des différentes valeurs de l'expérience

 lm(formula = z ~ I(x^2) + I(y^2) +I(x*y) + x + y ) -> lmpoly 

 cat( "Equation de la nappe dans lmpoly : \n" )

 print(lmpoly)

 if(length(lmpoly$coefficients[is.na(lmpoly$coefficients)==TRUE])>0){

 cat("ATTENTION ! Les paramètres x et y donnés sont incompatibles avec une régression 3D.\nARRET DE LA FONCTION PDE()...")

     plot(x,y,main="Répartition de vos valeurs",sub="Ces valeurs ne sont pas dispersées de façon assez irrégulière pour permettre une régression 3D.",cex=(z-min(z))/max(z-min(z))*4+1,col="blue",pch=16)

} else {

 # Fabrication de la grille de la nappe

 x.grid <- seq(xlim[1],xlim[2],length.out=dim)

 y.grid <- seq(ylim[1],ylim[2],length.out=dim)

 xy<-expand.grid(x=x.grid,y=y.grid)

 # Calcul des ordonnées des points de la nappe

 z_nappe<-with(xy,

     lmpoly$coefficients[2]*x^2+

     lmpoly$coefficients[3]*y^2+

     lmpoly$coefficients[4]*x*y+

     lmpoly$coefficients[5]*x+

     lmpoly$coefficients[6]*y+

     lmpoly$coefficients[1]

     )

 z.mat<- matrix(z_nappe,nrow=dim, ncol=dim,byrow=F)

 #Maximum de la nappe

 max <- max(z_nappe)

 #Détermination des coordonnées x et y du maximum de la nappe 

 xy.max  =c()

 xy.max$x <- xy[,1][z_nappe==max]

 xy.max$y <- xy[,2][z_nappe==max]

cat("Coordonnées du sommet de la nappe :\n")

print(xy.max)

# 

require("plotly") ; require("plot3D") ; require("rgl");require("plot3Drgl")

#############################################

#   CREATION DU GRADIENT DE COULEURS

#############################################

#colorscale = cbind(seq(0, 1, by=1/(length(z) - 1)), rev(rainbow(length(z))))

colfunc<-colorRampPalette(col)

colors <- (colfunc(dim))

#col2hex <- function(col, alpha=1) rgb(t(col2rgb(col)), alpha=alpha*255, maxColorValue=255) 

#colors <- col2hex(colors,alpha=alpha)

colors_tp <- colfunc(length(z)*2)

pas = (1/(length(z)*2-1))

colorscale <- cbind(seq(0, 1, by=pas),colors_tp)

#############################################

if (mode==1) {

print("mode 1")

 require(plotly)

 # Mise en forme

 font <- list(    family = "Courier New, monospace",    size = 12,    color = "#7f7f7f" )

 xlabel <- list(    title = paste(xlab," (x)"),    titlefont = font )

 ylabel <- list(    title = paste(ylab," (y)"),    titlefont = font )

 zlabel <- list(    title = paste(zlab," (z)"),    titlefont = font )

 scene = list(    xaxis = xlabel ,    yaxis = ylabel ,    zaxis = zlabel ) 

 plot.pde <- plot_ly(x = x.grid, y = y.grid, z = z.mat, colorscale = colorscale) %>%  

  layout(title = main, scene = scene)  %>% add_surface() %>% 

  add_trace(x = x, y = y, z = z,name="Data brutes",

        marker = list(

        size = cex.pt,

        color = col.pt,

        line = list( color = 'rgb(231, 99, 250)', width = 2 ) 

        )

  )

  #plot_ly(x = x, y = y, z = z, size = I(3))

 return(plot.pde)

 # Autres sources : https://plot.ly/r/line-and-scatter/

}

if (mode == 2) {

#############################################

#  GRAPHIQUE 2D

#############################################

#Représentation de ce maximum

require(c("plot3D","rgl"))

scatter3D(x=xy.max$x,y=xy.max$y,z=max, 

    col='red',    pch=16,

    cex=1.5,    surf=NULL,

    xlim=xlim,    ylim=ylim,

    zlim=zlim,    main =main,

    xlab=xlab,     ylab=ylab,

    zlab=zlab    )

# Projection du maximum sur les axes x et y

segments3D(x0=xy.max$x,y0=xy.max$y,z0=max,x1=xy.max$x,y1=xy.max$y,z1=0,lty='dashed',col='red',add=TRUE)

segments3D(x0=xy.max$x,y0=xy.max$y,z0=0,x1=xy.max$x,y1=0,z1=0,lty='dashed',col='red',add=TRUE)

segments3D(x0=xy.max$x,y0=xy.max$y,z0=0, x1=0,y1=xy.max$y,z1=0, lty='dashed',col='red',add=TRUE)

# Ajout des légendes

text3D(0,xy.max$y,0,paste('ymax=',round(xy.max$y,2)),col='red',add=TRUE)

text3D(xy.max$x,0,0,paste('ymax=',round(xy.max$x,2)),col='red',add=TRUE)

# Représentation des expériences  

scatter3D(x, y, z,pch=pch ,cex=1,

    #sphere.     #size=5     #surface.alpha=0.5,    #revolution=300,

   # colvar = NULL,     #col='blue',

col=colors,

    surf=list(

        x=x.grid,y=y.grid,z=z.mat,

        facets = NA, col=colors,

 #       colvar=z.mat,

        alpha=0.1    

        ),

    xlab="", ylab="",zlab="",add=TRUE)

scatter3D(x, y, z,

    pch=pch ,cex=1,col="black",

   #xlab="",ylab="",zlab="",

    add=TRUE

    )

###################################################

#  GRAPHIQUE 3D

###################################################

Unaccent <- function(text) {

 text <- gsub("Ç","C",text); text <- gsub("ç","c",text); text <- gsub("é","e",text); text <- gsub("È","e",text);

 text <- gsub("î","i",text); text <- gsub("ï","i",text); text <- gsub("è","e",text); text <- gsub("ê","e",text);

 text <- gsub(" ","_",text); text <- gsub("-","_",text);   text <- gsub("['`^~\"]", " ", text)

   text <- iconv(text, to="ASCII//TRANSLIT//IGNORE");   text <- gsub("['`^~\"]", "", text)

   return(text)

}

main <- Unaccent(main);xlab <- Unaccent(xlab);ylab <- Unaccent(ylab);zlab <- Unaccent(zlab);

require("plot3Drgl")

#Représentation de ce maximum

scatter3Drgl(x=xy.max$x,y=xy.max$y,z=max, 

    col='red',    pch=16,

    cex=1.5,    surf=NULL,

    xlim=xlim,    ylim=ylim,

    zlim=zlim,    main =main,

    xlab=xlab,     ylab=ylab,

    zlab=zlab    )

# Projection du maximum sur les axes x et y

segments3Drgl(x0=xy.max$x,y0=xy.max$y,z0=max,x1=xy.max$x,y1=xy.max$y,z1=0,lty='dashed',col='red',add=TRUE)

segments3Drgl(x0=xy.max$x,y0=xy.max$y,z0=0,x1=xy.max$x,y1=0,z1=0,lty='dashed',col='red',add=TRUE)

segments3Drgl(x0=xy.max$x,y0=xy.max$y,z0=0, x1=0,y1=xy.max$y,z1=0, lty='dashed',col='red',add=TRUE)

#print("Attente 01");Sys.sleep(5)

# Ajout des légendes

text3Drgl(0,xy.max$y,0, paste('ymax=',round(xy.max$y,2)),col='red',add=TRUE)

text3Drgl(xy.max$x,0,0,paste('ymax=',round(xy.max$x,2)),col='red',add=TRUE)

#print("Attente 02");Sys.sleep(5)

# Représentation des expériences  

persp3Drgl(x=x.grid,y=y.grid,z=z.mat,col=colors,alpha=alpha ,add=T)

scatter3Drgl(x, y, z,

    pch=pch ,cex=1,col="black",

   #xlab="",ylab="",zlab="",

    add=TRUE

    )

}

}}}

rr = function(Var=c(),optr=c(),app=c(),sigma=5.15) {

means = function(data=c(),parametres=c()) {

   temp = split(data,parametres,drop=TRUE) #drop permet de ne faire des calculs que sur les paramètres présentant des valeurs associées

   x = c() ; y = unique(parametres) ; z = c() ; w = c() ; resultat = list()

   for (i in (1:(length(y)))) {

x = c(x,mean(temp[[i]],na.rm = TRUE))

if (length(temp[[i]]) > 1) {

z = c(z,sd(temp[[i]],na.rm = TRUE)) ;  w = c(w,sd(temp[[i]],na.rm = TRUE)*1.96/sqrt(length(temp[[i]])))}

      else {z = c(z,0);w = c(w,0)}

   }   

    resultat$moyennes = x[order(y)] ; resultat$sd = z[order(y)] ; resultat$ic = w[order(y)] ; resultat$parametres = y[order(y)] #Intervalles de confiance

    return(resultat)     

}

d2_constants = c(

1.41,1.91,2.24,2.48,2.67,2.83,2.96,3.08,3.18,3.27,3.35,3.42,3.49,3.55 ,1.28,1.81,2.15,2.40,2.60,2.77,2.91,3.02,3.13,3.22,3.30,3.38,

3.45,3.51, 1.23,1.77,2.12,2.38,2.58,2.75,2.89,3.01,3.11,3.21,3.29,3.37,3.43,3.5, 1.21,1.75,2.11,2.37,2.57,2.74,2.88,3.00,3.1,3.2,3.28,3.36,

3.43,3.49, 1.19,1.74,2.1,2.36,2.56,2.78,2.87,2.99,3.1,3.19,3.28,3.36,3.42,3.49, 1.18,1.73,2.09,2.35,2.56,2.73,2.87,2.99,3.10,3.19,3.27,3.35,

3.42,3.49, 1.17,1.73,2.09,2.35,2.55,2.72,2.87,2.99,3.1,3.19,3.27,3.35,3.42,3.48, 1.17,1.72,2.08,2.35,2.55,2.72,2.87,2.98,3.09,3.19,3.27,3.35,

3.42,3.48, 1.16, 1.72, 2.08, 2.34, 2.55, 2.72, 2.86, 2.98, 3.09, 3.19, 3.27, 3.35, 3.42, 3.48, 1.16, 1.72, 2.08, 2.34, 2.55, 2.72, 2.86, 2.98, 3.09, 3.18, 3.27, 3.34, 3.42, 3.48, 1.15, 1.71, 2.08, 2.34, 2.55, 2.72, 2.86, 2.98, 3.09, 3.18, 3.27, 3.34, 3.41, 3.48, 1.15, 1.71, 2.07, 2.34, 2.55, 2.72, 2.85, 2.98, 3.09, 3.18, 3.27, 3.34, 3.41, 3.48, 1.15, 1.71, 2.07, 2.34, 2.55, 2.71, 2.85, 2.98, 3.09, 3.18, 3.27, 3.34, 3.41, 3.48, 1.15, 1.71, 

2.07, 2.34, 2.54, 2.71, 2.85, 2.98, 3.09, 3.18, 3.27, 3.34, 3.41, 3.48, 1.15, 1.71, 2.07, 2.34, 2.54, 2.71, 2.85, 2.98, 3.08, 3.18, 3.26, 3.34, 3.41, 3.48, 1.128, 1.693, 2.059, 2.326, 2.534, 2.704, 2.847, 2.97, 3.078, 3.173, 3.258, 3.336, 3.407, 3.472)

d2_constants = matrix(d2_constants,nc=14, nr=16, byrow=T)

# d2 se déduit d'une grille avec z et w en relation : 

# z : produit du nombre d'opérateurs fois le nombre de pièces

# w : nombre de répétition par test

colnames(d2_constants) = c(2:15) #W

rownames(d2_constants) = c(1:15,">15") #z

if ((length(Var)==0)|(length(optr)==0)|(length(app)==0)) {stop("Un des paramètres est un vecteur vide.\n")}

operateurs = unique(optr)

appareils = unique(app)

n = length(appareils)

r = length(Var[optr==operateurs[1]&app==appareils[1]]) #repetition

o = length(operateurs)

W = r-1 ; z = n*o

if (z > 15) {z=16}

d2 = d2_constants[z,W] #d2 = ss.cc.getd2(repetition) avec le package SixSigma

# Calcul de la répétatibilité

moyenne_par_operateur = c()

étendue = c()

for (op in operateurs) {

for (i in appareils) {

étendue = c(étendue,(max(Var[optr==op&app==i])-min(Var[optr==op&app==i])))

}

moyenne_par_operateur = c(moyenne_par_operateur,mean(Var[optr==op]))

}

#cat("étendues : ",étendue[1:10],"\n");cat("étendues : ",étendue[11:20],"\n"); cat("étendues : ",étendue[21:30],"\n")

R = mean(étendue)

répétabilité = sigma*R/d2

W = o-1 # W nombre d'opérateurs

z = 1

d2 = d2_constants[z,W] #d2 = ss.cc.getd2(repetition) avec le package SixSigma

# Calcul de la reproductibilité

X = max(moyenne_par_operateur)-min(moyenne_par_operateur) # étendue entre opérateurs

#reproductibilité = sqrt((sigma*X/d2)^2-((répétabilité^2)/(n*r)))

A = (sigma*X/d2)^2

B = ((répétabilité^2)/(n*r))

if (B>A) {stop("\nLa répétabilité trop importante ne permet pas le calcul de la reproductibilité.\n")}

reproductibilité = sqrt(A-B)

#cat("R : ",R,"\n");cat("Répétabilité : ",répétabilité ,"\n");cat("d2 : ",d2,"\n");cat("moyenne par opérateur : ",moyenne_par_operateur,"\n");cat("X : ",X,"\n");cat("Reproductibilité : ",reproductibilité,"\n")

R_R = sqrt(répétabilité^2+reproductibilité^2)

moyenne_par_appareil = means(Var,app)

Rp = max(moyenne_par_appareil$moyennes) - min(moyenne_par_appareil$moyennes) # Rp étendue entre l'appareil ayant une mesure moyenne max et celui ayant une mesure moyenne min

W = n-1 # W nombre d'appareils

z = 1

d2 = d2_constants[z,W] #d2 = ss.cc.getd2(repetition) avec le package SixSigma

Vp = sigma*Rp/d2

Vt = sqrt(R_R^2+Vp^2)

RR=c()

RR$répétabilité = répétabilité

RR$reproductibilité = reproductibilité

RR$RR = R_R

RR$Vp = Vp 

RR$Vt = Vt

RR$part_Equipement = répétabilité/Vt*100

RR$part_Opérateurs = reproductibilité/Vt*100

RR$part_R_R = R_R/Vt*100

RR$part_Pièces = Vp/Vt*100

cat("\nAvec Vp : variabilité pièce et Vt : variabilité totale :\n\n")

# Ajouter la variabilité totale --> capabilité

return(RR)

}

###################################

# Valideur de régression

library(lmtest) ; library(car)

valreg <- function(reg,analyse=F,nvar=5) {

  error <- "OK"

 nvar1 <- length(coef(reg))

 if (nvar1 > nvar) {if(analyse==T){cat("Plus de 6 variables\n")};error = "error"

 } else if ( (length(summary(reg)[[4]][,4])) != (length(coef(reg))) ) {error = "error"

 } else {

  raintest(reg)$p.value -> pval # test de rainbow Adéquation

  if (pval < 0.05) {if(analyse==T){cat("Mauvaise adéquation.\n")};error = "error"}

  dwtest(reg)$p.value -> pval # Indépendance des résidus DurbinWatson

  if (pval < 0.05) {if(analyse==T){cat("Mauvaise indépendance des résidus.\n")};error = "error"}

  shapiro.test(residuals(reg))$p.value->pval # Distribution normale des résidus

  if (pval < 0.05) {if(analyse==T){cat("Distribution non normale des résidus.\n")};error = "error"}

  if (length(coef(reg))>=2) {  

   bptest(reg)$p.value -> pval # Breush Pagan : variance constante des résidus

   if (pval < 0.05) {if(analyse==T){cat("Variance non constante des résidus.\n")};error = "error"}

  }

  cooks.distance(reg)->cooksd

  if (max(cooksd,na.rm=T) > 8/length(cooksd)) {if(analyse==T){cat("Effet de levier.\n")};error = "error"}

 }

 return(error)

}

valreg(reg,nvar=5,analyse=T)