[R] Getting HR from Cox model in R

Michael Dewey info at aghmed.fsnet.co.uk
Mon Dec 29 16:31:15 CET 2014


Comments in line

On 24/12/2014 22:09, Ruzan Udumyan wrote:
> Dear Michael,
>
> Thank you very much for your reply. The more complete information is as
> follows:
>
> I want to do a mediation analysis following the below-mentioned syntax
> from:
> http://www.biomedcentral.com/content/supplementary/1471-2288-14-9-s1.pdf
>
> I did not define categorical variables as logical variables. I modelled
> them as /*factor.X, factor.Xstar*/, etc. the X variable has 3 levels.
>
> It is all sorted but I am not sure about the last bit: to return the
> values. For example, I could not figure out what unname stands for, and
> whether it is correct to use when variables are modelled as factor.X.
>
> I wrote the syntax as:
>
>   TE2 = exp(sum(coef(cox)[c('factor(X)2', 'factor(Xstar)2')]))   # level
> 2 vs level 1(ref)
>   TE3 = exp(sum(coef(cox)[c('factor(X)3', 'factor(Xstar)3')]))   # level
> 3 vs level 1(ref)
>
>    DE2 = exp(unname(coef(cox)['factor(X)2']))
>    DE3 = exp(unname(coef(cox)['factor(X)3']))
>
>    IE2 = exp(sum(coef(cox)['factor(Xstar)2']))
>    IE3 = exp(sum(coef(cox)['factor(Xstar)3']))
>    PM2 = log(IE2) / log(TE2)
>    PM3 = log(IE3) / log(TE3)
>
>
> Thank you very much for your help.
>
> Wishing you happy holidays,
> Ruzan
>
>
> *The script from the link:*
> doEffectDecomp = function(d)
> {
>   # Step 1: Replicate exposure variable, predict mediator
>   d$TrialTemp = d$Trial
>   MOpti = glm(Opti ~ TrialTemp + Age5 + ECOG + Ascit + Comorb + Histo +
>   Grade, family=binomial(), data=d)
> # Step 2: Replicate data with different exposures for the mediator
>   d1 = d2 = d
>   d1$Med = d1$Trial
>   d2$Med = !d2$Trial
>   newd = rbind(d1, d2)
> # Step 3: Compute weights for the mediator
>   newd$TrialTemp = newd$Trial
>   w = predict(MOpti, newdata=newd, type='response')
>   direct = ifelse(newd$Opti, w, 1-w)
>   newd$TrialTemp = newd$Med
>   w = predict(MOpti, newdata=newd, type='response')
>   indirect = ifelse(newd$Opti, w, 1-w)
>   newd$W = indirect/direct
> # Step 4: Weighted Cox Model
>   cox = coxph(Surv(OS, Status) ~ Trial + Med + Age5 + ECOG + Ascit +
>   Comorb + Histo + Grade, weight=W, data=newd)
> # Return value: Estimates for total, direct, indirect effect
>   TE = exp(sum(coef(cox)[c('TrialTRUE', 'MedTRUE')]))
>   DE = exp(unname(coef(cox)['TrialTRUE']))
>   IE = exp(sum(coef(cox)['MedTRUE']))
>   PM = log(IE) / log(TE)
>   return(c(exp(coef(cox)), TE=TE, DE=DE, IE=IE, PM=PM))
> }
>
> On Tue, Dec 23, 2014 at 6:21 PM, Michael Dewey <info at aghmed.fsnet.co.uk
> <mailto:info at aghmed.fsnet.co.uk>> wrote:
>
>     Inline comments
>
>     On 23/12/2014 09:42, Ruzan Udumyan wrote:
>
>         Dear All,
>
>         I am not familiar with R language well. Could you please help me
>         interpret
>         these commands?:
>
>
>            TE = exp(sum(coef(cox)[c('aTRUE', 'bTRUE')]))   - does it
>         mean exp(coef(a
>         variable) + coef(b variable)) ?
>
>
>     You have not given us much to go on here.
>     I assume if you go
>     coef(cox)
>     you will find elements labelled aTRUE and bTRUE which implies the
>     existence of a logical covariate with values TRUE and FALSE.

You now tell me my assumption was wrong. Presumably you know what you 
are trying to do but we do not and you are not really helping us by 
giving us a load of code to read through with any details of the dataset.

  The
>     author of the code is trying to do what you suggest.
>
>            DE = exp(unname(coef(cox)['aTRUE'])__)  - what is unname for ?
>
>
>     ?unname
>
>         Thank you very much beforehand for your help.
>
>         Wishing you happy holidays,
>         Ruzan
>
>                  [[alternative HTML version deleted]]
>          > PLEASE do read the posting guide
>         http://www.R-project.org/__posting-guide.html
>         <http://www.R-project.org/posting-guide.html>
>         and provide commented, minimal, self-contained, reproducible code.
>
>
>     If you post again please do read the message above.
>

Commented, minimal, self-contained, reproducible code was asked for.

>
>
>         -----
>         No virus found in this message.
>         Checked by AVG - www.avg.com <http://www.avg.com>
>         Version: 2015.0.5577 / Virus Database: 4257/8792 - Release Date:
>         12/23/14
>
>
>
>     --
>     Michael
>     http://www.dewey.myzen.co.uk
>
>
> No virus found in this message.
> Checked by AVG - www.avg.com <http://www.avg.com>
> Version: 2015.0.5577 / Virus Database: 4257/8833 - Release Date: 12/29/14
>

-- 
Michael
http://www.dewey.myzen.co.uk



More information about the R-help mailing list