[R] Problem of data format using function of tmerger()

2019-03-20 Thread wong jane
We want to perform a survival analysis using time-dependent covariates in
the Cox
regression. In this analysis, ESRD1_END and ESRD1_TIME are the outcome and
follow-up time,respectively. AKI_END is a binary time-dependent variable,
while
AKI_TIME is the time point for AKI_END measurement. The Cox model is some-
thing like: Surv(tstart, tstop, ESRD1_END) ˜ AKI_END.
However, it seems the formatted data were incorrect using the functions of
event()
and tdc(), mainly including two problems:
1, The output of event(y, x) is inconsistent with the input data x.
2, Missing value NA was generated from tdc(y, x).


> library(survival)
> ori_data <- read.table("test_data.csv", sep=",", header=TRUE,
stringsAsFactors=F)
> str(ori_data)
' data.frame ' : 7 obs. of 5 variables:
$ data95_Obs: int 10112 15 11 12 13 14 16
$ ESRD1_TIME: num 6.27 9.35 16.77 10.84 18.33 ...
$ ESRD1_END : int 1 0 0 0 0 1 1
$ AKI_END : int 1 1 0 0 0 0 1
$ AKI_TIME : num 3.23 7.65 16.77 10.84 18.33 ...
> dim(ori_data)
[1] 7 5
> head(ori_data)
data95_Obs ESRD1_TIME ESRD1_END AKI_END AKI_TIME
 10112 6.267 1 1 3.233
 15 9.347 0 1 7.654
 11 16.775 0 0 16.775
 12 10.836 0 0 10.836
 13 18.330 0 0 18.330
 14 6.971 1 0 6.971

> #generate the variables "tstart", "tstop" and "status". However, the
"status" is different
> #from "ESRD1_END"
> ori_data <- tmerge(ori_data, ori_data, id=data95_Obs,
status=event(ESRD1_TIME, ESRD1_END))
> dim(ori_data)
[1] 7 9
> head(ori_data)
data95_Obs ESRD1_TIME ESRD1_END AKI_END AKI_TIME id tstart tstop status
 10112 6.267 1 1 3.233 10112 0 6.267 1
 15 9.347 0 1 7.654 15 0 9.347 0
 11 16.775 0 0 16.775 11 0 16.775 1
 12 10.836 0 0 10.836 12 0 10.836 0
 13 18.330 0 0 18.330 13 0 18.330 0
 14 6.971 1 0 6.971 14 0 6.971 0

> #generate the time-dependent variable "AKI_TC". However, why there are
some NAs.
> ori_data <- tmerge(ori_data, ori_data, id=data95_Obs,
AKI_TC=tdc(AKI_TIME, AKI_END))
> dim(ori_data)
[1] 9 10
> head(ori_data)
data95_Obs ESRD1_TIME ESRD1_END AKI_END AKI_TIME id tstart tstop status
 10112 6.267 1 1 3.233 10112 0.000 3.233 0
 10112 6.267 1 1 3.233 10112 3.233 6.267 1
 15 9.347 0 1 7.654 15 0.000 7.654 0
 15 9.347 0 1 7.654 15 7.654 9.347 0
 11 16.775 0 0 16.775 11 0.000 16.775 1
 12 10.836 0 0 10.836 12 0.000 10.836 0
AKI_TC
 NA
 1
 NA
 0
 NA
 NA

[[alternative HTML version deleted]]

__
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.


[R] Fitting problem for Cox model with Strata as interaction term

2018-05-08 Thread wong jane
Dear All,


I got a warning message "X matrix deemed to be singular" in Cox model with
a time dependent coefficient. In my analysis, the variable "SEX" is a
categorical variable which violate the PH assumption in Cox. I first used
the survSplit() function to break the data set into different time
intervals, and then fit the model. The procedures can be described as
follows:

> data <- survSplit(Surv(time, status)~., data=data, cut=c(5, 10), episode=
"tgroup", id="id")#data split
> data$SEX <- as.factor(data$SEX)
> data$SEX <- relevel(data$SEX, ref="1")
> str(data$SEX)
 Factor w/ 2 levels "1","2": 1 1 1 2 2 2 2 2 2 2 ...
> xtabs(~status+SEX, data)
  SEX
status12
 0 6764 5537
 1 1643 1468
> xtabs(~status+tgroup+SEX, data)
, , SEX = 1

  tgroup
status123
 0 3035 2403 1326
 1  759  458  426

