#StackBounty: #r #machine-learning #mathematical-statistics #survival Accelerated Failure Time Regression Performance (Survival Analysis)

Bounty: 100

Issue is around high intercept with AFT Regression. Let me explain below:

Suppose you are modelling the time to an event via an Accelerated Failure Time Regression i.e. given survival time $T$, suppose we have observed values of covariates $x_{i1}, …, x_{ip}$ and possibly censored survival time $t_i$, then:
$$ log(t_i) = beta_0 + beta_1 x_{i1} + … + beta_p x_{ip}+ sigma epsilon_i $$

Suppose we are looking at a Weibull AFT i.e. where $epsilon_i $ are IID according to a Gumbel Distribution (Extreme Value Type 1).

You are looking at the case of time varying covariates (assume just one for now) e.g. you have a dataset like the following example with a single time dependent covariate (TDC_1). Where Start is the enter time (period start) and End is the period end (exit time) and UNIT_ID is the ID for the entity in the study:

START END EVENT UNIT_ID TDC_1 
0     1   0     1       0.1 
1     2   0     1       0.2
2     3   0     1       0.3
...
19    20  1     1       1.9
0     1   0     2       0.1
1     2   0     2       0.2
2     3   0     2       0.3
...
19    20  1     2       1.9

With the aftreg function from the eha library in R you can construct a Weibull AFT e.g.

model <- aftreg(Surv(START, END, EVENT) ~ TDC_1, dist="weibull", data=df, id=UNIT_ID, param='lifeExp')

Calling model.coefficients gives:

             model.coefficients
TDC_1        -0.905
log(scale)    9.393
log(shape)    0.046

The expected time to event when $T$ follows a Weibull distrubtion is given by:
$$E(T|X_i) = exp left( beta_0 + x_i beta_1 right)Gamma(1 + sigma) = exp left( 9.393 – 0.905*TDC_1 right)*0.98 $$

As $beta_0 = log(scale)$ and $sigma = frac{1}{exp(log(shape))}$

