Dear R users,
Hopefully someone can help me, 
Maybe I just misunderstand the function in the package?
I am working with a logistic regression model.
Until now I always worked with the basic glm function, where for the model 
was:
бз glm( disease ~  test.value + cnct , family=binomial(link=бжlogitбж) би.

This works fine when test .value and concentration (cnct) are continuous 
vairables.
However, concentration is in fact a grouping variable over 5 experiments 
with 5 concentrations
( 25, 50, 100, 200 & 400).

Therefore I believe concentration to be an ordered factor ( in model : 
cnct_o).
To make this model I used the бзrmsби library (previously known as Design) 
and functions lrm (or Glm).
The lrm (or Glm) returns the odds for disease, the бзinv.logit (odds) би 
gives the probability of disease, but I have to do this with the Predict 
function of the бзrmsби package.

#####
The resulting model (with lrm or Glm) would be :
                Coef    S.E.    Wald Z  Pr(>|Z|)
Intercept       23.800  0.8891  2.68    0.0074
test.value      20.806  0.3409  6.10    <0.0001
cnct_o          -0.1127 0.0268  -4.21   <0.0001
cnct_o=100      77.393  17.542  4.41    <0.0001
cnct_o=200      204.291 45.080  4.53    <0.0001
cnct_o=400      427.829 98.180  4.36    <0.0001
#####
The results of the standard glm function  are very different :

#####
Standard glm
Deviance Residuals: 
    Min       1Q   Median       3Q      Max 
-2.7361  -0.2750   0.2177   0.5143   1.6897 

Coefficients:
            Estimate Std. Error z value Pr(>|z|) 
(Intercept)  -0.9022     0.3370  -2.677 0.007427 ** 
test.value    2.0806     0.3409   6.103 1.04e-09 ***
cnct_o.L      1.4363     0.4722   3.042 0.002352 ** 
cnct_o.Q      1.2208     0.4934   2.474 0.013359 * 
cnct_o.C     -2.0649     0.5610  -3.681 0.000232 ***
cnct_o^4      0.5599     0.4760   1.176 0.239485 
---
Signif. codes:  0 бе***бж 0.001 бе**бж 0.01 бе*бж 0.05 бе.бж 0.1 бе бж 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 252.32  on 220  degrees of freedom
Residual deviance: 167.68  on 215  degrees of freedom
AIC: 179.68

Number of Fisher Scoring iterations: 6
#####

As I use the model parameters of the standard glm model to calculate the 
odds and probability manualy, I believe cnct_o = 25 to be the reference 
category and that cnct_o.L = 50, cnct_o.Q=100, cnct_o.C = 200 and 
cnct_o^4= 400. But I am not sure of this. The formulas used are :
Odds <-  intercept + slope * test.value + cnct_o  , where cnct_o is the 
corresponding value for the given concentration.
Probability <- inv.logit ( Odds ),  the function inv.logit from package 
бзcarби.
The results of the glm are in the table below, which are first the odds 
and then the probabilityбжs ( inv.logit (odds)).

#####
glm odds        cntc: 
test.value      25              cnct_o.L        cnct_o.Q        cnct_o.C 
cnct_o^4 
0               -0.90218        0.534169        0.318649        -2.96706 
-0.34231
0.5             0.138132        1.574477        1.358957        -1.92675 
0.698001
1               1.17844         2.614785        2.399265        -0.88644 
1.738309
1.5             2.218747        3.655093        3.439572        0.153864 
2.778616
2               3.259055        4.6954          4.47988         1.194172 
3.818924
#####
glm prob.       cntc: 
test.value      25              cnct_o.L        cnct_o.Q        cnct_o.C 
cnct_o^4 
0               0.288604        0.630455        0.578995        0.048936 
0.415249
0.5             0.534478        0.828421        0.79559         0.127111 
0.667744
1               0.764667        0.931807        0.916771        0.291844 
0.850472
1.5             0.90192         0.974793        0.968919        0.53839 
0.941509
2               0.962997        0.990946        0.988792        0.767486 
0.97852
#####

If I compare this with the result of the Predict function in rms, the 
results seem very different, it can be because I misinterpret the glm 
model parameters for the ordered factor. How can I be sure which model 
parameter corresponds to which factor in the standard glm.

#####
Results of lrm:
Predict.lrm     cntc: 
test.value      25              50              100             200 400
0               -0,43815        -3,25628        -1,15323        0,264036 
0,072751
0,50580154      0,614227        -2,2039         -0,10085        1,316414 
1,125129
1,01160308      1,666606        -1,15153        0,951526        2,368793 
2,177508
1,505068        2,693316        -0,12482        1,978236        3,395503 
3,204218
2,01086954      3,745695        0,927563        3,030615        4,447882 
4,256597
#####
Prob. (inv.logit(odds)) cntc: 
test.value      25              50              100             200 400
0               0,392182        0,037102        0,239899        0,565628 
0,51818
0,50580154      0,648904        0,0994          0,474808        0,788585 
0,754939
1,01160308      0,841123        0,24021         0,721422        0,914416 
0,898211
1,505068        0,936631        0,468837        0,878493        0,967564 
0,960993
2,01086954      0,976926        0,716581        0,953938        0,988432 
0,986028
#####

One more problem occurs when I calculate the odds and probability 
manually, using the formulas as with the standard glm I get the following 
(very different) results: 

man.lrm.odds    cntc: 
test.value      25              50              100             200 400
0               2,37998         2,267255        10,11929        22,80909 
45,16286
0,5             3,420288        3,307563        11,1596         23,8494 
46,20316
1               4,460596        4,34787         12,19991        24,8897 
47,24347
1,5             5,500903        5,388178        13,24022        25,93001 
48,28378
2               6,541211        6,428486        14,28053        26,97032 
49,32409
man.lrm.prob 
test.value      25              50              100             200 400
0               0,915288        0,906129        0,99996 1       1
0,5             0,968333        0,964687        0,999986        1       1
1               0,988577        0,987231        0,999995        1       1
1,5             0,995934        0,995451        0,999998        1       1
2               0,998559        0,998388        0,999999        1       1
#####

It is clear to me that I am missing something essential about this lrm 
function.
Question: Are the model parameters for ordered factors converted in some 
way by from the lrm and Glm function?
Д╗      Which could expain why I get a different result for the Predict 
function, then when I calculate it manually. 
Question: Is it correct that the first level of the ordered factor is used 
as reference category by the standard glm function? It there a way of 
giving them a different name, as it does automatically for the nominal 
factor?

Thank you for reading my questions,
And please share any insight that you have.
Kind regards,
Tom 
 


Disclaimer: click here 

	[[alternative HTML version deleted]]

