Probit et régression logarithmique

**PhilLU** · 05/04/2016, 10h15

Bonjour,
Je souhaite calculer la dose minimale détectable à l'aide de la fonction Probit:
Voici un exemple de données à traiter:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
 
HBV	rep.	positifs	Hit Rate	Frct%0-50	SubRoutine	Prob 0-50	Probit
0	14	0,00					
2,5	14	8,00					
5	14	10,00					
10	14	13,00					
15	14	14,00					
20	14	14,00					
25	14	14,00

Je calcule le Probit à l'aide du code suivant:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
 
for I := 1 to ProbitGrid.rowcount-1 do
if ProbitGrid.Cells[0,i]<>'' then
begin
    if StrToFloat(ProbitGrid.Cells[1,i])<>0
    then
    ProbitGrid.Cells[3,i]:=   FloatToStr(StrToFloat(ProbitGrid.Cells[2,i])/  StrToFloat(ProbitGrid.Cells[1,i])*100)
    else
    ProbitGrid.Cells[3,i]:='0';
   //SI(C3>50;(100-C3)/100;C3/100)
    if StrToFloat(ProbitGrid.Cells[3,i])>50
    then
    ProbitGrid.Cells[4,i]:=  FloatToStr((100-StrToFloat(ProbitGrid.Cells[3,i]))/100)
    else
    ProbitGrid.Cells[4,i]:= FloatToStr(StrToFloat(ProbitGrid.Cells[3,i])/100);
 
  //RACINE(LN(1/PUISSANCE(C4;2)))
    if SQR(StrToFloat(ProbitGrid.Cells[4,i]))<>0 then   ProbitGrid.Cells[5,i]:=  FloatToStr(SQRT(ln(1/SQR(StrToFloat(ProbitGrid.Cells[4,i]))))) else ProbitGrid.Cells[5,i]:='0';
 
  //+C5-((2,515517+0,802853*C5+0,010328*PUISSANCE(C5;2))/(1+1,432788*C5+0,189269*(PUISSANCE(C5;2)+0,001308*PUISSANCE(C5;3))))
    if StrToFloat(ProbitGrid.Cells[5,i])<>0
    then
     ProbitGrid.Cells[6,i]:=FloatToStr(ABS(StrToFloat(ProbitGrid.Cells[5,i]))
     -((2.515517+0.802853*StrToFloat(ProbitGrid.Cells[5,i])+0.010328*SQR(StrToFloat(ProbitGrid.Cells[5,i])))
     /(1+1.432788*StrToFloat(ProbitGrid.Cells[5,i])+0.189269*(SQR(StrToFloat(ProbitGrid.Cells[5,i]))+0.001308*Power(StrToFloat(ProbitGrid.Cells[5,i]),3)))))
     else  ProbitGrid.Cells[6,i]:='0';
 
     //SI(L5>50;5+P5;5-P5)
      if ProbitGrid.Cells[6,i]='0' then  ProbitGrid.Cells[7,i]:='0'
      else
      begin
         if StrToFloat(ProbitGrid.Cells[3,i])>50 then   ProbitGrid.Cells[7,i]:= FloatToStr(5+StrToFloat(ProbitGrid.Cells[6,i])) else
         ProbitGrid.Cells[7,i]:= FloatToStr(5-StrToFloat(ProbitGrid.Cells[6,i]));
      end;
end;