My question is around these parameter estimates (in particular the the intercept term ($beta_0 = log(scale)$). No matter how I change the error term parameterisation e.g. if $epsilon_i$ are distributed normally (then $T$ lognormal) or if $epsilon_i$ ~ Logistic etc, the intercept is exceptionally high and appears not to be optimal in terms of minimising error on time to event.

For example if I manually subtract 2 from the intercept (9.393 – 2) I can reduce the root mean squared error on the time to event on the dataset fit:.

Intercept TIME_TO_EVENT_RMSE
9.393     776 days
7.393     97 days

Here TIME_TO_EVENT_RMSE is calculated as (with a dataset that only contains non-censored events):

$$ RMSE = sqrt{sum_{i}^{n} frac{(exp left( beta_0 + x_i beta_1 right)Gamma(1 + sigma) – t_i)^2}{n}} $$

For further illustration, suppose you model directly using exponential regression (i.e. linear regression and logging the target variable) with exactly the same dataset (only using non-censored events so the two are comparable). I know they are minimising different loss functions and aren’t directly comparable, but just for illustration purposes:

TIME_TO_EVENT UNIT_ID TDC_1 
19             1      0.1 
18             1      0.2
17             1      0.3
...

Here we have:

$$E(T|X_i) = exp left( beta_0 + x_i beta_1 right) = exp left( 8.03 – 0.5*x_i right) $$

I know that AFT Regression is not directly minimising RMSE, and that with the AFT regression the TDC_1 coefficient magnitude is larger in addition to a larger intercept, however with the intercept as high as it is, the model isn’t particularly useful (significantly over-predicting the time to event).

Questions:

  1. Has anyone experienced this before and have any advice on how to improve the AFT model?
  2. Is there anyway to fix the scale with time varying covariates in AFTRegression?


Get this bounty!!!

#StackBounty: #survival #hierarchical-bayesian #stan Survival analysis for different diseases on same patients

Bounty: 100

I want to apply survival analysis on UFC-fights. Each fighter represents a "disease" and each knock-out is a "death". Each UFC fight consists of a number of rounds and the number of rounds corresponds to the time alive. There are basically three different possible outcomes per fight:

  • Fighter KO’s opponent in round r – event is observed in round r
  • Opponent KO’s fighter in round r – event is right-censored in round r
  • Fighter or opponent wins by jury decision – event is right-censored in round r = 3 ( normal fight) or r = 5 ( title fight)

This looks for a single fighter as follows:
enter image description here

The goal is to determine the conditional probability of opponent o getting KO’ed when fighting against fighter f in round r give that he/she has survived until that round (a.k.a. the hazard rate):

enter image description here

The person period data looks as follows:

enter image description here

I assume that because fighters are fighting against each other within their weight classes, there is a lot of information that can borrowed between the opponents. Therefore, I intend to use (bayesian) discrete survival analysis. Singer and Willet (1993) show the likelihood function of this is the same as logistic regression / likelihood function for N independent Bernoulli trials with parameter λ (the hazard rate).

Is anyone aware of similar papers where the methodology consist of different diseases and more or less same patient pool? There are no competing risks between the fighters / diseases because there is enough rest between the fights (sometimes months). Although, I understand that it is difficult to die more than once, can please someone point me out how to determine the "deadliest" fighter?


Get this bounty!!!

#StackBounty: #mixed-model #multiple-regression #survival #panel-data #time-varying-covariate Include step function into Joint longitud…

Bounty: 50

I am fitting a joint longitudinal and time to event model on production data with the aim of making dynamic predictions of the time of assembly of a machine. I am using JMbayes R package.

Among the time-dependent variables in the longitudinal part of the model I have a dummy variable that tells whether the mechanical part of the assembly is finished or not, this variable is 0 up to a certain time point, thereafter it becomes one from the time the mechanical part is completed to the time of the event (the assembly has been completed). So basically it is a one step function.

Now I am fitting a binomial(link="logit") longitudinal model for this variable using time fixed effect and a random effect for the order number related to the machine being assembled.
I am linking this longitudinal model to the time event component using the "current value" association.

The idea is that if I complete the mechanical part of the assembly ahead of the historical average, it is likely that the total time of assembly will be lower, coversely if I spend much more time, it is likely that it will be higher.

I am sure this is not the right way to include this kind of variable in the model, both in modelling it in the longitudinal part of the model without taking into account that this is a binary step function over time monotone non-decreasing, and also for the way I linked it to the time component i.e. using the current value association, since the impact of the value of this variable depends on time: its value is always 0 during the first days of assembly and only from a certain time onwards it is informative.

Could you help me correct the model definition?


Get this bounty!!!

#StackBounty: #survival #cox-model #hazard #kaplan-meier Is censoring data necessary to calculate the hazard ratio between 2 KM curves?

Bounty: 50

I would like to know if censoring data is necessary to calculate the hazard ratio between 2 Kaplan Meier (KM) curves.

Censoring data is typically represented by small vertical bars atop of the curve on KM graphs.

Censoring data has an impact on the shape of a KM curve already (for instance, if 98/100 people drop from a trial, and 1 person dies, the curve goes down much more than if no one had dropped out). This makes me think that the curves are perhaps enough for calculating hazard ratios (since I understand hazard rates are akin to the slope of the curves). Is this true?

I am asking a practical question:
Is having the values of survival curves over time, without the censoring information, enough to estimate the hazard ratio between these curves?

Context:
I am right now digitizing KM curves. I would like to estimate the hazard ratio between them. I cannot retrieve the censoring data. Not having the censoring data is a limitation to my analysis, I know. But I would like to know: Can I nonetheless calculate the hazard ratio? Can I calculate a confidence interval for this hazard ratio?

Ideally, please let me know how the answer differs depending on whether or not we assume relative hazards between the 2 groups are constant across time.


Get this bounty!!!

#StackBounty: #r #survival #exponential-distribution rpart with survival data: what do the reported rates mean?

Bounty: 100

This is essentially the same question as this one but the accepted answer doesn’t answer it.

I made data consisting of three groups, each following an exponential distribution with rates 1, 3 and 6 respectively. Then I aksed rpart if it could recognize these groups, which it can, but then for each group it reports a rate that does not make any sense. I was told that rpart ‘scales’ the rates so that the ‘overall’ rate becomes 1, but how this scaling takes place is a mystery to me as I would expect that the scaling would preserve ratios, which clearly it doesn’t. Here is the code:

library(survival)
library(rpart)

#create data:

fakedata <- data.frame(patnr = seq(1, 100), A = rep(c(0, 1), each = 50), B = c(rep(0, 75), rep(1, 25)))
fakedata$lambda <- 1 + 2*fakedata$A + 3*fakedata$B                       
set.seed(2020)
fakedata$survivalyears <- rexp(100, rate = fakedata$lambda)

#to be clear: these are the rates and how often they appear:
table(fakedata$lambda)

#which produces:
#  1  3  6 
# 50 25 25 

#now let's see what rpart makes of it:

fit <- rpart(Surv(survivalyears, rep(1, 100)~ A + B, data = fakedata)

Now rpart gets the split into 50 vs 50 based on A and the further split of the 50 A == 1 patients into 25 + 25 based on B right, as indicated by the rpart output (not shown here) and also by the off-diagonal zeros in the table below. However, the rates it produces are a complete mystery to me:

fakedata$predictedrate <- predict(fit)
table(fakedata[, c("lambda", "predictedrate"))

#which produces:
#       predictedrate
# lambda 0.673787511221765 1.46471239153325 2.7203606776826
#      1                50                0               0
#      3                 0               25               0
#      6                 0                0              25 

What is the ratio behind (and between) those numbers? I would expect to find 2.7203606776826/0.673787511221765 equal to 6 and 1.46471239153325/0.673787511221765 equal to 3 but obviously this is false.

What is going on here?


Get this bounty!!!

#StackBounty: #r #logistic #survival #recurrent-events (Repeated) Event Model – How to model entry?

Bounty: 50

I am modeling the entry of 30 firms into 25 industries over time. The unit of analysis are industries.

I have cross-sectional panel data on the industries and know when a firm enters. I have 5 industry-level co-variates that change over time.

I would like to know which of these 5 co-variates best predicts the entry of a firm into the industry.

Note: In some industries multiple firms enter (at different points of time).

What type of model do I use?

  • A logistic panel model with fixed effects for industry (and year)?
  • An extended survival model? Perhaps a Cox regression that models repeated events?
  • Or something different altogether?

I am looking for a statistical approach to be performed in R.

Example data:

dt <- data.table(firm = c(NA, "F1", NA, "F2",NA, NA, NA, "F1", NA, NA, NA, "F3","F4", NA, NA, NA),
             industry = c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4),
             year = c(2000, 2001, 2002, 2003,2000, 2001, 2002, 2003,2000, 2001, 2002, 2003,2000, 2001, 2002, 2003),
             entry= c(0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0),
             x1 = c(runif(16)),
             x2 = c(runif(16)),
             x3 = c(runif(16)),
             x4 = c(runif(16)),
             x5 = c(runif(16)))

> dt
    firm industry year entry        x1         x2        x3         x4         x5
 1: <NA>        1 2000     0 0.8709432 0.53087423 0.8454785 0.25285903 0.05671375
 2:   F1        1 2001     1 0.2775248 0.57784857 0.3895359 0.29725940 0.17875968
 3: <NA>        1 2002     0 0.4863695 0.44049405 0.6658640 0.46859840 0.72456472
 4:   F2        1 2003     1 0.4604883 0.47147933 0.1241863 0.30133078 0.88424245
 5: <NA>        2 2000     0 0.8946535 0.41725531 0.7416106 0.47860660 0.15310571
 6: <NA>        2 2001     0 0.4675261 0.17393733 0.2680798 0.24721979 0.25179696
 7: <NA>        2 2002     0 0.1801262 0.76724741 0.3157204 0.20496050 0.21098053
 8:   F1        2 2003     1 0.9967714 0.54301656 0.9133624 0.41144085 0.39097097
 9: <NA>        3 2000     0 0.4584210 0.26767502 0.4698928 0.62344696 0.96677505
10: <NA>        3 2001     0 0.5241493 0.21622151 0.7802806 0.89008628 0.28688089
11: <NA>        3 2002     0 0.1472759 0.58341162 0.5209162 0.42375726 0.06895383
12:   F3        3 2003     1 0.4255832 0.09506578 0.4067654 0.25417083 0.32360013
13:   F4        4 2000     1 0.1922379 0.61931944 0.3609198 0.99608796 0.57967160
14: <NA>        4 2001     0 0.2681883 0.26763619 0.6279573 0.57704165 0.07834483
15: <NA>        4 2002     0 0.2431978 0.69426770 0.4675626 0.77913438 0.87128792
16: <NA>        4 2003     0 0.7301488 0.94025212 0.5586453 0.05350145 0.98182507


Get this bounty!!!

#StackBounty: #r #survival #cox-model #kaplan-meier Extended Kaplan-Meier for time-dependent covariates

Bounty: 50

I have read Snapinn et al. paper on "Illustrating the Impact of a Time-Varying Covariate With an Extended Kaplan-Meier Estimator" (https://doi.org/10.1198/000313005X70371). They describe an extended Kaplan-Meier survival analysis for working with time-dependent covariates. Please reach out if you need the full text article.

Figure from Snapinn paper

I have been trying to replicate their results but I am reaching out to you to validate my results. Since they haven’t made their data available, I made up some data and I already ran the following in R:

library(survival)
library(survminer)
library(tidyverse)

set.seed(99)

#arbitrary data
df1 <- data.frame(ID = rep(seq(1, 400, by = 1), 2),
                  score = factor(sample(1:4, 200, replace = TRUE)),
                  timetoFU = sample(1:200, 200, replace = TRUE),
                  status = sample(c(0, 1), 200, replace = TRUE, prob = c(0.9, 0.2))
                  )

df1 <- df1 %>% group_by(ID) %>% arrange(ID, timetoFU) %>% mutate(obs_n = row_number(), time_max = last(timetoFU)+sample(1:50, 1, replace = TRUE)) %>% ungroup()
df1$timetoFU[which(df1$obs_n == 1)] <- 0
#events more likely with higher scores
temp <- df1 %>% group_by(ID) %>% summarise(risk=sum(as.numeric(score))/10, status, time_max, obs_n) %>% ungroup() %>% mutate(status = ifelse(risk >0.4, sample(c(0,1), 100, replace = TRUE), status)) %>% mutate(status = ifelse(risk > 0.6, 1, status)) %>% filter(obs_n == 1)

#build time dependent variable data frame
td_df <- tmerge(temp, temp, id = ID, outcome = event(time_max, status))
td_df <- tmerge(td_df, df1, id = ID, td_score = tdc(timetoFU, score))

#survival analysis

s1 <- survfit(Surv(tstart, tstop, outcome) ~  td_score, data = td_df, id = ID)
ggsurvplot(s1, fun = "event", risk.table = TRUE, conf.int = TRUE, break.x.by = 10)

cox_fit_A <- coxph(formula = Surv(tstart, tstop, outcome) ~ td_score, data = td_df, id = ID)
summary(cox_fit_A)

Replication

Is this correct?


Get this bounty!!!

#StackBounty: #r #survival #cox-model #kaplan-meier Replicating extended Kaplan-Meier rates for time-dependent covariates

Bounty: 50

I have read Snapinn et al. paper on "Illustrating the Impact of a Time-Varying Covariate With an Extended Kaplan-Meier Estimator" (https://doi.org/10.1198/000313005X70371). They describe an extended Kaplan-Meier survival analysis for working with time-dependent covariates. Please reach out if you need the full text article.

Figure from Snapinn paper

I have been trying to replicate their results but I am reaching out to you because I am not sure whether I got it right. Since they haven’t made their data available, I made up some data and I already ran the following in R:

library(survival)
library(survminer)
library(tidyverse)

set.seed(99)

#arbitrary data
df1 <- data.frame(ID = rep(seq(1, 400, by = 1), 2),
                  score = factor(sample(1:4, 200, replace = TRUE)),
                  timetoFU = sample(1:200, 200, replace = TRUE),
                  status = sample(c(0, 1), 200, replace = TRUE, prob = c(0.9, 0.2))
                  )

df1 <- df1 %>% group_by(ID) %>% arrange(ID, timetoFU) %>% mutate(obs_n = row_number(), time_max = last(timetoFU)+sample(1:50, 1, replace = TRUE)) %>% ungroup()
df1$timetoFU[which(df1$obs_n == 1)] <- 0
#events more likely with higher scores
temp <- df1 %>% group_by(ID) %>% summarise(risk=sum(as.numeric(score))/10, status, time_max, obs_n) %>% ungroup() %>% mutate(status = ifelse(risk >0.4, sample(c(0,1), 100, replace = TRUE), status)) %>% mutate(status = ifelse(risk > 0.6, 1, status)) %>% filter(obs_n == 1)

#build time dependent variable data frame
td_df <- tmerge(temp, temp, id = ID, outcome = event(time_max, status))
td_df <- tmerge(td_df, df1, id = ID, td_score = tdc(timetoFU, score))

#survival analysis

s1 <- survfit(Surv(tstart, tstop, outcome) ~  td_score, data = td_df, id = ID)
ggsurvplot(s1, fun = "event", risk.table = TRUE, conf.int = TRUE, break.x.by = 10)

cox_fit_A <- coxph(formula = Surv(tstart, tstop, outcome) ~ td_score, data = td_df, id = ID)
summary(cox_fit_A)

Replication

Is this correct?


Get this bounty!!!

#StackBounty: #self-study #survival #interpretation #weibull Interpretation of Weibull Accelerated Failure Time Model Output

Bounty: 50

In this case study I have to assume a baseline Weibull distribution, and I’m fitting an Accelerated Failure Time model, which will be interpreted by me later on regarding both hazard ratio and survival time.

The data looks like this.

head(data1.1)

TimeSurv IndSurv Treat Age
1     6 days       1     D  27
2    33 days       1     D  43
3   361 days       1     I  36
4   488 days       1     I  54
5   350 days       1     D  49
6   721 days       1     I  49
7  1848 days       0     D  32
8   205 days       1     D  47
9   831 days       1     I  24
10  260 days       1     I  38

I’m fitting a model using the function Weibullreg() in R. The survival function is built reading TimeSurv as the time measures and IndSurv as the indicator of censoring. The covariates considered are Treat and Age.

My issue deals with understanding the output properly:

wei1 = WeibullReg(Surv(TimeSurv, IndSurv) ~ Treat + Age, data=data1.1)
wei1


$formula
Surv(TimeSurv, IndSurv) ~ Treat + Age

$coef
            Estimate           SE
lambda  0.0009219183 0.0006803664
gamma   0.9843411517 0.0931305471
TreatI -0.5042111027 0.2303038312
Age     0.0180225253 0.0089632209

$HR
              HR       LB       UB
TreatI 0.6039819 0.384582 0.948547
Age    1.0181859 1.000455 1.036231

$ETR
             ETR        LB        UB
TreatI 1.6690124 1.0574337 2.6343045
Age    0.9818574 0.9644488 0.9995801

$summary

Call:
survival::survreg(formula = formula, data = data, dist = "weibull")
               Value Std. Error     z      p
(Intercept)  7.10024    0.41283 17.20 <2e-16
TreatI       0.51223    0.23285  2.20  0.028
Age         -0.01831    0.00913 -2.01  0.045
Log(scale)   0.01578    0.09461  0.17  0.868

Scale= 1.02 

Weibull distribution
Loglik(model)= -599.1   Loglik(intercept only)= -604.1
    Chisq= 9.92 on 2 degrees of freedom, p= 0.007 
Number of Newton-Raphson Iterations: 5 
n= 120

I don’t really get how Scale = 1.02 and log(scale) = 0.015, and if the p-value of this log(scale) is a big non-signfificant one, from how the documentation of the function shows some conversions it makes, am I to assume that the values of the alphas are also not to be trusted (considering they were reached using the scale value)?


Get this bounty!!!

#StackBounty: #r #survival #data-transformation #cox-model rms/R: How to apply survSplit on two time-depedent covariates, one being a d…

Bounty: 100

I am doing a survival analysis of time p$os.neck to death p$mors using a Cox Regression.

Please, find my data sample p below.

Question: how can I apply a survSplit from the rms-package on this Cox Regression with two time-depedent covariates, one being a categorial covariate (p$uicc with four levels 1,2,3,4) and the other being a discrete covariate (p$n.sygdom currently ranging from 0 to 10 in p but theoretically could increase to higher values)?

First

library(rms)

p$sex <- factor(p$sex,levels=c("0","1"),labels=c("0","1"))
p$ecs <- factor(p$ecs,levels=c("0","1"),labels=c("0","1"))
p$uicc <- factor(p$uicc,levels=c("1","2","3","4"),labels=c("1","2","3","4"))
p$rt.kemo <- factor(p$rt.kemo,levels=c("0","1"),labels=c("0","1"))

And

d <- datadist(p)
options(datadist="d")

I then have

a < - cph(Surv(os.neck,mors)~alder+sex+n.fjernet+rcs(n.sygdom)+ecs+uicc+rt.kemo,data=p,surv=TRUE,x=TRUE,y=TRUE)

> cox.zph(a)
               chisq df      p
alder          0.539  1 0.4627
sex            0.593  1 0.4411
n.fjernet      1.052  1 0.3051
rcs(n.sygdom) 10.291  2 0.0058
ecs            0.646  1 0.4216
uicc          12.987  3 0.0047
rt.kemo        1.099  1 0.2945
GLOBAL        26.705 10 0.0029 

For the two time-depedent covariates:

> table(p$uicc)

  1   2   3   4 
126  99  59 146

And

> table(p$n.sygdom)

  0   1   2   3   4   5   6   7   9  10 
292  72  29  13  10   3   4   3   2   2 

Based on plot(cox.zph(a),var=..., I have found that one survival split at time=24 months may be adequate and should be investigated further.

However, I am not experienced in doing survSplit in case of (1) more than one time-dependent covariate and (2) other than categorial covariates with two levels, such as gender.

So, currently, I have something like

v <- survSplit(Surv(os.neck, mors) ~ ., cut=c(24), data=p, episode="time_group")

Please, how can I incorporate rcs(n.sygdom) and p$uicc in the abovementioned survSplit?

My data p

p <- structure(list(alder = structure(c(58.53, 51.43, 78.5, 48.44, 
68.61, 58.28, 55.06, 67.33, 86.51, 61.57, 76.98, 63.73, 63.72, 
55.29, 55.34, 60.85, 60.54, 56.13, 76.09, 71.54, 80.24, 81.67, 
59.49, 61.07, 58.28, 60.2, 58.57, 60, 71.95, 40.48, 81.41, 30.08, 
51.39, 62.44, 75.43, 69.68, 52.99, 34.77, 55.09, 57.18, 34.91, 
67.34, 68.6, 73.74, 52.82, 64.58, 59.18, 48.63, 73.14, 68.9, 
53.71, 58.13, 60.87, 55.65, 68.94, 61.49, 59.14, 89.1, 71.57, 
86.25, 59, 94.49, 46.5, 81.39, 57.28, 53.39, 60.37, 56.82, 73.79, 
62.41, 73.13, 48.68, 50.68, 65.01, 60.67, 71.99, 58.98, 50.76, 
64.04, 61.04, 65.57, 61, 67.92, 55.03, 54.33, 51.94, 82.55, 62.53, 
57.13, 65.87, 60.54, 60.93, 72.49, 61.87, 51.87, 63.94, 82.42, 
51.7, 76.35, 60.46, 65.49, 51.83, 61.07, 63.25, 74.82, 59.19, 
60.2, 52.85, 52.38, 53.64, 65.87, 59.94, 69.86, 60.91, 65.09, 
63.97, 67.49, 57.29, 50.1, 56.08, 76.79, 69.58, 58.48, 61.8, 
83.28, 66.18, 71.04, 45.58, 81.72, 52.92, 56.14, 56.2, 73.12, 
55.06, 63.84, 67.65, 45.81, 84.85, 65.72, 69.39, 63.69, 62.42, 
67.92, 44, 56.44, 87.48, 63.1, 54.79, 36.45, 28.08, 56.54, 52.56, 
59.92, 75.97, 47.35, 46.79, 29.12, 57.3, 66.9, 48.35, 49.7, 53.84, 
51.34, 53.83, 60.29, 72.79, 73.68, 73.63, 62.6, 32.78, 40.55, 
48.03, 67.11, 53.23, 70.34, 64.54, 87.24, 81.97, 55.27, 79.79, 
68.88, 53.22, 61.04, 63.91, 93.75, 58.33, 69.92, 63.66, 82.98, 
64.6, 74.47, 67.52, 65.67, 56.1, 71.71, 57.65, 83.1, 60.1, 49.07, 
59.52, 33.07, 49.69, 63.14, 40.61, 62.57, 78.63, 66.54, 55.35, 
55.43, 72.71, 65.31, 69.52, 69.03, 48.47, 56.74, 70.16, 56.94, 
95.7, 75.9, 67.49, 66.07, 78.65, 82.91, 63.76, 68.2, 54.28, 73.65, 
74.49, 76.37, 91.65, 66.31, 42.7, 68.14, 86.09, 38.79, 53.81, 
70.56, 63.36, 62.38, 77.92, 61.42, 50.07, 70.28, 63.85, 69.17, 
65.83, 58.17, 49.18, 50.27, 59.33, 53.08, 70.95, 62.99, 45.54, 
67.55, 57.72, 67.31, 59.91, 61.15, 69.92, 78.56, 68.9, 69.73, 
57.3, 51.94, 68.96, 60.58, 65.23, 67.02, 65.41, 64.12, 82.47, 
72.53, 58.44, 74.02, 75.52, 63.56, 66.73, 67.89, 60.17, 54.37, 
54.91, 58.34, 68.6, 60.02, 59.28, 48.95, 72.54, 54.16, 65.88, 
67.27, 45.78, 78.15, 36.62, 69.72, 61.72, 56.28, 69.47, 56.82, 
68.63, 73.13, 70.35, 55.47, 52.06, 87.93, 73.5, 66.1, 69.71, 
50.65, 62.57, 74.45, 63.75, 67.12, 79.28, 65.53, 63.38, 54.71, 
54.68, 68.66, 64.87, 94.64, 75.63, 88.05, 51.13, 66.58, 56.24, 
51.39, 52.47, 46.08, 59.73, 52.8, 64.19, 63.6, 68.64, 73.52, 
68.37, 57.05, 77.54, 70.7, 53.69, 68.34, 76.95, 51.52, 69.73, 
55.36, 56.26, 61.88, 60.64, 71.92, 69.59, 75.28, 71.66, 59.23, 
58.2, 61.8, 66.01, 56.3, 46.69, 45.61, 62.79, 59.76, 66.75, 73.65, 
48.46, 51.56, 79.86, 47.76, 58.45, 45.84, 64.38, 56.4, 63.02, 
49.47, 57.17, 68.35, 63.56, 61.11, 35.65, 61.18, 67.96, 75.21, 
62.62, 65.29, 74.27, 68.93, 61.2, 70.19, 51, 66.94, 53.47, 64.25, 
51.97, 67.07, 71.39, 58.03, 60.67, 73.35, 78.87, 75.14, 74.39, 
63.44, 79.67, 45.01, 58.78, 57.44, 67.86, 55.85, 65.79, 58.67, 
60.55, 76.89, 80.2, 62.94, 43.76, 65.12, 50.4, 67.4, 45.98, 23.17, 
30.57, 57.62, 70.49, 43.84, 77.53, 45.88, 63.86, 63.11, 68.27, 
83.6, 57.02), label = c(alder = "Age"), class = c("labelled", 
"numeric")), n.fjernet = structure(c(4L, 27L, 18L, 11L, 14L, 
15L, 9L, 6L, 3L, 16L, 4L, 6L, 10L, 13L, 33L, 16L, 6L, 9L, 15L, 
23L, 5L, 9L, 10L, 8L, 17L, 14L, 13L, 13L, 5L, 9L, 30L, 16L, 9L, 
25L, 3L, 19L, 10L, 8L, 9L, 9L, 10L, 12L, 7L, 38L, 21L, 24L, 5L, 
7L, 15L, 4L, 4L, 35L, 9L, 6L, 10L, 15L, 9L, 8L, 7L, 4L, 21L, 
6L, 10L, 6L, 3L, 8L, 4L, 9L, 10L, 14L, 14L, 3L, 4L, 6L, 6L, 20L, 
7L, 6L, 17L, 3L, 26L, 13L, 13L, 14L, 19L, 13L, 13L, 3L, 7L, 6L, 
8L, 18L, 23L, 6L, 5L, 6L, 5L, 4L, 10L, 7L, 15L, 29L, 13L, 18L, 
7L, 7L, 26L, 18L, 27L, 4L, 22L, 15L, 6L, 20L, 11L, 13L, 17L, 
17L, 26L, 8L, 5L, 14L, 17L, 17L, 9L, 12L, 56L, 16L, 18L, 35L, 
28L, 22L, 12L, 7L, 24L, 9L, 17L, 16L, 20L, 16L, 21L, 20L, 34L, 
7L, 9L, 8L, 4L, 8L, 6L, 8L, 16L, 6L, 11L, 3L, 15L, 3L, 10L, 4L, 
4L, 9L, 6L, 5L, 5L, 3L, 30L, 6L, 2L, 4L, 8L, 5L, 5L, 8L, 16L, 
18L, 7L, 12L, 9L, 9L, 13L, 9L, 22L, 20L, 24L, 8L, 18L, 8L, 15L, 
19L, 5L, 4L, 14L, 18L, 18L, 11L, 15L, 22L, 46L, 11L, 18L, 13L, 
9L, 12L, 13L, 26L, 8L, 30L, 11L, 14L, 22L, 23L, 26L, 5L, 4L, 
26L, 32L, 6L, 9L, 11L, 22L, 6L, 25L, 15L, 22L, 20L, 35L, 5L, 
5L, 20L, 8L, 18L, 7L, 15L, 22L, 13L, 7L, 20L, 11L, 4L, 2L, 7L, 
7L, 4L, 11L, 13L, 13L, 9L, 9L, 9L, 12L, 11L, 13L, 16L, 6L, 13L, 
8L, 17L, 5L, 8L, 22L, 12L, 19L, 3L, 15L, 14L, 7L, 18L, 24L, 9L, 
27L, 9L, 6L, 9L, 4L, 21L, 10L, 36L, 18L, 24L, 19L, 11L, 8L, 15L, 
37L, 7L, 7L, 6L, 18L, 9L, 4L, 22L, 5L, 2L, 24L, 2L, 23L, 30L, 
55L, 9L, 24L, 7L, 8L, 20L, 9L, 22L, 11L, 2L, 24L, 15L, 30L, 5L, 
10L, 8L, 11L, 11L, 11L, 15L, 6L, 16L, 7L, 9L, 16L, 11L, 33L, 
5L, 27L, 27L, 16L, 57L, 5L, 7L, 8L, 11L, 15L, 15L, 12L, 5L, 25L, 
9L, 21L, 13L, 3L, 55L, 27L, 28L, 33L, 23L, 49L, 49L, 11L, 7L, 
28L, 19L, 13L, 23L, 4L, 5L, 11L, 12L, 10L, 4L, 14L, 6L, 12L, 
7L, 32L, 13L, 5L, 12L, 10L, 4L, 4L, 11L, 8L, 17L, 25L, 10L, 8L, 
5L, 15L, 21L, 19L, 11L, 31L, 9L, 20L, 11L, 16L, 12L, 6L, 16L, 
27L, 30L, 18L, 18L, 10L, 7L, 23L, 16L, 15L, 4L, 12L, 9L, 10L, 
12L, 11L, 7L, 7L, 8L, 8L, 8L, 7L, 6L, 9L, 9L, 13L, 15L, 12L, 
35L, 12L, 5L, 5L, 19L, 13L, 27L, 34L, 10L, 16L, 18L, 6L, 22L), label = c(n.fjernet = "LNY"), class = c("labelled", 
"integer")), n.sygdom = structure(c(0L, 0L, 4L, 1L, 0L, 0L, 0L, 
0L, 0L, 4L, 0L, 0L, 0L, 0L, 4L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 
3L, 0L, 0L, 0L, 0L, 0L, 0L, 3L, 5L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 
0L, 1L, 0L, 3L, 0L, 0L, 1L, 1L, 0L, 0L, 5L, 1L, 1L, 0L, 0L, 0L, 
0L, 0L, 3L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 0L, 1L, 1L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 2L, 1L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 
1L, 0L, 10L, 3L, 0L, 0L, 0L, 0L, 0L, 6L, 1L, 2L, 0L, 0L, 0L, 
0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 2L, 
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 2L, 1L, 2L, 1L, 0L, 0L, 3L, 
0L, 0L, 1L, 0L, 0L, 4L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 2L, 0L, 0L, 
0L, 0L, 0L, 0L, 3L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 3L, 3L, 0L, 0L, 
2L, 0L, 0L, 0L, 3L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 
0L, 2L, 10L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 
0L, 0L, 0L, 0L, 0L, 7L, 4L, 0L, 2L, 1L, 0L, 4L, 0L, 2L, 0L, 7L, 
0L, 4L, 6L, 2L, 0L, 0L, 1L, 1L, 0L, 2L, 1L, 0L, 2L, 3L, 2L, 0L, 
0L, 0L, 0L, 4L, 0L, 0L, 1L, 1L, 1L, 0L, 2L, 3L, 2L, 0L, 1L, 3L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 
1L, 0L, 2L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 2L, 0L, 1L, 0L, 2L, 2L, 
0L, 0L, 0L, 0L, 9L, 0L, 2L, 6L, 0L, 9L, 0L, 1L, 0L, 7L, 0L, 0L, 
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 2L, 5L, 2L, 4L, 6L, 0L, 0L, 
1L, 0L, 4L, 0L, 0L, 1L, 1L, 2L, 1L), label = c(n.sygdom = "No. LN+"), class = c("labelled", 
"integer")), ecs = structure(c(1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 
1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 
1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 
1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 
2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 
2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 
1L, 1L, 2L, 2L, 1L, 1L), .Label = c("0", "1"), class = c("labelled", 
"factor"), label = c(ecs = "ECS")), uicc = structure(c(2L, 2L, 
4L, 3L, 3L, 2L, 2L, 2L, 2L, 4L, 1L, 1L, 2L, 1L, 4L, 2L, 1L, 2L, 
3L, 3L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 3L, 2L, 1L, 1L, 1L, 1L, 2L, 
3L, 2L, 3L, 1L, 2L, 4L, 1L, 1L, 1L, 2L, 1L, 4L, 4L, 4L, 1L, 3L, 
4L, 1L, 4L, 3L, 1L, 4L, 3L, 1L, 4L, 2L, 2L, 3L, 4L, 2L, 1L, 4L, 
4L, 3L, 2L, 4L, 1L, 4L, 2L, 4L, 4L, 2L, 1L, 1L, 1L, 4L, 4L, 1L, 
4L, 3L, 2L, 2L, 3L, 2L, 2L, 2L, 4L, 4L, 2L, 3L, 2L, 2L, 2L, 1L, 
4L, 4L, 4L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 4L, 1L, 4L, 2L, 3L, 1L, 
1L, 1L, 4L, 4L, 2L, 3L, 4L, 4L, 4L, 2L, 2L, 4L, 2L, 2L, 4L, 4L, 
4L, 4L, 2L, 1L, 1L, 4L, 3L, 4L, 2L, 4L, 3L, 3L, 2L, 3L, 2L, 2L, 
1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 
1L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, 4L, 2L, 1L, 2L, 4L, 1L, 1L, 1L, 
2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 3L, 1L, 2L, 3L, 
1L, 1L, 3L, 4L, 1L, 1L, 1L, 2L, 4L, 1L, 1L, 1L, 3L, 4L, 3L, 4L, 
4L, 1L, 2L, 4L, 2L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 4L, 4L, 
4L, 2L, 1L, 4L, 1L, 1L, 3L, 1L, 3L, 4L, 2L, 4L, 2L, 3L, 3L, 4L, 
1L, 1L, 3L, 1L, 4L, 2L, 1L, 3L, 4L, 1L, 2L, 1L, 1L, 4L, 1L, 1L, 
4L, 4L, 1L, 1L, 3L, 2L, 2L, 1L, 4L, 1L, 1L, 4L, 2L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 4L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 
1L, 1L, 4L, 1L, 4L, 4L, 4L, 4L, 1L, 1L, 2L, 3L, 4L, 2L, 4L, 1L, 
1L, 3L, 1L, 3L, 2L, 1L, 1L, 3L, 4L, 4L, 2L, 4L, 4L, 3L, 4L, 4L, 
4L, 1L, 4L, 1L, 4L, 4L, 3L, 2L, 2L, 4L, 3L, 1L, 4L, 3L, 3L, 4L, 
4L, 4L, 2L, 3L, 4L, 3L, 4L, 1L, 1L, 4L, 3L, 3L, 1L, 4L, 4L, 4L, 
2L, 3L, 4L, 2L, 2L, 4L, 4L, 1L, 4L, 2L, 4L, 2L, 1L, 4L, 3L, 1L, 
4L, 4L, 3L, 3L, 2L, 4L, 2L, 3L, 3L, 4L, 4L, 2L, 4L, 4L, 2L, 4L, 
4L, 4L, 4L, 1L, 4L, 4L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 1L, 2L, 4L, 3L, 2L, 1L, 2L, 1L, 2L, 2L, 4L, 4L, 4L, 4L, 4L, 
4L, 2L, 1L, 3L, 1L, 4L, 4L, 1L, 3L, 3L, 4L, 3L), .Label = c("1", 
"2", "3", "4"), class = c("labelled", "factor"), label = c(uicc = "UICC Stage")), 
    rt.kemo = structure(c(2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 
    2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 
    2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
    1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 
    2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 
    1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 
    1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 
    2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
    1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 
    2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 
    1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 
    2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 
    2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 
    1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 
    1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
    1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 
    1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 
    2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 
    2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 
    1L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 
    2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 
    1L), .Label = c("0", "1"), class = c("labelled", "factor"
    ), label = c(rt.kemo = "Radiochemotherapy")), sex = structure(c(2L, 
    2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 
    1L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 
    1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 
    1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 
    1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 
    1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 
    1L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 
    2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 
    1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 
    1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 
    2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 
    1L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 
    1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
    1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 
    1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 
    2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 
    2L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 
    1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 
    2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 
    2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 
    1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 
    1L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 
    1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 
    1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L), .Label = c("Female", 
    "Male"), class = c("labelled", "factor"), label = c(sex = "Sex")), 
    mors = structure(c(0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 
    0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 
    0L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 
    1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 
    1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 
    1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 
    0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 
    0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 
    1L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 
    0L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 
    0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 
    1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 
    1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 
    1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 
    0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 
    0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 
    0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 
    0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 
    0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 
    1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L
    ), label = c(os.neck = "os.neck"), class = c("labelled", 
    "integer")), os.neck = structure(c(77.01, 75.96, 11.5, 74.38, 
    17.02, 7.89, 96.03, 40.48, 17.74, 14.65, 62.46, 12.55, 9.92, 
    26.05, 45.47, 17.38, 39.72, 51.45, 119, 8.61, 117.39, 76.98, 
    115.78, 67.09, 113.74, 113.22, 111.64, 94.79, 72.15, 110.23, 
    93.93, 108.16, 106.91, 17.05, 12.48, 104.22, 103.69, 131.98, 
    91.6, 15.87, 101.85, 11.04, 67.22, 67.02, 120.28, 149.88, 
    8.94, 6.6, 5.09, 10.68, 150.21, 135.4, 128.69, 17.15, 122.78, 
    0.07, 5.19, 40.77, 0.2, 170.88, 164.7, 5.55, 1.61, 162.11, 
    167.53, 38.28, 10.58, 32.99, 110.98, 103.69, 122.32, 14.78, 
    42.74, 4.04, 8.28, 84.96, 144.04, 150.67, 145.05, 11.7, 49.97, 
    120.48, 52.6, 139.04, 137.83, 71.26, 16.3, 100.14, 55.03, 
    130.96, 123.44, 118.67, 114.04, 6.51, 119.1, 112.76, 89.89, 
    114.83, 51.71, 95.84, 24.97, 55.66, 85.39, 77.73, 83.42, 
    21.91, 88.41, 86.9, 85.92, 84.17, 71.56, 77.08, 81.48, 79.21, 
    30.92, 68.27, 1.58, 67.65, 64.53, 71.66, 61.47, 7.52, 61.21, 
    61.93, 61.14, 36.34, 35.71, 35.61, 30.75, 34.17, 32.3, 3.45, 
    32.89, 32.76, 31.93, 19.22, 31.74, 30.62, 28.72, 30, 29.64, 
    5.42, 17.68, 178.7, 45.54, 76.22, 151.07, 125.34, 146.96, 
    143.08, 142.36, 140.95, 83.62, 30.82, 137.92, 137.56, 136.41, 
    90.32, 1.84, 135.23, 134.34, 133.62, 19.98, 20.53, 130.47, 
    128.33, 32.59, 128.53, 54.77, 126.52, 2.3, 125.67, 125.64, 
    106.84, 22.28, 90.38, 82.99, 45.18, 4.47, 80.76, 80.46, 80, 
    78.23, 77.83, 39.66, 74.74, 71.33, 32.3, 70.41, 71.95, 16.23, 
    66.63, 64.13, 58.58, 57.92, 3.68, 3.88, 47.9, 47.02, 46.72, 
    46.69, 45.44, 44.55, 44.62, 40.87, 41.73, 40.84, 39.82, 37.98, 
    2.23, 31.38, 52.04, 23.59, 29.24, 28.32, 91.99, 74.09, 0.23, 
    62.39, 18.73, 56.31, 53.03, 45.37, 43.07, 43.37, 41.66, 36.63, 
    28.95, 29.24, 0.79, 27.07, 144.92, 33.61, 83.32, 180.34, 
    28.75, 29.83, 79.54, 14.46, 15.15, 54.97, 48.59, 34.83, 58.42, 
    35.29, 45.73, 57.53, 63.11, 65.05, 29.54, 132.57, 77.21, 
    63.48, 83.35, 34.3, 64.49, 29.54, 62.69, 21.62, 67.52, 49.35, 
    99.02, 15.8, 41.89, 12.98, 13.8, 35.19, 163.78, 44.81, 43.6, 
    90.48, 81.68, 36.14, 137.96, 57.23, 94.33, 31.38, 70.74, 
    59.34, 39.46, 32.07, 20.76, 49.94, 67.22, 91.11, 127.15, 
    121.56, 89.6, 74.12, 31.8, 77.31, 159.35, 1.97, 40.38, 7.39, 
    40.54, 40.02, 38.9, 38.41, 37.49, 25.17, 28.22, 14, 36.53, 
    20.83, 19.55, 40.77, 27.76, 62.56, 45.31, 42.32, 34.46, 35.55, 
    26.94, 9.43, 10.51, 6.8, 8.18, 8.02, 14.29, 6.11, 13.8, 4.9, 
    141.21, 4.04, 40.94, 14.82, 11.66, 73.07, 92.91, 99.98, 10.64, 
    10.05, 95.8, 7.23, 12.81, 114.93, 43.99, 61.93, 66.2, 34, 
    32.99, 30.39, 48.69, 29.31, 27.34, 33.18, 13.9, 10.25, 45.04, 
    16.36, 18.2, 18.76, 12.32, 145.12, 173.7, 8.64, 11.79, 112.04, 
    70.97, 31.28, 28.85, 21.49, 138.68, 19.94, 22.14, 148.31, 
    29.44, 175.61, 164.08, 67.62, 11.01, 84.17, 45.24, 46.82, 
    110.72, 154.71, 20.24, 14.06, 12.88, 31.51, 8.08, 13.08, 
    21.45, 24.28, 21.98, 32.89, 23.26, 15.41, 15.41, 13.8, 40.12, 
    8.02, 15.77, 49.81, 18.17, 24.21, 47.08, 6.6, 37.16, 13.01, 
    8.38, 14.36, 91.86, 18.27, 80.43, 17.28, 66.76, 73.76, 68.21, 
    22.83, 2.66, 69.06, 17.05, 8.61, 23.33, 13.34, 12.65, 8.77, 
    152.45, 128.92, 16.1, 42.28, 4.99, 11.73, 22.97, 40.12, 20.37, 
    2.04, 45.73), label = c(mors = "mors"), class = c("labelled", 
    "numeric"))), row.names = c(NA, 430L), class = "data.frame")


Get this bounty!!!