14 replies to this topic

[sinyinlim]

Hi,

 

This is an example from Phoenix NLME 1.3 for three category multinomial model with a covariate (pg231): 

 

Prob(y=0) = ilogit(th1+ th4*PER + th5*DOSE)
Prob(y=1) = ilogit(th1+ th4*PER + th5*DOSE + th2) – Prob(category=0)
Prob(y=2) =1-(Prob(y=1) + Prob(y=0) )

 

PER and DOSE are covariates

 

how should I write the "multi" statement? 

multi(EObs, ilogit, ...)?

 

Thanks.

[serge guzy]

Dear colleague

The best way to do it is to start the other way with first modeling the probability to be in the highest level , 2 here

 

You write 

prob(2)=ilogit   (th1+ th4*PER + th5*DOSE)

or logit(prob(2))=(th1+ th4*PER + th5*DOSE)

​prob(2 or 1)=ilogit(th1+ th4*PER + th5*DOSE+th2))   th2 must be positive

 

​logit(prob(2 or 1)=(th1+ th4*PER + th5*DOSE+th2))

 

​now that you start witgh the highest level first you can easily write the multi statement

 

​multi(DV,(th1+ th4*PER + th5*DOSE+th2),th1+ th4*PER + th5*DOSE)

​do not forget to constraint th2 to be positive

 

​Note that I changed the order of the probability that it changes the meaning of the parameters.

 

​If you want to have the same meaning as in your original model, just change the levels

 

​2                to    0

1               to      1

​0                to     2

 

​Then now you use the same model as you had and write

multi(DV,(th1+ th4*PER + th5*DOSE+th2),th1+ th4*PER + th5*DOSE)​    th2 positive

 

Try it and if you still have problems, send me the project for me to review.

best Regards

Serge

[smouksassi1]

please note that there is several ways to code the multiple categories models depending on software and statement. 

we can do it in cumulative logits that can be done in ascending or descending ways. I have a project that compares R outputs with the multi and ordered phoenix statement event then the meaning of the "intercept" can be different depending on how you define the logits


As Serge pointed out you need to make sure that the probabilities are ordered appropriately for the model to work.

 

can you share a test data that I can use to illustrate on it if not will use a generic R free dataset?

 

[sinyinlim]

Thanks for the response. In my project, my PD response are in 6 ordered categories (2,1,0,-1,-2,-3), 2 is the lowest (correlates with low concentration) and -3 is the highest (correlates with high concentration). they are sedation scores, negative values represent sedated and calm while positive values are agitated and awake.

 

I have a PK model and I fixed the parameters, now I am adding the PD data. I attempted to write the code for the logistic model as below: 

 

stparm(th1=tvth1+nth1)

stparm(th2 = tvth2)

stparm(th3 = tvth3)

stparm(th4 = tvth4)

stparm(th5 = tvth5)

stparm(th6 = tvth6)

fixef(tvth1= c(, 1, )) 

fixef(tvth2 = c(0, 1, ))

fixef(tvth3=c(0,1,))

fixef(tvth4 = c(0,1,))

fixef(tvth5 = c(0,1,))

fixef(tvth6 = c(0,1,))

 

#ilogit

Prob(-3) = ilogit(th1+th2*C)

Prob(-3 or -2) = ilogit(th1+th2*C+th3)

Prob(-3 or -2 or -1) = ilogit(th1+th2*C+th3+th4)

Prob(-3 or -2 or -1 or 0) = ilogit(th1+th2*C+th3+th4+th5)

Prob(-3 or -2 or -1 or 0 or 1) = ilogit(th1+th2*C+th3+th4+th5+th6)

 

multi(EObs,ilogit,(th1+th2*C+th3+th4+th5+th6),(th1+th2*C+th3+th4+th5),(th1+th2*C+th3+th4),(th1+th2*C+th3),(th1+th2*C))

 

model wouldn't run, i'm not sure if my code is correct, could you help? 

 

thank you!

[smouksassi1]

Dear Sinyinlim,
It is hard to help without a fully reproducible example and a complete code that I can run on my end. You can always simualte fake data that illustrate the problem and share it.