, , SEX = 2

  tgroup
status123
 0 2596 1940 1001
 1  701  476  291


> fit1 <- coxph(Surv(tstart, time, status) ~ *SEX:strata(tgroup)*,
data=data)
Warning message:
In coxph(Surv(tstart, time, status) ~ SEX:strata(tgroup), data = data) :
  X matrix deemed to be singular; variable 2 4 6
> fit1
Call:
coxph(formula = Surv(tstart, time, status) ~ SEX:strata(tgroup),
data = data)

   coef exp(coef) se(coef) z   p
SEX1:strata(tgroup)tgroup=1 -0.07150.9310   0.0524 -1.360.17
SEX2:strata(tgroup)tgroup=1  NANA   0.NA  NA
SEX1:strata(tgroup)tgroup=2 -0.26240.7692   0.0655 -4.01 6.1e-05
SEX2:strata(tgroup)tgroup=2  NANA   0.NA  NA
SEX1:strata(tgroup)tgroup=3 -0.03090.9696   0.0761 -0.410.68
SEX2:strata(tgroup)tgroup=3  NANA   0.NA  NA

Likelihood ratio test=18.1  on 3 df, p=0.000425
n= 15412, number of events= 3111
>

#It works well if the main effect was included
> fit2 <- coxph(Surv(tstart, time, status) ~* SEX*strata(tgroup)*,
data=data)
> fit2
Call:
coxph(formula = Surv(tstart, time, status) ~ SEX * strata(tgroup),
data = data)

   coef exp(coef) se(coef) z p
SEX2 0.07151.0741   0.0524  1.36 0.172
SEX2:strata(tgroup)tgroup=2  0.19091.2104   0.0838  2.28 0.023
SEX2:strata(tgroup)tgroup=3 -0.04060.9602   0.0924 -0.44 0.660

Likelihood ratio test=18.1  on 3 df, p=0.000425
n= 15412, number of events= 3111
>

Much appreciated for any suggestions. Many thanks in advance.

[[alternative HTML version deleted]]

__
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.


[R] How to pass a variable to a function which use variable name as a parameter

2015-05-26 Thread wong jane
There are functions which use variable names as parameters in some R
packages. However, if the variable name is stored in another variable, how
can I pass this variable to the function. Taking the rms package as an
example:

library(rms)
n - 1000
age - rnorm(n, 50, 10)
sex - factor(sample(c('female','male'), n,TRUE))

y - rnorm(n, 200, 25)
ddist - datadist(age, sex)
options(datadist='ddist')
fit - lrm(y ~ age)
Predict(fit, age, np=4)
options(datadist=NULL)

Here age was a variable name passed to Predict() function, but if age
was stored in variable var, that is, var - age, how can I pass var
to Predict() function? The purpose is that I want to change the parameter
of Predict()  in a loop.

[[alternative HTML version deleted]]

__
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.


[R] Cox regression model for matched data with replacement

2014-06-23 Thread wong jane
My problem was how to build a Cox model for the matched data (1:n) with
replacement. Usually, we can use stratified Cox regression model when the
data were matched without replacement. However, if the data were matched
with replacement, due to the re-use of subjects, we should give a weight
for each pair, then how to incorporate this weight into a Cox model. I also
checked the clogit function in survival package, it seems suitable to the
logistic model for the matched data with replacement, rather than Cox
model. Because it sets the time to a constant. Anyone can give me some
suggestions?

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
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.


[R] Restricted Cubic Spline using rms package

2014-02-25 Thread wong jane
I am trying to use restricted cubic spline to examine the association
between an independent variable and outcome. When I plot the relationship
between HR and the independent variable, I found the HR is not equal to 1
for the specified reference value. For example:

dat - read.table(mydata.csv)
ddist - datadist(dat[HBA1C])
ddist$limits$HBA1C[2] -  6.5;  # specify the reference
value for HBA1C
options(datadist='ddist');

fit - cph(Surv(ftime,outcome)~rcs(HBA1C,3), dat);   # fit the cox model
using rcs with 3 knots
pred_HR - Predict(fit, HBA1C, fun=exp);  # obtain the HR
plot(pred_HR,ylab=HR,abline=list(h=1, v=6.5));

The figure presents that the HR is not equal to 1 for the reference HBA1C
value (say, 6.5). I am not quite familiar with rcs, can somebody help me
explain it? If I want to obtain the HR=1 for the specified value, how
should I do? Thanks.

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
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.