Correlaties

Een correlatie is een associatiemaat. Het geeft de sterkte van een relatie tussen twee variabelen. De parametrische vorm is de Pearson methode. Deze methode berekent de p-waarde en correlatie exact en geeft daarom ook een confidence interval. De Kendall en Spearman methodes* zijn non-parametrisch en benaderen de p-waarde en geven daarom geen confidence interval. Er zijn meer associatiematen beschikbaar maar vallen buiten het spectrum van deze handleiding.

Correlatie coefficient: De output van een correlatie test is een correlatie coefficient. Dit is een waarde die ligt tussen-1 en 1. 0 = geen correlatie. Toenemend of afnemend vanaf de 0 is er toenemende overeenstemming (-1 is inverse correlatie) tussen de variabelen.

Pearson’s product moment coefficient: Dit is een parametrische test gebaseerd op de sterkte van de lineare samenhang tussen twee (normaal verdeelde continue/interval) variabelen. Deze test is niet geschikt voor ordinale data. Lees meer

Spearman’s rank correlatie coefficient: Dit is een non parametrische test gebaseerd op de rangnummers van de data in plaats van de data zelf. Deze test is ook geschikt voor ordinale data. Lees meer

Kendall tau rank correlatie coefficient: Dit is een non parametrische test gebaseerd op de rangnummers van de data in plaats van de data zelf. Deze test is ook geschikt voor ordinale data. Lees meer

Correlatie testen

Voorbeeld van correlatie coefficiënten:

# gebruik de functie cor() om de correlatie te berekenen middels de verschillende methoden
cor(mydata$age,mydata$mmse, use="complete.obs", method="spearman")
## [1] -0.3709608
cor(mydata$age,mydata$mmse, use="complete.obs", method="pearson")
## [1] -0.4402338
cor(mydata$age,mydata$mmse, use="complete.obs", method="kendall")
## [1] -0.2732419

Voorbeeld van een correlatie test:

# gebruik de cor.test() functie
cor.test(mydata$age,mydata$visus, method = "kendall", 
         alternative = "two.sided")                     # 2 zijdige p-waarde
## 
##  Kendall's rank correlation tau
## 
## data:  mydata$age and mydata$visus
## z = -5.8548, p-value = 4.776e-09
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##        tau 
## -0.1781067
cor.test(mydata$age,mydata$visus, method = "pearson",
         alternative = "two.sided")                     # 2 zijdige p-waarde
## 
##  Pearson's product-moment correlation
## 
## data:  mydata$age and mydata$visus
## t = -6.848, df = 591, p-value = 1.882e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3441412 -0.1948770
## sample estimates:
##        cor 
## -0.2711384
cor.test(mydata$age,mydata$visus, method = "spearman", 
         alternative = "two.sided")                     # 2 zijdige p-waarde
## 
##  Spearman's rank correlation rho
## 
## data:  mydata$age and mydata$visus
## S = 43281292, p-value = 1.407e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.245342

Correlatie matrix:

library(Hmisc)
## Correlatie matrix tussen numerieke variabelen
# tussenstap 1: selecteer numerieke variabelen
mydata_num <- subset(mydata, select=c("age", "mmse", "visus", "apache", "creaurer"))
head(mydata_num)
##   age mmse visus apache creaurer
## 1  87   19  0.30     18    11.20
## 2  92   28  0.40     19    12.70
## 3  83   27  0.50     10     9.40
## 4  82   24  0.30     14     9.51
## 5  95   17  0.20     12     7.72
## 6  94   16  0.05     15     6.36
# correlatie matrix
cor(mydata_num)
##                 age        mmse       visus      apache    creaurer
## age       1.0000000 -0.44023384 -0.27113838  0.39112972 -0.19080041
## mmse     -0.4402338  1.00000000  0.36285760 -0.35545868  0.07243213
## visus    -0.2711384  0.36285760  1.00000000 -0.28421350  0.06874053
## apache    0.3911297 -0.35545868 -0.28421350  1.00000000  0.03286541
## creaurer -0.1908004  0.07243213  0.06874053  0.03286541  1.00000000
## correlatie matrix met p-waarde. Zie http://goo.gl/nahmV voor bron van deze functie
cor.prob <- function (X, dfr = nrow(X) - 2) {
R <- cor(X, use="pairwise.complete.obs")
above <- row(R) < col(R)
r2 <- R[above]^2
Fstat <- r2 * dfr/(1 - r2)
R[above] <- 1 - pf(Fstat, 1, dfr)
R[row(R) == col(R)] <- NA
print(R)
}
# correlatie matrix met p-waarde
cor.prob(mydata_num)
##                 age        mmse         visus       apache     creaurer
## age              NA  0.00000000  1.881872e-11 0.000000e+00 2.875095e-06
## mmse     -0.4402338          NA  0.000000e+00 0.000000e+00 7.799643e-02
## visus    -0.2711384  0.36285760            NA 1.758593e-12 9.444783e-02
## apache    0.3911297 -0.35545868 -2.842135e-01           NA 4.243761e-01
## creaurer -0.1908004  0.07243213  6.874053e-02 3.286541e-02           NA