The nlme user guide has examples using the multi statement and the software comes with the Multinomial_3cat_with_covariates.phxproj and the phoenix modeling language guide has a section regarding:
Multi statement for Multinomial models 
Check also the meaning and differences between ordinal versus multi statement and how to code it your own way using LL.

The built in multi / ordered statement start by category 0  up to highest category e.g. 0,1,2,3 I am not sure 
it is reading your EObs correctly try changing your categories to something like 0,1,2,3,4,5 with the order that makes sense. 

or try to code the LL yourself by testing for the categories you have in your data 

 

Attached is how the multi/ordinal/LL fit the same model and compares with standard R functions like :
MASS polr, ordinal clm, or rms lrm.

Bests,
 
Samer

 

Attached Files

•  ordinalmultitutorial.phxproj   1011.5KB   345 downloads

Edited by smouksassi1, 15 January 2018 - 05:23 PM.

[sinyinlim]

Hi Samer,

 

thank you for your help. I meant to attach my project in my last post. I have made changes based on your suggestions and according to the PML guide. 

 

I was trying to look at your attachment and learn about the multi/ordinal/LL fitting, however, it's not downloading correctly, could you re-upload the file?

 

Thanks!

 

SYL

Attached Files

•  Logistic model PD 1_11_18.phxproj   1.14MB   385 downloads

[Simon Davis]

SYL what error message are you getting when you try to download? I just downloaded it without issue in Firefox, what browser are you using?

 

If you are getting a message that you need a password, be aware you do not.  it is not a Zip file, rename the extension to .phxproj and it will open no trouble.

 

Internet explorer (and perhaps other browsers) sometimes do this to PHX files when attached to the forum, I guess since it's a compressed an encrypted format.

 

 Simon.

  Simon.

[sinyinlim]

SYL what error message are you getting when you try to download? I just downloaded it without issue in Firefox, what browser are you using?

 

If you are getting a message that you need a password, be aware you do not.  it is not a Zip file, rename the extension to .phxproj and it will open no trouble.

 

Internet explorer (and perhaps other browsers) sometimes do this to PHX files when attached to the forum, I guess since it's a compressed an encrypted format.

 

 Simon.

  Simon.

Hello Simon,

 

Thanks. sorry it's actually an error from loading the project. It says "There was an error loading the project file: C:\Users\slim2\Downloads\ordinalmultitutorial.phxproj"

 

Not sure why..

[smouksassi1]

do you have phoenix 8.0 ?

[sinyinlim]

do you have phoenix 8.0 ?

I have Phoenix 64 (Build 6.4.0.768) WinNonlin 6.4

[smouksassi1]

Hi sinyinlim,

Your version is much older than the current released 8.0 and that explains why you can't open it.

I would try coding you model like below

 

 

stparm(th1=tvth1+nth1)

stparm(th2 = tvth2)

stparm(th3 = tvth3)

stparm(th4 = tvth4)

stparm(th5 = tvth5)

stparm(th6 = tvth6)

 

fixef(tvth1 = c( , 1, )) 

fixef(tvth2 = c(0, 1, ))

fixef(tvth3 = c(0, 5,))

fixef(tvth4 = c(0,10,))

fixef(tvth5 = c(0,15,))

fixef(tvth6 = c(0,20,))

 

conceffects= th2*C

 

#ilogit  Y>=2 ( =1) , Y>=1, Y>=0,  Y>=-1,  Y>=-2,  Y>=-3 

PGE1 =ilogit(conceffects - th1)

PGE0 =ilogit(conceffects - th3)

PGEn1=ilogit(conceffects - th4)

PGEn2=ilogit(conceffects - th5)

PGEn3=ilogit(conceffects - th6)

 

#Probabilities for Y=2, Y=1, Y=0, Y=-1, Y=-2 

P2  =    1  - PGE1 

P1  = PGE1  - PGE0

P0  = PGE0  - PGEn1

Pn1 = PGEn1 - PGEn2

Pn2 = PGEn2 - PGEn3

Pn3 = PGEn3 - 0

 

LL(EObs,EObs==-3 ? Pn3:

EObs==-2 ? Pn2:

EObs==-1 ? Pn1:

EObs== 0 ? P0 :

EObs== 1 ? P1 :

           P2)

