Données disponibles';
$nb_station = count(array_keys($hype_data));
//$nb_parametre = count(array_unique (array_column($hype_data, 'CODE_PARAMETRE')));
//$nb_chronique = count(array_unique (array_column($hype_data, 'CHRONIQUE')));
//$output.=''.$nb_chronique.' chronique(s) correspondant à '.$nb_station.' point(s) de prélèvement et '.$nb_parametre.' paramètre(s) ont été lues.
';
$output.=''.$nb_station.' point(s) de prélèvement
';
//$output.= theme(
// 'table',
// [
// 'header'=>array_keys($hype_data[0]),
// 'rows'=>$hype_data
// ]
//);
$output.='Caractérisation';
//Tri par date du jeu de données
//array_multisort (array_column($hype_data, 'DATE_DEBUT_PRELEVEMENT'), SORT_ASC, $hype_data);
//$carac = hype_caracterisation($hype_data);
//dsm($carac);
$output.='Tendances/Ruptures';
//dsm($hype_data);
$tmp =[];
$tmp['06594X0019/S'] = $hype_data['06594X0019/S'];
//$tmp['06363X0017/SOURCE'] = $hype_data['06363X0017/SOURCE'];
//krumo($tmp);
$tendance = hype_tendance($tmp);
//$tendance = hype_tendance($hype_data);
krumo($tendance);
}*/
//============================================ Tendances
function hype_tendance($data){
$tendance = [];
//krumo($data);
//exit();
foreach($data as $code=>$stq){
foreach($stq as $cdpar=>$items){
//stat caracterisation generale
$chronique = $items['CHRONIQUE'];
//Tri par date du jeu de données
array_multisort (array_column($chronique , 'DATE_PREL'), SORT_ASC, $chronique );
//array_multisort (array_column($chronique , 'RESULTAT_TH'), SORT_ASC, $chronique );
$carac = hype_caracterisation_chronique($chronique);
//krumo($carac);
//Calcul des dates de rupture avec la P-value pour le test de Pettitt et pour le test de Buishand 1 si c'est significatif 0 si ça ne l'est pas
//test de Pettitt si données non normalement distribuée
//if(pval_shapi[i]<0.05&&is.na(pval_shapi[i])==FALSE){
//Si la pop est <= 3 on ne va même pas plus loin
if($carac['population']>3 && $carac['pval_shapi']){
//test de Pettitt si données non normalement distribuée
if($carac['pval_shapi'] < 0.05){
$carac['testrupt_graph']='Pettitt';
$carac['test_rupt']='test non paramétrique de Pettitt';
$Pettitt = hype_tendance_pettitt(array_column($chronique, 'RESULTAT_TH'));
$carac['rupt_pvalue'] = $Pettitt['pvalue'];
//Si pas rupture
if($Pettitt['pvalue']<=0.05){
$carac['pres_rupt'] = "non";
$carac['rupt_graph1'] = "Pas de rupture significative détectée";
$carac['pres_rupt'] = "";
$carac['rupt'] = null;
}
//SI rupture
else{
$carac['pres_rupt'] = "oui";
$carac['rupt'] = $chronique[$Pettitt['key']]['DATE_DEBUT_PRELEVEMENT']; //DATE_PREL
}
}
//test de Buishand si données normalement distribuées et population >=10
elseif($carac['pval_shapi']>=0.05 && $carac['population']>=10){
$carac['testrupt_graph']='Buishand';
$carac['test_rupt']='test paramétrique de Buishand';
$bu = hype_tendance_bu($data);
if($bu['pvalue']<=0.05){
$carac['pres_rupt'] = "oui";
$carac['rupt'] = $chronique[$bu['key']]['DATE_DEBUT_PRELEVEMENT'];
}else{
$carac['pres_rupt'] = "non";
$carac['rupt_graph1'] = "Pas de rupture significative détectée";
$carac['pres_rupt'] = "";
$carac['rupt'] = null;
}
}
else{
$carac['sig_ruptB']=null;
$carac['pval_rupt']=null;
$carac['rupt']=null;
$carac['pres_rupt'] = "trop peu de données";
$carac['rupt_graph1'] = "Non effectué (pas assez de données)";
$carac['rupt_graph2']='';
$carac['testrupt_graph']='';
$carac['test_rupt']='';
}
//krumo($carac);
}
//"chronique stationnaire"
else{
$carac['sig_ruptB']=null;
$carac['pval_rupt']=null;
$carac['rupt']=null;
$carac['pres_rupt'] = "chronique stationnaire";
$carac['rupt_graph1'] = "Non effectué";
$carac['rupt_graph2']='';
$carac['testrupt_graph']='';
$carac['test_rupt']='';
}
//Si les données sont autocorrélées
if($carac['autocorr_test']){
if($carac['population'] > 40){
$mk = hype_tendance_mk_modif(array_column($chronique , 'RESULTAT_TH'), 0.95);
$carac['pval_mkmodif'] = $mk['cpvalue'];
$carac['test_toto'] = $mk['pvalue'];
$carac['tau_mkmodif'] = $mk['tau'];
if($mk['cpvalue'] && $mk['cpvalue'] <= 0.05){
$carac['text_mkmodif'] = "Un test de Mann-Kendall modifié a été appliqué. Une tendance significative est détectée.";
}
else{
$carac['text_mkmodif'] ="Un test de Mann-Kendall modifié a été appliqué. Aucune tendance significative n'est détectée.";
}
$carac['tendMKmodif_graph1'] = "";
$carac['tendMKmodif_graph2'] = $mk['cpvalue'];
}
else{
$carac['text_mkmodif'] = "Les données sont autocorrélées mais il n'y a pas suffisamment de données pour appliquer un test de Mann-Kendall modifié.";
$carac['tendMKmodif_graph1'] = "Non effectué (pas assez de données)";
$carac['tendMKmodif_graph2'] = "";
$carac['pval_mkmodif'] = null;
}
}
else{
$carac['tendMKmodif_graph1'] = "Non effectué (pas assez de données)";
$carac['tendMKmodif_graph2'] = "";
$carac['text_mkmodif'] = "Les données ne sont pas autocorrélées. Le test de Mann-Kendall modifié n'a donc pas été appliqué.";
$carac['pval_mkmodif'] = null;
}
#Pour toutes les chroniques : Mann-Kendall
if(count(array_unique(array_column($chronique, 'RESULTAT_TH'))) > 1){
$carac['tend'] = "La chronique est stationnaire";
$carac['tendMK_graph1'] = "Chronique stationnaire";
$carac['tendL_graph1'] = "Chronique stationnaire";
$carac['tendMK_graph2'] = "";
$carac['tendL_graph2'] = "";
$carac['interc'] = null;
$carac['slope_Sen_an'] = null;
$carac['rcarre'] = null;
$carac['pval_lin'] = null;
$carac['pente_lin_an'] = null;
$carac['interc_lin'] = null;
$carac['tend_lin'] = "La chronique est stationnaire";
}
else{
#Mann-Kendall
//si au moins 10 analyses
if (count($chronique) >=10) {
$cor = hype_tendance_cor(array_column($chronique, 'RESULTAT_TH'), array_column($chronique, 'DATE_PREL'), 'everything','kendall');
$carac['pval_mk'] = $cor['pvalue'];
if($carac['pval_mk'] == 0 )$carac['pval_graph'] ="<2.2e-16";
else{$carac['pval_graph'] =$cor['pvalue'];}
$carac['estimate_mk'] = $cor['estimate'];
#Calcul pente de Sen
if ($carac['pval_mk']<=0.05){
//============================ TODO ========================================================================
/*
donneesx<-outer(donnees,donnees,"-")
datex<-outer(date_an,date_an,"-")
z<-donneesx/matrix(as.numeric(datex),ncol=nbr_anal[i])
s<-z[lower.tri(z)]
slope_Sen[i]<-median(s,na.rm=TRUE)
if (slope_Sen[i]==0) {
if(estimate_mk[i]<0) {slope_Sen[i]<-max(s[which(s<0)])}
else {slope_Sen[i]<-min(s[which(s>0)])}
}
slope_Sen_an[i]<-slope_Sen[i]*365
#Calcul intercept (Conover,1980)
interc[i]<-median(donnees,na.rm=TRUE)-slope_Sen[i]*as.numeric(median(date_an,na.rm=TRUE))
tend[i]<-paste("Un test de tendance de Mann-Kendall a été appliqué. Une tendance significative de pente",format(slope_Sen[i]*365,digits=3),unite,"/an est détectée.")
tendMK_graph1<-paste(format(slope_Sen[i]*365,digits=3,scientific=TRUE),unite,"/an")
*/
$carac['tendMK_graph2']=$carac['pval_graph'];
}
else{
$carac['tend']="Un test de tendance de Mann-Kendall a été appliqué. Aucune tendance significative n'est détectée";
$carac['tendMK_graph1']="Aucune tendance significative détectée";
$carac['tendMK_graph2']=$carac['pval_graph'];
$carac['interc']=null;
$carac['slope_Sen_an']=null;
}
}
else{
$carac['tend'] = "Il n'y a pas assez de données pour effectuer un test de tendance de Mann-Kendall.";
$carac['tendMK_graph1'] = "Non effectué (pas assez de données)";
$carac['tendMK_graph2'] = "";
$carac['interc'] = null;
$carac['slope_Sen_an'] = null;
$carac['pval_mk'] = null;
$carac['estimate_mk'] = null;
}
}
$tendance[$code.'_'.$cdpar] = $carac;
}
}
return $tendance;
}
function hype_tendance_cor($x, $y, $use, $method){
$mk_modif = db_query("
WITH cor AS (SELECT r_cor(ARRAY[".implode(', ', $x)."], ARRAY[".implode(', ', $y)."], '".$use."', '".$method."') as mk_modif)
SELECT cor[1] as pvalue, cor[2] as estimate
FROM test
")->fetchAssoc();
return $cor;
}
//Test de Mann Kendal modifié
function hype_tendance_mk_modif($data, $ci = 0.95){
$mk_modif = db_query("
WITH test AS (SELECT r_mk_modif(ARRAY[".implode(', ', $data)."], ".$ci.") as mk_modif)
SELECT mk_modif[1] as z, mk_modif[2] as pvalue, mk_modif[3] as zc, mk_modif[4] as cpvalue, mk_modif[5] as tau, mk_modif[6] as nns, mk_modif[7] as sen_slope,
FROM test
")->fetchAssoc();
return $mk_modif;
}
/*-------------------------------------------------------------------------------
* Breakpoint detection in a time series according to Pettitt 1979
* R script by Pascal Haenggi, v20090819
*-------------------------------------------------------------------------------
* x = time series with 1. column: year and 2. column: value
* alpha = significance level for test, e.g. 0.05
*/
function hype_tendance_pettitt($data){
$pettitt = db_query("
WITH test AS (SELECT r_pettitt(ARRAY[".implode(', ', $data)."]) as pettitt)
SELECT pettitt[1] as statistic, pettitt[2] as pvalue, pettitt[3] as nobs, pettitt[4] as key
FROM test
")->fetchAssoc();
return $pettitt;
}
//Buishand U Test (Buishand, 1984), Breakpoint detection in a time series
function hype_tendance_bu($data){
$bu = db_query("
WITH test AS (SELECT r_bu(ARRAY[".implode(', ', $data)."]) as bu)
SELECT bu[1] as statistic, bu[2] as pvalue, bu[3] as key
FROM test
")->fetchAssoc();
return $bu;
}
//============================================ Caracterisation
function hype_caracterisation($data){
$carac = [];
//krumo($data);
foreach($data as $code=>$stq){
foreach($data as $code=>$stq){
foreach($stq as $cdpar=>$items){
//id_chronique = $code.'_'.$cdpar
$carac[$code.'_'.$cdpar] = hype_caracterisation_chronique($items['CHRONIQUE']);
}
}
}
return $carac;
}
function hype_caracterisation_chronique($chronique){
//Stat générale ==> nécessite R
$stats = hype_caracterisation_general(array_column($chronique, 'RESULTAT_TH'));
//==============Taux de quantification - LQs
$lqs = hype_caracterisation_lqs($chronique, $stats);
//Calcul de la moyenne et de l'écart-type des écarts entre deux dates d'analyses
$prel = array_column($chronique, 'DATE_PREL');
sort($prel);
$ecarts = hype_caracterisation_get_ecarts($prel, (3600*24));
$ecarts_stat = hype_caracterisation_general($ecarts);
//identification d'outliers dans la fréquence de prélèvements
$ecarts_stat2 = hype_caracterisation_ecarts2($ecarts, $ecarts_stat);
//test de Shapiro-Wilk
$shapiro = hype_caracterisation_shapiro(array_column($chronique, 'RESULTAT_TH'));
//Autocorrélation des données
//FIX ME : le correlation doit-elle se faire sur les résultats d'analyse ou sur les dates de prélévement
$correlation = hype_caracterisation_correlation($chronique);
//dsm($stats);
return [
"date_min"=>date('Y-m-d', min(array_column($chronique, 'DATE_PREL'))),
"date_max"=>date('Y-m-d', max(array_column($chronique, 'DATE_PREL'))),
"duree"=>(max(array_column($chronique, 'DATE_PREL')) - min(array_column($chronique, 'DATE_PREL')))/(3600*24), //en jour
"population"=>count($chronique),
"moyenne"=>$stats['avg'],
"mediane"=>$stats['p50'],
"percentil1"=>$stats['p25'],
"percentil2"=>$stats['p75'],
"taux_quanti"=>$lqs["taux_quanti"],
"ecarttype"=>$stats['stddev'],
"pval_shapi"=>$shapiro['p.value'],
"normalite"=>$shapiro['normalite'],
"normalite_chart"=>$shapiro['chart'],
"LQmin"=>$lqs["LQmin"],
"LQmax"=>$lqs["LQmax"],
"LQmin2"=>$lqs["LQmin2"],
"LQmax2"=>$lqs["LQmax2"],
"text_med"=>$lqs["LQmessage_p50"],
"text_perc1"=>$lqs["LQmessage_p25"],
"text_perc2"=>$lqs["LQmessage_p75"],
"moy_ecarts"=>isset($ecarts_stat['avg'])?$ecarts_stat['avg']:null,
"med_ecarts"=>isset($ecarts_stat['p50'])?$ecarts_stat['p50']:null,
"ec_typ_ecarts"=>isset($ecarts_stat['stddev'])?$ecarts_stat['stddev']:null,
"moy_ecarts2"=>isset($ecarts_stat2['avg'])?$ecarts_stat2['avg']:null,
"ec_typ_ecarts2"=>isset($ecarts_stat2['stddev'])?$ecarts_stat2['stddev']:null,
"autocorr_r1"=>$correlation['autocorr_r1'],
"autocorr_test"=>$correlation['result'],
"signif_autocorr"=>$correlation['analyse'],
"signif_autocorr_short"=>$correlation['short'],
"test_unite"=>count(array_unique(array_column($chronique, 'UNITE_GRAPH')))>1 ? false:true, //Si une seule unité tout va bien (valeur = true), sinon mettte false qui permettra d'alerter l'utilisateur sur la fait que son jeu de données contient plusieurs unités
];
}
#identification d'outliers dans la fréquence de prélèvements
function hype_caracterisation_ecarts2($ecarts, $stat){
$e2 = [];
if(!empty($ecarts)){
$cap = ($stat['p50']+2) * $stat['stddev'];
$ecarts2 = array_filter($ecarts, create_function('$v', 'return $v>='. $cap .';'));
if(!empty($ecarts2)){
$e2 = hype_caracterisation_general($ecarts2);
}else{
$e2 = $stat;
}
}
return $e2;
}
function hype_caracterisation_get_ecarts($src, $divide = 1){
$ecarts = [];
if(!empty($src)){
$nb = count($src);
for($i=0;$inull,
"LQmax"=>null,
"LQmin2"=>null,
"LQmax2"=>null,
"taux_quanti"=>null,
"LQmessage_p25"=>'',
"LQmessage_p50"=>'',
"LQmessage_p75"=>''
];
if(!empty($chronique)){
$nb_ana = count($chronique);
$nb_quanti = 0;
foreach($chronique as $ana){
if($ana['CODE_SIGNE']==1){
$nb_quanti++;
}
else{
//LQmin
if(is_null($lqs['LQmin']) || $ana['RESULTAT']<$lqs['LQmin']){
$lqs['LQmin'] = $ana['RESULTAT'];
}
//LQmax
if(is_null($lqs['LQmax']) || $ana['RESULTAT']>$lqs['LQmax']){
$lqs['LQmax'] = $ana['RESULTAT'];
}
if($ana['CODE_SIGNE']==10){
//LQmin2
if(is_null($lqs['LQmin2']) || $ana['RESULTAT']<$lqs['LQmin2']){
$lqs['LQmin2'] = $ana['RESULTAT'];
}
//LQmax2
if(is_null($lqs['LQmax2']) || $ana['RESULTAT']>$lqs['LQmax2']){
$lqs['LQmax2'] = $ana['RESULTAT'];
}
}
}
}
#Remarque dans le cas où mediane/decile$stats['p25']){
$lqs['LQmessage_p25'] ="Attention : la valeur du premier quartile est inférieure à la limite de quantification maximale
";
}
if($lqs['LQmax']>$stats['p50']){
$lqs['LQmessage_p50'] ="Attention : la valeur de la médiane est inférieure à la limite de quantification maximale
";
}
if($lqs['LQmax']>$stats['p75']){
$lqs['LQmessage_p75'] ="Attention : la valeur du dernier quartile est inférieure à la limite de quantification maximale
";
}
$lqs['taux_quanti'] = ($nb_quanti/$nb_ana)*100; //on veux des %
}
return $lqs;
}
//Permet la caractérisation d'un jeu de données en utilisant les fonctions PLR (Postgres + R)
function hype_caracterisation_general($item){
//krumo($item);
$r =[];
if(!empty($item)){
//=================== stat de base
$r = db_query("
WITH data AS (SELECT unnest(ARRAY[".implode(', ', $item)."]) as src)
SELECT
min(data.src) as min,
max(data.src) as max,
avg(data.src) as avg,
variance(data.src) as variance,
stddev(data.src) as stddev,
r_quantile(array_accum(data.src), 0.25) as p25,
r_quantile(array_accum(data.src), 0.50) as p50,
r_quantile(array_accum(data.src), 0.75) as p75,
r_quantile(array_accum(data.src), 0.90) as p90
FROM data;
")->fetchAssoc();
$r['pop'] = count($item);
}
return $r;
}
#Autocorrélation des données
function hype_caracterisation_correlation($item){
$correlation = [];
//Tri par date du jeu de données
array_multisort (array_column($item, 'RESULTAT_TH'), SORT_ASC, $item);
#Autocorrélation au rang 1
$rank = 3; //n° de rang + 1 (dans acf l'index 0 vaut 1) + 1 dans pg, les index commencent à 1 pas à 0
$autocorr_r1 = db_query("
WITH acf AS (SELECT r_acf(ARRAY[".implode(', ', array_column($item, 'RESULTAT'))."], 'correlation', FALSE) as data)
SELECT acf.data[".$rank."][1][1] as rank
FROM acf
")->fetchColumn(0);
$correlation['autocorr_r1'] = $autocorr_r1;
if($autocorr_r1){
#Significativité de l'autocorrélation au seuil de 5%
$qnorm = db_query("SELECT r_qnorm((1 + 0.95)/2) as qnorm")->fetchColumn(0);
$seuil = $qnorm / sqrt(count($item));
if(abs($autocorr_r1) > abs($seuil)){
$correlation['analyse'] = "La probabilité que les données soient autocorrélées au rang 1 est supérieure à 95%";
$correlation['short'] = "Données autocorrélées";
$correlation['result'] = true;
}
else {
$correlation['analyse'] = "Il est peu probable que les données soient autocorrélées";
$correlation['short'] ="Données non autocorrélées";
$correlation['result'] = false;
}
}else{
$correlation['analyse'] = "L'autocorrélation des données n'a pas pu être testée.";
$correlation['short'] ="Données non autocorrélées";
$correlation['result'] = null;
}
return $correlation;
}
function hype_caracterisation_shapiro($item){
//krumo($item);
if(!empty($item)){
//===============Shapiro
//Si pop <= 3 on ne peut pas faire le test de shapiro
//Si toutes les valeurs sont identiques on ne pas peux pas faire le shapiro non plus
if(count($item)<=3){
$shapiro=[
'p.value'=>false,
'normalite'=>"Il n'y a pas assez de données pour estimer la normalité de la distribution",
'chart'=>"Test de normalité non effectué"
];
}
//Si toutes les valeurs sont identiques on ne pas peux pas faire le shapiro non plus
elseif(count(array_unique($item))==1){
$shapiro=[
'p.value'=>false,
'normalite'=>"Toutes les valeurs sont identiques. Le test de Shapiro-Wilk n'a pas été appliqué.",
'chart'=>"Chronique stationnaire"
];
}
//on tente le Shapiro
else{
$sh = db_query("SELECT r_shapiro(ARRAY[".implode(', ', $item)."]) as pvalue")->fetchAssoc();
if($sh['pvalue']>=0.05){
$shapiro=[
'p.value'=>$sh['pvalue'],
'normalite'=>"Un test de Shapiro-Wilk a été appliqué. Les données de cette chronique sont normalement distribuées",
'chart'=>"Données normalement distribuées"
];
}
else{
$shapiro=[
'p.value'=>$sh['pvalue'],
'normalite'=>"Un test de Shapiro-Wilk a été appliqué. Les données de cette chronique ne sont pas normalement distribuées",
'chart'=>"Données non normalement distribuées"
];
}
}
}
return $shapiro;
}