Correlatie tabel:

## Gebruik bovenstaande functie om een correlatie-matrix te maken.
## Gebruik onderstaande functie om van een correlatie-matrix een tabel te maken.
# zie StackOverflow discussie: http://goo.gl/fCUcQ voor de bron.
cortable <- function(m) {
if( (class(m) != "matrix") | (nrow(m) != ncol(m))) stop("De data moet een matrix zijn, gebruik cor.prob()")
if(!identical(rownames(m), colnames(m))) stop("rij en kolomnamen moeten hetzelfde zijn.")
ut <- upper.tri(m)
data.frame(i = rownames(m)[row(m)[ut]],
j = rownames(m)[col(m)[ut]],
cor=t(m)[ut],
p=m[ut])
}
# Correlatie tabel met de functies corprob() en cortable()
cortable(cor.prob(mydata_num))
##                 age        mmse         visus       apache     creaurer
## age              NA  0.00000000  1.881872e-11 0.000000e+00 2.875095e-06
## mmse     -0.4402338          NA  0.000000e+00 0.000000e+00 7.799643e-02
## visus    -0.2711384  0.36285760            NA 1.758593e-12 9.444783e-02
## apache    0.3911297 -0.35545868 -2.842135e-01           NA 4.243761e-01
## creaurer -0.1908004  0.07243213  6.874053e-02 3.286541e-02           NA
##         i        j         cor            p
## 1     age     mmse -0.44023384 0.000000e+00
## 2     age    visus -0.27113838 1.881872e-11
## 3    mmse    visus  0.36285760 0.000000e+00
## 4     age   apache  0.39112972 0.000000e+00
## 5    mmse   apache -0.35545868 0.000000e+00
## 6   visus   apache -0.28421350 1.758593e-12
## 7     age creaurer -0.19080041 2.875095e-06
## 8    mmse creaurer  0.07243213 7.799643e-02
## 9   visus creaurer  0.06874053 9.444783e-02
## 10 apache creaurer  0.03286541 4.243761e-01

Lees meer over de matrix en tabellen van correlaties.

Ordinale correlaties

Voor het analyseren van de sterkte van een ordinale factor variabele zijn de functies cor() en cortest() niet geschikt. Gebruik hiervoor de functie polyserial() en polychoric() uit de package polycor. Kies de correcte functie op basis van de type data en type relatie. Zie ‘??polyserial’ en ‘??polychoric’ voor meer informatie.

  • polyserial (numeriek - ordinaal) [parametrische test: zie ?polyserial]
  • polychoric (ordinaal - ordinaal) [parametrische test: zie ?polyserial]