[sinyinlim]

 

Hi sinyinlim,

Your version is much older than the current released 8.0 and that explains why you can't open it.

I would try coding you model like below

 

 

stparm(th1=tvth1+nth1)

stparm(th2 = tvth2)

stparm(th3 = tvth3)

stparm(th4 = tvth4)

stparm(th5 = tvth5)

stparm(th6 = tvth6)

 

fixef(tvth1 = c( , 1, )) 

fixef(tvth2 = c(0, 1, ))

fixef(tvth3 = c(0, 5,))

fixef(tvth4 = c(0,10,))

fixef(tvth5 = c(0,15,))

fixef(tvth6 = c(0,20,))

 

conceffects= th2*C

 

#ilogit  Y>=2 ( =1) , Y>=1, Y>=0,  Y>=-1,  Y>=-2,  Y>=-3 

PGE1 =ilogit(conceffects - th1)

PGE0 =ilogit(conceffects - th3)

PGEn1=ilogit(conceffects - th4)

PGEn2=ilogit(conceffects - th5)

PGEn3=ilogit(conceffects - th6)

 

#Probabilities for Y=2, Y=1, Y=0, Y=-1, Y=-2 

P2  =    1  - PGE1 

P1  = PGE1  - PGE0

P0  = PGE0  - PGEn1

Pn1 = PGEn1 - PGEn2

Pn2 = PGEn2 - PGEn3

Pn3 = PGEn3 - 0

 

LL(EObs,EObs==-3 ? Pn3:

EObs==-2 ? Pn2:

EObs==-1 ? Pn1:

EObs== 0 ? P0 :

EObs== 1 ? P1 :

           P2)

 

Thanks for your help!

[sinyinlim]

 

Hi sinyinlim,

Your version is much older than the current released 8.0 and that explains why you can't open it.

I would try coding you model like below

 

 

stparm(th1=tvth1+nth1)

stparm(th2 = tvth2)

stparm(th3 = tvth3)

stparm(th4 = tvth4)

stparm(th5 = tvth5)

stparm(th6 = tvth6)

 

fixef(tvth1 = c( , 1, )) 

fixef(tvth2 = c(0, 1, ))

fixef(tvth3 = c(0, 5,))

fixef(tvth4 = c(0,10,))

fixef(tvth5 = c(0,15,))

fixef(tvth6 = c(0,20,))

 

conceffects= th2*C

 

#ilogit  Y>=2 ( =1) , Y>=1, Y>=0,  Y>=-1,  Y>=-2,  Y>=-3 

PGE1 =ilogit(conceffects - th1)

PGE0 =ilogit(conceffects - th3)

PGEn1=ilogit(conceffects - th4)

PGEn2=ilogit(conceffects - th5)

PGEn3=ilogit(conceffects - th6)

 

#Probabilities for Y=2, Y=1, Y=0, Y=-1, Y=-2 

P2  =    1  - PGE1 

P1  = PGE1  - PGE0

P0  = PGE0  - PGEn1

Pn1 = PGEn1 - PGEn2

Pn2 = PGEn2 - PGEn3

Pn3 = PGEn3 - 0

 

LL(EObs,EObs==-3 ? Pn3:

EObs==-2 ? Pn2:

EObs==-1 ? Pn1:

EObs== 0 ? P0 :

EObs== 1 ? P1 :

           P2)

 

Hello,

 

In this case, concentration © from the PK model is a covariate of the response. is this model connecting the observed response (EObs) to the corresponding predicted concentration (from the PK model, the predicted concentration at the time of the response recorded)? 

 

thanks.

[smouksassi1]

Note that I only included the ordinal logistic part.
It depends on what is the meaning of C and how you compute it.

Typically it is the individual predicted concentration using your PK model and then it is the prediction of PK across time at any time point you have a PD observation at.

Bests,

Samer

[sinyinlim]

Note that I only included the ordinal logistic part.
It depends on what is the meaning of C and how you compute it.

Typically it is the individual predicted concentration using your PK model and then it is the prediction of PK across time at any time point you have a PD observation at.

Bests,

Samer

I see. Thank you very much!