###Assignement 12
setwd("~/Desktop/ACTU 5846/R files")
library(glmnet)
## Warning: package 'glmnet' was built under R version 3.4.2
## Loading required package: Matrix
## Loading required package: foreach
## Warning: package 'foreach' was built under R version 3.4.3
## Loaded glmnet 2.0-13
y = scan('Swiss.txt')
x = read.table('Trendfullregression.txt', header = FALSE)
x = as.matrix(x)
N = length(y)
fit1 = glmnet(x, y, standardize = FALSE)
answer = print(fit1)
##
## Call: glmnet(x = x, y = y, standardize = FALSE)
##
## Df %Dev Lambda
## [1,] 0 0.0000 5.6980000
## [2,] 1 0.1477 5.1920000
## [3,] 1 0.2703 4.7310000
## [4,] 1 0.3721 4.3110000
## [5,] 1 0.4566 3.9280000
## [6,] 1 0.5267 3.5790000
## [7,] 1 0.5850 3.2610000
## [8,] 1 0.6333 2.9710000
## [9,] 1 0.6735 2.7070000
## [10,] 1 0.7068 2.4670000
## [11,] 1 0.7345 2.2480000
## [12,] 1 0.7575 2.0480000
## [13,] 1 0.7765 1.8660000
## [14,] 1 0.7924 1.7000000
## [15,] 1 0.8055 1.5490000
## [16,] 1 0.8164 1.4120000
## [17,] 1 0.8255 1.2860000
## [18,] 1 0.8330 1.1720000
## [19,] 1 0.8393 1.0680000
## [20,] 1 0.8444 0.9729000
## [21,] 1 0.8487 0.8865000
## [22,] 1 0.8523 0.8077000
## [23,] 1 0.8553 0.7360000
## [24,] 1 0.8577 0.6706000
## [25,] 1 0.8598 0.6110000
## [26,] 1 0.8615 0.5567000
## [27,] 1 0.8629 0.5073000
## [28,] 1 0.8641 0.4622000
## [29,] 1 0.8650 0.4212000
## [30,] 1 0.8658 0.3837000
## [31,] 1 0.8665 0.3496000
## [32,] 1 0.8671 0.3186000
## [33,] 1 0.8675 0.2903000
## [34,] 1 0.8679 0.2645000
## [35,] 2 0.8767 0.2410000
## [36,] 2 0.8850 0.2196000
## [37,] 2 0.8921 0.2001000
## [38,] 3 0.8981 0.1823000
## [39,] 3 0.9039 0.1661000
## [40,] 3 0.9095 0.1514000
## [41,] 3 0.9141 0.1379000
## [42,] 4 0.9179 0.1257000
## [43,] 4 0.9213 0.1145000
## [44,] 4 0.9240 0.1043000
## [45,] 4 0.9263 0.0950500
## [46,] 5 0.9284 0.0866100
## [47,] 5 0.9301 0.0789200
## [48,] 4 0.9314 0.0719100
## [49,] 4 0.9325 0.0655200
## [50,] 4 0.9334 0.0597000
## [51,] 4 0.9341 0.0543900
## [52,] 5 0.9349 0.0495600
## [53,] 5 0.9360 0.0451600
## [54,] 5 0.9369 0.0411500
## [55,] 6 0.9376 0.0374900
## [56,] 6 0.9385 0.0341600
## [57,] 6 0.9393 0.0311300
## [58,] 6 0.9399 0.0283600
## [59,] 6 0.9405 0.0258400
## [60,] 8 0.9410 0.0235500
## [61,] 8 0.9420 0.0214500
## [62,] 9 0.9433 0.0195500
## [63,] 9 0.9446 0.0178100
## [64,] 8 0.9456 0.0162300
## [65,] 8 0.9465 0.0147900
## [66,] 8 0.9472 0.0134700
## [67,] 9 0.9479 0.0122800
## [68,] 9 0.9486 0.0111900
## [69,] 9 0.9492 0.0101900
## [70,] 9 0.9497 0.0092870
## [71,] 9 0.9501 0.0084620
## [72,] 9 0.9505 0.0077100
## [73,] 9 0.9507 0.0070250
## [74,] 9 0.9510 0.0064010
## [75,] 11 0.9513 0.0058320
## [76,] 12 0.9519 0.0053140
## [77,] 12 0.9523 0.0048420
## [78,] 12 0.9526 0.0044120
## [79,] 13 0.9530 0.0040200
## [80,] 15 0.9532 0.0036630
## [81,] 15 0.9537 0.0033380
## [82,] 16 0.9541 0.0030410
## [83,] 16 0.9544 0.0027710
## [84,] 17 0.9548 0.0025250
## [85,] 17 0.9551 0.0023000
## [86,] 17 0.9553 0.0020960
## [87,] 17 0.9555 0.0019100
## [88,] 17 0.9556 0.0017400
## [89,] 17 0.9558 0.0015860
## [90,] 17 0.9559 0.0014450
## [91,] 17 0.9560 0.0013160
## [92,] 17 0.9560 0.0011990
## [93,] 17 0.9561 0.0010930
## [94,] 17 0.9562 0.0009958
## [95,] 17 0.9562 0.0009073
## [96,] 17 0.9562 0.0008267
## [97,] 17 0.9563 0.0007533
## [98,] 18 0.9563 0.0006864
## [99,] 18 0.9563 0.0006254
## [100,] 18 0.9564 0.0005698
answer
## Df %Dev Lambda
## [1,] 0 0.0000 5.6980000
## [2,] 1 0.1477 5.1920000
## [3,] 1 0.2703 4.7310000
## [4,] 1 0.3721 4.3110000
## [5,] 1 0.4566 3.9280000
## [6,] 1 0.5267 3.5790000
## [7,] 1 0.5850 3.2610000
## [8,] 1 0.6333 2.9710000
## [9,] 1 0.6735 2.7070000
## [10,] 1 0.7068 2.4670000
## [11,] 1 0.7345 2.2480000
## [12,] 1 0.7575 2.0480000
## [13,] 1 0.7765 1.8660000
## [14,] 1 0.7924 1.7000000
## [15,] 1 0.8055 1.5490000
## [16,] 1 0.8164 1.4120000
## [17,] 1 0.8255 1.2860000
## [18,] 1 0.8330 1.1720000
## [19,] 1 0.8393 1.0680000
## [20,] 1 0.8444 0.9729000
## [21,] 1 0.8487 0.8865000
## [22,] 1 0.8523 0.8077000
## [23,] 1 0.8553 0.7360000
## [24,] 1 0.8577 0.6706000
## [25,] 1 0.8598 0.6110000
## [26,] 1 0.8615 0.5567000
## [27,] 1 0.8629 0.5073000
## [28,] 1 0.8641 0.4622000
## [29,] 1 0.8650 0.4212000
## [30,] 1 0.8658 0.3837000
## [31,] 1 0.8665 0.3496000
## [32,] 1 0.8671 0.3186000
## [33,] 1 0.8675 0.2903000
## [34,] 1 0.8679 0.2645000
## [35,] 2 0.8767 0.2410000
## [36,] 2 0.8850 0.2196000
## [37,] 2 0.8921 0.2001000
## [38,] 3 0.8981 0.1823000
## [39,] 3 0.9039 0.1661000
## [40,] 3 0.9095 0.1514000
## [41,] 3 0.9141 0.1379000
## [42,] 4 0.9179 0.1257000
## [43,] 4 0.9213 0.1145000
## [44,] 4 0.9240 0.1043000
## [45,] 4 0.9263 0.0950500
## [46,] 5 0.9284 0.0866100
## [47,] 5 0.9301 0.0789200
## [48,] 4 0.9314 0.0719100
## [49,] 4 0.9325 0.0655200
## [50,] 4 0.9334 0.0597000
## [51,] 4 0.9341 0.0543900
## [52,] 5 0.9349 0.0495600
## [53,] 5 0.9360 0.0451600
## [54,] 5 0.9369 0.0411500
## [55,] 6 0.9376 0.0374900
## [56,] 6 0.9385 0.0341600
## [57,] 6 0.9393 0.0311300
## [58,] 6 0.9399 0.0283600
## [59,] 6 0.9405 0.0258400
## [60,] 8 0.9410 0.0235500
## [61,] 8 0.9420 0.0214500
## [62,] 9 0.9433 0.0195500
## [63,] 9 0.9446 0.0178100
## [64,] 8 0.9456 0.0162300
## [65,] 8 0.9465 0.0147900
## [66,] 8 0.9472 0.0134700
## [67,] 9 0.9479 0.0122800
## [68,] 9 0.9486 0.0111900
## [69,] 9 0.9492 0.0101900
## [70,] 9 0.9497 0.0092870
## [71,] 9 0.9501 0.0084620
## [72,] 9 0.9505 0.0077100
## [73,] 9 0.9507 0.0070250
## [74,] 9 0.9510 0.0064010
## [75,] 11 0.9513 0.0058320
## [76,] 12 0.9519 0.0053140
## [77,] 12 0.9523 0.0048420
## [78,] 12 0.9526 0.0044120
## [79,] 13 0.9530 0.0040200
## [80,] 15 0.9532 0.0036630
## [81,] 15 0.9537 0.0033380
## [82,] 16 0.9541 0.0030410
## [83,] 16 0.9544 0.0027710
## [84,] 17 0.9548 0.0025250
## [85,] 17 0.9551 0.0023000
## [86,] 17 0.9553 0.0020960
## [87,] 17 0.9555 0.0019100
## [88,] 17 0.9556 0.0017400
## [89,] 17 0.9558 0.0015860
## [90,] 17 0.9559 0.0014450
## [91,] 17 0.9560 0.0013160
## [92,] 17 0.9560 0.0011990
## [93,] 17 0.9561 0.0010930
## [94,] 17 0.9562 0.0009958
## [95,] 17 0.9562 0.0009073
## [96,] 17 0.9562 0.0008267
## [97,] 17 0.9563 0.0007533
## [98,] 18 0.9563 0.0006864
## [99,] 18 0.9563 0.0006254
## [100,] 18 0.9564 0.0005698
lambda = answer[,3]
rsq = answer[,2]
dof = answer[,1]
sst = sum((y-mean(y))^2)
ssr = (1-rsq)*sst
sigsq = ssr/(N-dof)
NLL = N*log(sigsq)+ssr/2/sigsq
beta = fit1[[2]]
h=dim(answer)[1]
k=dim(x)[2]
L1 = c(1:h)
for(i in 1:h) L1[i] = sum(abs(beta[2:k,i])) #sum of absolute value of non-intercept coefs
min = NLL-dof*log(lambda)+lambda*L1
plot(fit1, label=TRUE)
min
## [1] 139.1546578 128.5099482 119.2727479 110.2654986 101.5464556
## [6] 93.1729916 85.1716841 77.6168526 70.5099937 63.9140876
## [11] 57.8211853 52.2532368 47.2379003 42.7012551 38.6963721
## [16] 35.1559199 32.0395586 29.3520893 27.0057140 25.0460777
## [21] 23.3507511 21.9024199 20.6774497 19.6905216 18.8214565
## [26] 18.1204360 17.5483675 17.0637192 16.7143903 16.4110022
## [31] 16.1538759 15.9435496 15.8276049 15.7128830 13.4839918
## [36] 9.3518397 5.5858747 4.4655629 1.1074878 -2.3413246
## [41] -5.3021017 -5.2037415 -7.4618867 -9.2628828 -10.8067730
## [46] -9.2122464 -10.2469277 -14.1664652 -14.8062671 -15.2758039
## [51] -15.5672680 -12.3776308 -12.9763013 -13.3960749 -9.7473662
## [56] -10.0955550 -10.3551425 -10.4168976 -10.4854673 -1.7043340
## [61] -2.0199045 1.9086412 1.3079381 -3.7673173 -4.0607240
## [66] -4.1319612 0.8366980 0.8321761 0.9443275 1.1637626
## [71] 1.5035898 1.8398307 2.4239278 2.8805064 15.0096378
## [76] 21.2284912 21.8258324 22.5492985 29.4083223 43.1541567
## [81] 43.8813588 51.3694641 52.4497302 60.2344237 61.4072652
## [86] 62.7084885 64.0092954 65.4535972 66.7481645 68.1896481
## [91] 69.6377607 71.2197972 72.6515078 74.0926753 75.6742375
## [96] 77.2551261 78.6934431 88.4513041 90.1260388 91.6594330
lambda
## [1] 5.6980000 5.1920000 4.7310000 4.3110000 3.9280000 3.5790000 3.2610000
## [8] 2.9710000 2.7070000 2.4670000 2.2480000 2.0480000 1.8660000 1.7000000
## [15] 1.5490000 1.4120000 1.2860000 1.1720000 1.0680000 0.9729000 0.8865000
## [22] 0.8077000 0.7360000 0.6706000 0.6110000 0.5567000 0.5073000 0.4622000
## [29] 0.4212000 0.3837000 0.3496000 0.3186000 0.2903000 0.2645000 0.2410000
## [36] 0.2196000 0.2001000 0.1823000 0.1661000 0.1514000 0.1379000 0.1257000
## [43] 0.1145000 0.1043000 0.0950500 0.0866100 0.0789200 0.0719100 0.0655200
## [50] 0.0597000 0.0543900 0.0495600 0.0451600 0.0411500 0.0374900 0.0341600
## [57] 0.0311300 0.0283600 0.0258400 0.0235500 0.0214500 0.0195500 0.0178100
## [64] 0.0162300 0.0147900 0.0134700 0.0122800 0.0111900 0.0101900 0.0092870
## [71] 0.0084620 0.0077100 0.0070250 0.0064010 0.0058320 0.0053140 0.0048420
## [78] 0.0044120 0.0040200 0.0036630 0.0033380 0.0030410 0.0027710 0.0025250
## [85] 0.0023000 0.0020960 0.0019100 0.0017400 0.0015860 0.0014450 0.0013160
## [92] 0.0011990 0.0010930 0.0009958 0.0009073 0.0008267 0.0007533 0.0006864
## [99] 0.0006254 0.0005698
# the lowest min is #51 with later lamda value equals to 0.0543900.
# 63rd is before increse starts,with lambda 0.0178100.
coef(fit1, s=0.017)
## 19 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 4.999046230
## V1 0.094253475
## V2 .
## V3 .
## V4 .
## V5 -0.155108366
## V6 .
## V7 .
## V8 .
## V9 -1.373967820
## V10 -0.256331867
## V11 0.415489107
## V12 0.585430198
## V13 .
## V14 0.124508436
## V15 0.001808786
## V16 .
## V17 0.200832902
## V18 .
cvfit = cv.glmnet(x, y, standardize = FALSE)
cvfit$lambda.1se
## [1] 0.1379082
cvfit$lambda.min
## [1] 0.01019241
#Coeffcients for 0.0543900 and 0.0178100 & lambda.1se and lambda.min.
coef(cvfit, s = "lambda.1se")
## 19 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 4.6914905
## V1 .
## V2 .
## V3 .
## V4 .
## V5 .
## V6 .
## V7 .
## V8 .
## V9 -1.0897684
## V10 .
## V11 .
## V12 0.2382458
## V13 0.1801836
## V14 .
## V15 .
## V16 .
## V17 .
## V18 .
coef(cvfit, s = "lambda.min")
## 19 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 4.95058031
## V1 0.16819969
## V2 -0.07381200
## V3 .
## V4 .
## V5 -0.16635734
## V6 .
## V7 .
## V8 .
## V9 -1.29119121
## V10 -0.52531170
## V11 0.70022564
## V12 0.49026123
## V13 .
## V14 0.09235081
## V15 .
## V16 .
## V17 0.54212511
## V18 .
coef(fit1, s = c(0.0543900, 0.0178100))
## 19 x 2 sparse Matrix of class "dgCMatrix"
## 1 2
## (Intercept) 4.960018e+00 5.001415563
## V1 1.940272e-02 0.092453880
## V2 . .
## V3 . .
## V4 . .
## V5 -5.443126e-06 -0.151370158
## V6 . .
## V7 . .
## V8 . .
## V9 -1.321330e+00 -1.383386398
## V10 . -0.226120101
## V11 . 0.384747373
## V12 6.936972e-01 0.594702547
## V13 . .
## V14 6.540153e-02 0.127350933
## V15 . 0.003709564
## V16 . .
## V17 . 0.159812155
## V18 . .
#Based on the coeffcients above, we choose the column of 0.01781 for our redueced model.
#Although the lambda.min and the column of 0.05439 have same variables but the previous one shrinks more, constant higher.
library(rstan)
## Warning: package 'rstan' was built under R version 3.4.3
## Loading required package: ggplot2
## Loading required package: StanHeaders
## Warning: package 'StanHeaders' was built under R version 3.4.3
## rstan (Version 2.17.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
library("loo")
## This is loo version 1.1.0
set.seed(9)
x1 = read.table('reduced.txt', header = FALSE)
x1 = as.matrix(x1)
U = ncol(x1)
N = length(y)
fit2 = stan(file = 'logregression.stan', verbose = FALSE, chains = 4, iter = 2000)
## Warning in strptime(xx, f <- "%Y-%m-%d %H:%M:%OS", tz = tz): unknown
## timezone 'zone/tz/2017c.1.0/zoneinfo/America/New_York'
## recompiling to avoid crashing R session
## In file included from file369b6ed61f.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/var.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/math/tools/config.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/config.hpp:39:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/config/compiler/clang.hpp:200:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined]
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command line>:6:9: note: previous definition is here
## #define BOOST_NO_CXX11_RVALUE_REFERENCES 1
## ^
## In file included from file369b6ed61f.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:44:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: unused function 'set_zero_all_adjoints' [-Wunused-function]
## static void set_zero_all_adjoints() {
## ^
## In file included from file369b6ed61f.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:45:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints_nested.hpp:17:17: warning: 'static' function 'set_zero_all_adjoints_nested' declared in header file should be declared 'static inline' [-Wunneeded-internal-declaration]
## static void set_zero_all_adjoints_nested() {
## ^
## In file included from file369b6ed61f.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:58:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/autocorrelation.hpp:17:14: warning: function 'fft_next_good_size' is not needed and will not be emitted [-Wunneeded-internal-declaration]
## size_t fft_next_good_size(size_t N) {
## ^
## In file included from file369b6ed61f.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:298:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/arr.hpp:38:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:17:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/numeric/odeint.hpp:61:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/numeric/odeint/util/multi_array_adaption.hpp:29:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/multi_array.hpp:21:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/multi_array/base.hpp:28:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:42:43: warning: unused typedef 'index_range' [-Wunused-local-typedef]
## typedef typename Array::index_range index_range;
## ^
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:43:37: warning: unused typedef 'index' [-Wunused-local-typedef]
## typedef typename Array::index index;
## ^
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:53:43: warning: unused typedef 'index_range' [-Wunused-local-typedef]
## typedef typename Array::index_range index_range;
## ^
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:54:37: warning: unused typedef 'index' [-Wunused-local-typedef]
## typedef typename Array::index index;
## ^
## 8 warnings generated.
print(fit2, pars=c("cn", "v", "sig", "s"), probs=c(.05, 0.2, 0.5, 0.8, 0.95))
## Inference for Stan model: logregression.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 5% 20% 50% 80% 95% n_eff Rhat
## cn 4.95 0.00 0.25 4.55 4.75 4.95 5.16 5.36 2684 1
## v[1] 0.11 0.00 0.06 0.01 0.06 0.11 0.16 0.21 2208 1
## v[2] -0.19 0.00 0.12 -0.40 -0.29 -0.19 -0.08 0.01 2397 1
## v[3] -1.24 0.01 0.26 -1.64 -1.45 -1.25 -1.02 -0.80 1943 1
## v[4] -0.61 0.01 0.44 -1.36 -0.98 -0.58 -0.21 0.04 1698 1
## v[5] 0.75 0.01 0.43 0.07 0.37 0.74 1.11 1.47 1981 1
## v[6] 0.45 0.01 0.32 -0.05 0.17 0.43 0.72 0.99 2619 1
## v[7] 0.12 0.01 0.27 -0.31 -0.09 0.10 0.33 0.55 2696 1
## v[8] 0.02 0.01 0.28 -0.44 -0.20 0.03 0.24 0.48 2163 1
## v[9] 0.56 0.01 0.51 -0.14 0.11 0.50 0.99 1.51 2289 1
## sig 0.59 0.00 0.06 0.50 0.54 0.59 0.64 0.69 2867 1
## s 0.55 0.01 0.23 0.27 0.37 0.51 0.71 1.00 1461 1
##
## Samples were drawn using NUTS(diag_e) at Thu Dec 21 17:23:13 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
plot(fit2, pars = "v")
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)
out <- get_posterior_mean(fit2)
write.csv(out, file="out_fit2.csv")
log_LL <- extract_log_lik(fit2)
loo_LL <- loo(log_LL)
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-
## diagnostic') for details.
loo_LL
## Computed from 4000 by 62 log-likelihood matrix
##
## Estimate SE
## elpd_loo -62.3 13.9
## p_loo 12.2 6.4
## looic 124.6 27.8
##
## Pareto k diagnostic values:
## Count Pct
## (-Inf, 0.5] (good) 61 98.4%
## (0.5, 0.7] (ok) 0 0.0%
## (0.7, 1] (bad) 0 0.0%
## (1, Inf) (very bad) 1 1.6%
## See help('pareto-k-diagnostic') for details.
loo(log_LL[1:1000,])
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-
## diagnostic') for details.
## Computed from 1000 by 62 log-likelihood matrix
##
## Estimate SE
## elpd_loo -61.5 13.3
## p_loo 11.3 5.5
## looic 123.1 26.5
##
## Pareto k diagnostic values:
## Count Pct
## (-Inf, 0.5] (good) 58 93.5%
## (0.5, 0.7] (ok) 3 4.8%
## (0.7, 1] (bad) 0 0.0%
## (1, Inf) (very bad) 1 1.6%
## See help('pareto-k-diagnostic') for details.
loo(log_LL[1001:2000,])
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-
## diagnostic') for details.
## Computed from 1000 by 62 log-likelihood matrix
##
## Estimate SE
## elpd_loo -61.8 13.4
## p_loo 11.7 5.9
## looic 123.5 26.8
##
## Pareto k diagnostic values:
## Count Pct
## (-Inf, 0.5] (good) 60 96.8%
## (0.5, 0.7] (ok) 1 1.6%
## (0.7, 1] (bad) 0 0.0%
## (1, Inf) (very bad) 1 1.6%
## See help('pareto-k-diagnostic') for details.
loo(log_LL[2001:3000,])
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-
## diagnostic') for details.
## Computed from 1000 by 62 log-likelihood matrix
##
## Estimate SE
## elpd_loo -62.3 14.0
## p_loo 12.3 6.6
## looic 124.6 28.1
##
## Pareto k diagnostic values:
## Count Pct
## (-Inf, 0.5] (good) 61 98.4%
## (0.5, 0.7] (ok) 0 0.0%
## (0.7, 1] (bad) 0 0.0%
## (1, Inf) (very bad) 1 1.6%
## See help('pareto-k-diagnostic') for details.
loo(log_LL[3001:4000,])
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-
## diagnostic') for details.
## Computed from 1000 by 62 log-likelihood matrix
##
## Estimate SE
## elpd_loo -61.9 13.3
## p_loo 11.7 6.0
## looic 123.8 26.7
##
## Pareto k diagnostic values:
## Count Pct
## (-Inf, 0.5] (good) 61 98.4%
## (0.5, 0.7] (ok) 0 0.0%
## (0.7, 1] (bad) 0 0.0%
## (1, Inf) (very bad) 1 1.6%
## See help('pareto-k-diagnostic') for details.
fit2_ss = extract(fit2, permuted = FALSE)
dim(fit2_ss)
## [1] 1000 4 139
fit2_ss[1,,]
## parameters
## chains cn v[1] v[2] v[3] v[4]
## chain:1 4.728812 0.270724047 -0.50821083 -1.2703896 -0.5043485
## chain:2 4.634634 -0.003660029 0.02493424 -0.6305092 -1.1218700
## chain:3 4.853581 0.132238212 -0.14602855 -1.1488543 -0.6079608
## chain:4 5.458185 -0.005011595 0.08089612 -1.3986261 -0.4279217
## parameters
## chains v[5] v[6] v[7] v[8] v[9]
## chain:1 0.4772071 0.7076945 0.65087836 -0.84628817 0.9137740
## chain:2 0.6647887 0.2835790 0.48596011 -0.25146738 0.4842604
## chain:3 0.5148438 0.5442275 0.40832317 -0.51551025 1.0317544
## chain:4 0.1240131 1.2521386 0.03771278 -0.06554464 0.3996669
## parameters
## chains logs logsig sig s mu[1] mu[2]
## chain:1 0.004379172 -0.4243227 0.6542127 1.0043888 5.270261 5.540985
## chain:2 -0.491501611 -0.5723744 0.5641842 0.6117072 4.627314 4.623654
## chain:3 -0.779278930 -0.4052164 0.6668325 0.4587367 5.118058 5.250296
## chain:4 -0.855652137 -0.5296998 0.5887817 0.4250059 5.448162 5.443150
## parameters
## chains mu[3] mu[4] mu[5] mu[6] mu[7] mu[8] mu[9]
## chain:1 4.861762 5.574222 5.811709 5.336735 5.099248 4.999536 4.728812
## chain:2 4.705091 4.641268 4.619994 4.662543 4.683817 4.630974 4.634634
## chain:3 5.327373 5.368744 5.382534 5.354953 5.341163 4.985819 4.853581
## chain:4 5.741677 5.514023 5.438139 5.589908 5.665792 5.453173 5.458185
## parameters
## chains mu[10] mu[11] mu[12] mu[13] mu[14] mu[15] mu[16]
## chain:1 3.458423 4.066345 3.729147 3.999871 4.541319 4.270595 4.303832
## chain:2 4.004125 4.032033 4.000465 3.996805 3.989485 3.993145 4.010759
## chain:3 3.704727 4.206099 3.836965 3.969203 4.233680 4.101441 4.219889
## chain:4 4.059559 4.191282 4.054547 4.049536 4.039513 4.044524 4.115397
## parameters
## chains mu[17] mu[18] mu[19] mu[20] mu[21] mu[22] mu[23]
## chain:1 3.591372 3.828859 2.225133 1.954409 1.683685 2.495857 2.766581
## chain:2 4.074582 4.053308 2.244426 2.248086 2.251746 2.240766 2.237106
## chain:3 4.178518 4.192309 2.212388 2.080150 1.947912 2.344626 2.476865
## chain:4 4.343051 4.267166 2.222988 2.228000 2.233011 2.217976 2.212965
## parameters
## chains mu[24] mu[25] mu[26] mu[27] mu[28] mu[29]
## chain:1 2.054121 2.529094 2.291607 1.816634 0.4042399 0.3861537
## chain:2 2.300928 2.258380 2.279654 2.322203 0.3880521 1.1641553
## chain:3 2.435494 2.463074 2.449284 2.421703 0.5095687 0.7059404
## chain:4 2.440618 2.288849 2.364734 2.516503 0.2118030 0.5304765
## parameters
## chains mu[30] mu[31] mu[32] mu[33] mu[34] mu[35]
## chain:1 0.87921346 1.1983259 1.4690499 0.7565896 0.9276018 1.2315631
## chain:2 0.34550368 1.1531752 1.1495152 1.2133378 1.1568352 1.1707894
## chain:3 0.53714942 1.1026550 1.2348932 1.1935222 0.9704168 1.2211029
## chain:4 0.06003396 0.5154417 0.5104301 0.7380836 0.5204533 0.5863146
## parameters
## chains mu[36] mu[37] mu[38] mu[39] mu[40] mu[41]
## chain:1 0.9940764 -0.20368273 0.33776536 0.60848941 0.6568778 0.3504189
## chain:2 1.1920636 0.36014380 0.35282374 0.34916371 1.1604952 -0.7765591
## chain:3 1.2073126 0.00819657 0.27267300 0.40491121 0.8381786 -0.4500151
## chain:4 0.6621991 0.08008034 0.07005715 0.06504556 0.5254649 -0.8030455
## parameters
## chains mu[42] mu[43] mu[44] mu[45] mu[46]
## chain:1 0.2893770 0.6417267 -1.247000 0.01865294 0.05189021
## chain:2 -0.4585078 0.3667779 -1.335098 -0.45484776 -0.43723357
## chain:3 -0.1605944 0.5233591 -1.651661 -0.29283262 -0.17438475
## chain:4 -0.3903622 0.1359185 -1.266239 -0.38535056 -0.31447763
## parameters
## chains mu[47] mu[48] mu[49] mu[50] mu[51]
## chain:1 -0.2520711 0.06704132 -0.7055514 -0.7935192 -0.1910292
## chain:2 -0.4511877 0.35648377 -1.3424179 -0.4438677 -0.7692391
## chain:3 -0.4250708 0.14043478 -1.3871842 -0.6895473 -0.7144915
## chain:4 -0.3803390 0.07506875 -1.2762619 -0.3703158 -0.7930223
## parameters
## chains mu[52] mu[53] mu[54] mu[55] mu[56] mu[57]
## chain:1 -0.4617533 -2.032246 -0.5227952 0.07969483 -1.903718 -1.517724
## chain:2 -0.7655791 -1.904617 -0.4475277 -0.77289912 -1.989875 -1.331438
## chain:3 -0.8467297 -2.456592 -0.5573090 -0.58225328 -2.229768 -1.783899
## chain:4 -0.7880107 -1.744467 -0.3753274 -0.79803389 -1.823028 -1.261227
## parameters
## chains mu[58] mu[59] mu[60] mu[61] mu[62] log_lik[1]
## chain:1 -1.761522 -0.7324773 -2.174442 -0.9762755 -2.302970 -0.4953445
## chain:2 -1.908277 -0.7619190 -1.986215 -1.3387578 -1.900957 -0.9464451
## chain:3 -2.324353 -0.9789679 -2.362006 -1.5194224 -2.588830 -0.5319235
## chain:4 -1.749478 -0.7829991 -1.818016 -1.2712503 -1.739455 -0.4486028
## parameters
## chains log_lik[2] log_lik[3] log_lik[4] log_lik[5] log_lik[6]
## chain:1 -0.5986698 -0.6586663 -0.6381510 -0.9066694 -0.5168871
## chain:2 -0.9482222 -0.7901503 -0.8794331 -0.9079558 -0.7980586
## chain:3 -0.5137898 -0.5230086 -0.5373773 -0.5442320 -0.5411886
## chain:4 -0.4472839 -0.7573308 -0.5108070 -0.4592544 -0.6100183
## parameters
## chains log_lik[7] log_lik[8] log_lik[9] log_lik[10] log_lik[11]
## chain:1 -0.5048259 -0.5307301 -0.6245237 -1.1164272 -0.4997694
## chain:2 -0.7534040 -0.8120851 -0.6338362 -0.3996661 -0.3625027
## chain:3 -0.5384947 -0.5541178 -0.5626967 -0.7763241 -0.5197695
## chain:4 -0.7120060 -0.5005586 -0.6153107 -0.4130279 -0.3941778
## parameters
## chains log_lik[12] log_lik[13] log_lik[14] log_lik[15] log_lik[16]
## chain:1 -0.6548515 -0.5051533 -0.7563370 -0.5487004 -0.5785235
## chain:2 -0.3619700 -0.3616625 -0.3562478 -0.3526581 -0.3475516
## chain:3 -0.5912226 -0.5314723 -0.5445876 -0.5161026 -0.5518149
## chain:4 -0.3921530 -0.3921997 -0.3904093 -0.3894103 -0.3983695
## parameters
## chains log_lik[17] log_lik[18] log_lik[19] log_lik[20] log_lik[21]
## chain:1 -0.7121095 -0.5345983 -0.5115823 -0.5364263 -0.6695533
## chain:2 -0.3507682 -0.3490088 -0.3626576 -0.3637167 -0.3980787
## chain:3 -0.5409708 -0.5495299 -0.5336892 -0.5182476 -0.5306624
## chain:4 -0.5371189 -0.4817857 -0.4109379 -0.3995154 -0.4272585
## parameters
## chains log_lik[22] log_lik[23] log_lik[24] log_lik[25] log_lik[26]
## chain:1 -0.7265670 -1.2924528 -0.5760555 -1.2115437 -0.8870873
## chain:2 -0.4035671 -0.4850563 -0.7564780 -0.7594161 -0.8527431
## chain:3 -0.6111500 -0.8375958 -0.9821010 -1.0923632 -1.1249631
## chain:4 -0.4298045 -0.4965644 -0.9996071 -0.8147147 -1.0037658
## parameters
## chains log_lik[27] log_lik[28] log_lik[29] log_lik[30] log_lik[31]
## chain:1 -0.5112108 -1.426693 -1.4357536 -0.6477533 -0.5040289
## chain:2 -0.9596846 -1.645682 -0.3690156 -1.6069867 -0.3496905
## chain:3 -1.1035412 -1.211750 -0.8890774 -1.0712015 -0.5137614
## chain:4 -1.3568180 -2.089240 -1.2075492 -2.4017287 -0.8966358
## parameters
## chains log_lik[32] log_lik[33] log_lik[34] log_lik[35] log_lik[36]
## chain:1 -0.6823733 -0.5684323 -0.4971919 -0.6260937 -0.5077881
## chain:2 -0.3569626 -0.4128230 -0.3987538 -0.4651000 -0.4918982
## chain:3 -0.5449842 -0.5524413 -0.5137415 -0.6325021 -0.6284488
## chain:4 -0.8378793 -0.4942859 -0.6866631 -0.5276431 -0.4627061
## parameters
## chains log_lik[37] log_lik[38] log_lik[39] log_lik[40] log_lik[41]
## chain:1 -1.2573151 -0.5049422 -0.5681077 -0.6092779 -0.5383190
## chain:2 -0.4402168 -0.3563574 -0.3466779 -1.3948308 -1.7156014
## chain:3 -0.9133015 -0.5421883 -0.5162312 -0.7887808 -0.9280468
## chain:4 -0.7856212 -0.5779593 -0.5127470 -0.4369486 -1.7186165
## parameters
## chains log_lik[42] log_lik[43] log_lik[44] log_lik[45] log_lik[46]
## chain:1 -0.5308231 -0.9053383 -2.340345 -0.5003314 -0.5127125
## chain:2 -0.8602195 -0.5054001 -3.188431 -0.6023829 -0.5554510
## chain:3 -0.5980933 -0.7669638 -3.618247 -0.5793234 -0.5253782
## chain:4 -0.7551616 -0.4001879 -2.738284 -0.5501933 -0.4736418
## parameters
## chains log_lik[47] log_lik[48] log_lik[49] log_lik[50] log_lik[51]
## chain:1 -0.5258243 -0.8216321 -0.5467529 -0.5393495 -0.7057850
## chain:2 -0.3485630 -1.3989763 -1.4764753 -0.3838030 -0.3833610
## chain:3 -0.5138248 -0.9218608 -1.4101809 -0.5231796 -0.5245887
## chain:4 -0.3910235 -0.8053203 -1.2711738 -0.4639017 -0.4343416
## parameters
## chains log_lik[52] log_lik[53] log_lik[54] log_lik[55] log_lik[56]
## chain:1 -0.5685664 -1.818821 -0.8557797 -2.7486168 -0.7147035
## chain:2 -0.3508491 -1.725796 -0.9725649 -0.7985876 -0.7716405
## chain:3 -0.5337262 -3.006768 -0.8195283 -1.1081524 -1.1633564
## chain:4 -0.3972786 -1.259746 -1.1030243 -0.7663001 -0.5693230
## parameters
## chains log_lik[57] log_lik[58] log_lik[59] log_lik[60] log_lik[61]
## chain:1 -0.4946310 -0.5064837 -2.272509 -0.5453187 -2.889133
## chain:2 -0.3989914 -0.4428223 -2.624392 -0.3471989 -2.142274
## chain:3 -0.5955550 -1.0089194 -1.609439 -0.6899575 -1.401438
## chain:4 -0.4814877 -0.4005985 -2.408137 -0.4208727 -2.252825
## parameters
## chains log_lik[62] lp__
## chain:1 -6.686434 -5.159025018
## chain:2 -11.833657 -2.318758974
## chain:3 -5.085280 -0.004473056
## chain:4 -12.234030 -2.253576189
fit2_ss = fit2_ss[,,1:13]
dim(fit2_ss) = c(4000,13)
fit2_ss[1:4,]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 4.728812 0.2707240 -0.50821083 -1.270390 -0.5043485 0.4772071
## [2,] 4.868296 0.1542473 -0.25314313 -1.172641 -0.3557149 0.3259970
## [3,] 5.109792 0.0954108 -0.13878352 -1.331257 -0.7408038 1.1960905
## [4,] 5.230506 -0.0299422 0.02545317 -1.112368 -0.9128487 1.1123810
## [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0.7076945 0.6508784 -0.84628817 0.91377400 0.004379172 -0.4243227
## [2,] 0.3644840 0.6156502 -0.30986422 0.28476890 -1.053084617 -0.6084608
## [3,] 0.3800718 -0.4384711 0.58015988 -0.02970696 -0.667796885 -0.4061983
## [4,] 0.1841625 0.2062065 -0.02041861 0.64158126 -0.599702177 -0.7222139
## [,13]
## [1,] 0.6542127
## [2,] 0.5441878
## [3,] 0.6661780
## [4,] 0.4856758
write.csv(fit2_ss, file="samples_fit2.csv")