[R] Coxme, package updated on May 2018, Similar case to Re: CoxME: Family relatedness

[Sarahí Rueda Salazar]

Dear Terry,

I can imagine that you are overwhelming with a quite large demand of
request emails, so I hope (always optimistic) that you have time to see
this email.

It is related to the recently updated version of coxme package.  Currently,
I work with multistate models with changes in health status on elder
population.  I use coxPH by stratifying the hazard by different transition
type ( deterioration, health improvements, death from healthy status and
death from not healthy status). I use mstate and Biograph packages to
reshape my data disentangling different events by transition type (1:4)

So, I wondered if I can apply mixed effect on my multistate model to see
the differences by countries (shared frailty).  I´ve read the material you
post on May 11, 2018 and other several materials in this subject :Steele et
al:2004, Austin : 2017,Putter: 2014, Willekens: 2014, Brostön:2011, Mills,
2011 and finally Allison which provide me a quite understandable approach
in the line of my humble knowledge in mathematical demography.

I thought it is not possible to use to multilevel in
multistate with coxme package. By stratifying the hazard  (reference risk
by interested covariates )  by my transitions type, I make hazard be
free(no proportional by transition type) and it means that the model
estimates separate baseline hazard for the different values of transition
type ( following this material by Putter( 2018:7
<https://cran.r-project.org/web/packages/mstate/vignettes/Tutorial.pdf>)

It is a simple model using sex , reference category "male" with my data:

> modelSex.0 <- coxph(Surv(Tstarta,Tstopa,status) ~
+                        SexF.1+ SexF.2+SexF.3+SexF.4+
+                       strata(trans),
+                     data=d0)

SexF.1= related to covariate sex  in transition 1(healthy to not healthy)
SexF.2 = transition healthy to death
SexF.3= transition  Not healthy to Healthy
SexF.4= transition  Not healthy to Death

But, I came across  this post
<https://stat.ethz.ch/pipermail/r-help/2014-September/421690.html>  where
you replied a request ( I copy a chunk that I´m interested in):

2. The model above is the correct covariance structure for a set of
> families. There is a
> single intercept per subject, with a complex correlation matrix. The
> simpler "per family"
> frailty model would be



> model4 <- coxme(Surv(Survival, Event) ~ Sex + strata(cohort) + SNP1 + SNP2
> + SNP3 +
> (1|famid), death.dat)



> This model lets each family have a separate risk, with everyone in the
> same family sharing
> the exact same risk. It is less general than model3 above which lets a
> family have higher
> risk plus has variation between family members. A model with both
> per-subject and per family terms is identical to one with a covariance
> matrix of s1 K + s2 B, where K is the kinship matrix, B is a block
> diagonal matrix which
> has a solid block of "1" for each family, and s1 s2 are the fitted
> variance coefficients.


 So, my strata would be my transition type (trans) instead "cohort"  and my
groups (instead famid) would be country, as follows


> modelSex.1 <- coxme(Surv(Tstarta,Tstopa,status) ~
+                        SexF.1+ SexF.2+SexF.3+SexF.4+
+                       (1|Country),
+                     data=d0)

> anova(modelSex.0,modelSex.1)

Analysis of Deviance Table
 Cox model: response is  Surv(Tstarta, Tstopa, status)
 Model 1: ~ SexF.1+ SexF.2+SexF.3+SexF.4+ strata(trans)
 Model 2: ~ SexF.1+ SexF.2+SexF.3+SexF.4+ strata(trans) + (1 | Country)
    loglik  Chisq Df P(>|Chi|)
1 -1342082
2 -1339605 4953.8  1 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> stem(exp(ranef(modelSex.1)[[1]]))

  The decimal point is 1 digit(s) to the left of the |

   6 | 589
   8 | 1448477
  10 | 115669
  12 | 0846

> exp(ranef(modelSex.1)[[1]])
       AT        BE        BG        CY        CZ        EE
1.1981346 0.9712875 0.8092424 1.3614033 0.8382588 1.0117071
       EL        ES        HU        IE        IT        LT
0.7514988 1.2761863 1.0085404 0.8760514 1.1615717 0.9708070
       LU        LV        MT        PL        PT        RO
1.1893096 1.3381012 1.0522161 0.7843473 1.1642434 0.7911292
       SK        UK
0.9440166 0.8428222

> fixed.effects(modelSex.1)
     SexF.1      SexF.2      SexF.3      SexF.4
 0.16349517 -0.63184370 -0.08578787 -0.61023260


The thing is that I´m not sure on how to interpret the frailty values by
countries (random effect described by the variance within groups)  because
I have four different effects (each related with my transition types). I
know that by using strata for my transition type is the same as I applied 4
different cox model (related for each type of transition). If I
apply separated model  I would obtain frailty random effect by groups
(countries) regarding specific transition but, doing this model with the
strata  ( modelSex.1 ) with my 4 transition type at once, I do not know
what those frailty values are telling me about the different type of hazard.

For other side, I have other question related to the recent article of
Mixed Effect (May, 2018) . Might be it is a very very silly question but I
need to understand .
In pag 9, second parraph , you describe the simple cox model : when you
describe the standard deviation (excess of risk for each group)  , you
state that "... 15% of the families to be 1 std dev or more above the
mean..."  Im working with that data and code but I couldnt find where you
got the value of 15%.


I hope this makes any sense to you,

Best wishes

 Sarahí


-- 
*Sarahí Rueda Salazar*
*Investigadora en Formación (FPI/CED)*












*Centre d'Estudis DemogràficsCarrer de Ca n'Altayó, Edifici E2Universitat
Autònoma de Barcelona,08193
Bellaterra, Barcelona/SPAINPhone: 34/93.581.30.60
<34%2F93.581.30.60>Fax: 34/93.581.30.61
<34%2F93.581.30.61>e-mail: srueda using ced.uab.es
<ssancho using ced.uab.es>http://www.ced.uab.es <http://www.ced.uab.es/>*