Probit et régression logarithmique

Version imprimable

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

Code:

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:

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:

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;

Code:

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:

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 :calim2:
Ici par exemple je trouve 13.31 au lieu de 10.2 :pleure:

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

05/04/2016, 10h26
anapurna

salut

les arrondi peut être non ? ... quand je vois plus de 5 décimal j'ai des doute sur la suite du calcul
05/04/2016, 10h29
PhilLU

Citation:

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...
05/04/2016, 10h35
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.
05/04/2016, 16h33
PhilLU

Citation:

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
05/04/2016, 17h30
PhilLU

1 pièce(s) jointe(s)

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)

Pièce jointe 206029

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...
:help:
05/04/2016, 18h28
anapurna
salut

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

Code:

=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:

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
05/04/2016, 18h51
PhilLU

Citation:

Envoyé par anapurna

salut

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

Code:

=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:

=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 :D
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...

:merci:
05/04/2016, 19h07
anapurna

salut

je ne sais pas c'est une éventualité
05/04/2016, 20h17
PhilLU

Citation:

Envoyé par anapurna

salut

je ne sais pas c'est une éventualité

Ben non, je n'obtiens plus de droites :?