Answers.Body: <p>Unlike <a href="https://en.wikipedia.org/wiki/Proportional_hazards_model" rel="nofollow noreferrer">wikipedia</a>, the linear predictor (LP) in <code>coxph</code> is made dependent from the baseline hazard function <span class="math-container">$\lambda$</span>. For example, if a model has two predictors <code>x1</code>, <code>x2</code> stratified by <code>sex</code>, the linear predictor would be like:
  <span class="math-container">$$LP = log (  \lambda _{sex}) + \beta _{1}x _{1}+\beta _{2}x _{2}$$</span></p>
  <p>Where the <span class="math-container">$log(\lambda _{sex})$</span> can be obtained by <code>predict()</code> with <code>x1=x2=0</code>.</p>
  <p>We could calculate the LP without strata information (sex) by hand using the <span class="math-container">$\beta _{1} \, \beta _{2}$</span> provided by <code>coef(fit)</code> and a customised <span class="math-container">$\lambda _{unkown \, sex}$</span>.</p>
  <p>Below is a simulation comparing the results from <code>predict()</code> and by hand <code>coef(fit)</code>.</p>
  <pre><code>## Toy data
  df &lt;- list(time=c(4,3,1,1,2,2,3), 
            status=c(1,1,1,0,1,1,0), 
            x1=c(0,2,1,1,1,0,0), 
            x2=1:7,
            sex=c(0,0,0,0,1,1,1)) 
  
  ## Fit a stratified model 
  fit &lt;- coxph(Surv(time, status) ~ x1 + x2 + strata(sex), df) 
  coef(fit)
  #         x1          x2 
  # 0.79451881 -0.01917633
  
  ## Baseline linear predictor
  lambda_sex1 &lt;- predict(fit, newdata=list(x1=0, x2=0, sex=0))
  # -0.746578
  lambda_sex2 &lt;- predict(fit, newdata=list(x1=0, x2=0, sex=1))
  # -0.1497816
  
  ## by predict function
  predict(fit, newdata=list(x1=1, x2=2, sex=0))
  # 0.009588167
  predict(fit, newdata=list(x1=1, x2=2, sex=1))
  # 0.6063845
  
  ## by hand
  lambda_sex1 + coef(fit)[1]*1 + coef(fit)[2]*2
  # 0.009588167
  lambda_sex2 + coef(fit)[1]*1 + coef(fit)[2]*2
  # 0.6063845
  </code></pre>
  <p>Finally, to calculate &quot;a&quot; linear predictor for an unknown <code>sex</code>:</p>
  <pre><code># prevented by predict function
  predict(fit, newdata=list(x1=1, x2=2, sex=2))
  #Error in model.frame.default(data = list(x1 = 1, x2 = 2, sex = 2), formula = ~x1 +  : 
  # factor strata(sex) has new level sex=2
  
  # by hand, assuming lambda for sex 2 is 0
  0 + coef(fit)[1]*1 + coef(fit)[2]*2
  # 0.7561661
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-23T13:12:40.710Z
Answers.LastActivityDate: 2020-07-23T13:38:52.647Z
Answers.LastEditDate: 2020-07-23T13:38:52.647Z
Answers.UserDisplayName: Ryan SY Kwan
Body: <p>I have fitted a cox model on a pooled dataset of multiple studies, say studies A, B, and C. As these studies have different baseline hazard functions, I stratified for 'study' in the model.</p>
  <p>Like this:</p>
  <pre><code>df$SurvObj &lt;- with(df, Surv(event_rd, event == 1))
  
  fit &lt;- coxph(SurvObj ~ cov1 + cov2 + cov3 + cov4 + strata(study),data=df)
  </code></pre>
  <p>Now, I want to assess whether I can use the beta coefficients from the above model to predict the event probability in a new study, study D.
  The baseline hazard in study D is different from those in studies A, B, and C.</p>
  <p>When I have</p>
  <pre><code>predict(fit,type=&quot;lp&quot;,newdata=studyD,reference=&quot;strata&quot;)
  </code></pre>
  <p>I get this error:</p>
  <blockquote>
  <p>Error in model.frame.default(data = studyD, formula = ~hba1c + sbp +
  : factor strata(study) has new levels study=4</p>
  </blockquote>
  <p>Why does R require that I match the strata of study D with those in study A, B, and/or C when I only want to output the linear predictor (<code>type=&quot;lp&quot;</code>)? As the linear predictor is independent of the baseline hazard function in this case.</p>
  <p>I suppose it should be possible to extract the linear predictor for individuals in study D and then manually calculate the event probability using the baseline hazard of study D.</p>
  <p>Does anybody know how to do this?<br />
  I could not find any question like this so maybe I am thinking the wrong way.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-02-04T12:21:55.597Z
Comments.Text: I would recommend trying out Stats.SE (aka CrossValidated) as some folks there might be of more help.
Comments.UserDisplayName: posdef
Tags: <r>
Title: Output the linear predictor from a stratified cox model?
ThreadText: # Output the linear predictor from a stratified cox model?
  
  (Asked by: username_0 on 2019-01-26)
  
  I have fitted a cox model on a pooled dataset of multiple studies, say studies A, B, and C. As these studies have different baseline hazard functions, I stratified for 'study' in the model.
  
  Like this:
  
  ```
  df$SurvObj <- with(df, Surv(event_rd, event == 1))
  
  fit <- coxph(SurvObj ~ cov1 + cov2 + cov3 + cov4 + strata(study),data=df)
  
  ```
  
  Now, I want to assess whether I can use the beta coefficients from the above model to predict the event probability in a new study, study D.
  The baseline hazard in study D is different from those in studies A, B, and C.
  
  When I have
  
  ```
  predict(fit,type="lp",newdata=studyD,reference="strata")
  
  ```
  
  I get this error:
  
  > Error in model.frame.default(data = studyD, formula = ~hba1c + sbp +
  > : factor strata(study) has new levels study=4
  
  Why does R require that I match the strata of study D with those in study A, B, and/or C when I only want to output the linear predictor (`type="lp"`)? As the linear predictor is independent of the baseline hazard function in this case.
  
  I suppose it should be possible to extract the linear predictor for individuals in study D and then manually calculate the event probability using the baseline hazard of study D.
  
  Does anybody know how to do this?  
  I could not find any question like this so maybe I am thinking the wrong way.
  
  **Tags:** r
  
  > **Comments:**
  > **username_1**: I would recommend trying out Stats.SE (aka CrossValidated) as some folks there might be of more help.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2020-07-23. Score: 3
  
  Unlike [wikipedia](https://en.wikipedia.org/wiki/Proportional_hazards_model), the linear predictor (LP) in `coxph` is made dependent from the baseline hazard function $\lambda$. For example, if a model has two predictors `x1`, `x2` stratified by `sex`, the linear predictor would be like:
  $$LP = log ( \lambda \_{sex}) + \beta \_{1}x \_{1}+\beta \_{2}x \_{2}$$
  
  Where the $log(\lambda \_{sex})$ can be obtained by `predict()` with `x1=x2=0`.
  
  We could calculate the LP without strata information (sex) by hand using the $\beta \_{1} \, \beta \_{2}$ provided by `coef(fit)` and a customised $\lambda \_{unkown \, sex}$.
  
  Below is a simulation comparing the results from `predict()` and by hand `coef(fit)`.
  
  ```
  ## Toy data
  df <- list(time=c(4,3,1,1,2,2,3), 
            status=c(1,1,1,0,1,1,0), 
            x1=c(0,2,1,1,1,0,0), 
            x2=1:7,
            sex=c(0,0,0,0,1,1,1)) 
  
  ## Fit a stratified model 
  fit <- coxph(Surv(time, status) ~ x1 + x2 + strata(sex), df) 
  coef(fit)
  #         x1          x2 
  # 0.79451881 -0.01917633
  
  ## Baseline linear predictor
  lambda_sex1 <- predict(fit, newdata=list(x1=0, x2=0, sex=0))
  # -0.746578
  lambda_sex2 <- predict(fit, newdata=list(x1=0, x2=0, sex=1))
  # -0.1497816
  
  ## by predict function
  predict(fit, newdata=list(x1=1, x2=2, sex=0))
  # 0.009588167
  predict(fit, newdata=list(x1=1, x2=2, sex=1))
  # 0.6063845
  
  ## by hand
  lambda_sex1 + coef(fit)[1]*1 + coef(fit)[2]*2
  # 0.009588167
  lambda_sex2 + coef(fit)[1]*1 + coef(fit)[2]*2
  # 0.6063845
  
  ```
  
  Finally, to calculate "a" linear predictor for an unknown `sex`:
  
  ```
  # prevented by predict function
  predict(fit, newdata=list(x1=1, x2=2, sex=2))
  #Error in model.frame.default(data = list(x1 = 1, x2 = 2, sex = 2), formula = ~x1 +  : 
  # factor strata(sex) has new level sex=2
  
  # by hand, assuming lambda for sex 2 is 0
  0 + coef(fit)[1]*1 + coef(fit)[2]*2
  # 0.7561661
  
  ```
  
  
  ---
  
  