framingham_initial-analysis
Shimon
2024-02-27
#library
pacman::p_load(
tidyverse,data.table,tableone,tibble,stringi,readxl,dplyr,scales,ggsci,
skimr,Epi,janitor,summarytools,broom,easystats,car,rms)
#conflict解消
select<-dplyr::select
library(pROC);library(glmnetUtils);library(glmnet);library(ISLR)## Type 'citation("pROC")' for a citation.##
## 次のパッケージを付け加えます: 'pROC'## 以下のオブジェクトは 'package:parameters' からマスクされています:
##
## ci## 以下のオブジェクトは 'package:bayestestR' からマスクされています:
##
## auc, ci## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## cov, smooth, var## 要求されたパッケージ Matrix をロード中です##
## 次のパッケージを付け加えます: 'Matrix'## 以下のオブジェクトは 'package:tidyr' からマスクされています:
##
## expand, pack, unpack## Loaded glmnet 4.1-8##
## 次のパッケージを付け加えます: 'glmnet'## 以下のオブジェクトは 'package:glmnetUtils' からマスクされています:
##
## cv.glmnet, glmnet#factor var
fac_var<-c("RANDID","SEX", "CURSMOKE", "DIABETES", "BPMEDS","educ","PREVCHD", "PREVAP", "PREVMI", "PREVSTRK", "PREVHYP", "PERIOD", "DEATH",
"ANGINA", "HOSPMI", "MI_FCHD", "ANYCHD", "STROKE", "CVD", "HYPERTEN")
cont_var<-c("TOTCHOL", "AGE", "SYSBP", "DIABP", "CIGPDAY", "BMI", "HEARTRTE", "GLUCOSE")
full_var<-c(fac_var[-c(1,6,12:20)],cont_var)
var<-c("SEX", "CURSMOKE", "DIABETES", "BPMEDS", "PREVCHD", "PREVAP", "PREVMI", "PREVHYP", cont_var) #PREVSTRは除外している
inter_var<-c("SEX:PREVCHD", "SEX:PREVAP", "SEX:PREVMI", "SEX:PREVHYP", "SEX:DIABP", "CURSMOKE:HEARTRTE","DIABETES:TOTCHOL", "DIABETES:SYSBP",
"BPMEDS:PREVCHD", "BPMEDS:TOTCHOL", "BPMEDS:AGE", "PREVCHD:GLUCOSE", "PREVHYP:BMI", "TOTCHOL:AGE","HEARTRTE:GLUCOSE")
nonlin_term<-c("BMI","HEARTRTE","GLUCOSE")
#data import
fr<-read_csv("frmgham2.csv") %>%
mutate(across(fac_var,as.factor))## Rows: 11627 Columns: 39
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (39): RANDID, SEX, TOTCHOL, AGE, SYSBP, DIABP, CURSMOKE, CIGPDAY, BMI, D...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(fac_var, as.factor)`.
## Caused by warning:
## ! Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(fac_var)
##
## # Now:
## data %>% select(all_of(fac_var))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.#欠測補完後データのimport
load("fr.missf.RData")
fr.missf<-fr.missf$ximp
#plot in baseline variable
fr_p1<-fr %>%
filter(PERIOD==1)
fr.missf_p1<-
fr.missf %>%
filter(PERIOD==1)full model formula
# full model formula
fullmodel<-
paste0("STROKE~",
str_c(var[!var %in% nonlin_term],collapse = "+"), #非線形以外のvarを＋でつなぐ
"+",
str_c(paste0("rcs(",nonlin_term,",4)"),collapse = "+"), #非線形のvarを＋でつなぐ
"+",
str_c(inter_var,collapse = "+")) #交互作用項を＋でつなぐfull model fitting results
# Complete case
full_results<-
glm(as.formula(fullmodel),data=fr_p1,
family = binomial(link = "logit")) %>%
summary() %>%
.$coefficients %>%
as.data.frame() %>%
rownames_to_column() %>%
rename("var"="rowname") %>%
mutate(OR=exp(Estimate),
Lower=exp(Estimate-1.96*`Std. Error`),
Upper=exp(Estimate+1.96*`Std. Error`)) %>%
select(var,OR,Lower,Upper,"p"="Pr(>|z|)") %>%
#すべての変数を丸める
mutate_if(is.numeric,~round(.,3))
full_results## var OR Lower Upper p
## 1 (Intercept) 0.031 0.000 21.641 0.298
## 2 SEX2 0.346 0.060 2.005 0.237
## 3 CURSMOKE1 0.226 0.050 1.013 0.052
## 4 DIABETES1 52.994 1.393 2015.673 0.032
## 5 BPMEDS1 127.580 2.055 7920.611 0.021
## 6 PREVCHD1 3.288 0.288 37.541 0.338
## 7 PREVAP1 0.637 0.163 2.487 0.516
## 8 PREVMI1 1.392 0.394 4.915 0.607
## 9 PREVHYP1 5.301 1.053 26.683 0.043
## 10 TOTCHOL 1.002 0.984 1.020 0.825
## 11 AGE 1.085 1.003 1.174 0.043
## 12 SYSBP 1.001 0.992 1.010 0.816
## 13 DIABP 1.012 0.993 1.031 0.204
## 14 CIGPDAY 1.003 0.988 1.018 0.677
## 15 rcs(BMI, 4)BMI 0.866 0.755 0.994 0.040
## 16 rcs(BMI, 4)BMI' 1.469 0.940 2.298 0.091
## 17 rcs(BMI, 4)BMI'' 0.361 0.089 1.463 0.154
## 18 rcs(HEARTRTE, 4)HEARTRTE 1.013 0.967 1.062 0.585
## 19 rcs(HEARTRTE, 4)HEARTRTE' 0.787 0.658 0.940 0.008
## 20 rcs(HEARTRTE, 4)HEARTRTE'' 1.856 1.179 2.923 0.008
## 21 rcs(GLUCOSE, 4)GLUCOSE 0.947 0.904 0.991 0.020
## 22 rcs(GLUCOSE, 4)GLUCOSE' 1.102 0.934 1.300 0.249
## 23 rcs(GLUCOSE, 4)GLUCOSE'' 0.768 0.492 1.198 0.244
## 24 SEX2:PREVCHD1 0.720 0.020 26.534 0.858
## 25 SEX2:PREVAP1 2.493 0.082 75.953 0.600
## 26 SEX2:PREVMI1 2.786 0.106 73.500 0.539
## 27 SEX2:PREVHYP1 1.048 0.589 1.863 0.874
## 28 SEX2:DIABP 1.011 0.989 1.033 0.329
## 29 CURSMOKE0:HEARTRTE 0.976 0.957 0.996 0.017
## 30 DIABETES1:TOTCHOL 0.992 0.981 1.003 0.144
## 31 DIABETES1:SYSBP 0.992 0.975 1.011 0.413
## 32 BPMEDS1:PREVCHD1 0.463 0.105 2.038 0.309
## 33 BPMEDS1:TOTCHOL 0.993 0.983 1.003 0.153
## 34 BPMEDS1:AGE 0.957 0.904 1.013 0.133
## 35 PREVCHD0:GLUCOSE 1.013 0.991 1.036 0.249
## 36 PREVHYP0:BMI 1.044 0.985 1.107 0.144
## 37 TOTCHOL:AGE 1.000 1.000 1.000 0.816
## 38 HEARTRTE:GLUCOSE 1.000 1.000 1.001 0.069#imputation data
full_results_mi<-
glm(as.formula(fullmodel),data=fr.missf_p1,
family = binomial(link = "logit")) %>%
summary() %>%
.$coefficients %>%
as.data.frame() %>%
rownames_to_column() %>%
rename("var"="rowname") %>%
mutate(OR=exp(Estimate),
Lower=exp(Estimate-1.96*`Std. Error`),
Upper=exp(Estimate+1.96*`Std. Error`)) %>%
select(var,OR,Lower,Upper,"p"="Pr(>|z|)") %>%
#すべての変数を丸める
mutate_if(is.numeric,~round(.,3))
full_results_mi## var OR Lower Upper p
## 1 (Intercept) 0.005 0.000 2.849 0.102
## 2 SEX2 0.303 0.058 1.583 0.157
## 3 CURSMOKE1 0.290 0.071 1.188 0.085
## 4 DIABETES1 55.743 1.632 1903.573 0.026
## 5 BPMEDS1 198.320 4.319 9106.958 0.007
## 6 PREVCHD1 3.015 0.328 27.704 0.329
## 7 PREVAP1 0.689 0.186 2.558 0.578
## 8 PREVMI1 1.057 0.331 3.368 0.926
## 9 PREVHYP1 4.549 0.974 21.246 0.054
## 10 TOTCHOL 1.007 0.990 1.024 0.447
## 11 AGE 1.106 1.027 1.192 0.008
## 12 SYSBP 1.001 0.993 1.010 0.748
## 13 DIABP 1.013 0.995 1.031 0.159
## 14 CIGPDAY 1.003 0.988 1.017 0.718
## 15 rcs(BMI, 4)BMI 0.917 0.803 1.048 0.204
## 16 rcs(BMI, 4)BMI' 1.256 0.817 1.930 0.299
## 17 rcs(BMI, 4)BMI'' 0.552 0.144 2.112 0.385
## 18 rcs(HEARTRTE, 4)HEARTRTE 1.010 0.965 1.056 0.681
## 19 rcs(HEARTRTE, 4)HEARTRTE' 0.809 0.683 0.958 0.014
## 20 rcs(HEARTRTE, 4)HEARTRTE'' 1.721 1.120 2.645 0.013
## 21 rcs(GLUCOSE, 4)GLUCOSE 0.941 0.901 0.982 0.006
## 22 rcs(GLUCOSE, 4)GLUCOSE' 1.134 0.990 1.300 0.070
## 23 rcs(GLUCOSE, 4)GLUCOSE'' 0.672 0.439 1.028 0.067
## 24 SEX2:PREVCHD1 0.638 0.021 19.464 0.797
## 25 SEX2:PREVAP1 2.654 0.100 70.232 0.559
## 26 SEX2:PREVMI1 3.705 0.165 83.060 0.409
## 27 SEX2:PREVHYP1 1.073 0.623 1.845 0.800
## 28 SEX2:DIABP 1.012 0.992 1.033 0.240
## 29 CURSMOKE0:HEARTRTE 0.979 0.962 0.998 0.026
## 30 DIABETES1:TOTCHOL 0.993 0.983 1.004 0.217
## 31 DIABETES1:SYSBP 0.989 0.972 1.006 0.192
## 32 BPMEDS1:PREVCHD1 0.434 0.108 1.743 0.239
## 33 BPMEDS1:TOTCHOL 0.991 0.982 1.001 0.070
## 34 BPMEDS1:AGE 0.957 0.907 1.009 0.107
## 35 PREVCHD0:GLUCOSE 1.010 0.991 1.030 0.299
## 36 PREVHYP0:BMI 1.038 0.982 1.098 0.187
## 37 TOTCHOL:AGE 1.000 1.000 1.000 0.428
## 38 HEARTRTE:GLUCOSE 1.000 1.000 1.001 0.045full_results_mi %>%
ggplot(aes(y=var,x=OR))+
geom_point()+
geom_errorbar(aes(xmin=Lower,xmax=Upper),width=0.2)+
geom_vline(xintercept=1,color="red")
ORが異常に幅広くなっている因子があるため、sparse data biasがかかっている可能性が高い →penalized regressionを行うか、モデル式の変更をすべき状態
Backwards stepwiseとElastic Netによるパラメータ推定を行う
Elastic Netによるパラメータ推定
下準備
# install.packages("glmnet")
# install.packages("ISLR")
# install.packages("glmnetUtils")
# library(glmnetUtils);library(glmnet);library(ISLR)Ridge regression
#to detect optimal lambda
ridge.cv.res<-
cv.glmnet(
x=fr.missf_p1[,var] %>% data.matrix(),
y=fr.missf_p1[,"STROKE"] %>% data.matrix(),
family="binomial",
alpha=0) #ridge
plot(ridge.cv.res) 
ridge.cv.res$lambda.min #min lambda## [1] 0.009523244# model coefficient with min lambda
glmnet(
x=fr.missf_p1[,var] %>% data.matrix(),
y=fr.missf_p1[,"STROKE"] %>% data.matrix(),
family="binomial",
lambda=ridge.cv.res$lambda.min,
alpha=0) %>%
.$beta## 16 x 1 sparse Matrix of class "dgCMatrix"
## s0
## SEX -0.0813962421
## CURSMOKE 0.2468760397
## DIABETES 0.7155327927
## BPMEDS 0.5592037443
## PREVCHD 0.3002169372
## PREVAP -0.2447784453
## PREVMI 0.0955726098
## PREVHYP 0.5155415106
## TOTCHOL -0.0011747152
## AGE 0.0604815428
## SYSBP 0.0041215475
## DIABP 0.0149010624
## CIGPDAY 0.0040394071
## BMI 0.0159035694
## HEARTRTE -0.0102841712
## GLUCOSE -0.0009616904Lasso regression
#to detect optimal lambda
lasso.cv.res<-
cv.glmnet(
x=fr.missf_p1[,var] %>% data.matrix(),
y=fr.missf_p1[,"STROKE"] %>% data.matrix(),
family="binomial",
alpha=1) #lasso
#log lambdaとdevianceの関連性plot
plot(lasso.cv.res)
#devianceが最小のときのlambda
lasso.cv.res$lambda.min #min lambda## [1] 0.0021#min lambdaのもとでのlasso regression coefficient；二通りのcodeがあるが、結果は同じ
coef(lasso.cv.res,s="lambda.min") #interceptの値も表示される## 17 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) -9.3732919829
## SEX -0.0405473179
## CURSMOKE 0.2688147316
## DIABETES 0.6136939222
## BPMEDS 0.5142181504
## PREVCHD 0.0298690706
## PREVAP .
## PREVMI 0.1240355798
## PREVHYP 0.5511057332
## TOTCHOL -0.0008671472
## AGE 0.0671689119
## SYSBP 0.0012846234
## DIABP 0.0178463734
## CIGPDAY 0.0025925936
## BMI 0.0112364738
## HEARTRTE -0.0092393943
## GLUCOSE .glmnet(
x=fr.missf_p1[,var] %>% data.matrix(),
y=fr.missf_p1[,"STROKE"] %>% data.matrix(),
family="binomial",
lambda=lasso.cv.res$lambda.min,
alpha=1) %>%
.$beta ## 16 x 1 sparse Matrix of class "dgCMatrix"
## s0
## SEX -0.040476551
## CURSMOKE 0.268983024
## DIABETES 0.613718984
## BPMEDS 0.514118152
## PREVCHD 0.029927580
## PREVAP .
## PREVMI 0.123893960
## PREVHYP 0.551330129
## TOTCHOL -0.000867187
## AGE 0.067179689
## SYSBP 0.001277099
## DIABP 0.017852439
## CIGPDAY 0.002589202
## BMI 0.011237120
## HEARTRTE -0.009239599
## GLUCOSE .Elastic Net
alphaを定める
Elastic Netでは、Ridge,Lassoの罰則項を割合alphaで混ぜ合わせる事ができる 当てはまりが最も良い時の最適なalpha求める
#cv.glmnetではなく、cv"a".glmnetを使う; to examine optimal alpha and lambda
elastic.cv.res<-
cva.glmnet(
x=fr.missf_p1[,var] %>% data.matrix(),
y=fr.missf_p1[,"STROKE"] %>% data.matrix(),
family="binomial"
)
plot(elastic.cv.res)
#alphaをデフォルトでは0-1の間で11個に振り分け、それぞれでのdevianceとlambdaの関係性を求める
#各alphaのうち、最もDevianceが小さいときのalphaを探し、その中で最も小さいlambdaを探す
##alpha=1はlassoと同じ；alpha=0はRidgeと同じということになる
sapply(elastic.cv.res$modlist,"[[","cvm") %>% #各alphaにおけるcvmを取り出す
sapply(min) %>% #各alphaの中で最小の値を取り出す
min() #すべてのalphaの中で最小のcvmを取り出す## [1] 0.55716#結果を踏まえると、alpha=1のとき最小の値を示しており、すなわちlassoが最適と考えられる
#alpha=1のときの最小のlamdaを求める
cv.glmnet(
x=fr.missf_p1[,var] %>% data.matrix(),
y=fr.missf_p1[,"STROKE"] %>% data.matrix(),
family="binomial",
alpha=1
) %>%
.$lambda.min## [1] 0.001588573 #これは一個前のchankのlamda minと同じ値になる
attributes(elastic.cv.res)## $names
## [1] "alpha" "nfolds" "modlist" "call"
##
## $class
## [1] "cva.glmnet"#minlossplot(elastic.cv.res)最適な(α,λ)は(1,0.001913442)であり、lasso regressionが良いということになる →これは間違っていそう。binomial devianceが最小のときのα、λのセットを選べていないかも →
ということで、更に再検討 Grid search ref: https://asmquantmacro.com/2016/04/26/fitting-elastic-net-model-in-r/
#install.packages("caret")
library(caret)## 要求されたパッケージ lattice をロード中です##
## 次のパッケージを付け加えます: 'caret'## 以下のオブジェクトは 'package:parameters' からマスクされています:
##
## compare_models## 以下のオブジェクトは 'package:purrr' からマスクされています:
##
## liftGrid Searchの準備
#develop grid search frame work
lambda.grid<-10^seq(2, -2, length=100)
alpha.grid<-seq(0, 1, len = 11)^3 #alphaの分割回数(cva.glmnetと同じに合わせた)
srchGrd<-expand.grid(alpha=alpha.grid,lambda=lambda.grid) #alpha 11個、lambda 100個の組み合わせを作る;計1100行
#Grid search with Cross validation
trnCntrl<-
trainControl(method = "repeatedcv", #bootstrapも可能
number=10, #fold数
repeats=5,
search="grid") #random searchも可能
#perform Cross validation
set.seed(2024)
grid.res<-
train(x=fr.missf_p1[,var] %>% data.matrix(),
y=fr.missf_p1[,"STROKE"] %>% data.matrix(),
data=fr.missf_p1,
method="glmnet", #method list is in http://topepo.github.io/caret/train-models-by-tag.html
#metric="RMSE",#value to determine optimal parameter
tuneGrid=srchGrd, #grid search framework
trControl=trnCntrl #searching method condition.
)
plot(grid.res)
#得られた結果から任意の項目を抽出する項目名の指定
attributes(grid.res) ## $names
## [1] "method" "modelInfo" "modelType" "results" "pred"
## [6] "bestTune" "call" "dots" "metric" "control"
## [11] "finalModel" "preProcess" "trainingData" "ptype" "resample"
## [16] "resampledCM" "perfNames" "maximize" "yLimits" "times"
## [21] "levels"
##
## $class
## [1] "train"#accuracy最適なときのalpha, lambda
grid.res$bestTune## alpha lambda
## 28 0 0.1232847grid.res$results## alpha lambda Accuracy Kappa AccuracySD KappaSD
## 1 0.000 0.01000000 0.9061354 1.088097e-02 0.002272464 0.0231936150
## 2 0.000 0.01097499 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 3 0.000 0.01204504 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 4 0.000 0.01321941 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 5 0.000 0.01450829 0.9062256 1.181208e-02 0.002236353 0.0234894941
## 6 0.000 0.01592283 0.9062256 1.181208e-02 0.002236353 0.0234894941
## 7 0.000 0.01747528 0.9062256 1.181208e-02 0.002236353 0.0234894941
## 8 0.000 0.01917910 0.9062707 1.189718e-02 0.002147844 0.0234213703
## 9 0.000 0.02104904 0.9062255 1.034582e-02 0.002146289 0.0216204550
## 10 0.000 0.02310130 0.9062255 1.034582e-02 0.002146289 0.0216204550
## 11 0.000 0.02535364 0.9062708 9.670396e-03 0.001992461 0.0210364255
## 12 0.000 0.02782559 0.9063158 8.994152e-03 0.001945444 0.0206199572
## 13 0.000 0.03053856 0.9063158 8.236576e-03 0.001891483 0.0200397462
## 14 0.000 0.03351603 0.9063158 8.236576e-03 0.001891483 0.0200397462
## 15 0.000 0.03678380 0.9063158 8.236576e-03 0.001891483 0.0200397462
## 16 0.000 0.04037017 0.9063158 7.453662e-03 0.001778129 0.0193530744
## 17 0.000 0.04430621 0.9063158 7.453662e-03 0.001778129 0.0193530744
## 18 0.000 0.04862602 0.9063158 6.674133e-03 0.001718651 0.0186261584
## 19 0.000 0.05336699 0.9064060 6.851051e-03 0.001660926 0.0185389761
## 20 0.000 0.05857021 0.9064962 7.038989e-03 0.001595930 0.0186518966
## 21 0.000 0.06428073 0.9064512 6.193030e-03 0.001594263 0.0179667139
## 22 0.000 0.07054802 0.9064963 6.281564e-03 0.001589830 0.0179241718
## 23 0.000 0.07742637 0.9064512 5.467220e-03 0.001527643 0.0150640003
## 24 0.000 0.08497534 0.9064513 3.923264e-03 0.001451245 0.0131545349
## 25 0.000 0.09326033 0.9064513 3.923264e-03 0.001451245 0.0131545349
## 26 0.000 0.10235310 0.9065415 4.096896e-03 0.001295782 0.0130416236
## 27 0.000 0.11233240 0.9065415 3.339320e-03 0.001213262 0.0117864848
## 28 0.000 0.12328467 0.9065866 3.427854e-03 0.001203981 0.0117435073
## 29 0.000 0.13530478 0.9065415 2.559791e-03 0.001124287 0.0102356292
## 30 0.000 0.14849683 0.9064963 1.714023e-03 0.001132386 0.0084829470
## 31 0.000 0.16297508 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 32 0.000 0.17886495 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 33 0.000 0.19630407 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 34 0.000 0.21544347 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 35 0.000 0.23644894 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 36 0.000 0.25950242 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 37 0.000 0.28480359 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 38 0.000 0.31257158 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 39 0.000 0.34304693 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 40 0.000 0.37649358 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 41 0.000 0.41320124 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 42 0.000 0.45348785 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 43 0.000 0.49770236 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 44 0.000 0.54622772 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 45 0.000 0.59948425 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 46 0.000 0.65793322 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 47 0.000 0.72208090 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 48 0.000 0.79248290 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 49 0.000 0.86974900 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 50 0.000 0.95454846 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 51 0.000 1.04761575 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 52 0.000 1.14975700 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 53 0.000 1.26185688 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 54 0.000 1.38488637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 55 0.000 1.51991108 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 56 0.000 1.66810054 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 57 0.000 1.83073828 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 58 0.000 2.00923300 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 59 0.000 2.20513074 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 60 0.000 2.42012826 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 61 0.000 2.65608778 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 62 0.000 2.91505306 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 63 0.000 3.19926714 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 64 0.000 3.51119173 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 65 0.000 3.85352859 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 66 0.000 4.22924287 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 67 0.000 4.64158883 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 68 0.000 5.09413801 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 69 0.000 5.59081018 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 70 0.000 6.13590727 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 71 0.000 6.73415066 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 72 0.000 7.39072203 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 73 0.000 8.11130831 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 74 0.000 8.90215085 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 75 0.000 9.77009957 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 76 0.000 10.72267222 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 77 0.000 11.76811952 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 78 0.000 12.91549665 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 79 0.000 14.17474163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 80 0.000 15.55676144 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 81 0.000 17.07352647 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 82 0.000 18.73817423 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 83 0.000 20.56512308 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 84 0.000 22.57019720 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 85 0.000 24.77076356 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 86 0.000 27.18588243 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 87 0.000 29.83647240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 88 0.000 32.74549163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 89 0.000 35.93813664 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 90 0.000 39.44206059 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 91 0.000 43.28761281 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 92 0.000 47.50810162 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 93 0.000 52.14008288 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 94 0.000 57.22367659 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 95 0.000 62.80291442 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 96 0.000 68.92612104 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 97 0.000 75.64633276 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 98 0.000 83.02175681 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 99 0.000 91.11627561 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 100 0.000 100.00000000 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 101 0.001 0.01000000 0.9061354 1.088097e-02 0.002272464 0.0231936150
## 102 0.001 0.01097499 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 103 0.001 0.01204504 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 104 0.001 0.01321941 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 105 0.001 0.01450829 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 106 0.001 0.01592283 0.9062256 1.181208e-02 0.002236353 0.0234894941
## 107 0.001 0.01747528 0.9062256 1.181208e-02 0.002236353 0.0234894941
## 108 0.001 0.01917910 0.9062256 1.105122e-02 0.002141771 0.0230620373
## 109 0.001 0.02104904 0.9062255 1.034582e-02 0.002146289 0.0216204550
## 110 0.001 0.02310130 0.9062255 1.034582e-02 0.002146289 0.0216204550
## 111 0.001 0.02535364 0.9062708 9.670396e-03 0.001992461 0.0210364255
## 112 0.001 0.02782559 0.9062708 8.148193e-03 0.001939808 0.0200865050
## 113 0.001 0.03053856 0.9063158 8.236576e-03 0.001891483 0.0200397462
## 114 0.001 0.03351603 0.9063158 8.236576e-03 0.001891483 0.0200397462
## 115 0.001 0.03678380 0.9063610 8.321725e-03 0.001835233 0.0199753241
## 116 0.001 0.04037017 0.9063158 7.453662e-03 0.001778129 0.0193530744
## 117 0.001 0.04430621 0.9063158 7.453662e-03 0.001778129 0.0193530744
## 118 0.001 0.04862602 0.9063610 6.762667e-03 0.001718172 0.0185827872
## 119 0.001 0.05336699 0.9064060 6.851051e-03 0.001660926 0.0185389761
## 120 0.001 0.05857021 0.9064962 7.038989e-03 0.001595930 0.0186518966
## 121 0.001 0.06428073 0.9064512 6.193030e-03 0.001594263 0.0179667139
## 122 0.001 0.07054802 0.9064963 6.281564e-03 0.001589830 0.0179241718
## 123 0.001 0.07742637 0.9064060 4.599157e-03 0.001461392 0.0140497368
## 124 0.001 0.08497534 0.9064513 3.923264e-03 0.001451245 0.0131545349
## 125 0.001 0.09326033 0.9064513 3.923264e-03 0.001451245 0.0131545349
## 126 0.001 0.10235310 0.9065415 4.096896e-03 0.001295782 0.0130416236
## 127 0.001 0.11233240 0.9065415 3.339320e-03 0.001213262 0.0117864848
## 128 0.001 0.12328467 0.9065866 3.427854e-03 0.001203981 0.0117435073
## 129 0.001 0.13530478 0.9065415 2.559791e-03 0.001124287 0.0102356292
## 130 0.001 0.14849683 0.9064963 1.714023e-03 0.001132386 0.0084829470
## 131 0.001 0.16297508 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 132 0.001 0.17886495 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 133 0.001 0.19630407 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 134 0.001 0.21544347 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 135 0.001 0.23644894 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 136 0.001 0.25950242 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 137 0.001 0.28480359 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 138 0.001 0.31257158 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 139 0.001 0.34304693 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 140 0.001 0.37649358 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 141 0.001 0.41320124 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 142 0.001 0.45348785 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 143 0.001 0.49770236 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 144 0.001 0.54622772 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 145 0.001 0.59948425 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 146 0.001 0.65793322 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 147 0.001 0.72208090 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 148 0.001 0.79248290 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 149 0.001 0.86974900 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 150 0.001 0.95454846 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 151 0.001 1.04761575 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 152 0.001 1.14975700 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 153 0.001 1.26185688 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 154 0.001 1.38488637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 155 0.001 1.51991108 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 156 0.001 1.66810054 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 157 0.001 1.83073828 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 158 0.001 2.00923300 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 159 0.001 2.20513074 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 160 0.001 2.42012826 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 161 0.001 2.65608778 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 162 0.001 2.91505306 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 163 0.001 3.19926714 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 164 0.001 3.51119173 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 165 0.001 3.85352859 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 166 0.001 4.22924287 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 167 0.001 4.64158883 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 168 0.001 5.09413801 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 169 0.001 5.59081018 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 170 0.001 6.13590727 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 171 0.001 6.73415066 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 172 0.001 7.39072203 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 173 0.001 8.11130831 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 174 0.001 8.90215085 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 175 0.001 9.77009957 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 176 0.001 10.72267222 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 177 0.001 11.76811952 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 178 0.001 12.91549665 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 179 0.001 14.17474163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 180 0.001 15.55676144 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 181 0.001 17.07352647 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 182 0.001 18.73817423 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 183 0.001 20.56512308 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 184 0.001 22.57019720 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 185 0.001 24.77076356 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 186 0.001 27.18588243 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 187 0.001 29.83647240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 188 0.001 32.74549163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 189 0.001 35.93813664 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 190 0.001 39.44206059 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 191 0.001 43.28761281 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 192 0.001 47.50810162 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 193 0.001 52.14008288 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 194 0.001 57.22367659 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 195 0.001 62.80291442 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 196 0.001 68.92612104 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 197 0.001 75.64633276 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 198 0.001 83.02175681 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 199 0.001 91.11627561 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 200 0.001 100.00000000 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 201 0.008 0.01000000 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 202 0.008 0.01097499 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 203 0.008 0.01204504 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 204 0.008 0.01321941 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 205 0.008 0.01450829 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 206 0.008 0.01592283 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 207 0.008 0.01747528 0.9062256 1.105122e-02 0.002141771 0.0230620373
## 208 0.008 0.01917910 0.9062255 1.034582e-02 0.002146289 0.0216204550
## 209 0.008 0.02104904 0.9062707 1.043440e-02 0.002098869 0.0215676921
## 210 0.008 0.02310130 0.9062256 9.588442e-03 0.002092654 0.0211172183
## 211 0.008 0.02535364 0.9063610 9.862215e-03 0.001997771 0.0211406162
## 212 0.008 0.02782559 0.9063610 9.104639e-03 0.001945262 0.0206076814
## 213 0.008 0.03053856 0.9063610 9.104639e-03 0.001945262 0.0206076814
## 214 0.008 0.03351603 0.9063610 9.104639e-03 0.001945262 0.0206076814
## 215 0.008 0.03678380 0.9063610 8.321725e-03 0.001835233 0.0199753241
## 216 0.008 0.04037017 0.9064061 8.410259e-03 0.001833650 0.0199274204
## 217 0.008 0.04430621 0.9064060 7.630580e-03 0.001722399 0.0192618766
## 218 0.008 0.04862602 0.9064060 6.851051e-03 0.001660926 0.0185389761
## 219 0.008 0.05336699 0.9064512 6.950605e-03 0.001657927 0.0186963496
## 220 0.008 0.05857021 0.9064512 6.193030e-03 0.001594263 0.0179667139
## 221 0.008 0.06428073 0.9064512 6.193030e-03 0.001594263 0.0179667139
## 222 0.008 0.07054802 0.9064512 5.467220e-03 0.001527643 0.0150640003
## 223 0.008 0.07742637 0.9063609 3.731094e-03 0.001390496 0.0128968321
## 224 0.008 0.08497534 0.9064513 3.923264e-03 0.001451245 0.0131545349
## 225 0.008 0.09326033 0.9065415 4.096896e-03 0.001295782 0.0130416236
## 226 0.008 0.10235310 0.9065866 3.427854e-03 0.001203981 0.0117435073
## 227 0.008 0.11233240 0.9065866 3.427854e-03 0.001203981 0.0117435073
## 228 0.008 0.12328467 0.9065415 2.559791e-03 0.001124287 0.0102356292
## 229 0.008 0.13530478 0.9064963 1.714023e-03 0.001132386 0.0084829470
## 230 0.008 0.14849683 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 231 0.008 0.16297508 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 232 0.008 0.17886495 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 233 0.008 0.19630407 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 234 0.008 0.21544347 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 235 0.008 0.23644894 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 236 0.008 0.25950242 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 237 0.008 0.28480359 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 238 0.008 0.31257158 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 239 0.008 0.34304693 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 240 0.008 0.37649358 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 241 0.008 0.41320124 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 242 0.008 0.45348785 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 243 0.008 0.49770236 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 244 0.008 0.54622772 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 245 0.008 0.59948425 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 246 0.008 0.65793322 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 247 0.008 0.72208090 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 248 0.008 0.79248290 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 249 0.008 0.86974900 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 250 0.008 0.95454846 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 251 0.008 1.04761575 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 252 0.008 1.14975700 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 253 0.008 1.26185688 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 254 0.008 1.38488637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 255 0.008 1.51991108 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 256 0.008 1.66810054 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 257 0.008 1.83073828 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 258 0.008 2.00923300 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 259 0.008 2.20513074 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 260 0.008 2.42012826 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 261 0.008 2.65608778 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 262 0.008 2.91505306 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 263 0.008 3.19926714 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 264 0.008 3.51119173 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 265 0.008 3.85352859 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 266 0.008 4.22924287 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 267 0.008 4.64158883 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 268 0.008 5.09413801 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 269 0.008 5.59081018 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 270 0.008 6.13590727 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 271 0.008 6.73415066 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 272 0.008 7.39072203 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 273 0.008 8.11130831 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 274 0.008 8.90215085 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 275 0.008 9.77009957 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 276 0.008 10.72267222 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 277 0.008 11.76811952 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 278 0.008 12.91549665 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 279 0.008 14.17474163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 280 0.008 15.55676144 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 281 0.008 17.07352647 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 282 0.008 18.73817423 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 283 0.008 20.56512308 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 284 0.008 22.57019720 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 285 0.008 24.77076356 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 286 0.008 27.18588243 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 287 0.008 29.83647240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 288 0.008 32.74549163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 289 0.008 35.93813664 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 290 0.008 39.44206059 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 291 0.008 43.28761281 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 292 0.008 47.50810162 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 293 0.008 52.14008288 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 294 0.008 57.22367659 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 295 0.008 62.80291442 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 296 0.008 68.92612104 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 297 0.008 75.64633276 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 298 0.008 83.02175681 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 299 0.008 91.11627561 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 300 0.008 100.00000000 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 301 0.027 0.01000000 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 302 0.027 0.01097499 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 303 0.027 0.01204504 0.9061806 1.096612e-02 0.002229593 0.0231280435
## 304 0.027 0.01321941 0.9062256 1.105122e-02 0.002141771 0.0230620373
## 305 0.027 0.01450829 0.9062256 1.105122e-02 0.002141771 0.0230620373
## 306 0.027 0.01592283 0.9063159 1.200787e-02 0.002143910 0.0233964351
## 307 0.027 0.01747528 0.9063158 1.130246e-02 0.002148429 0.0220081093
## 308 0.027 0.01917910 0.9063158 1.130246e-02 0.002148429 0.0220081093
## 309 0.027 0.02104904 0.9063610 1.064832e-02 0.002148264 0.0216959403
## 310 0.027 0.02310130 0.9063611 9.869417e-03 0.001992909 0.0209689648
## 311 0.027 0.02535364 0.9064061 9.957801e-03 0.001943769 0.0209167578
## 312 0.027 0.02782559 0.9064061 9.957801e-03 0.001943769 0.0209167578
## 313 0.027 0.03053856 0.9064513 1.004295e-02 0.001886873 0.0208478728
## 314 0.027 0.03351603 0.9064512 9.263271e-03 0.001778956 0.0202768221
## 315 0.027 0.03678380 0.9064512 8.483742e-03 0.001719507 0.0196575856
## 316 0.027 0.04037017 0.9064963 8.572276e-03 0.001715397 0.0196081593
## 317 0.027 0.04430621 0.9064513 6.957808e-03 0.001652059 0.0185031921
## 318 0.027 0.04862602 0.9064963 7.046192e-03 0.001589830 0.0184582386
## 319 0.027 0.05336699 0.9064512 6.193030e-03 0.001594263 0.0179667139
## 320 0.027 0.05857021 0.9064512 6.193030e-03 0.001594263 0.0179667139
## 321 0.027 0.06428073 0.9063609 3.731094e-03 0.001390496 0.0128968321
## 322 0.027 0.07054802 0.9064963 4.008362e-03 0.001302815 0.0130849567
## 323 0.027 0.07742637 0.9065866 4.185430e-03 0.001287096 0.0129975308
## 324 0.027 0.08497534 0.9065416 3.339470e-03 0.001288258 0.0117863840
## 325 0.027 0.09326033 0.9065415 2.559791e-03 0.001124287 0.0102356292
## 326 0.027 0.10235310 0.9064963 1.714023e-03 0.001132386 0.0084829470
## 327 0.027 0.11233240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 328 0.027 0.12328467 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 329 0.027 0.13530478 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 330 0.027 0.14849683 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 331 0.027 0.16297508 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 332 0.027 0.17886495 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 333 0.027 0.19630407 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 334 0.027 0.21544347 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 335 0.027 0.23644894 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 336 0.027 0.25950242 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 337 0.027 0.28480359 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 338 0.027 0.31257158 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 339 0.027 0.34304693 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 340 0.027 0.37649358 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 341 0.027 0.41320124 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 342 0.027 0.45348785 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 343 0.027 0.49770236 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 344 0.027 0.54622772 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 345 0.027 0.59948425 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 346 0.027 0.65793322 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 347 0.027 0.72208090 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 348 0.027 0.79248290 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 349 0.027 0.86974900 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 350 0.027 0.95454846 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 351 0.027 1.04761575 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 352 0.027 1.14975700 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 353 0.027 1.26185688 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 354 0.027 1.38488637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 355 0.027 1.51991108 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 356 0.027 1.66810054 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 357 0.027 1.83073828 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 358 0.027 2.00923300 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 359 0.027 2.20513074 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 360 0.027 2.42012826 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 361 0.027 2.65608778 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 362 0.027 2.91505306 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 363 0.027 3.19926714 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 364 0.027 3.51119173 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 365 0.027 3.85352859 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 366 0.027 4.22924287 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 367 0.027 4.64158883 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 368 0.027 5.09413801 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 369 0.027 5.59081018 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 370 0.027 6.13590727 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 371 0.027 6.73415066 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 372 0.027 7.39072203 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 373 0.027 8.11130831 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 374 0.027 8.90215085 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 375 0.027 9.77009957 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 376 0.027 10.72267222 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 377 0.027 11.76811952 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 378 0.027 12.91549665 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 379 0.027 14.17474163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 380 0.027 15.55676144 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 381 0.027 17.07352647 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 382 0.027 18.73817423 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 383 0.027 20.56512308 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 384 0.027 22.57019720 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 385 0.027 24.77076356 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 386 0.027 27.18588243 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 387 0.027 29.83647240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 388 0.027 32.74549163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 389 0.027 35.93813664 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 390 0.027 39.44206059 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 391 0.027 43.28761281 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 392 0.027 47.50810162 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 393 0.027 52.14008288 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 394 0.027 57.22367659 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 395 0.027 62.80291442 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 396 0.027 68.92612104 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 397 0.027 75.64633276 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 398 0.027 83.02175681 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 399 0.027 91.11627561 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 400 0.027 100.00000000 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 401 0.064 0.01000000 0.9063610 1.285383e-02 0.002148046 0.0237159144
## 402 0.064 0.01097499 0.9063610 1.285383e-02 0.002148046 0.0237159144
## 403 0.064 0.01204504 0.9063610 1.285383e-02 0.002148046 0.0237159144
## 404 0.064 0.01321941 0.9064060 1.295629e-02 0.002151211 0.0238339119
## 405 0.064 0.01450829 0.9064060 1.300159e-02 0.002151211 0.0225424198
## 406 0.064 0.01592283 0.9063611 1.141310e-02 0.002191502 0.0219846637
## 407 0.064 0.01747528 0.9064061 1.150149e-02 0.002146912 0.0219285264
## 408 0.064 0.01917910 0.9063611 9.869417e-03 0.001992909 0.0209689648
## 409 0.064 0.02104904 0.9063610 9.858247e-03 0.001945262 0.0207909087
## 410 0.064 0.02310130 0.9063610 9.858247e-03 0.001945262 0.0207909087
## 411 0.064 0.02535364 0.9064061 9.943396e-03 0.001889512 0.0207220228
## 412 0.064 0.02782559 0.9064512 9.252250e-03 0.001778956 0.0200921854
## 413 0.064 0.03053856 0.9064963 8.572276e-03 0.001715397 0.0196081593
## 414 0.064 0.03351603 0.9064964 7.825871e-03 0.001709724 0.0191759911
## 415 0.064 0.03678380 0.9064964 7.825871e-03 0.001709724 0.0191759911
## 416 0.064 0.04037017 0.9064963 7.046192e-03 0.001589830 0.0184582386
## 417 0.064 0.04430621 0.9064962 6.278128e-03 0.001460445 0.0179048296
## 418 0.064 0.04862602 0.9064962 5.552318e-03 0.001387414 0.0149943413
## 419 0.064 0.05336699 0.9064962 4.008161e-03 0.001310318 0.0130850891
## 420 0.064 0.05857021 0.9064513 2.382874e-03 0.001218210 0.0103178904
## 421 0.064 0.06428073 0.9064964 2.471407e-03 0.001212399 0.0102771793
## 422 0.064 0.07054802 0.9064963 1.714023e-03 0.001132386 0.0084829470
## 423 0.064 0.07742637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 424 0.064 0.08497534 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 425 0.064 0.09326033 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 426 0.064 0.10235310 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 427 0.064 0.11233240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 428 0.064 0.12328467 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 429 0.064 0.13530478 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 430 0.064 0.14849683 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 431 0.064 0.16297508 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 432 0.064 0.17886495 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 433 0.064 0.19630407 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 434 0.064 0.21544347 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 435 0.064 0.23644894 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 436 0.064 0.25950242 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 437 0.064 0.28480359 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 438 0.064 0.31257158 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 439 0.064 0.34304693 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 440 0.064 0.37649358 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 441 0.064 0.41320124 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 442 0.064 0.45348785 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 443 0.064 0.49770236 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 444 0.064 0.54622772 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 445 0.064 0.59948425 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 446 0.064 0.65793322 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 447 0.064 0.72208090 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 448 0.064 0.79248290 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 449 0.064 0.86974900 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 450 0.064 0.95454846 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 451 0.064 1.04761575 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 452 0.064 1.14975700 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 453 0.064 1.26185688 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 454 0.064 1.38488637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 455 0.064 1.51991108 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 456 0.064 1.66810054 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 457 0.064 1.83073828 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 458 0.064 2.00923300 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 459 0.064 2.20513074 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 460 0.064 2.42012826 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 461 0.064 2.65608778 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 462 0.064 2.91505306 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 463 0.064 3.19926714 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 464 0.064 3.51119173 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 465 0.064 3.85352859 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 466 0.064 4.22924287 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 467 0.064 4.64158883 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 468 0.064 5.09413801 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 469 0.064 5.59081018 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 470 0.064 6.13590727 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 471 0.064 6.73415066 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 472 0.064 7.39072203 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 473 0.064 8.11130831 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 474 0.064 8.90215085 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 475 0.064 9.77009957 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 476 0.064 10.72267222 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 477 0.064 11.76811952 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 478 0.064 12.91549665 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 479 0.064 14.17474163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 480 0.064 15.55676144 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 481 0.064 17.07352647 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 482 0.064 18.73817423 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 483 0.064 20.56512308 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 484 0.064 22.57019720 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 485 0.064 24.77076356 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 486 0.064 27.18588243 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 487 0.064 29.83647240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 488 0.064 32.74549163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 489 0.064 35.93813664 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 490 0.064 39.44206059 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 491 0.064 43.28761281 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 492 0.064 47.50810162 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 493 0.064 52.14008288 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 494 0.064 57.22367659 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 495 0.064 62.80291442 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 496 0.064 68.92612104 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 497 0.064 75.64633276 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 498 0.064 83.02175681 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 499 0.064 91.11627561 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 500 0.064 100.00000000 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 501 0.125 0.01000000 0.9064512 1.380946e-02 0.002148896 0.0239621219
## 502 0.125 0.01097499 0.9064061 1.225906e-02 0.002194602 0.0223473420
## 503 0.125 0.01204504 0.9063610 1.140193e-02 0.002148264 0.0218157061
## 504 0.125 0.01321941 0.9063610 1.140193e-02 0.002148264 0.0218157061
## 505 0.125 0.01450829 0.9062708 9.687908e-03 0.002094469 0.0209249372
## 506 0.125 0.01592283 0.9063610 9.858247e-03 0.001945262 0.0207909087
## 507 0.125 0.01747528 0.9063610 9.858247e-03 0.001945262 0.0207909087
## 508 0.125 0.01917910 0.9064964 1.013148e-02 0.001883127 0.0207945022
## 509 0.125 0.02104904 0.9064061 8.395358e-03 0.001776033 0.0197063825
## 510 0.125 0.02310130 0.9064512 8.483742e-03 0.001719507 0.0196575856
## 511 0.125 0.02535364 0.9064964 7.825871e-03 0.001709724 0.0191759911
## 512 0.125 0.02782559 0.9065415 7.914255e-03 0.001648414 0.0191285262
## 513 0.125 0.03053856 0.9064512 6.193030e-03 0.001594263 0.0179667139
## 514 0.125 0.03351603 0.9064962 6.278128e-03 0.001460445 0.0179048296
## 515 0.125 0.03678380 0.9064060 3.073668e-03 0.001312087 0.0119104917
## 516 0.125 0.04037017 0.9063609 2.205605e-03 0.001232634 0.0103966579
## 517 0.125 0.04430621 0.9064061 1.514810e-03 0.001134027 0.0084537425
## 518 0.125 0.04862602 0.9063611 -8.838384e-05 0.001128015 0.0006249681
## 519 0.125 0.05336699 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 520 0.125 0.05857021 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 521 0.125 0.06428073 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 522 0.125 0.07054802 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 523 0.125 0.07742637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 524 0.125 0.08497534 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 525 0.125 0.09326033 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 526 0.125 0.10235310 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 527 0.125 0.11233240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 528 0.125 0.12328467 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 529 0.125 0.13530478 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 530 0.125 0.14849683 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 531 0.125 0.16297508 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 532 0.125 0.17886495 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 533 0.125 0.19630407 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 534 0.125 0.21544347 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 535 0.125 0.23644894 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 536 0.125 0.25950242 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 537 0.125 0.28480359 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 538 0.125 0.31257158 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 539 0.125 0.34304693 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 540 0.125 0.37649358 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 541 0.125 0.41320124 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 542 0.125 0.45348785 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 543 0.125 0.49770236 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 544 0.125 0.54622772 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 545 0.125 0.59948425 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 546 0.125 0.65793322 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 547 0.125 0.72208090 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 548 0.125 0.79248290 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 549 0.125 0.86974900 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 550 0.125 0.95454846 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 551 0.125 1.04761575 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 552 0.125 1.14975700 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 553 0.125 1.26185688 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 554 0.125 1.38488637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 555 0.125 1.51991108 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 556 0.125 1.66810054 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 557 0.125 1.83073828 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 558 0.125 2.00923300 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 559 0.125 2.20513074 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 560 0.125 2.42012826 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 561 0.125 2.65608778 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 562 0.125 2.91505306 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 563 0.125 3.19926714 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 564 0.125 3.51119173 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 565 0.125 3.85352859 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 566 0.125 4.22924287 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 567 0.125 4.64158883 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 568 0.125 5.09413801 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 569 0.125 5.59081018 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 570 0.125 6.13590727 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 571 0.125 6.73415066 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 572 0.125 7.39072203 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 573 0.125 8.11130831 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 574 0.125 8.90215085 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 575 0.125 9.77009957 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 576 0.125 10.72267222 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 577 0.125 11.76811952 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 578 0.125 12.91549665 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 579 0.125 14.17474163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 580 0.125 15.55676144 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 581 0.125 17.07352647 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 582 0.125 18.73817423 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 583 0.125 20.56512308 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 584 0.125 22.57019720 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 585 0.125 24.77076356 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 586 0.125 27.18588243 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 587 0.125 29.83647240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 588 0.125 32.74549163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 589 0.125 35.93813664 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 590 0.125 39.44206059 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 591 0.125 43.28761281 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 592 0.125 47.50810162 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 593 0.125 52.14008288 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 594 0.125 57.22367659 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 595 0.125 62.80291442 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 596 0.125 68.92612104 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 597 0.125 75.64633276 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 598 0.125 83.02175681 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 599 0.125 91.11627561 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 600 0.125 100.00000000 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 601 0.216 0.01000000 0.9063159 1.055597e-02 0.002144129 0.0214095069
## 602 0.216 0.01097499 0.9062708 9.687908e-03 0.002094469 0.0209249372
## 603 0.216 0.01204504 0.9063611 9.869417e-03 0.001992909 0.0209689648
## 604 0.216 0.01321941 0.9064061 9.174887e-03 0.001833650 0.0203275912
## 605 0.216 0.01450829 0.9064513 9.263421e-03 0.001830931 0.0202767189
## 606 0.216 0.01592283 0.9064512 8.483742e-03 0.001719507 0.0196575856
## 607 0.216 0.01747528 0.9064512 8.483742e-03 0.001719507 0.0196575856
## 608 0.216 0.01917910 0.9065415 7.914255e-03 0.001648414 0.0191285262
## 609 0.216 0.02104904 0.9064963 7.046192e-03 0.001589830 0.0184582386
## 610 0.216 0.02310130 0.9063610 4.479007e-03 0.001526446 0.0163251343
## 611 0.216 0.02535364 0.9064963 4.756074e-03 0.001453777 0.0164614212
## 612 0.216 0.02782559 0.9063609 2.205605e-03 0.001232634 0.0103966579
## 613 0.216 0.03053856 0.9062708 -2.654514e-04 0.001127638 0.0010613545
## 614 0.216 0.03351603 0.9063159 -1.769176e-04 0.001128748 0.0008755156
## 615 0.216 0.03678380 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 616 0.216 0.04037017 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 617 0.216 0.04430621 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 618 0.216 0.04862602 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 619 0.216 0.05336699 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 620 0.216 0.05857021 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 621 0.216 0.06428073 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 622 0.216 0.07054802 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 623 0.216 0.07742637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 624 0.216 0.08497534 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 625 0.216 0.09326033 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 626 0.216 0.10235310 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 627 0.216 0.11233240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 628 0.216 0.12328467 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 629 0.216 0.13530478 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 630 0.216 0.14849683 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 631 0.216 0.16297508 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 632 0.216 0.17886495 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 633 0.216 0.19630407 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 634 0.216 0.21544347 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 635 0.216 0.23644894 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 636 0.216 0.25950242 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 637 0.216 0.28480359 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 638 0.216 0.31257158 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 639 0.216 0.34304693 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 640 0.216 0.37649358 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 641 0.216 0.41320124 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 642 0.216 0.45348785 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 643 0.216 0.49770236 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 644 0.216 0.54622772 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 645 0.216 0.59948425 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 646 0.216 0.65793322 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 647 0.216 0.72208090 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 648 0.216 0.79248290 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 649 0.216 0.86974900 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 650 0.216 0.95454846 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 651 0.216 1.04761575 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 652 0.216 1.14975700 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 653 0.216 1.26185688 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 654 0.216 1.38488637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 655 0.216 1.51991108 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 656 0.216 1.66810054 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 657 0.216 1.83073828 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 658 0.216 2.00923300 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 659 0.216 2.20513074 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 660 0.216 2.42012826 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 661 0.216 2.65608778 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 662 0.216 2.91505306 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 663 0.216 3.19926714 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 664 0.216 3.51119173 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 665 0.216 3.85352859 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 666 0.216 4.22924287 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 667 0.216 4.64158883 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 668 0.216 5.09413801 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 669 0.216 5.59081018 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 670 0.216 6.13590727 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 671 0.216 6.73415066 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 672 0.216 7.39072203 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 673 0.216 8.11130831 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 674 0.216 8.90215085 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 675 0.216 9.77009957 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 676 0.216 10.72267222 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 677 0.216 11.76811952 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 678 0.216 12.91549665 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 679 0.216 14.17474163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 680 0.216 15.55676144 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 681 0.216 17.07352647 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 682 0.216 18.73817423 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 683 0.216 20.56512308 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 684 0.216 22.57019720 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 685 0.216 24.77076356 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 686 0.216 27.18588243 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 687 0.216 29.83647240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 688 0.216 32.74549163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 689 0.216 35.93813664 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 690 0.216 39.44206059 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 691 0.216 43.28761281 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 692 0.216 47.50810162 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 693 0.216 52.14008288 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 694 0.216 57.22367659 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 695 0.216 62.80291442 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 696 0.216 68.92612104 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 697 0.216 75.64633276 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 698 0.216 83.02175681 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 699 0.216 91.11627561 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 700 0.216 100.00000000 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 701 0.343 0.01000000 0.9064513 9.263421e-03 0.001830931 0.0202767189
## 702 0.343 0.01097499 0.9064061 8.395358e-03 0.001776033 0.0197063825
## 703 0.343 0.01204504 0.9064962 8.572126e-03 0.001659812 0.0196082607
## 704 0.343 0.01321941 0.9064512 7.726166e-03 0.001658209 0.0190273574
## 705 0.343 0.01450829 0.9064513 6.200232e-03 0.001588160 0.0177659385
## 706 0.343 0.01592283 0.9063610 4.479007e-03 0.001526446 0.0163251343
## 707 0.343 0.01747528 0.9063610 4.479007e-03 0.001526446 0.0163251343
## 708 0.343 0.01917910 0.9063609 2.205605e-03 0.001232634 0.0103966579
## 709 0.343 0.02104904 0.9062708 -2.654514e-04 0.001127638 0.0010613545
## 710 0.343 0.02310130 0.9063611 -8.838384e-05 0.001128015 0.0006249681
## 711 0.343 0.02535364 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 712 0.343 0.02782559 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 713 0.343 0.03053856 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 714 0.343 0.03351603 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 715 0.343 0.03678380 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 716 0.343 0.04037017 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 717 0.343 0.04430621 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 718 0.343 0.04862602 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 719 0.343 0.05336699 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 720 0.343 0.05857021 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 721 0.343 0.06428073 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 722 0.343 0.07054802 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 723 0.343 0.07742637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 724 0.343 0.08497534 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 725 0.343 0.09326033 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 726 0.343 0.10235310 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 727 0.343 0.11233240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 728 0.343 0.12328467 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 729 0.343 0.13530478 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 730 0.343 0.14849683 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 731 0.343 0.16297508 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 732 0.343 0.17886495 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 733 0.343 0.19630407 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 734 0.343 0.21544347 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 735 0.343 0.23644894 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 736 0.343 0.25950242 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 737 0.343 0.28480359 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 738 0.343 0.31257158 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 739 0.343 0.34304693 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 740 0.343 0.37649358 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 741 0.343 0.41320124 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 742 0.343 0.45348785 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 743 0.343 0.49770236 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 744 0.343 0.54622772 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 745 0.343 0.59948425 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 746 0.343 0.65793322 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 747 0.343 0.72208090 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 748 0.343 0.79248290 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 749 0.343 0.86974900 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 750 0.343 0.95454846 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 751 0.343 1.04761575 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 752 0.343 1.14975700 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 753 0.343 1.26185688 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 754 0.343 1.38488637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 755 0.343 1.51991108 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 756 0.343 1.66810054 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 757 0.343 1.83073828 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 758 0.343 2.00923300 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 759 0.343 2.20513074 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 760 0.343 2.42012826 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 761 0.343 2.65608778 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 762 0.343 2.91505306 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 763 0.343 3.19926714 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 764 0.343 3.51119173 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 765 0.343 3.85352859 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 766 0.343 4.22924287 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 767 0.343 4.64158883 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 768 0.343 5.09413801 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 769 0.343 5.59081018 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 770 0.343 6.13590727 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 771 0.343 6.73415066 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 772 0.343 7.39072203 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 773 0.343 8.11130831 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 774 0.343 8.90215085 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 775 0.343 9.77009957 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 776 0.343 10.72267222 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 777 0.343 11.76811952 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 778 0.343 12.91549665 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 779 0.343 14.17474163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 780 0.343 15.55676144 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 781 0.343 17.07352647 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 782 0.343 18.73817423 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 783 0.343 20.56512308 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 784 0.343 22.57019720 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 785 0.343 24.77076356 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 786 0.343 27.18588243 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 787 0.343 29.83647240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 788 0.343 32.74549163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 789 0.343 35.93813664 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 790 0.343 39.44206059 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 791 0.343 43.28761281 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 792 0.343 47.50810162 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 793 0.343 52.14008288 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 794 0.343 57.22367659 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 795 0.343 62.80291442 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 796 0.343 68.92612104 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 797 0.343 75.64633276 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 798 0.343 83.02175681 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 799 0.343 91.11627561 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 800 0.343 100.00000000 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 801 0.512 0.01000000 0.9064061 6.111698e-03 0.001591294 0.0178084460
## 802 0.512 0.01097499 0.9063610 4.479007e-03 0.001526446 0.0163251343
## 803 0.512 0.01204504 0.9063610 4.479007e-03 0.001526446 0.0163251343
## 804 0.512 0.01321941 0.9064512 3.941731e-03 0.001385496 0.0131943665
## 805 0.512 0.01450829 0.9062708 -2.654514e-04 0.001127638 0.0010613545
## 806 0.512 0.01592283 0.9062708 -2.654514e-04 0.001127638 0.0010613545
## 807 0.512 0.01747528 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 808 0.512 0.01917910 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 809 0.512 0.02104904 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 810 0.512 0.02310130 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 811 0.512 0.02535364 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 812 0.512 0.02782559 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 813 0.512 0.03053856 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 814 0.512 0.03351603 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 815 0.512 0.03678380 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 816 0.512 0.04037017 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 817 0.512 0.04430621 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 818 0.512 0.04862602 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 819 0.512 0.05336699 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 820 0.512 0.05857021 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 821 0.512 0.06428073 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 822 0.512 0.07054802 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 823 0.512 0.07742637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 824 0.512 0.08497534 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 825 0.512 0.09326033 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 826 0.512 0.10235310 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 827 0.512 0.11233240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 828 0.512 0.12328467 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 829 0.512 0.13530478 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 830 0.512 0.14849683 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 831 0.512 0.16297508 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 832 0.512 0.17886495 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 833 0.512 0.19630407 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 834 0.512 0.21544347 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 835 0.512 0.23644894 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 836 0.512 0.25950242 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 837 0.512 0.28480359 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 838 0.512 0.31257158 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 839 0.512 0.34304693 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 840 0.512 0.37649358 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 841 0.512 0.41320124 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 842 0.512 0.45348785 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 843 0.512 0.49770236 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 844 0.512 0.54622772 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 845 0.512 0.59948425 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 846 0.512 0.65793322 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 847 0.512 0.72208090 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 848 0.512 0.79248290 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 849 0.512 0.86974900 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 850 0.512 0.95454846 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 851 0.512 1.04761575 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 852 0.512 1.14975700 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 853 0.512 1.26185688 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 854 0.512 1.38488637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 855 0.512 1.51991108 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 856 0.512 1.66810054 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 857 0.512 1.83073828 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 858 0.512 2.00923300 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 859 0.512 2.20513074 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 860 0.512 2.42012826 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 861 0.512 2.65608778 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 862 0.512 2.91505306 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 863 0.512 3.19926714 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 864 0.512 3.51119173 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 865 0.512 3.85352859 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 866 0.512 4.22924287 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 867 0.512 4.64158883 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 868 0.512 5.09413801 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 869 0.512 5.59081018 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 870 0.512 6.13590727 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 871 0.512 6.73415066 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 872 0.512 7.39072203 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 873 0.512 8.11130831 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 874 0.512 8.90215085 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 875 0.512 9.77009957 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 876 0.512 10.72267222 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 877 0.512 11.76811952 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 878 0.512 12.91549665 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 879 0.512 14.17474163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 880 0.512 15.55676144 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 881 0.512 17.07352647 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 882 0.512 18.73817423 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 883 0.512 20.56512308 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 884 0.512 22.57019720 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 885 0.512 24.77076356 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 886 0.512 27.18588243 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 887 0.512 29.83647240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 888 0.512 32.74549163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 889 0.512 35.93813664 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 890 0.512 39.44206059 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 891 0.512 43.28761281 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 892 0.512 47.50810162 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 893 0.512 52.14008288 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 894 0.512 57.22367659 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 895 0.512 62.80291442 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 896 0.512 68.92612104 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 897 0.512 75.64633276 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 898 0.512 83.02175681 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 899 0.512 91.11627561 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 900 0.512 100.00000000 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 901 0.729 0.01000000 0.9062707 5.138766e-04 0.001224382 0.0063068401
## 902 0.729 0.01097499 0.9062708 -2.654514e-04 0.001127638 0.0010613545
## 903 0.729 0.01204504 0.9063159 -1.769176e-04 0.001128748 0.0008755156
## 904 0.729 0.01321941 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 905 0.729 0.01450829 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 906 0.729 0.01592283 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 907 0.729 0.01747528 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 908 0.729 0.01917910 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 909 0.729 0.02104904 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 910 0.729 0.02310130 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 911 0.729 0.02535364 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 912 0.729 0.02782559 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 913 0.729 0.03053856 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 914 0.729 0.03351603 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 915 0.729 0.03678380 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 916 0.729 0.04037017 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 917 0.729 0.04430621 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 918 0.729 0.04862602 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 919 0.729 0.05336699 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 920 0.729 0.05857021 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 921 0.729 0.06428073 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 922 0.729 0.07054802 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 923 0.729 0.07742637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 924 0.729 0.08497534 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 925 0.729 0.09326033 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 926 0.729 0.10235310 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 927 0.729 0.11233240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 928 0.729 0.12328467 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 929 0.729 0.13530478 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 930 0.729 0.14849683 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 931 0.729 0.16297508 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 932 0.729 0.17886495 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 933 0.729 0.19630407 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 934 0.729 0.21544347 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 935 0.729 0.23644894 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 936 0.729 0.25950242 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 937 0.729 0.28480359 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 938 0.729 0.31257158 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 939 0.729 0.34304693 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 940 0.729 0.37649358 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 941 0.729 0.41320124 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 942 0.729 0.45348785 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 943 0.729 0.49770236 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 944 0.729 0.54622772 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 945 0.729 0.59948425 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 946 0.729 0.65793322 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 947 0.729 0.72208090 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 948 0.729 0.79248290 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 949 0.729 0.86974900 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 950 0.729 0.95454846 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 951 0.729 1.04761575 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 952 0.729 1.14975700 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 953 0.729 1.26185688 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 954 0.729 1.38488637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 955 0.729 1.51991108 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 956 0.729 1.66810054 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 957 0.729 1.83073828 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 958 0.729 2.00923300 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 959 0.729 2.20513074 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 960 0.729 2.42012826 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 961 0.729 2.65608778 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 962 0.729 2.91505306 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 963 0.729 3.19926714 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 964 0.729 3.51119173 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 965 0.729 3.85352859 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 966 0.729 4.22924287 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 967 0.729 4.64158883 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 968 0.729 5.09413801 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 969 0.729 5.59081018 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 970 0.729 6.13590727 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 971 0.729 6.73415066 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 972 0.729 7.39072203 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 973 0.729 8.11130831 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 974 0.729 8.90215085 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 975 0.729 9.77009957 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 976 0.729 10.72267222 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 977 0.729 11.76811952 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 978 0.729 12.91549665 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 979 0.729 14.17474163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 980 0.729 15.55676144 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 981 0.729 17.07352647 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 982 0.729 18.73817423 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 983 0.729 20.56512308 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 984 0.729 22.57019720 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 985 0.729 24.77076356 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 986 0.729 27.18588243 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 987 0.729 29.83647240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 988 0.729 32.74549163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 989 0.729 35.93813664 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 990 0.729 39.44206059 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 991 0.729 43.28761281 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 992 0.729 47.50810162 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 993 0.729 52.14008288 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 994 0.729 57.22367659 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 995 0.729 62.80291442 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 996 0.729 68.92612104 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 997 0.729 75.64633276 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 998 0.729 83.02175681 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 999 0.729 91.11627561 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1000 0.729 100.00000000 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1001 1.000 0.01000000 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1002 1.000 0.01097499 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1003 1.000 0.01204504 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1004 1.000 0.01321941 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1005 1.000 0.01450829 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1006 1.000 0.01592283 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1007 1.000 0.01747528 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1008 1.000 0.01917910 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1009 1.000 0.02104904 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1010 1.000 0.02310130 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1011 1.000 0.02535364 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1012 1.000 0.02782559 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1013 1.000 0.03053856 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1014 1.000 0.03351603 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1015 1.000 0.03678380 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1016 1.000 0.04037017 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1017 1.000 0.04430621 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1018 1.000 0.04862602 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1019 1.000 0.05336699 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1020 1.000 0.05857021 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1021 1.000 0.06428073 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1022 1.000 0.07054802 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1023 1.000 0.07742637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1024 1.000 0.08497534 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1025 1.000 0.09326033 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1026 1.000 0.10235310 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1027 1.000 0.11233240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1028 1.000 0.12328467 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1029 1.000 0.13530478 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1030 1.000 0.14849683 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1031 1.000 0.16297508 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1032 1.000 0.17886495 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1033 1.000 0.19630407 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1034 1.000 0.21544347 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1035 1.000 0.23644894 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1036 1.000 0.25950242 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1037 1.000 0.28480359 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1038 1.000 0.31257158 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1039 1.000 0.34304693 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1040 1.000 0.37649358 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1041 1.000 0.41320124 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1042 1.000 0.45348785 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1043 1.000 0.49770236 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1044 1.000 0.54622772 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1045 1.000 0.59948425 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1046 1.000 0.65793322 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1047 1.000 0.72208090 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1048 1.000 0.79248290 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1049 1.000 0.86974900 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1050 1.000 0.95454846 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1051 1.000 1.04761575 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1052 1.000 1.14975700 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1053 1.000 1.26185688 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1054 1.000 1.38488637 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1055 1.000 1.51991108 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1056 1.000 1.66810054 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1057 1.000 1.83073828 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1058 1.000 2.00923300 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1059 1.000 2.20513074 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1060 1.000 2.42012826 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1061 1.000 2.65608778 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1062 1.000 2.91505306 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1063 1.000 3.19926714 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1064 1.000 3.51119173 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1065 1.000 3.85352859 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1066 1.000 4.22924287 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1067 1.000 4.64158883 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1068 1.000 5.09413801 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1069 1.000 5.59081018 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1070 1.000 6.13590727 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1071 1.000 6.73415066 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1072 1.000 7.39072203 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1073 1.000 8.11130831 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1074 1.000 8.90215085 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1075 1.000 9.77009957 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1076 1.000 10.72267222 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1077 1.000 11.76811952 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1078 1.000 12.91549665 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1079 1.000 14.17474163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1080 1.000 15.55676144 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1081 1.000 17.07352647 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1082 1.000 18.73817423 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1083 1.000 20.56512308 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1084 1.000 22.57019720 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1085 1.000 24.77076356 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1086 1.000 27.18588243 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1087 1.000 29.83647240 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1088 1.000 32.74549163 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1089 1.000 35.93813664 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1090 1.000 39.44206059 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1091 1.000 43.28761281 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1092 1.000 47.50810162 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1093 1.000 52.14008288 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1094 1.000 57.22367659 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1095 1.000 62.80291442 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1096 1.000 68.92612104 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1097 1.000 75.64633276 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1098 1.000 83.02175681 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1099 1.000 91.11627561 0.9064061 0.000000e+00 0.001038735 0.0000000000
## 1100 1.000 100.00000000 0.9064061 0.000000e+00 0.001038735 0.0000000000coef(grid.res$finalModel,
s=grid.res$bestTune$lambda) #最適なalpha, lambdaのもとでの回帰係数## 17 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) -6.976642e+00
## SEX -5.213194e-02
## CURSMOKE 6.053955e-02
## DIABETES 3.924487e-01
## BPMEDS 4.669276e-01
## PREVCHD 1.658391e-01
## PREVAP 6.559511e-02
## PREVMI 1.711201e-01
## PREVHYP 3.117723e-01
## TOTCHOL -7.454905e-05
## AGE 2.562016e-02
## SYSBP 5.841279e-03
## DIABP 8.863818e-03
## CIGPDAY 1.484293e-03
## BMI 1.200504e-02
## HEARTRTE -3.768294e-03
## GLUCOSE 6.241628e-04Lasso model (再掲)
cva.glmnetで求めた最適なalpha, lamdaのモデル
lasso.cv.res<-
cv.glmnet(
x=fr.missf_p1[,var] %>% data.matrix(),
y=fr.missf_p1[,"STROKE"] %>% data.matrix(),
family="binomial",
alpha=1)
coef(lasso.cv.res,s="lambda.min") ## 17 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) -9.4090042597
## SEX -0.0421977679
## CURSMOKE 0.2735297679
## DIABETES 0.6209602182
## BPMEDS 0.5182370318
## PREVCHD 0.0312762345
## PREVAP .
## PREVMI 0.1297005059
## PREVHYP 0.5534816855
## TOTCHOL -0.0009247672
## AGE 0.0675262370
## SYSBP 0.0012615861
## DIABP 0.0179817879
## CIGPDAY 0.0027332234
## BMI 0.0116862895
## HEARTRTE -0.0094909742
## GLUCOSE .cf<-coef(lasso.cv.res,s="lambda.min") %>%
as.matrix() %>%
as.data.frame() %>%
mutate(exp=exp(s1)) %>%
add_rownames(var="var") ## Warning: `add_rownames()` was deprecated in dplyr 1.0.0.
## ℹ Please use `tibble::rownames_to_column()` instead.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.#plot
cf %>%
ggplot(aes(y=var,x=exp))+
geom_point()+
geom_vline(xintercept=1,color="red")
- Lasso regressionでは一部の変数の回帰係数が0にshrinkageしており、実質的に変数選択が行われたことがわかる
- Shrinkageしたのは、PREVAP, GLUCOSEであり、これらは予測寄与が低いことを意味する
Grid searchによるalpha, labdaの最適化を図ったときのモデル
cf2<-
coef(grid.res$finalModel,
s=grid.res$bestTune$lambda) %>%
as.matrix() %>%
as.data.frame() %>%
mutate(exp=exp(s1)) %>%
add_rownames(var="var") ## Warning: `add_rownames()` was deprecated in dplyr 1.0.0.
## ℹ Please use `tibble::rownames_to_column()` instead.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.#plot
cf2 %>%
ggplot(aes(y=var,x=exp))+
geom_point()+
geom_vline(xintercept=1,color="red")
- cva.glmnetと大きくは変わらない結果となった
- 予測寄与度の高い因子も変わらず
- Glucoseなど予測寄与度の低い因子は完全に0にshrinkはしなかったが実質効果なし(OR=1)に近似している
Lasso modelを使った予測性能評価
fr.missf_p1$lasso<-
predict(lasso.cv.res,
s="lambda.min",
newx=fr.missf_p1[,var] %>% data.matrix(),
type="response"
)
fr.missf_p1<-
fr.missf_p1 %>%
mutate(STROKEn=as.numeric(STROKE)-1)ROC curve
#ROC curve
roc.lasso<-
roc(STROKEn~lasso,
data=fr.missf_p1,
ci=T)## Setting levels: control = 0, case = 1## Warning in roc.default(response, predictors[, 1], ...): Deprecated use a matrix
## as predictor. Unexpected results may be produced, please pass a numeric vector.## Setting direction: controls < casesroc.lasso##
## Call:
## roc.formula(formula = STROKEn ~ lasso, data = fr.missf_p1, ci = T)
##
## Data: lasso in 4019 controls (STROKEn 0) < 415 cases (STROKEn 1).
## Area under the curve: 0.7515
## 95% CI: 0.7276-0.7753 (DeLong)roc.lasso %>%
ggroc(legacy.axes=T)+
geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="grey", linetype="dashed")
- AUC=0.7514547と中等度の予測性能
Calibration plot
val.prob(
p=fr.missf_p1$lasso,
y=fr.missf_p1$STROKEn,
g=10,
cex=0.5
)
## Dxy C (ROC) R2 D D:Chi-sq
## 0.502857211 0.751428606 0.148173033 0.070828019 315.051435537
## D:p U U:Chi-sq U:p Q
## NA -0.000268482 0.809550747 0.667126648 0.071096501
## Brier Intercept Slope Emax E90
## 0.077920351 0.111304815 1.056365459 0.009301043 0.008167332
## Eavg S:z S:p
## 0.004960714 -0.128408818 0.897825453- Calibration性能はかなり高い事がわかる
BootstrapによるLasso modelの内的妥当性検証
optimismの算出
B<-100 #number of resampling
N<-dim(fr.missf_p1)[1] #dataのdimentionを返す→4434行x39列→そのうち一番目（4434)を取得しNとする
AUC.i1<-AUC.i2<-numeric(B)
for(i in 1:B){
#bootstrap random sampling
bs.i<-sample(1:N,N,replace=TRUE) #重複を許して1~4434までの数を4434回サンプリングする
#bootstrap samplingに基づき、dataから集団を再構成 →data.iとする
fr.missf_p1.i<-fr.missf_p1[bs.i,]
#bootstrapしたデータセットごとで回帰モデル構築→各bootstrap dataからリスク予測確率の算出→data2.1のprob_1列に格納
lasso.cv.res.i<-
cv.glmnet(
x=fr.missf_p1.i[,var] %>% data.matrix(),
y=fr.missf_p1.i[,"STROKE"] %>% data.matrix(),
family="binomial",
alpha=1)
fr.missf_p1.i$lasso1<- #lasso1では、bootstrap dataから構築した新しいLasso回帰モデルを使い、そのbootstrap dataに対するリスク予測確率を格納
predict(lasso.cv.res.i,
s="lambda.min",
newx=fr.missf_p1.i[,var] %>% data.matrix(),
type="response"
)
#AUCの算出；Bootstrap回数（B=100)分のAUCを算出して結果を格納；各bootstrap dataに対するAUC
AUC.i1[i]<-roc(STROKEn~lasso1,data=fr.missf_p1.i)$auc
#bootstrapで作成した各回帰モデルを使って、オリジナルデータ(fr.missf_p1)に対するリスク予測確率を算出→fr.missf_p1のlasso2列に格納;AUC算出
fr.missf_p1$lasso2<-
predict(lasso.cv.res.i,
s="lambda.min",
newx=fr.missf_p1[,var] %>% data.matrix(),
type="response"
)
AUC.i2[i]<-roc(STROKEn~lasso2,data=fr.missf_p1)$auc #ここで得られるAUCはbootstrap dataを用いてoriginal dataに外挿したときのAUC
#print(paste(i,"th bootstrap iteration is completed.",sep=""))
}
opt1<-AUC.i1-AUC.i2
#(bootstrap samplingした各data setを使い、それぞれ導いた回帰モデルにおけるAUC)ー(bootstrap samplingした各data setを使い、original dataに外挿したときのAUC)
summary(opt1);hist(opt1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.031467 -0.002995 0.006353 0.006422 0.015910 0.034586
optimismを用いたbias corrected AUCの算出
lam1<-mean(opt1) #estimate of the optimism
cor.AUC<-roc.lasso$auc-lam1 #bias corrected AUC estimate
cor.AUC## [1] 0.7450327- Lasso modelによる予測モデル構築を行い、Bootstrapによるbias corrected AUCを算出したところ、0.7450327であった
**疑問 - 今回のモデルでは、事前に考えていた非線形性や交互作用項がモデリングに組み込まれなかったが、実際にはそれらの柔軟なモデリングをLassoに適用することはできるのか？交互作用項に関しては、ペアの変数の掛け算項を列結合すれば良さそうだが、非線形性はどのように表現すればよいのかわからない。
Cross validationによる内的妥当性の検証
#k-fold CV
k <-5
AUC_CV <-data.frame(matrix(ncol=2,nrow=k))
#Randomly shuffle the data
fr.missf_p1.r<-fr.missf_p1[sample(nrow(fr.missf_p1)),]
#Create k equally size folds
folds <- cut(seq(1,nrow(fr.missf_p1.r)),breaks=k,labels=FALSE)
#Perform k fold cross validation
for(i in 1:k){
#Segment data2r by fold using the which() function
testIndexes <- which(folds==i,arr.ind=TRUE)
testData <- fr.missf_p1.r[testIndexes, ]
trainData <- fr.missf_p1.r[-testIndexes, ]
fit1<-
cv.glmnet(
x=trainData[,var] %>% data.matrix(),
y=trainData[,"STROKE"] %>% data.matrix(),
family="binomial",
alpha=1)
testData$fitted<-
predict(lasso.cv.res,
s="lambda.min",
newx=testData[,var] %>% data.matrix(),
type="response"
)
ROC <- roc(testData$STROKE, testData$fitted)
AUC_CV[i,2] <- ROC$auc
AUC_CV[i,1] <- i
#print(paste(i,"th cross-validation iteration is completed.", sep=""))
} ## Setting levels: control = 0, case = 1## Warning in roc.default(testData$STROKE, testData$fitted): Deprecated use a
## matrix as predictor. Unexpected results may be produced, please pass a numeric
## vector.## Setting direction: controls < cases
## Setting levels: control = 0, case = 1## Warning in roc.default(testData$STROKE, testData$fitted): Deprecated use a
## matrix as predictor. Unexpected results may be produced, please pass a numeric
## vector.## Setting direction: controls < cases
## Setting levels: control = 0, case = 1## Warning in roc.default(testData$STROKE, testData$fitted): Deprecated use a
## matrix as predictor. Unexpected results may be produced, please pass a numeric
## vector.## Setting direction: controls < cases
## Setting levels: control = 0, case = 1## Warning in roc.default(testData$STROKE, testData$fitted): Deprecated use a
## matrix as predictor. Unexpected results may be produced, please pass a numeric
## vector.## Setting direction: controls < cases
## Setting levels: control = 0, case = 1## Warning in roc.default(testData$STROKE, testData$fitted): Deprecated use a
## matrix as predictor. Unexpected results may be produced, please pass a numeric
## vector.## Setting direction: controls < casesnames(AUC_CV)[2] <- "AUC_cv"
names(AUC_CV)[1] <- "k"
print(summary(AUC_CV$AUC_cv)) ## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.7144 0.7307 0.7632 0.7497 0.7695 0.7705- CVによる内的妥当性検証ではAUCは0.7496755であり、Bootstrapによるbias corrected AUC=0.7450327と遜色ない。
Backward stepwiseによる変数選択とモデリング
ここでは、Lasso modelとの比較可能性を高めるために、2つの回帰モデルを考える 1. 交互作用項、非線形項を含まない線形の全変数モデル 2. 交互作用項、非線形項を含む全変数モデル(full model)
1. 交互作用項、非線形項を含まない線形の全変数モデルを考える(lin_model)
用いる変数はである
lin_form<-paste("STROKE~",paste(full_var,collapse="+"))
lin_model<-glm(as.formula(lin_form),
data=fr.missf_p1,
family=binomial(link="logit"))VIFの検討
rms::vif(lin_model)## SEX2 CURSMOKE1 DIABETES1 BPMEDS1 PREVCHD1 PREVAP1 PREVMI1 PREVSTRK1
## 1.245566 2.439088 1.694071 1.107002 9.598226 6.717491 2.379642 1.000000
## PREVHYP1 TOTCHOL AGE SYSBP DIABP CIGPDAY BMI HEARTRTE
## 1.937207 1.067710 1.317387 3.700669 3.017425 2.568418 1.177581 1.100995
## GLUCOSE
## 1.723184- PREVCHD、PREVAPのVIFが高く、多重共線性が疑われるため、PREVCHDを除外して再度VIFを算出
glm(STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVAP +
PREVMI + PREVSTRK + PREVHYP + TOTCHOL + AGE + SYSBP + DIABP +
CIGPDAY + BMI + HEARTRTE + GLUCOSE,
data=fr.missf_p1,
family=binomial(link="logit")) %>%
rms::vif()## SEX2 CURSMOKE1 DIABETES1 BPMEDS1 PREVAP1 PREVMI1 PREVSTRK1 PREVHYP1
## 1.244364 2.439228 1.693826 1.106665 1.187531 1.164856 1.000000 1.935830
## TOTCHOL AGE SYSBP DIABP CIGPDAY BMI HEARTRTE GLUCOSE
## 1.067423 1.316767 3.695881 3.010815 2.568423 1.177247 1.100259 1.722415- 見事に各変数のVIFが低下したため、PREVCHDを除外したモデルを採用する(lin_model2)
lin_form2<-
as.formula(STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVAP +
PREVMI + PREVSTRK + PREVHYP + TOTCHOL + AGE + SYSBP + DIABP +
CIGPDAY + BMI + HEARTRTE + GLUCOSE)
lin_model2<-glm(lin_form2,
data=fr.missf_p1,
family=binomial(link="logit")) Backwards stepwiseによる変数選択
step(lin_model2,
direction="both") %>%
summary()## Start: AIC=2356.87
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVAP + PREVMI +
## PREVSTRK + PREVHYP + TOTCHOL + AGE + SYSBP + DIABP + CIGPDAY +
## BMI + HEARTRTE + GLUCOSE
##
## Df Deviance AIC
## - PREVAP 1 2322.9 2354.9
## - PREVMI 1 2322.9 2354.9
## - SYSBP 1 2322.9 2354.9
## - SEX 1 2323.1 2355.1
## - CIGPDAY 1 2323.4 2355.4
## - GLUCOSE 1 2323.5 2355.5
## - TOTCHOL 1 2323.6 2355.6
## - BMI 1 2324.6 2356.6
## <none> 2322.9 2356.9
## - BPMEDS 1 2326.4 2358.4
## - CURSMOKE 1 2326.5 2358.5
## - HEARTRTE 1 2327.9 2359.9
## - DIABETES 1 2329.9 2361.9
## - DIABP 1 2331.5 2363.5
## - PREVHYP 1 2333.0 2365.0
## - AGE 1 2412.6 2444.6
## - PREVSTRK 1 2438.7 2470.7
##
## Step: AIC=2354.87
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI + PREVSTRK +
## PREVHYP + TOTCHOL + AGE + SYSBP + DIABP + CIGPDAY + BMI +
## HEARTRTE + GLUCOSE
##
## Df Deviance AIC
## - PREVMI 1 2322.9 2352.9
## - SYSBP 1 2322.9 2352.9
## - SEX 1 2323.1 2353.1
## - CIGPDAY 1 2323.4 2353.4
## - GLUCOSE 1 2323.5 2353.5
## - TOTCHOL 1 2323.6 2353.6
## - BMI 1 2324.6 2354.6
## <none> 2322.9 2354.9
## - BPMEDS 1 2326.4 2356.4
## - CURSMOKE 1 2326.5 2356.5
## + PREVAP 1 2322.9 2356.9
## - HEARTRTE 1 2328.0 2358.0
## - DIABETES 1 2329.9 2359.9
## - DIABP 1 2331.5 2361.5
## - PREVHYP 1 2333.0 2363.0
## - AGE 1 2413.8 2443.8
## - PREVSTRK 1 2438.8 2468.8
##
## Step: AIC=2352.88
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVSTRK + PREVHYP +
## TOTCHOL + AGE + SYSBP + DIABP + CIGPDAY + BMI + HEARTRTE +
## GLUCOSE
##
## Df Deviance AIC
## - SYSBP 1 2322.9 2350.9
## - SEX 1 2323.1 2351.1
## - CIGPDAY 1 2323.4 2351.4
## - GLUCOSE 1 2323.5 2351.5
## - TOTCHOL 1 2323.6 2351.6
## - BMI 1 2324.6 2352.6
## <none> 2322.9 2352.9
## - BPMEDS 1 2326.4 2354.4
## - CURSMOKE 1 2326.5 2354.5
## + PREVMI 1 2322.9 2354.9
## + PREVAP 1 2322.9 2354.9
## - HEARTRTE 1 2328.0 2356.0
## - DIABETES 1 2329.9 2357.9
## - DIABP 1 2331.5 2359.5
## - PREVHYP 1 2333.0 2361.0
## - AGE 1 2414.8 2442.8
## - PREVSTRK 1 2439.4 2467.4
##
## Step: AIC=2350.94
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVSTRK + PREVHYP +
## TOTCHOL + AGE + DIABP + CIGPDAY + BMI + HEARTRTE + GLUCOSE
##
## Df Deviance AIC
## - SEX 1 2323.2 2349.2
## - CIGPDAY 1 2323.5 2349.5
## - GLUCOSE 1 2323.5 2349.5
## - TOTCHOL 1 2323.7 2349.7
## - BMI 1 2324.7 2350.7
## <none> 2322.9 2350.9
## - CURSMOKE 1 2326.6 2352.6
## - BPMEDS 1 2326.6 2352.6
## + SYSBP 1 2322.9 2352.9
## + PREVMI 1 2322.9 2352.9
## + PREVAP 1 2322.9 2352.9
## - HEARTRTE 1 2328.0 2354.0
## - DIABETES 1 2330.0 2356.0
## - PREVHYP 1 2334.9 2360.9
## - DIABP 1 2339.6 2365.6
## - AGE 1 2427.5 2453.5
## - PREVSTRK 1 2439.5 2465.5
##
## Step: AIC=2349.15
## STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVSTRK + PREVHYP +
## TOTCHOL + AGE + DIABP + CIGPDAY + BMI + HEARTRTE + GLUCOSE
##
## Df Deviance AIC
## - GLUCOSE 1 2323.7 2347.7
## - CIGPDAY 1 2323.9 2347.9
## - TOTCHOL 1 2324.1 2348.1
## - BMI 1 2324.9 2348.9
## <none> 2323.2 2349.2
## - BPMEDS 1 2326.7 2350.7
## - CURSMOKE 1 2326.8 2350.8
## + SEX 1 2322.9 2350.9
## + SYSBP 1 2323.1 2351.1
## + PREVMI 1 2323.1 2351.1
## + PREVAP 1 2323.2 2351.2
## - HEARTRTE 1 2328.5 2352.5
## - DIABETES 1 2330.3 2354.3
## - PREVHYP 1 2335.1 2359.1
## - DIABP 1 2339.9 2363.9
## - AGE 1 2427.9 2451.9
## - PREVSTRK 1 2439.8 2463.8
##
## Step: AIC=2347.72
## STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVSTRK + PREVHYP +
## TOTCHOL + AGE + DIABP + CIGPDAY + BMI + HEARTRTE
##
## Df Deviance AIC
## - CIGPDAY 1 2324.5 2346.5
## - TOTCHOL 1 2324.6 2346.6
## - BMI 1 2325.4 2347.4
## <none> 2323.7 2347.7
## + GLUCOSE 1 2323.2 2349.2
## - BPMEDS 1 2327.2 2349.2
## - CURSMOKE 1 2327.3 2349.3
## + SEX 1 2323.5 2349.5
## + PREVMI 1 2323.7 2349.7
## + SYSBP 1 2323.7 2349.7
## + PREVAP 1 2323.7 2349.7
## - HEARTRTE 1 2329.4 2351.4
## - DIABETES 1 2331.7 2353.7
## - PREVHYP 1 2335.6 2357.6
## - DIABP 1 2340.7 2362.7
## - AGE 1 2428.0 2450.0
## - PREVSTRK 1 2440.1 2462.1
##
## Step: AIC=2346.51
## STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVSTRK + PREVHYP +
## TOTCHOL + AGE + DIABP + BMI + HEARTRTE
##
## Df Deviance AIC
## - TOTCHOL 1 2325.4 2345.4
## - BMI 1 2326.2 2346.2
## <none> 2324.5 2346.5
## + CIGPDAY 1 2323.7 2347.7
## - BPMEDS 1 2327.9 2347.9
## + GLUCOSE 1 2323.9 2347.9
## + SEX 1 2324.1 2348.1
## + PREVMI 1 2324.5 2348.5
## + SYSBP 1 2324.5 2348.5
## + PREVAP 1 2324.5 2348.5
## - HEARTRTE 1 2330.0 2350.0
## - DIABETES 1 2332.4 2352.4
## - PREVHYP 1 2336.3 2356.3
## - CURSMOKE 1 2339.1 2359.1
## - DIABP 1 2341.6 2361.6
## - AGE 1 2428.0 2448.0
## - PREVSTRK 1 2440.6 2460.6
##
## Step: AIC=2345.44
## STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVSTRK + PREVHYP +
## AGE + DIABP + BMI + HEARTRTE
##
## Df Deviance AIC
## - BMI 1 2327.1 2345.1
## <none> 2325.4 2345.4
## + TOTCHOL 1 2324.5 2346.5
## + CIGPDAY 1 2324.6 2346.6
## - BPMEDS 1 2328.7 2346.7
## + SEX 1 2324.8 2346.8
## + GLUCOSE 1 2324.8 2346.8
## + PREVMI 1 2325.4 2347.4
## + PREVAP 1 2325.4 2347.4
## + SYSBP 1 2325.4 2347.4
## - HEARTRTE 1 2331.4 2349.4
## - DIABETES 1 2333.3 2351.3
## - PREVHYP 1 2337.1 2355.1
## - CURSMOKE 1 2340.3 2358.3
## - DIABP 1 2342.4 2360.4
## - AGE 1 2428.2 2446.2
## - PREVSTRK 1 2442.3 2460.3
##
## Step: AIC=2345.13
## STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVSTRK + PREVHYP +
## AGE + DIABP + HEARTRTE
##
## Df Deviance AIC
## <none> 2327.1 2345.1
## + BMI 1 2325.4 2345.4
## + TOTCHOL 1 2326.2 2346.2
## + CIGPDAY 1 2326.3 2346.3
## - BPMEDS 1 2330.3 2346.3
## + SEX 1 2326.5 2346.5
## + GLUCOSE 1 2326.6 2346.6
## + PREVMI 1 2327.1 2347.1
## + PREVAP 1 2327.1 2347.1
## + SYSBP 1 2327.1 2347.1
## - HEARTRTE 1 2333.0 2349.0
## - DIABETES 1 2335.8 2351.8
## - PREVHYP 1 2339.5 2355.5
## - CURSMOKE 1 2340.9 2356.9
## - DIABP 1 2347.3 2363.3
## - AGE 1 2429.3 2445.3
## - PREVSTRK 1 2443.8 2459.8##
## Call:
## glm(formula = STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVSTRK +
## PREVHYP + AGE + DIABP + HEARTRTE, family = binomial(link = "logit"),
## data = fr.missf_p1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3340 -0.4512 -0.3184 -0.2316 2.9566
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.767033 0.667549 -11.635 < 2e-16 ***
## CURSMOKE1 0.427568 0.115490 3.702 0.000214 ***
## DIABETES1 0.750248 0.241140 3.111 0.001863 **
## BPMEDS1 0.407887 0.222002 1.837 0.066163 .
## PREVSTRK1 17.571864 238.184327 0.074 0.941190
## PREVHYP1 0.516553 0.146211 3.533 0.000411 ***
## AGE 0.071498 0.007278 9.823 < 2e-16 ***
## DIABP 0.023741 0.005264 4.510 6.49e-06 ***
## HEARTRTE -0.011267 0.004715 -2.389 0.016872 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2756.0 on 4433 degrees of freedom
## Residual deviance: 2327.1 on 4425 degrees of freedom
## AIC: 2345.1
##
## Number of Fisher Scoring iterations: 14- Backwardsの結果、glm(formula = STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVSTRK + PREVHYP + AGE + DIABP + HEARTRTE, family = binomial(link = “logit”), data = fr.missf_p1) が最良のモデルとして選ばれた
stepwise_results<-
glm(formula = STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVSTRK +
PREVHYP + AGE + DIABP + HEARTRTE, family = binomial(link = "logit"),
data = fr.missf_p1) %>%
summary() %>%
.$coefficients %>%
as.data.frame() %>%
rownames_to_column() %>%
rename("var"="rowname") %>%
mutate(OR=exp(Estimate),
Lower=exp(Estimate-1.96*`Std. Error`),
Upper=exp(Estimate+1.96*`Std. Error`)) %>%
select(var,OR,Lower,Upper,"p"="Pr(>|z|)") %>%
#すべての変数を丸める
mutate_if(is.numeric,~round(.,3))
#結果を表示
stepwise_results## var OR Lower Upper p
## 1 (Intercept) 0.00000e+00 0.000 2.000000e-03 0.000
## 2 CURSMOKE1 1.53400e+00 1.223 1.923000e+00 0.000
## 3 DIABETES1 2.11800e+00 1.320 3.397000e+00 0.002
## 4 BPMEDS1 1.50400e+00 0.973 2.323000e+00 0.066
## 5 PREVSTRK1 4.27921e+07 0.000 2.387569e+210 0.941
## 6 PREVHYP1 1.67600e+00 1.259 2.233000e+00 0.000
## 7 AGE 1.07400e+00 1.059 1.090000e+00 0.000
## 8 DIABP 1.02400e+00 1.014 1.035000e+00 0.000
## 9 HEARTRTE 9.89000e-01 0.980 9.980000e-01 0.017#plot
stepwise_results %>%
#filter(!var=="PREVSTRK1") %>% #これだけ非常にCIが大きくグラフ化すると視認性が悪いので一時的に削除
ggplot(aes(y=var,x=OR))+
geom_point()+
geom_errorbar(aes(xmin=Lower,xmax=Upper),width=0.2)+
geom_vline(xintercept=1,color="red")
-
変数選択の結果残った変数のうち、PREVSTRKのORが非常に大きくなり、推定が不安定になっている
- アウトカムとPREVSTROKEの関係を詳しく見てみると、分割表において、 4019,
383, 0, 32
であり、STROKE=0のとき、PREVSTRK=1に該当する患者が0であることがわかった→これがsparse
biasの原因である
xtabs(~STROKE+PREVSTRK,
data=fr.missf_p1)## PREVSTRK
## STROKE 0 1
## 0 4019 0
## 1 383 32ここで、すべてのカテゴリカル変数に対して同様にアウトカムとの分割表を考えてみる
f<-
function(x){
form<-as.formula(paste("~STROKE+",x))
xtabs(form,data=fr.missf_p1)
}
for(i in fac_var[-1]){
print(f(i))
}## SEX
## STROKE 1 2
## 0 1751 2268
## 1 193 222
## CURSMOKE
## STROKE 0 1
## 0 2037 1982
## 1 216 199
## DIABETES
## STROKE 0 1
## 0 3926 93
## 1 387 28
## BPMEDS
## STROKE 0 1
## 0 3917 102
## 1 373 42
## educ
## STROKE 1 2 3 4
## 0 1659 1206 676 478
## 1 224 104 52 35
## PREVCHD
## STROKE 0 1
## 0 3862 157
## 1 378 37
## PREVAP
## STROKE 0 1
## 0 3899 120
## 1 388 27
## PREVMI
## STROKE 0 1
## 0 3949 70
## 1 399 16
## PREVSTRK
## STROKE 0 1
## 0 4019 0
## 1 383 32
## PREVHYP
## STROKE 0 1
## 0 2840 1179
## 1 164 251
## PERIOD
## STROKE 1 2 3
## 0 4019 0 0
## 1 415 0 0
## DEATH
## STROKE 0 1
## 0 2746 1273
## 1 138 277
## ANGINA
## STROKE 0 1
## 0 3394 625
## 1 315 100
## HOSPMI
## STROKE 0 1
## 0 3636 383
## 1 344 71
## MI_FCHD
## STROKE 0 1
## 0 3398 621
## 1 305 110
## ANYCHD
## STROKE 0 1
## 0 2958 1061
## 1 236 179
## STROKE
## 0 1
## 4019 415
## CVD
## STROKE 0 1
## 0 3277 742
## 1 0 415
## HYPERTEN
## STROKE 0 1
## 0 1142 2877
## 1 40 375- カテゴリカル変数に対してアウトカムとの分割表を考えてみると、sparseな部分が可視化できるため、事前の変数選択を考える上で非常に有用であることがわかった
- 今回は、PREVSTRKのみが完全にsparse(cell=0がある)であるため、この変数を除外して再度変数選択を行うこととする（本当はこれを解析前にできればよかった）
PREVSTROKEを除外したモデルにおいて、VIFから再度検討する
glm(
STROKE~ SEX+CURSMOKE+DIABETES+BPMEDS+PREVCHD+PREVAP+PREVMI+PREVHYP+TOTCHOL+AGE+SYSBP+DIABP+CIGPDAY+BMI+HEARTRTE+GLUCOSE,
data=fr.missf_p1,
family=binomial(link="logit")) %>%
rms::vif()## SEX2 CURSMOKE1 DIABETES1 BPMEDS1 PREVCHD1 PREVAP1 PREVMI1 PREVHYP1
## 1.240065 2.456395 1.705096 1.112819 8.979822 5.835302 2.735766 1.922086
## TOTCHOL AGE SYSBP DIABP CIGPDAY BMI HEARTRTE GLUCOSE
## 1.066603 1.317287 3.664247 2.980477 2.583968 1.182116 1.101576 1.730470#ココでも同様にPREVCHDのVIFが高いため、除外して再度VIFを算出
glm(
STROKE~ SEX+CURSMOKE+DIABETES+BPMEDS+PREVAP+PREVMI+PREVHYP+TOTCHOL+AGE+SYSBP+DIABP+CIGPDAY+BMI+HEARTRTE+GLUCOSE,
data=fr.missf_p1,
family=binomial(link="logit")) %>%
rms::vif()## SEX2 CURSMOKE1 DIABETES1 BPMEDS1 PREVAP1 PREVMI1 PREVHYP1 TOTCHOL
## 1.240237 2.456599 1.705139 1.112127 1.179417 1.159379 1.919544 1.065887
## AGE SYSBP DIABP CIGPDAY BMI HEARTRTE GLUCOSE
## 1.316273 3.658322 2.971481 2.583932 1.181654 1.100465 1.729751#いい感じ→これを初期モデルとする
lin_model3<-
glm(STROKE~ SEX+CURSMOKE+DIABETES+BPMEDS+PREVAP+PREVMI+PREVHYP+TOTCHOL+AGE+SYSBP+DIABP+CIGPDAY+BMI+HEARTRTE+GLUCOSE,
data=fr.missf_p1,
family=binomial(link="logit"))step(lin_model3,
direction="both") %>%
summary()## Start: AIC=2470.74
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVAP + PREVMI +
## PREVHYP + TOTCHOL + AGE + SYSBP + DIABP + CIGPDAY + BMI +
## HEARTRTE + GLUCOSE
##
## Df Deviance AIC
## - PREVAP 1 2438.8 2468.8
## - SYSBP 1 2438.8 2468.8
## - SEX 1 2439.0 2469.0
## - CIGPDAY 1 2439.1 2469.1
## - GLUCOSE 1 2439.1 2469.1
## - PREVMI 1 2439.4 2469.4
## - TOTCHOL 1 2440.2 2470.2
## - BMI 1 2440.3 2470.3
## <none> 2438.7 2470.7
## - CURSMOKE 1 2442.3 2472.3
## - DIABETES 1 2445.2 2475.2
## - HEARTRTE 1 2445.4 2475.4
## - BPMEDS 1 2445.7 2475.7
## - DIABP 1 2446.6 2476.6
## - PREVHYP 1 2453.2 2483.2
## - AGE 1 2532.4 2562.4
##
## Step: AIC=2468.82
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI + PREVHYP +
## TOTCHOL + AGE + SYSBP + DIABP + CIGPDAY + BMI + HEARTRTE +
## GLUCOSE
##
## Df Deviance AIC
## - SYSBP 1 2438.9 2466.9
## - SEX 1 2439.1 2467.1
## - CIGPDAY 1 2439.1 2467.1
## - GLUCOSE 1 2439.1 2467.1
## - PREVMI 1 2439.4 2467.4
## - TOTCHOL 1 2440.3 2468.3
## - BMI 1 2440.4 2468.4
## <none> 2438.8 2468.8
## - CURSMOKE 1 2442.4 2470.4
## + PREVAP 1 2438.7 2470.7
## - DIABETES 1 2445.2 2473.2
## - HEARTRTE 1 2445.5 2473.5
## - BPMEDS 1 2445.7 2473.7
## - DIABP 1 2446.7 2474.7
## - PREVHYP 1 2453.3 2481.3
## - AGE 1 2533.2 2561.2
##
## Step: AIC=2466.91
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI + PREVHYP +
## TOTCHOL + AGE + DIABP + CIGPDAY + BMI + HEARTRTE + GLUCOSE
##
## Df Deviance AIC
## - SEX 1 2439.1 2465.1
## - GLUCOSE 1 2439.2 2465.2
## - CIGPDAY 1 2439.2 2465.2
## - PREVMI 1 2439.5 2465.5
## - TOTCHOL 1 2440.4 2466.4
## - BMI 1 2440.5 2466.5
## <none> 2438.9 2466.9
## - CURSMOKE 1 2442.5 2468.5
## + SYSBP 1 2438.8 2468.8
## + PREVAP 1 2438.8 2468.8
## - DIABETES 1 2445.4 2471.4
## - HEARTRTE 1 2445.6 2471.6
## - BPMEDS 1 2446.0 2472.0
## - DIABP 1 2454.4 2480.4
## - PREVHYP 1 2456.1 2482.1
## - AGE 1 2546.7 2572.7
##
## Step: AIC=2465.12
## STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVMI + PREVHYP + TOTCHOL +
## AGE + DIABP + CIGPDAY + BMI + HEARTRTE + GLUCOSE
##
## Df Deviance AIC
## - GLUCOSE 1 2439.4 2463.4
## - CIGPDAY 1 2439.6 2463.6
## - PREVMI 1 2439.8 2463.8
## - BMI 1 2440.7 2464.7
## - TOTCHOL 1 2440.8 2464.8
## <none> 2439.1 2465.1
## - CURSMOKE 1 2442.7 2466.7
## + SEX 1 2438.9 2466.9
## + PREVAP 1 2439.1 2467.1
## + SYSBP 1 2439.1 2467.1
## - DIABETES 1 2445.7 2469.7
## - BPMEDS 1 2446.1 2470.1
## - HEARTRTE 1 2446.2 2470.2
## - DIABP 1 2454.8 2478.8
## - PREVHYP 1 2456.2 2480.2
## - AGE 1 2546.9 2570.9
##
## Step: AIC=2463.41
## STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVMI + PREVHYP + TOTCHOL +
## AGE + DIABP + CIGPDAY + BMI + HEARTRTE
##
## Df Deviance AIC
## - CIGPDAY 1 2439.9 2461.9
## - PREVMI 1 2440.1 2462.1
## - BMI 1 2441.0 2463.0
## - TOTCHOL 1 2441.1 2463.1
## <none> 2439.4 2463.4
## - CURSMOKE 1 2443.0 2465.0
## + GLUCOSE 1 2439.1 2465.1
## + SEX 1 2439.2 2465.2
## + PREVAP 1 2439.4 2465.4
## + SYSBP 1 2439.4 2465.4
## - BPMEDS 1 2446.3 2468.3
## - HEARTRTE 1 2446.7 2468.7
## - DIABETES 1 2447.4 2469.4
## - DIABP 1 2455.2 2477.2
## - PREVHYP 1 2456.5 2478.5
## - AGE 1 2547.0 2569.0
##
## Step: AIC=2461.93
## STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVMI + PREVHYP + TOTCHOL +
## AGE + DIABP + BMI + HEARTRTE
##
## Df Deviance AIC
## - PREVMI 1 2440.6 2460.6
## - BMI 1 2441.5 2461.5
## - TOTCHOL 1 2441.6 2461.6
## <none> 2439.9 2461.9
## + CIGPDAY 1 2439.4 2463.4
## + SEX 1 2439.5 2463.5
## + GLUCOSE 1 2439.6 2463.6
## + PREVAP 1 2439.9 2463.9
## + SYSBP 1 2439.9 2463.9
## - BPMEDS 1 2446.7 2466.7
## - HEARTRTE 1 2447.1 2467.1
## - DIABETES 1 2447.9 2467.9
## - CURSMOKE 1 2453.0 2473.0
## - DIABP 1 2455.8 2475.8
## - PREVHYP 1 2456.9 2476.9
## - AGE 1 2547.0 2567.0
##
## Step: AIC=2460.57
## STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVHYP + TOTCHOL + AGE +
## DIABP + BMI + HEARTRTE
##
## Df Deviance AIC
## - BMI 1 2442.1 2460.1
## - TOTCHOL 1 2442.3 2460.3
## <none> 2440.6 2460.6
## + PREVMI 1 2439.9 2461.9
## + CIGPDAY 1 2440.1 2462.1
## + SEX 1 2440.1 2462.1
## + GLUCOSE 1 2440.3 2462.3
## + SYSBP 1 2440.6 2462.6
## + PREVAP 1 2440.6 2462.6
## - BPMEDS 1 2447.4 2465.4
## - HEARTRTE 1 2447.8 2465.8
## - DIABETES 1 2448.8 2466.8
## - CURSMOKE 1 2454.1 2472.1
## - DIABP 1 2456.3 2474.3
## - PREVHYP 1 2457.8 2475.8
## - AGE 1 2550.6 2568.6
##
## Step: AIC=2460.12
## STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVHYP + TOTCHOL + AGE +
## DIABP + HEARTRTE
##
## Df Deviance AIC
## - TOTCHOL 1 2443.8 2459.8
## <none> 2442.1 2460.1
## + BMI 1 2440.6 2460.6
## + PREVMI 1 2441.5 2461.5
## + CIGPDAY 1 2441.6 2461.6
## + SEX 1 2441.6 2461.6
## + GLUCOSE 1 2441.9 2461.9
## + SYSBP 1 2442.1 2462.1
## + PREVAP 1 2442.1 2462.1
## - BPMEDS 1 2449.0 2465.0
## - HEARTRTE 1 2449.2 2465.2
## - DIABETES 1 2451.1 2467.1
## - CURSMOKE 1 2454.6 2470.6
## - PREVHYP 1 2460.2 2476.2
## - DIABP 1 2461.0 2477.0
## - AGE 1 2551.7 2567.7
##
## Step: AIC=2459.81
## STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVHYP + AGE + DIABP +
## HEARTRTE
##
## Df Deviance AIC
## <none> 2443.8 2459.8
## + TOTCHOL 1 2442.1 2460.1
## + BMI 1 2442.3 2460.3
## + SEX 1 2443.0 2461.0
## + PREVMI 1 2443.2 2461.2
## + CIGPDAY 1 2443.2 2461.2
## + GLUCOSE 1 2443.6 2461.6
## + SYSBP 1 2443.8 2461.8
## + PREVAP 1 2443.8 2461.8
## - BPMEDS 1 2450.4 2464.4
## - HEARTRTE 1 2451.6 2465.6
## - DIABETES 1 2452.7 2466.7
## - CURSMOKE 1 2456.7 2470.7
## - PREVHYP 1 2461.6 2475.6
## - DIABP 1 2462.5 2476.5
## - AGE 1 2551.7 2565.7##
## Call:
## glm(formula = STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVHYP +
## AGE + DIABP + HEARTRTE, family = binomial(link = "logit"),
## data = fr.missf_p1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3955 -0.4649 -0.3255 -0.2367 2.9230
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.459839 0.645143 -11.563 < 2e-16 ***
## CURSMOKE1 0.401439 0.112012 3.584 0.000339 ***
## DIABETES1 0.739673 0.235375 3.143 0.001675 **
## BPMEDS1 0.547228 0.207449 2.638 0.008342 **
## PREVHYP1 0.597373 0.141426 4.224 2.40e-05 ***
## AGE 0.071244 0.007065 10.085 < 2e-16 ***
## DIABP 0.021968 0.005071 4.332 1.48e-05 ***
## HEARTRTE -0.012629 0.004590 -2.751 0.005933 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2756.0 on 4433 degrees of freedom
## Residual deviance: 2443.8 on 4426 degrees of freedom
## AIC: 2459.8
##
## Number of Fisher Scoring iterations: 6-Backwardsで最良のモデルとして選ばれたのは glm(formula = STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVHYP + AGE + DIABP + HEARTRTE, family = binomial(link = “logit”), data = fr.missf_p1)
stepwise_results<-
glm(formula = STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVHYP +
AGE + DIABP + HEARTRTE, family = binomial(link = "logit"),
data = fr.missf_p1) %>%
summary() %>%
.$coefficients %>%
as.data.frame() %>%
rownames_to_column() %>%
rename("var"="rowname") %>%
mutate(OR=exp(Estimate),
Lower=exp(Estimate-1.96*`Std. Error`),
Upper=exp(Estimate+1.96*`Std. Error`)) %>%
select(var,OR,Lower,Upper,"p"="Pr(>|z|)") %>%
#すべての変数を丸める
mutate_if(is.numeric,~round(.,3))
#結果を表示
stepwise_results## var OR Lower Upper p
## 1 (Intercept) 0.001 0.000 0.002 0.000
## 2 CURSMOKE1 1.494 1.199 1.861 0.000
## 3 DIABETES1 2.095 1.321 3.323 0.002
## 4 BPMEDS1 1.728 1.151 2.596 0.008
## 5 PREVHYP1 1.817 1.377 2.398 0.000
## 6 AGE 1.074 1.059 1.089 0.000
## 7 DIABP 1.022 1.012 1.032 0.000
## 8 HEARTRTE 0.987 0.979 0.996 0.006#plot
stepwise_results %>%
ggplot(aes(y=var,x=OR))+
geom_point()+
geom_errorbar(aes(xmin=Lower,xmax=Upper),width=0.2)+
geom_vline(xintercept=1,color="red")
- ORが発散することなくとてもいい感じの推定ができたので、このモデルを 1.
交互作用項、非線形項を含まない線形の全変数モデル(lin_model)における最終モデル(lin_final_model)とする
lin_final_formula<-as.formula(STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVHYP + AGE + DIABP + HEARTRTE)
lin_final_model<-
glm(formula = lin_final_formula,
family = binomial(link = "logit"),
#x=T,y=T,
data = fr.missf_p1)予測性能の評価
ROC AUC
# risk prediction fitting
fr.missf_p1$fitted.lin.fin<-
predict(lin_final_model,
type="response")
# ROC AUC
roc.lin.final<-
roc(STROKE~fitted.lin.fin,
data=fr.missf_p1,
ci=T)## Setting levels: control = 0, case = 1## Setting direction: controls < casesroc.lin.final##
## Call:
## roc.formula(formula = STROKE ~ fitted.lin.fin, data = fr.missf_p1, ci = T)
##
## Data: fitted.lin.fin in 4019 controls (STROKE 0) < 415 cases (STROKE 1).
## Area under the curve: 0.7504
## 95% CI: 0.7265-0.7742 (DeLong)roc.lin.final %>%
ggroc(legacy.axes=T)+
geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="grey", linetype="dashed")
- ROCは0.7503518であり、中等度の予測性能と考えられる
#calibration plot
cal.lin.final<-
val.prob(
p=fr.missf_p1$fitted.lin.fin,
y=fr.missf_p1$STROKEn,
g=10,
cex=0.5) 
- calibration plotからみても予測性能は高い
Bootstrapによるliner final modelの内的妥当性検証
optimismの算出
B<-100 #number of resampling
N<-dim(fr.missf_p1)[1] #dataのdimentionを返す→4434行x39列→そのうち一番目（4434)を取得しNとする
AUC.i1<-AUC.i2<-numeric(B)
for(i in 1:B){
#bootstrap random sampling
bs.i<-sample(1:N,N,replace=TRUE) #重複を許して1~4434までの数を4434回サンプリングする
#bootstrap samplingに基づき、dataから集団を再構成 →data.iとする
fr.missf_p1.i<-fr.missf_p1[bs.i,]
#bootstrapしたデータセットごとで回帰モデル構築→各bootstrap dataからリスク予測確率の算出→data2.1のprob_1列に格納
lin.fin.cv.res.i<-
glm(formula = lin_final_formula,
family = binomial(link = "logit"),
data = fr.missf_p1.i)
fr.missf_p1.i$lin.fin.1<- #lin.fin.1では、bootstrap dataから構築した新しいロジスティク回帰モデルを使い、そのbootstrap dataに対するリスク予測確率を格納
predict(lin.fin.cv.res.i,
data=fr.missf_p1.i,
type="response")
#AUCの算出；Bootstrap回数（B=100)分のAUCを算出して結果を格納；各bootstrap dataに対するAUC
AUC.i1[i]<-roc(STROKE~lin.fin.1,data=fr.missf_p1.i)$auc
#bootstrapで作成した各回帰モデルを使って、オリジナルデータ(fr.missf_p1)に対するリスク予測確率を算出→fr.missf_p1のlin.fin.2列に格納;AUC算出
fr.missf_p1$lin.fin.2<-
predict(lin.fin.cv.res.i,
newdata=fr.missf_p1, ##dataではなくnewdataにする必要がある！！！！
type="response")
AUC.i2[i]<-roc(STROKE~lin.fin.2,data=fr.missf_p1)$auc #ここで得られるAUCはbootstrap dataを用いてoriginal dataに外挿したときのAUC
#print(paste(i,"th bootstrap iteration is completed.",sep=""))
}
opt2<-AUC.i1-AUC.i2
#(bootstrap samplingした各data setを使い、それぞれ導いた回帰モデルにおけるAUC)ー(bootstrap samplingした各data setを使い、original dataに外挿したときのAUC)
summary(opt2);hist(opt2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.024696 -0.005136 0.006434 0.004972 0.014323 0.033388
optimismを用いたbias corrected AUCの算出
lam2<-mean(opt2) #estimate of the optimism
lin.cor.AUC<-roc.lin.final$auc-lam2 #bias corrected AUC estimate
lin.cor.AUC## [1] 0.74538- Linerity final modelによる予測モデル構築を行い、Bootstrapによるbias corrected AUCを算出したところ、0.74538であった
- 当初、２個目のpredictのdata部分をnewdataにしていなかったためoptimismが非常に大きく出てしまっていた。回帰モデルの内部におけるdataとpredictにおける外挿したいdataが異なるときは、newdataに設定する必要がある
Cross validationによる内的妥当性の検証
#k-fold CV
k <-5
AUC_CV <-data.frame(matrix(ncol=2,nrow=k))
lin_final_formula<-as.formula(STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVHYP + AGE + DIABP + HEARTRTE)
lin_final_model<-
glm(formula = lin_final_formula,
family = binomial(link = "logit"),
data = fr.missf_p1)
#Randomly shuffle the data
fr.missf_p1.r<-fr.missf_p1[sample(nrow(fr.missf_p1)),]
#Create k equally size folds
folds <- cut(seq(1,nrow(fr.missf_p1.r)),breaks=k,labels=FALSE)
#Perform k fold cross validation
for(i in 1:k){
#Segment data2r by fold using the which() function
testIndexes <- which(folds==i,arr.ind=TRUE)
testData <- fr.missf_p1.r[testIndexes, ]
trainData <- fr.missf_p1.r[-testIndexes, ]
fit1<-glm(formula = lin_final_formula, family = binomial(link = "logit"), data = trainData)
testData$fitted<-predict(fit1,testData,type="response")
ROC <- roc(testData$STROKE, testData$fitted)
AUC_CV[i,2] <- ROC$auc
AUC_CV[i,1] <- i
#print(paste(i,"th cross-validation iteration is completed.", sep=""))
} ## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < casesnames(AUC_CV)[2] <- "AUC_cv"
names(AUC_CV)[1] <- "k"
print(summary(AUC_CV$AUC_cv)) ## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.7312 0.7332 0.7342 0.7463 0.7659 0.7671- CVではAUCの平均値は0.7463201
2. 交互作用項、非線形項を含む線形の全変数モデルを考える(full_model)
用いる変数はである
- full modelのうち、前述で多重共線性やsparseの問題から削除した変数（PREVSTRK,PREVCHD)を除いたモデルを考える
fullmodel## [1] "STROKE~SEX+CURSMOKE+DIABETES+BPMEDS+PREVCHD+PREVAP+PREVMI+PREVHYP+TOTCHOL+AGE+SYSBP+DIABP+CIGPDAY+rcs(BMI,4)+rcs(HEARTRTE,4)+rcs(GLUCOSE,4)+SEX:PREVCHD+SEX:PREVAP+SEX:PREVMI+SEX:PREVHYP+SEX:DIABP+CURSMOKE:HEARTRTE+DIABETES:TOTCHOL+DIABETES:SYSBP+BPMEDS:PREVCHD+BPMEDS:TOTCHOL+BPMEDS:AGE+PREVCHD:GLUCOSE+PREVHYP:BMI+TOTCHOL:AGE+HEARTRTE:GLUCOSE"glm(as.formula(fullmodel),data=fr.missf_p1,
family = binomial(link = "logit")) ##
## Call: glm(formula = as.formula(fullmodel), family = binomial(link = "logit"),
## data = fr.missf_p1)
##
## Coefficients:
## (Intercept) SEX2
## -5.2535723 -1.1953373
## CURSMOKE1 DIABETES1
## -1.2365073 4.0207549
## BPMEDS1 PREVCHD1
## 5.2898809 1.1035782
## PREVAP1 PREVMI1
## -0.3720497 0.0550895
## PREVHYP1 TOTCHOL
## 1.5148234 0.0066005
## AGE SYSBP
## 0.1010951 0.0013414
## DIABP CIGPDAY
## 0.0127927 0.0026360
## rcs(BMI, 4)BMI rcs(BMI, 4)BMI'
## -0.0862948 0.2278817
## rcs(BMI, 4)BMI'' rcs(HEARTRTE, 4)HEARTRTE
## -0.5949000 0.0094925
## rcs(HEARTRTE, 4)HEARTRTE' rcs(HEARTRTE, 4)HEARTRTE''
## -0.2123066 0.5429251
## rcs(GLUCOSE, 4)GLUCOSE rcs(GLUCOSE, 4)GLUCOSE'
## -0.0611513 0.1259855
## rcs(GLUCOSE, 4)GLUCOSE'' SEX2:PREVCHD1
## -0.3977397 -0.4496324
## SEX2:PREVAP1 SEX2:PREVMI1
## 0.9761805 1.3096815
## SEX2:PREVHYP1 SEX2:DIABP
## 0.0700994 0.0122176
## CURSMOKE0:HEARTRTE CURSMOKE1:HEARTRTE
## -0.0207419 NA
## DIABETES1:TOTCHOL DIABETES1:SYSBP
## -0.0066810 -0.0112544
## BPMEDS1:PREVCHD1 BPMEDS1:TOTCHOL
## -0.8348873 -0.0085763
## BPMEDS1:AGE PREVCHD0:GLUCOSE
## -0.0439211 0.0103560
## PREVCHD1:GLUCOSE PREVHYP0:BMI
## NA 0.0374607
## PREVHYP1:BMI TOTCHOL:AGE
## NA -0.0001260
## HEARTRTE:GLUCOSE
## 0.0003235
##
## Degrees of Freedom: 4433 Total (i.e. Null); 4396 Residual
## Null Deviance: 2756
## Residual Deviance: 2395 AIC: 2471fullmodel2<-
STROKE~SEX+CURSMOKE+DIABETES+BPMEDS+PREVAP+PREVMI+PREVHYP+TOTCHOL+AGE+SYSBP+DIABP+CIGPDAY+
rcs(BMI,4)+rcs(HEARTRTE,4)+rcs(GLUCOSE,4)+
SEX:PREVAP+SEX:PREVMI+SEX:PREVHYP+SEX:DIABP+CURSMOKE:HEARTRTE+DIABETES:TOTCHOL+DIABETES:SYSBP+BPMEDS:TOTCHOL+
BPMEDS:AGE+PREVHYP:BMI+TOTCHOL:AGE+HEARTRTE:GLUCOSE
fullmodel3<-STROKE~SEX+CURSMOKE+DIABETES+BPMEDS+PREVAP+PREVMI+PREVHYP+TOTCHOL+AGE+SYSBP+DIABP+CIGPDAY+
rcs(BMI,4)+rcs(HEARTRTE,4)+rcs(GLUCOSE,4)+SEX:PREVHYP+SEX:PREVMIVIFの評価
#dfが大きいのでGVIFを算出するcar::vifを用いる
glm(as.formula(fullmodel3),data=fr.missf_p1,
family = binomial(link = "logit")) %>%
car::vif(type="predictor")## Warning in vif.lm(., type = "predictor"): type = 'predictor' is available only for unweighted linear models;
## type = 'terms' will be used## GVIF Df GVIF^(1/(2*Df))
## SEX 2.570312 1 1.603219
## CURSMOKE 2.463038 1 1.569407
## DIABETES 1.751968 1 1.323619
## BPMEDS 1.134620 1 1.065186
## PREVAP 1.204167 1 1.097345
## PREVMI 1.687747 1 1.299133
## PREVHYP 2.854402 1 1.689498
## TOTCHOL 1.089812 1 1.043941
## AGE 1.329260 1 1.152935
## SYSBP 3.727491 1 1.930671
## DIABP 2.981403 1 1.726674
## CIGPDAY 2.594084 1 1.610616
## rcs(BMI, 4) 1.323785 3 1.047859
## rcs(HEARTRTE, 4) 1.150514 3 1.023643
## rcs(GLUCOSE, 4) 1.866674 3 1.109630
## SEX:PREVHYP 3.811174 1 1.952223
## SEX:PREVMI 1.557143 1 1.247855- fullmodel2では交互作用項が多すぎるためか計算不能になるので、臨床的な（感覚的に）視点から交互作用項を減らしたmodel3を考える
- fullmodel3ではGVIFの問題は気にならないためこのまま進める
Backwards stepwiseによる変数選択
glm(as.formula(fullmodel3),data=fr.missf_p1,
family = binomial(link = "logit")) %>%
step(direction="both") %>%
summary()## Start: AIC=2472.26
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVAP + PREVMI +
## PREVHYP + TOTCHOL + AGE + SYSBP + DIABP + CIGPDAY + rcs(BMI,
## 4) + rcs(HEARTRTE, 4) + rcs(GLUCOSE, 4) + SEX:PREVHYP + SEX:PREVMI
##
## Df Deviance AIC
## - rcs(BMI, 4) 3 2426.7 2468.7
## - PREVAP 1 2424.3 2470.3
## - SYSBP 1 2424.3 2470.3
## - CIGPDAY 1 2424.5 2470.5
## - rcs(GLUCOSE, 4) 3 2428.9 2470.9
## - SEX:PREVHYP 1 2425.2 2471.2
## - TOTCHOL 1 2425.4 2471.4
## <none> 2424.3 2472.3
## - SEX:PREVMI 1 2426.5 2472.5
## - CURSMOKE 1 2428.0 2474.0
## - DIABETES 1 2429.0 2475.0
## - BPMEDS 1 2430.7 2476.7
## - rcs(HEARTRTE, 4) 3 2436.6 2478.6
## - DIABP 1 2433.4 2479.4
## - AGE 1 2515.0 2561.0
##
## Step: AIC=2468.74
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVAP + PREVMI +
## PREVHYP + TOTCHOL + AGE + SYSBP + DIABP + CIGPDAY + rcs(HEARTRTE,
## 4) + rcs(GLUCOSE, 4) + SEX:PREVHYP + SEX:PREVMI
##
## Df Deviance AIC
## - SYSBP 1 2426.7 2466.7
## - PREVAP 1 2426.7 2466.7
## - CIGPDAY 1 2427.0 2467.0
## - rcs(GLUCOSE, 4) 3 2431.6 2467.6
## - SEX:PREVHYP 1 2427.9 2467.9
## - TOTCHOL 1 2428.1 2468.1
## <none> 2426.7 2468.7
## - SEX:PREVMI 1 2429.0 2469.0
## - CURSMOKE 1 2430.2 2470.2
## - DIABETES 1 2432.0 2472.0
## + rcs(BMI, 4) 3 2424.3 2472.3
## - BPMEDS 1 2433.3 2473.3
## - rcs(HEARTRTE, 4) 3 2438.9 2474.9
## - DIABP 1 2437.6 2477.6
## - AGE 1 2517.8 2557.8
##
## Step: AIC=2466.74
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVAP + PREVMI +
## PREVHYP + TOTCHOL + AGE + DIABP + CIGPDAY + rcs(HEARTRTE,
## 4) + rcs(GLUCOSE, 4) + SEX:PREVHYP + SEX:PREVMI
##
## Df Deviance AIC
## - PREVAP 1 2426.7 2464.7
## - CIGPDAY 1 2427.0 2465.0
## - rcs(GLUCOSE, 4) 3 2431.6 2465.6
## - SEX:PREVHYP 1 2427.9 2465.9
## - TOTCHOL 1 2428.1 2466.1
## <none> 2426.7 2466.7
## - SEX:PREVMI 1 2429.0 2467.0
## - CURSMOKE 1 2430.2 2468.2
## + SYSBP 1 2426.7 2468.7
## - DIABETES 1 2432.0 2470.0
## + rcs(BMI, 4) 3 2424.3 2470.3
## - BPMEDS 1 2433.3 2471.3
## - rcs(HEARTRTE, 4) 3 2438.9 2472.9
## - DIABP 1 2445.8 2483.8
## - AGE 1 2527.5 2565.5
##
## Step: AIC=2464.74
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI + PREVHYP +
## TOTCHOL + AGE + DIABP + CIGPDAY + rcs(HEARTRTE, 4) + rcs(GLUCOSE,
## 4) + SEX:PREVHYP + SEX:PREVMI
##
## Df Deviance AIC
## - CIGPDAY 1 2427.0 2463.0
## - rcs(GLUCOSE, 4) 3 2431.6 2463.6
## - SEX:PREVHYP 1 2427.9 2463.9
## - TOTCHOL 1 2428.1 2464.1
## <none> 2426.7 2464.7
## - SEX:PREVMI 1 2429.0 2465.0
## - CURSMOKE 1 2430.2 2466.2
## + PREVAP 1 2426.7 2466.7
## + SYSBP 1 2426.7 2466.7
## - DIABETES 1 2432.0 2468.0
## + rcs(BMI, 4) 3 2424.3 2468.3
## - BPMEDS 1 2433.4 2469.4
## - rcs(HEARTRTE, 4) 3 2438.9 2470.9
## - DIABP 1 2445.8 2481.8
## - AGE 1 2529.1 2565.1
##
## Step: AIC=2462.99
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI + PREVHYP +
## TOTCHOL + AGE + DIABP + rcs(HEARTRTE, 4) + rcs(GLUCOSE, 4) +
## SEX:PREVHYP + SEX:PREVMI
##
## Df Deviance AIC
## - rcs(GLUCOSE, 4) 3 2431.9 2461.9
## - SEX:PREVHYP 1 2428.2 2462.2
## - TOTCHOL 1 2428.3 2462.3
## <none> 2427.0 2463.0
## - SEX:PREVMI 1 2429.3 2463.3
## + CIGPDAY 1 2426.7 2464.7
## + PREVAP 1 2427.0 2465.0
## + SYSBP 1 2427.0 2465.0
## - DIABETES 1 2432.2 2466.2
## + rcs(BMI, 4) 3 2424.5 2466.5
## - BPMEDS 1 2433.6 2467.6
## - rcs(HEARTRTE, 4) 3 2439.0 2469.0
## - CURSMOKE 1 2437.6 2471.6
## - DIABP 1 2446.1 2480.1
## - AGE 1 2529.2 2563.2
##
## Step: AIC=2461.88
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI + PREVHYP +
## TOTCHOL + AGE + DIABP + rcs(HEARTRTE, 4) + SEX:PREVHYP +
## SEX:PREVMI
##
## Df Deviance AIC
## - TOTCHOL 1 2433.4 2461.4
## - SEX:PREVHYP 1 2433.4 2461.4
## <none> 2431.9 2461.9
## - SEX:PREVMI 1 2433.9 2461.9
## + rcs(GLUCOSE, 4) 3 2427.0 2463.0
## + CIGPDAY 1 2431.6 2463.6
## + PREVAP 1 2431.9 2463.9
## + SYSBP 1 2431.9 2463.9
## + rcs(BMI, 4) 3 2429.1 2465.1
## - BPMEDS 1 2438.7 2466.7
## - DIABETES 1 2438.9 2466.9
## - rcs(HEARTRTE, 4) 3 2443.6 2467.6
## - CURSMOKE 1 2442.6 2470.6
## - DIABP 1 2450.7 2478.7
## - AGE 1 2536.6 2564.6
##
## Step: AIC=2461.38
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI + PREVHYP +
## AGE + DIABP + rcs(HEARTRTE, 4) + SEX:PREVHYP + SEX:PREVMI
##
## Df Deviance AIC
## - SEX:PREVHYP 1 2434.8 2460.8
## <none> 2433.4 2461.4
## - SEX:PREVMI 1 2435.4 2461.4
## + TOTCHOL 1 2431.9 2461.9
## + rcs(GLUCOSE, 4) 3 2428.3 2462.3
## + CIGPDAY 1 2433.2 2463.2
## + PREVAP 1 2433.4 2463.4
## + SYSBP 1 2433.4 2463.4
## + rcs(BMI, 4) 3 2430.3 2464.3
## - BPMEDS 1 2439.9 2465.9
## - DIABETES 1 2440.3 2466.3
## - rcs(HEARTRTE, 4) 3 2445.5 2467.5
## - CURSMOKE 1 2444.2 2470.2
## - DIABP 1 2451.9 2477.9
## - AGE 1 2536.6 2562.6
##
## Step: AIC=2460.77
## STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI + PREVHYP +
## AGE + DIABP + rcs(HEARTRTE, 4) + SEX:PREVMI
##
## Df Deviance AIC
## <none> 2434.8 2460.8
## - SEX:PREVMI 1 2437.0 2461.0
## + SEX:PREVHYP 1 2433.4 2461.4
## + TOTCHOL 1 2433.4 2461.4
## + rcs(GLUCOSE, 4) 3 2429.4 2461.4
## + CIGPDAY 1 2434.5 2462.5
## + SYSBP 1 2434.8 2462.8
## + PREVAP 1 2434.8 2462.8
## + rcs(BMI, 4) 3 2431.4 2463.4
## - BPMEDS 1 2441.8 2465.8
## - DIABETES 1 2441.9 2465.9
## - rcs(HEARTRTE, 4) 3 2446.9 2466.9
## - CURSMOKE 1 2445.4 2469.4
## - PREVHYP 1 2451.8 2475.8
## - DIABP 1 2453.6 2477.6
## - AGE 1 2542.1 2566.1##
## Call:
## glm(formula = STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI +
## PREVHYP + AGE + DIABP + rcs(HEARTRTE, 4) + SEX:PREVMI, family = binomial(link = "logit"),
## data = fr.missf_p1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3035 -0.4660 -0.3275 -0.2321 2.9086
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.297986 1.317164 -7.059 1.68e-12 ***
## SEX2 -0.118785 0.115042 -1.033 0.30182
## CURSMOKE1 0.376133 0.115624 3.253 0.00114 **
## DIABETES1 0.672881 0.240467 2.798 0.00514 **
## BPMEDS1 0.569736 0.209547 2.719 0.00655 **
## PREVMI1 -0.052790 0.369570 -0.143 0.88641
## PREVHYP1 0.586447 0.141663 4.140 3.48e-05 ***
## AGE 0.071503 0.007106 10.062 < 2e-16 ***
## DIABP 0.022128 0.005091 4.347 1.38e-05 ***
## rcs(HEARTRTE, 4)HEARTRTE 0.018041 0.018255 0.988 0.32303
## rcs(HEARTRTE, 4)HEARTRTE' -0.184369 0.084305 -2.187 0.02875 *
## rcs(HEARTRTE, 4)HEARTRTE'' 0.488790 0.213925 2.285 0.02232 *
## SEX2:PREVMI1 1.017895 0.670499 1.518 0.12899
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2756.0 on 4433 degrees of freedom
## Residual deviance: 2434.8 on 4421 degrees of freedom
## AIC: 2460.8
##
## Number of Fisher Scoring iterations: 6- glm(formula = STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI + PREVHYP + AGE + DIABP + rcs(HEARTRTE, 4) + SEX:PREVMI, family = binomial(link = “logit”), data = fr.missf_p1) がモデルとして選ばれた
stepwise_results<-
glm(formula = STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI +
PREVHYP + AGE + DIABP + rcs(HEARTRTE, 4) + SEX:PREVMI, family = binomial(link = "logit"),
data = fr.missf_p1) %>%
summary() %>%
.$coefficients %>%
as.data.frame() %>%
rownames_to_column() %>%
rename("var"="rowname") %>%
mutate(OR=exp(Estimate),
Lower=exp(Estimate-1.96*`Std. Error`),
Upper=exp(Estimate+1.96*`Std. Error`)) %>%
select(var,OR,Lower,Upper,"p"="Pr(>|z|)") %>%
#すべての変数を丸める
mutate_if(is.numeric,~round(.,3))
#結果を表示
stepwise_results## var OR Lower Upper p
## 1 (Intercept) 0.000 0.000 0.001 0.000
## 2 SEX2 0.888 0.709 1.113 0.302
## 3 CURSMOKE1 1.457 1.161 1.827 0.001
## 4 DIABETES1 1.960 1.223 3.140 0.005
## 5 BPMEDS1 1.768 1.172 2.666 0.007
## 6 PREVMI1 0.949 0.460 1.957 0.886
## 7 PREVHYP1 1.798 1.362 2.373 0.000
## 8 AGE 1.074 1.059 1.089 0.000
## 9 DIABP 1.022 1.012 1.033 0.000
## 10 rcs(HEARTRTE, 4)HEARTRTE 1.018 0.982 1.055 0.323
## 11 rcs(HEARTRTE, 4)HEARTRTE' 0.832 0.705 0.981 0.029
## 12 rcs(HEARTRTE, 4)HEARTRTE'' 1.630 1.072 2.480 0.022
## 13 SEX2:PREVMI1 2.767 0.744 10.299 0.129#plot
stepwise_results %>%
ggplot(aes(y=var,x=OR))+
geom_point()+
geom_errorbar(aes(xmin=Lower,xmax=Upper),width=0.2)+
geom_vline(xintercept=1,color="red")
- 交互作用項のORが増大しており発散気味だが、このまま計算を進める
#複雑なロジスティク回帰モデル
step_final_model<-
glm(formula = STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI +
PREVHYP + AGE + DIABP + rcs(HEARTRTE, 4) + SEX:PREVMI, family = binomial(link = "logit"),
data = fr.missf_p1)予測性能の評価
ROC AUC
# risk prediction fitting
fr.missf_p1$fitted.step.fin<-
predict(step_final_model,
type="response")
# ROC AUC
roc.step.final<-
roc(STROKE~fitted.step.fin,
data=fr.missf_p1,
ci=T)## Setting levels: control = 0, case = 1## Setting direction: controls < casesroc.step.final##
## Call:
## roc.formula(formula = STROKE ~ fitted.step.fin, data = fr.missf_p1, ci = T)
##
## Data: fitted.step.fin in 4019 controls (STROKE 0) < 415 cases (STROKE 1).
## Area under the curve: 0.7533
## 95% CI: 0.7292-0.7774 (DeLong)roc.step.final %>%
ggroc(legacy.axes=T)+
geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="grey", linetype="dashed")
- 複雑なロジスティク回帰モデル(step_final_model)のROC＿AUCは0.7532998と中等度の予測性能である
Calibration plot
#calibration plot
cal.step.final<-
val.prob(
p=fr.missf_p1$fitted.step.fin,
y=fr.missf_p1$STROKEn,
g=10,
cex=0.5) 
- calibration能は良好だが、liner modelやlassoのほうが当てはまりは視覚的に良さそうに見える
Bootstrapによるliner final modelの内的妥当性検証
optimismの算出
B<-100 #number of resampling
N<-dim(fr.missf_p1)[1] #dataのdimentionを返す→4434行x39列→そのうち一番目（4434)を取得しNとする
AUC.i1<-AUC.i2<-numeric(B)
for(i in 1:B){
#bootstrap random sampling
bs.i<-sample(1:N,N,replace=TRUE) #重複を許して1~4434までの数を4434回サンプリングする
#bootstrap samplingに基づき、dataから集団を再構成 →data.iとする
fr.missf_p1.i<-fr.missf_p1[bs.i,]
#bootstrapしたデータセットごとで回帰モデル構築→各bootstrap dataからリスク予測確率の算出→data2.1のprob_1列に格納
step.fin.cv.res.i<-
glm(formula = STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI +
PREVHYP + AGE + DIABP + rcs(HEARTRTE, 4) + SEX:PREVMI, family = binomial(link = "logit"),
data = fr.missf_p1.i)
fr.missf_p1.i$step.fin.1<- #step.fin.1では、bootstrap dataから構築した新しいロジスティク回帰モデルを使い、そのbootstrap dataに対するリスク予測確率を格納
predict(step.fin.cv.res.i,
data=fr.missf_p1.i,
type="response")
#AUCの算出；Bootstrap回数（B=100)分のAUCを算出して結果を格納；各bootstrap dataに対するAUC
AUC.i1[i]<-roc(STROKE~step.fin.1,data=fr.missf_p1.i)$auc
#bootstrapで作成した各回帰モデルを使って、オリジナルデータ(fr.missf_p1)に対するリスク予測確率を算出→fr.missf_p1のlin.fin.2列に格納;AUC算出
fr.missf_p1$step.fin.2<-
predict(step.fin.cv.res.i,
newdata=fr.missf_p1,
type="response")
AUC.i2[i]<-roc(STROKE~step.fin.2,data=fr.missf_p1)$auc #ここで得られるAUCはbootstrap dataを用いてoriginal dataに外挿したときのAUC
#print(paste(i,"th bootstrap iteration is completed.",sep=""))
}
opt3<-AUC.i1-AUC.i2
#(bootstrap samplingした各data setを使い、それぞれ導いた回帰モデルにおけるAUC)ー(bootstrap samplingした各data setを使い、original dataに外挿したときのAUC)
summary(opt3);hist(opt3)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.026013 -0.001946 0.006097 0.005969 0.015236 0.036204
lam3<-mean(opt3) #estimate of the optimism
step.cor.AUC<-roc.step.final$auc-lam3 #bias corrected AUC estimate
step.cor.AUC## [1] 0.7473308- liner modelと同じく、複雑なロジスティク回帰モデルでもoptimismが大きく、非常にoverfittingが生じていることが示唆される
Cross validationによる内的妥当性の検証
#k-fold CV
k <-5
AUC_CV <-data.frame(matrix(ncol=2,nrow=k))
#Randomly shuffle the data
fr.missf_p1.r<-fr.missf_p1[sample(nrow(fr.missf_p1)),]
#Create k equally size folds
folds <- cut(seq(1,nrow(fr.missf_p1.r)),breaks=k,labels=FALSE)
#Perform k fold cross validation
for(i in 1:k){
#Segment data2r by fold using the which() function
testIndexes <- which(folds==i,arr.ind=TRUE)
testData <- fr.missf_p1.r[testIndexes, ]
trainData <- fr.missf_p1.r[-testIndexes, ]
fit1<-
glm(formula = STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI +
PREVHYP + AGE + DIABP + rcs(HEARTRTE, 4) + SEX:PREVMI, family = binomial(link = "logit"),
data = trainData)
testData$fitted<-predict(fit1,testData,type="response")
ROC <- roc(testData$STROKE, testData$fitted)
AUC_CV[i,2] <- ROC$auc
AUC_CV[i,1] <- i
#print(paste(i,"th cross-validation iteration is completed.", sep=""))
} ## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Setting levels: control = 0, case = 1## Setting direction: controls < casesnames(AUC_CV)[2] <- "AUC_cv"
names(AUC_CV)[1] <- "k"
print(summary(AUC_CV$AUC_cv)) ## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.7133 0.7226 0.7453 0.7458 0.7583 0.7894- CVではAUCはオリジナルデータにおけるAUC: 0.7532998よりも多少小さく（0.7457838）なった
Bootstrapを別の関数を用いて再度検証
使うモデルはSimple model(lin_final_formula)
STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVHYP + AGE + DIABP + HEARTRTE
lin_final_formula<-as.formula(STROKE ~ CURSMOKE + DIABETES + BPMEDS + PREVHYP + AGE + DIABP + HEARTRTE)
lin_final_model.lrm<-
lrm(lin_final_formula,
x=T,
y=T,
data=fr.missf_p1)
boot<-validate(lin_final_model.lrm,
bw=FALSE,
B=100,
method="boot")
boot## index.orig training test optimism index.corrected n
## Dxy 0.5007 0.5085 0.4966 0.0119 0.4888 100
## R2 0.1469 0.1532 0.1438 0.0094 0.1375 100
## Intercept 0.0000 0.0000 -0.0611 0.0611 -0.0611 100
## Slope 1.0000 1.0000 0.9687 0.0313 0.9687 100
## Emax 0.0000 0.0000 0.0187 0.0187 0.0187 100
## D 0.0702 0.0734 0.0687 0.0047 0.0654 100
## U -0.0005 -0.0005 0.0001 -0.0005 0.0001 100
## Q 0.0706 0.0738 0.0686 0.0053 0.0653 100
## B 0.0780 0.0777 0.0783 -0.0006 0.0786 100
## g 1.0429 1.0670 1.0295 0.0375 1.0055 100
## gp 0.0837 0.0852 0.0828 0.0025 0.0812 100corrected.AUC<-boot[1,5]*0.5 + 0.5 #corrected AUC=Somer's D/2 + 0.5- rmsによるBootstrapでは、corrected AUC=0.7443924であった。
- こちらのコードのほうが非常にシンプルで、様々な結果が出てくる
- bias corrected AUCはSomer’s D/2 + 0.5で算出する必要があるのが厄介ではあるが。
Complex modelでも同様に実施
formula = STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI + PREVHYP + AGE + DIABP + rcs(HEARTRTE, 4) + SEX:PREVMI
step_final_model.lrm<-
lrm(STROKE ~ SEX + CURSMOKE + DIABETES + BPMEDS + PREVMI +
PREVHYP + AGE + DIABP + rcs(HEARTRTE, 4) + SEX:PREVMI,
x=T,
y=T,
data=fr.missf_p1
)
boot2<-validate(step_final_model.lrm,
bw=FALSE,
B=100,
method="boot")##
## Divergence or singularity in 4 samplesboot2## index.orig training test optimism index.corrected n
## Dxy 0.5066 0.5141 0.5000 0.0141 0.4925 96
## R2 0.1510 0.1577 0.1454 0.0123 0.1387 96
## Intercept 0.0000 0.0000 -0.0874 0.0874 -0.0874 96
## Slope 1.0000 1.0000 0.9540 0.0460 0.9540 96
## Emax 0.0000 0.0000 0.0272 0.0272 0.0272 96
## D 0.0722 0.0756 0.0695 0.0061 0.0661 96
## U -0.0005 -0.0005 0.0002 -0.0006 0.0002 96
## Q 0.0727 0.0761 0.0693 0.0068 0.0659 96
## B 0.0778 0.0773 0.0782 -0.0010 0.0788 96
## g 1.0623 1.0869 1.0344 0.0525 1.0098 96
## gp 0.0846 0.0861 0.0829 0.0032 0.0814 96corrected.AUC2<-boot2[1,5]*0.5 + 0.5 #AUC=Somer's D/2 + 0.5
corrected.AUC2## [1] 0.7462438