ensuite je dois calculer la régression logarithmique des données X (colonne 0) et Y (colonne 7) à l'aide du code suivant:
(proposé par http://www.developpez.net/forums/u75679/gilbert-geyer/)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 
procedure RegressLog(Nuage: TNuage; var B0, B1: Extended; var StrRes: string);
var SigX, SigY, SigLnX, SigYLnX, SigLnX2: Extended; i: integer;
begin
  SigX := 0.0; SigY := 0.0; SigLnX := 0.0; SigYLnX := 0.0; SigLnX2 := 0.0;
  for i := 0 to High(Nuage) do
  begin
    if Nuage[i].x <= 0.0 then
    begin
      StrRes := 'Non calculable car présence d''un point X <= 0 : EXIT';
      EXIT;
    end;
    SigX := SigX + Nuage[i].x;
    SigY := SigY + Nuage[i].y;
    SigLnX := SigLnX + ln(Nuage[i].x);
    SigYLnX := SigYLnX + Nuage[i].y * ln(Nuage[i].x);
    SigLnX2 := SigLnX2 + sqr(ln(Nuage[i].x));
  end;
  b1 := (SigYLnX - SigLnX * SigY / length(Nuage)) / (SigLnX2 - sqr(SigLnX) / length(Nuage));
  b0 := (SigY / length(Nuage)) - (b1 * SigLnX / length(Nuage));
  StrRes := 'y = ' + FloatToStr(b0) + ' + ' + FloatToStr(b1) + ' * Ln(x)';
end;

et

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
 
  j:=0;
for I := 1 to ProbitGrid.rowcount-1 do
If ProbitGrid.Cells[3,i]<>'' then
if (StrToFloat(ProbitGrid.Cells[3,i])>0) AND (StrToFloat(ProbitGrid.Cells[3,i])<100) then
  begin
   SetLength(Nuage, j+1);
   Nuage[j].X := StrToFloat(ProbitGrid.Cells[0,i]);
   Nuage[j].Y := StrToFloat(ProbitGrid.Cells[7,i]);
   //RichEdit1.Lines.Add(ProbitGrid.Cells[0,i]+'|'+ProbitGrid.Cells[7,i]);
   j:=j+1;
  end;
 
  RegressLog(Nuage, B0, B1, StrRes);
  equation.Caption:=StrRes;
  equation.Caption:='y = '+FormatFloat('0.0000',(B0))+' + '+FormatFloat('0.0000',(B1))+' * Ln(x)';

ensuite la limite de détection est renvoyée par la formule:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
LOD95.Caption:='LOD 95%= '+ FormatFloat('0.00',(exp((6.645-b0)/b1)));

Malheureusement, j'obtiens des résultats différents que ceux trouvés dans la littérature

Ici par exemple je trouve 13.31 au lieu de 10.2

Si qq a une idée...
Merci pour vos lumières

**anapurna** · 05/04/2016, 10h26

salut

les arrondi peut être non ? ... quand je vois plus de 5 décimal j'ai des doute sur la suite du calcul

**PhilLU** · 05/04/2016, 10h29

Envoyé par anapurna

salut

les arrondi peut être non ? ... quand je vois plus de 5 décimal j'ai des doute sur la suite du calcul

Disons que je n'arrondi que l'affichage de la régression... pour le calcul du LOD95 je n'arrondi rien...

**tourlourou** · 05/04/2016, 10h35

Au passage, tu devrais pê nommer des constantes pour remplacer les entiers désignant des colonnes du tableau par des intitulés parlants, facilitant la compréhension du code et des calculs.

**PhilLU** · 05/04/2016, 16h33

Envoyé par tourlourou

Au passage, tu devrais pê nommer des constantes pour remplacer les entiers désignant des colonnes du tableau par des intitulés parlants, facilitant la compréhension du code et des calculs.

Oui, de fait.
Voici les colonnes dans l'ordre pour ProbitGrid:
HBV =c0
rep. =c1
positifs =c2
Hit Rate =c3
Frct%0-50 =c4
SubRoutine =c5
Prob 0-50 =c6
Probit =c7

**PhilLU** · 05/04/2016, 17h30

J'ai refais le tout sur excel et j'obtiens la même chose.
Mon calcul d'après la formule de Gilbert Geyer comparée à la formule loi.normale.standard.inverse de Excel donne aussi la même chose.
(colonnes I et M)

Nom : Probit comparaison Excel.jpg
Affichages : 249
Taille : 278,1 Ko

Nom : Probit comparaison Excel.jpg
Affichages : 249
Taille : 278,1 Ko

Le but étant de calculer la valeur de X pour un Y (Probit) de 6.645
J'utilise donc la formule:
exp((6.645-b0)/b1)

En arrondissant les Hit Rate comme dans la littérature, je passe de 13.31 à 12 (mais bon, je dois trouver 10.2)
(voir CAP/CTM dans la publication: http://www.ncbi.nlm.nih.gov/pubmed/22535983 )

Pour info, j'ai le même type d'erreur pour d'autres évaluations...

**anapurna** · 05/04/2016, 18h28

salut

dans subroutine ne devrais tu pas faire un test pour éviter cette erreur de division par zéro
genre :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

  =SI(PUISSANCE(F5;2)=0;0;RACINE(LN(1/PUISSANCE(F5;2))))

et si tu fais ce test la colonne Probit va changer car toutes tes valeur de probit seront renseigné

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
=SI(E5>50;5+H5;5-H5)

car la valeur E seras a zéro ce qui implique que H aussi seras a zero donc 5-0 = 5

**PhilLU** · 05/04/2016, 18h51

Envoyé par anapurna

salut

dans subroutine ne devrais tu pas faire un test pour éviter cette erreur de division par zéro
genre :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

  =SI(PUISSANCE(F5;2)=0;0;RACINE(LN(1/PUISSANCE(F5;2))))

et sit ut fait ce test la colonne Probit va changer car toutes tes valeur de probit seront renseigné

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

=SI(E5>50;5+H5;5-H5)

car la valeur E seras a zéro ce qui implique que H aussi seras a zero donc 5-0 = 5

Merci Anapurna

Tu penses que je dois inclure l'ensemble des points y compris le 0 et le(s) 100 dans le calcul de la droite???

Je teste et te dis quoi...

**anapurna** · 05/04/2016, 19h07

salut

je ne sais pas c'est une éventualité

**PhilLU** · 05/04/2016, 20h17

Envoyé par anapurna

salut

je ne sais pas c'est une éventualité

Ben non, je n'obtiens plus de droites

Probit et régression logarithmique

Delphi

Discussions similaires

Partager

Partager