# polyserial correlatie
library(polycor)
ordcor <- polyserial(mydata$age,mydata$delier)
summary(ordcor)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.4216 -0.4216 -0.4216 -0.4216 -0.4216 -0.4216
# polychoric correlatie
library(polycor)
ordcor2 <- polychor(mydata$medicijn,mydata$delier)
summary(ordcor2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.2778 -0.2778 -0.2778 -0.2778 -0.2778 -0.2778

Visualiseren van correlaties

Zie deze discussie op Stackoverflow voor meer manieren van het visualiseren van een correlatie matrix.

Scatterplot matrix:

## scatterplot matrix
# Gebruik de pairs() functie voor het maken van een scatterplot matrix
pairs(~ age + mmse + apache, data = mydata,
      main = "Scatterplot Matrix")

Visuele inspectie van correlatie per groep:

## scatterplot matrix PER GROEP 
# Voorbeeld met het 'iris' dataset
require(lattice)
  data(iris)
  super.sym <- trellis.par.get("superpose.symbol")
  splom(~iris[1:4], groups = Species, data = iris,
        panel = panel.superpose,
        key = list(title = "Three Varieties of Iris",
                   columns = 3, 
                   points = list(pch = super.sym$pch[1:3],
                                 col = super.sym$col[1:3]),
                   text = list(c("Setosa", "Versicolor", "Virginica"))))

Lees meer over visuele inspectie van correlatie’s per groep.

correlatie-matrix plot:

# maak een nieuwe data.frame met alleen numerieke variabelen
mydata_num <- mydata[c(4,5,6,7,14)] # apache + mmse + visus + creaurer + age

cormatrixplot(mydata_num)

Bron voor de functie cormatrixplot().


Regressie analyse

De belangrijkste vormen van regressie analyse zijn lineare regressie en logistische regressie. Dit zijn vormen van een Generalised Linear Model (GLM).

Een linear regressie model kan worden beschreven met de formule: \[\begin{equation} Y = \beta_0+ \beta_1X_1 + \beta_2X_2 + ... + \beta_nX_n + \varepsilon \end{equation}\] In deze formule is X de onafhankelijke variabele (IV) en Y is de afhankelijke variabele (DV). \(\beta_0\) is de intercept. \(\beta_1\) is de slope. \(\varepsilon\) is de noise / model error

Een logistische regressie analyse is te formuleren als: \[Z = \beta_1X_1 + \beta_2X_2 + ... + \beta_kX_k + \varepsilon\] De wiskundige formule: \[Y = \int (Z) = \frac{1}{1+e^{-Z}}\]

Type regressie analyses:

Uitkomstvariabele Type analyse Voorbeeld Formule

1 continue DV & 1 IV
1 continue DV & 2+ IV

Univariate lineare regressie
Multivariate lineare regressie

\(\begin{equation} \mathbf y=\beta_0+\beta_1 x + \varepsilon\end{equation}\)
\(\begin{equation} \mathbf y=\beta_0+\beta_1 x_1 +\beta_2 x_2 + \varepsilon\end{equation}\)

1 dichotome DV & 1 IV
1 dichotome DV & 2+ IV

Univariate logistische regressie
Multivariate logistische regressie

\(\begin{equation} \mathbf Z=\beta_0+\beta_1 x + \varepsilon\end{equation}\)
\(\begin{equation} \mathbf Z=\beta_0+\beta_1 x_1 +\beta_2 x_2 + \varepsilon\end{equation}\)


Lineare regressie

Univariate lineare regressie

De formule van een univariate lineare regressie: \(\begin{equation} Y = \beta_0+ \beta_1X_1 + \varepsilon \end{equation}\)

lineare regressie kan met de funtie glm() en de famliy atribute gaussian. Ook kan gebruik worden gemaakt van de functie lm(). De functie lm() heeft de voorkeur.

# Maak een lineare regressie model van 1 DV en 1 IV 
model <- lm(mydata$age ~ mydata$visus)  # formule = Y ~ X. (Y = Uitkomstvariabele)

# Wat bevat dit model
attributes(model)
## $names
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"        
## 
## $class
## [1] "lm"

Lees hier een uitleg van de verschillende atributen en interpretatie ervan.

# Diagnostische functies
coefficients(model)               # Model coefficienten
confint(model, level=0.95)        # Confidence Interval
anova(model)                      # anova tabel
vcov(model)                       # covariantie matrix voor de model parameters
fitted(model)                     # Voorspelde waardes
residuals(model)                  # Residuals

## output wordt niet gegenereerd
# print een statistische samenvatting van het model
summary(model)
## 
## Call:
## lm(formula = mydata$age ~ mydata$visus)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.5655  -4.5921  -0.5655   3.4345  17.3549 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   82.6982     0.6685 123.710  < 2e-16 ***
## mydata$visus -10.2654     1.4990  -6.848 1.88e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.728 on 591 degrees of freedom
## Multiple R-squared:  0.07352,    Adjusted R-squared:  0.07195 
## F-statistic:  46.9 on 1 and 591 DF,  p-value: 1.882e-11

Lees meer over univariate lineare regressie analyse in R. Voor meer informatie over de interpretatie van een univariate regressie analyse zie AMC Wiki.

Multivariate lineare regressie

De formule van een multivariate lineare regressie: \(\begin{equation} Y = \beta_0+ \beta_1X_1 + \beta_2X_2 + ... + \beta_nX_n + \varepsilon \end{equation}\)

# Multiple linear regressie model
model <- lm(mydata$age ~ mydata$visus + mydata$mmse + mydata$apache, data=mydata)
summary(model) # Statistische samenvatting van het model
## 
## Call:
## lm(formula = mydata$age ~ mydata$visus + mydata$mmse + mydata$apache, 
##     data = mydata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.3289  -3.8143  -0.0353   3.2146  16.5762 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   84.81759    1.97839  42.872  < 2e-16 ***
## mydata$visus  -3.13887    1.46098  -2.148   0.0321 *  
## mydata$mmse   -0.45714    0.05656  -8.082 3.62e-15 ***
## mydata$apache  0.49562    0.07511   6.598 9.27e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.119 on 589 degrees of freedom
## Multiple R-squared:  0.2626, Adjusted R-squared:  0.2588 
## F-statistic: 69.92 on 3 and 589 DF,  p-value: < 2.2e-16

Colineariteit analyse Bij gebruik van een multilineare regressie moet rekening worden gehouden met Colineariteit. Of er sprake is van colineariteit kan worden getest met de Variance inflatie factor. Lees meer.

# Variance Inflation Factors
model2 <- lm(mydata$age ~ mydata$visus + mydata$mmse + mydata$apache, data=mydata)

library(car)
vif(model2)                   # variance inflation factors
##  mydata$visus   mydata$mmse mydata$apache 
##      1.189411      1.251456      1.182175
sqrt(vif(model2)) > 2         # Grofweg is er sprake van colineariteit bij een wortel(vif) > 2.
##  mydata$visus   mydata$mmse mydata$apache 
##         FALSE         FALSE         FALSE

Lees meer over Variance inflation factor.


Logistische regressie

Een logistische regressie analyse is te formuleren als: \(Z = \beta_1X_1 + \beta_2X_2 + ... + \beta_kX_k + \varepsilon\) met als wiskundige formule: \(Y = \int (Z) = \frac{1}{1+e^{-Z}}\)

Univariate logistische regressie

# Logistische regressie met glm()
model <- glm(delier ~ age, data=mydata, family=binomial)  
# Wat bevat dit model
attributes(model)
## $names
##  [1] "coefficients"      "residuals"         "fitted.values"    
##  [4] "effects"           "R"                 "rank"             
##  [7] "qr"                "family"            "linear.predictors"
## [10] "deviance"          "aic"               "null.deviance"    
## [13] "iter"              "weights"           "prior.weights"    
## [16] "df.residual"       "df.null"           "y"                
## [19] "converged"         "boundary"          "model"            
## [22] "call"              "formula"           "terms"            
## [25] "data"              "offset"            "control"          
## [28] "method"            "contrasts"         "xlevels"          
## 
## $class
## [1] "glm" "lm"
# print een statistische samenvatting van het model
summary(model)
## 
## Call:
## glm(formula = delier ~ age, family = binomial, data = mydata)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5148   0.3125   0.3971   0.5324   1.2721  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 11.89918    1.70086   6.996 2.63e-12 ***
## age         -0.12367    0.02077  -5.955 2.60e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 434.54  on 592  degrees of freedom
## Residual deviance: 397.53  on 591  degrees of freedom
## AIC: 401.53
## 
## Number of Fisher Scoring iterations: 5

Multivariate logistische regressie

# multivariate logistisch regressie model
model <- glm(delier ~ age + mmse + medicijn, data=mydata, family=binomial)  
# print een statistische samenvatting van het model
summary(model)
## 
## Call:
## glm(formula = delier ~ age + mmse + medicijn, family = binomial, 
##     data = mydata)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8398   0.2421   0.3764   0.4915   2.0785  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      4.87536    2.29320   2.126  0.03350 *  
## age             -0.06359    0.02367  -2.687  0.00721 ** 
## mmse             0.12393    0.03076   4.029 5.61e-05 ***
## medicijnhaldol  -0.93519    0.48030  -1.947  0.05152 .  
## medicijnplacebo -0.96382    0.48193  -2.000  0.04551 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 434.54  on 592  degrees of freedom
## Residual deviance: 372.71  on 588  degrees of freedom
## AIC: 382.71
## 
## Number of Fisher Scoring iterations: 6

Visuele inspectie

De makkelijkste manier van een visuele inspectie is door een scatterplot te maken inclusief de regressielijn. Dit is echter alleen mogelijk voor een univariate lineare regressie.

# Maak een scatterplot van twee variabelen (Y~X)
plot(mydata$mmse~mydata$visus)        

# Univaraite linear regressie model
model <- lm(mydata$mmse ~ mydata$visus) # formule = Y ~ X 

# Voeg de coefficiënt lijn toe (slope)
abline(model, col="red")

De plot() functie kan een regressie model visualiseren. Deze functie genereert 4 grafieken en is geschikt voor zowel univariate als multivariate regressie. Lees meer over de betekenis van deze grafieken.

# Univariate linear regressie model
model <- lm(mydata$mmse ~ mydata$visus) # formule = Y ~ X 

# Visuele inspectie met de functie plot()
layout(matrix(1:4,2,2))
plot(model)                  # ook geschikt voor multivariate lineare regressie

Lees meer over de interpretatie van de plot(lm() ) functie.

Een andere methode die ook geschikt is voor multivariate regressie is met een residuals vs components plot per onafhankelijke variabele (X). Dit kan met de functie crPlots() uit de package car.

# multivariate linear regressie model
multimodel <- lm(mydata$mmse ~ mydata$visus+mydata$age+mydata$apache+mydata$creaurer) # formule = Y ~ X1+X2+X3+X4. 

# component + residual plot per onafhankelijke variabele (X)
library(car)
crPlots(multimodel)

Validatie

Sensitiviteit en Specificiteit

Een Confusionmatrix is een frequentietabel met 2 rijen en 2 kollomen: false positive, false negative, true postive en true negative. Deze wordt gebruikt voor het testen van de kwaliteit van een logistisch voorspelmodel. Daarnaast kan uit een confusionmatrix de sensitiviteit en specificiteit worden berekend.

confusion

## Een voorbeeld van een confusiontable op basis van een logistisch predictiemodel inclusief sensitiviteit en specificiteit
# maak een GLM
glm <- glm(delier~apache+mmse+visus+creaurer, data=mydata, family="binomial")

# voorspel GLM
prop <- predict(glm, mydata, type="response")

# maak er 0/1 predictie van
# bij een voorspelde kans >50% op delier; voorspel delier.
prediction <- ifelse (prop >= 0.5, 1, 0)   

# maak een confusionmatrix
mytable <- table(prediction, mydata$delier)
print(mytable)
##           
## prediction delier geen delier
##          0      4          10
##          1     67         512
# Sensitiviteit
SE <- mytable[2,1] / (mytable[2,1] + mytable[1,1])
SE*100   # sensitiviteit in procenten
## [1] 94.3662
# Specificiteit
SP <- mytable[1,2] / (mytable[1,2] + mytable[2,2])
SP*100  # specificiteit in procenten
## [1] 1.915709
# Area Under the Curve (descriminatie-index)
library(caTools)
colAUC(prediction, mydata$delier)
##                             [,1]
## delier vs. geen delier 0.5185905

Visuele Validatie

De visuele validatie van een predictiemodel (werkt dus ook op andere predictie technieken dan alleen regressie-analyses) kan met de package ‘ROCR’.

Een ROC-curve is een grafiek van de gevoeligheid (sensitiviteit) in functie van de aspecificiteit (1 - specificiteit) voor een binaire classifier als zijn discriminatiedrempel wordt gevarieerd. De ROC kan ook worden weergegeven door de fractie van true positives (TPR = true positive rate) uit te zetten tegen de fractie van fout-positieven (FPR = false positive rate). Lees meer.

# maak een GLM
glm <- glm(delier~apache+mmse+visus+creaurer, data=mydata, family="binomial")

# voorspel GLM
predict <- predict(glm, mydata, type="response")

# laad package ROCR
library(ROCR)

# voorspel delier
pred <- prediction(predict, mydata$delier)

# 2 bij 2 multigrafiek
par(mfrow=c(2,2))


# maak een ROC-curve
perf <- performance(pred, measure = "tpr", x.measure = "fpr") 
plot(perf, avg= "threshold", colorize=T, lwd= 3,
     coloraxis.at=seq(0,1,by=0.2),
     main= "ROC curve")

# plot de accuracy
perf <- performance(pred, "acc")
plot(perf, avg= "vertical", spread.estimate="boxplot", lwd=3,col='blue',
     show.spread.at= seq(0.1, 0.9, by=0.1),
     main= "Accuracy across the range of possible cutoffs")

# plot de calibratie
plot(performance(pred, "cal", window.size= 10),
     avg="vertical",
     main= "How well are the probability predictions calibrated?")

# plot de predicties van de dichotome uitkomstvariabele los van elkaar
plot(0,0,type="n", xlim= c(0,1), ylim=c(0,15),
     xlab="Cutoff", ylab="Density",
     main="How well do the predictions separate the classes?")
for (runi in 1:length(pred@predictions)) {
    lines(density(pred@predictions[[runi]][pred@labels[[runi]]=="delier"]), col= "red")
    lines(density(pred@predictions[[runi]][pred@labels[[runi]]=="geen delier"]), col="green")
}

Lees meer over ROC-curves. Voor meer informatie over de package ‘ROCR’ zie deze Blog.


Survival analyse

Voor het verrichten van een survival analyse moeten eerst een aantal nieuwe variabelen worden gemaakt uit de data.

## Data management
# begin met een nieuwe data.frame met de geselecteerde variabelen
# kies de variabelen om te analyseren
myvars <- c("opndat", "datdood", "dood", "medicijn","age","visus")
data <- mydata[c(myvars)]

# Maak EERST nieuwe variabele met tijd tot doodgaan
data$overleving <- as.numeric(data$datdood - data$opndat)
# verander DAARNA alle NA's in een maximale tijd op basis van de maximale tijd die al voorkomt
library(plyr)
time <- mapvalues(data$overleving, from = c(NA), to = max(na.omit(data$overleving)) + 1 )

# Maak een nieuwe event variabele (logical 0/1) op basis van variabele data$dood
event <- as.logical(revalue(data$dood, c("overleden"=TRUE, "niet overleden"=FALSE)))

# Voeg de nieuwe variabelen toe aan het data.frame
data$event = event
data$time = time

Kaplan-Meier

Voor het uitvoeren van een survivalanalyse is de package survival geschikt. Er zijn 3 functies van belang voor het uitvoeren van de analyse: Surv(), survreg() en survfit(). Zie ‘?survival’ voor meer informatie. Voor het visualiseren van een survivalplot wordt is combinatie plot(survfit(Surv(X,Y)~1)) geschikt. Lees meer over survival curves met de ggsurv functie.

# Gebruik de package survival voor een kaplan meier analyse
library(survival)

# maak een survival-object
data$SurvObj <- with(data, Surv(time, event == 1))

# Kaplan-Meier estimator met "log-log" confidence interval.
km.as.one <- survfit(SurvObj ~ 1, data = data)
km.as.one
## Call: survfit(formula = SurvObj ~ 1, data = data)
## 
## records   n.max n.start  events  median 0.95LCL 0.95UCL 
##     593     593     593     192      NA      NA      NA
# Kaplan-Meier estimator per groep (medicijn) met "log-log" confidence interval
km.by.med <- survfit(SurvObj ~ medicijn, data = data)
km.by.med
## Call: survfit(formula = SurvObj ~ medicijn, data = data)
## 
##                    records n.max n.start events median 0.95LCL 0.95UCL
## medicijn=geen med.     172   172     172     30     NA      NA      NA
## medicijn=haldol        208   208     208     83     NA      NA      NA
## medicijn=placebo       213   213     213     79     NA      NA      NA
# verricht de analyse met survreg() en gebruik summary() voor de resultaten
kmr.as.one <- survreg( Surv(data$time, data$event)~1, data=data) # ~1 = niet vergelijken van groepen
summary(kmr.as.one)
## 
## Call:
## survreg(formula = Surv(data$time, data$event) ~ 1, data = data)
##             Value Std. Error     z        p
## (Intercept) 9.047     0.1425 63.51 0.00e+00
## Log(scale)  0.334     0.0693  4.82 1.42e-06
## 
## Scale= 1.4 
## 
## Weibull distribution
## Loglik(model)= -1837.7   Loglik(intercept only)= -1837.7
## Number of Newton-Raphson Iterations: 5 
## n= 593
# Verricht de analyse van survivaldata per groep
kmr.by.med <- survreg( Surv(data$time, data$event)~medicijn, data=data) # ~medicijn = vergelijken van groepen
summary(kmr.by.med)
## 
## Call:
## survreg(formula = Surv(data$time, data$event) ~ medicijn, data = data)
##                  Value Std. Error     z         p
## (Intercept)     10.039     0.3005 33.40 1.22e-244
## medicijnhaldol  -1.394     0.3070 -4.54  5.60e-06
## medicijnplacebo -1.239     0.3065 -4.04  5.30e-05
## Log(scale)       0.322     0.0688  4.68  2.81e-06
## 
## Scale= 1.38 
## 
## Weibull distribution
## Loglik(model)= -1823.4   Loglik(intercept only)= -1837.7
##  Chisq= 28.71 on 2 degrees of freedom, p= 5.8e-07 
## Number of Newton-Raphson Iterations: 5 
## n= 593
# Kaplan Meier plot
survival <- survfit( Surv(data$time, data$event)~1, data=data) # ~1 = niet vergelijken van groepen
plot(survival)

# Kaplan Meier plot per groep
survival <- survfit( Surv(data$time, data$event)~medicijn, data=data) # ~X = vergelijken van groepen
plot(survival, 
     col = c(4,2,3),              # volgorde van de kleuren = level volgorde van data$medicijn (blauw,rood,groen)
    xlim = c(0,2000),             # x-as van 0 tot 2000
    ylim = c(0.6,1.0),            # y-as van 1.0 tot 0.6
    xlab = "Tijd in dagen",       # x-as label
    ylab = "overlevingskans",     # y-as label
    main = "Kaplan meier plot",   # grafiek titel
    sub  = "subtitel",            # per medicatie
    )
legend("bottomleft", inset=.05, title="legenda", c(levels(data$medicijn)), fill = c(4,2,3), cex=0.8)

Cox-hazard ratio

## Gebruik de package survival
library(survival)

## Verricht de analyse met coxph() en gebruik summary() voor de resultaten
cox.reg <- coxph( Surv(data$time, data$event)~1, data=data) # ~1 = niet vergelijken van groepen
summary(cox.reg)
## Call:  coxph(formula = Surv(data$time, data$event) ~ 1, data = data)
## 
## Null model
##   log likelihood= -1191.037 
##   n= 593
## Verricht de analyse met coxph() en gebruik summary() voor de resultaten
cox.by.covariates <- coxph( Surv(data$time, data$event)~age+visus, data=data) # ~X = vergelijken van groepen
summary(cox.by.covariates)
## Call:
## coxph(formula = Surv(data$time, data$event) ~ age + visus, data = data)
## 
##   n= 593, number of events= 192 
## 
##           coef exp(coef) se(coef)      z Pr(>|z|)    
## age    0.09152   1.09584  0.01195  7.658 1.89e-14 ***
## visus -2.24947   0.10545  0.49437 -4.550 5.36e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## age      1.0958     0.9125   1.07047    1.1218
## visus    0.1055     9.4827   0.04002    0.2779
## 
## Concordance= 0.693  (se = 0.021 )
## Rsquare= 0.154   (max possible= 0.982 )
## Likelihood ratio test= 99.17  on 2 df,   p=0
## Wald test            = 105.9  on 2 df,   p=0
## Score (logrank) test = 111.9  on 2 df,   p=0
## Check voor overtreding van proportionele hazard (constante Hazard Ratio over tijd)
res.zph1 <- cox.zph(cox.by.covariates)
attributes(res.zph1)
## $names
## [1] "table"     "x"         "y"         "var"       "call"      "transform"
## 
## $class
## [1] "cox.zph"
res.zph1
##            rho chisq       p
## age    -0.0895  1.63 0.20106
## visus   0.1877  7.13 0.00760
## GLOBAL      NA 11.44 0.00328
## Plot de 'scaled Schoenfeld residuals' inclusief een 'smooth curve'.
par(mfrow=c(1,2))
plot(res.zph1)

Lees meer over Cox regressie.


Creative Commons-Licentie Dit werk valt onder een Creative Commons Naamsvermelding - Gelijk Delen 3.0 Internationaal-licentie .