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 <i...@aghmed.fsnet.co.uk
<mailto:i...@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
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
