20221102 Bài tập số 3

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

library(BayesFactor)

## Loading required package: coda

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## ************
## Welcome to BayesFactor 0.9.12-4.4. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
## 
## Type BFManual() to open the manual.
## ************

bf = ttestBF(formula = uwes_tol  ~ gioi, data=thuan)
mcmc = posterior(bf, iter =10000)
summary(mcmc)

## 
## Iterations = 1:10000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                  Mean      SD Naive SE Time-series SE
## mu            37.5719  0.7232 0.007232       0.008258
## beta (1 - 2)   3.0545  1.4153 0.014153       0.017048
## sig2         155.6661 11.9600 0.119600       0.119600
## delta          0.2453  0.1138 0.001138       0.001367
## g              3.2200 86.0171 0.860171       0.860171
## 
## 2. Quantiles for each variable:
## 
##                   2.5%      25%      50%      75%    97.5%
## mu            36.17116  37.0878  37.5598  38.0520  39.0027
## beta (1 - 2)   0.21792   2.1280   3.0624   3.9791   5.8857
## sig2         134.06199 147.3031 154.9149 163.2828 180.6724
## delta          0.01741   0.1705   0.2452   0.3203   0.4695
## g              0.07589   0.2046   0.4071   0.9656  11.4939

table(mcmc [, 2] <0) # nam nho hon nu

## 
## FALSE  TRUE 
##  9826   174

plot(mcmc [, 2], col="blue")

# hoi qui tuyen tinh
library(rstanarm)

## Loading required package: Rcpp

## This is rstanarm version 2.21.3

## - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!

## - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.

## - For execution on a local, multicore CPU with excess RAM we recommend calling

##   options(mc.cores = parallel::detectCores())

library(BAS)
library(broom)

 m1 = stan_glm(uwes_tol ~ congviec, data=thuan)

## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.081 seconds (Warm-up)
## Chain 1:                0.085 seconds (Sampling)
## Chain 1:                0.166 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.001 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.073 seconds (Warm-up)
## Chain 2:                0.085 seconds (Sampling)
## Chain 2:                0.158 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.075 seconds (Warm-up)
## Chain 3:                0.097 seconds (Sampling)
## Chain 3:                0.172 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.084 seconds (Warm-up)
## Chain 4:                0.094 seconds (Sampling)
## Chain 4:                0.178 seconds (Total)
## Chain 4:

 summary(m1)

## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      uwes_tol ~ congviec
##  algorithm:    sampling
##  sample:       4000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 345
##  predictors:   5
## 
## Estimates:
##               mean   sd    10%   50%   90%
## (Intercept)  41.2    1.2  39.7  41.2  42.8
## congviecBB   -2.2    1.7  -4.3  -2.2   0.0
## congviecCC  -13.0    1.8 -15.3 -13.0 -10.7
## congviecDD   -2.7    2.2  -5.6  -2.8   0.2
## congviecEE   -3.5    2.7  -6.8  -3.5  -0.1
## sigma        11.6    0.4  11.1  11.6  12.2
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 37.0    0.9 35.9  37.0  38.2 
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.0  1.0  2301 
## congviecBB    0.0  1.0  2746 
## congviecCC    0.0  1.0  2764 
## congviecDD    0.0  1.0  2805 
## congviecEE    0.0  1.0  3378 
## sigma         0.0  1.0  3718 
## mean_PPD      0.0  1.0  4430 
## log-posterior 0.0  1.0  1722 
## 
## For each parameter, mcse is Monte Carlo standard error, 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).

 prior_summary(m1)

## Priors for model 'm1' 
## ------
## Intercept (after predictors centered)
##   Specified prior:
##     ~ normal(location = 37, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 37, scale = 31)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [0,0,0,...], scale = [66.78,74.75,96.63,...])
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.08)
## ------
## See help('prior_summary.stanreg') for more details

 post.r2 = bayes_R2(m1)
 summary(post.r2)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0549  0.1348  0.1574  0.1583  0.1817  0.2758

# Tương tự có thể dùng cách khác

 library(MCMCpack)

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)

## ## Copyright (C) 2003-2022 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park

## ##
## ## Support provided by the U.S. National Science Foundation

## ## (Grants SES-0350646 and SES-0350613)
## ##

 library(broom)
 
 posterior = MCMCregress(formula= uwes_tol ~ congviec, b0=0, B0=0.01,
                         data=thuan, verbose=1000)

## 
## 
## MCMCregress iteration 1 of 11000 
## beta = 
##   42.03532
##   -2.17307
##  -12.93560
##   -3.35589
##   -6.26265
## sigma2 =  138.62058
## 
## 
## MCMCregress iteration 1001 of 11000 
## beta = 
##   41.88114
##   -3.40952
##  -13.06426
##   -7.59840
##   -5.27424
## sigma2 =  121.17754
## 
## 
## MCMCregress iteration 2001 of 11000 
## beta = 
##   40.08398
##   -0.85947
##  -12.54126
##   -1.90569
##    1.18606
## sigma2 =  132.29392
## 
## 
## MCMCregress iteration 3001 of 11000 
## beta = 
##   40.27066
##   -1.10404
##   -9.75080
##   -1.06093
##   -3.70308
## sigma2 =  131.26334
## 
## 
## MCMCregress iteration 4001 of 11000 
## beta = 
##   39.08653
##    1.80259
##  -12.40554
##   -3.07740
##   -3.54418
## sigma2 =  135.02175
## 
## 
## MCMCregress iteration 5001 of 11000 
## beta = 
##   39.67919
##   -0.82773
##   -9.67508
##    0.69031
##    3.34752
## sigma2 =  135.31327
## 
## 
## MCMCregress iteration 6001 of 11000 
## beta = 
##   39.94254
##   -3.99331
##  -10.43431
##   -1.47548
##    0.31092
## sigma2 =  145.00269
## 
## 
## MCMCregress iteration 7001 of 11000 
## beta = 
##   41.27269
##   -0.73909
##  -13.54514
##   -1.16271
##   -3.52284
## sigma2 =  124.55027
## 
## 
## MCMCregress iteration 8001 of 11000 
## beta = 
##   40.68559
##   -3.38986
##  -10.94732
##   -4.20216
##   -4.20274
## sigma2 =  137.37726
## 
## 
## MCMCregress iteration 9001 of 11000 
## beta = 
##   41.43884
##   -2.14746
##  -12.25018
##   -5.67377
##   -1.65582
## sigma2 =  123.66740
## 
## 
## MCMCregress iteration 10001 of 11000 
## beta = 
##   42.57927
##   -4.76734
##  -12.63398
##   -2.00652
##   -6.87933
## sigma2 =  132.84830

 plot(posterior, col= "red")

 raftery.diag(posterior)

## 
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95 
##                                                    
##              Burn-in  Total Lower bound  Dependence
##              (M)      (N)   (Nmin)       factor (I)
##  (Intercept) 2        3710  3746         0.990     
##  congviecBB  2        3834  3746         1.020     
##  congviecCC  2        3741  3746         0.999     
##  congviecDD  2        3819  3746         1.020     
##  congviecEE  2        3710  3746         0.990     
##  sigma2      2        3741  3746         0.999

 summary(posterior)

## 
## Iterations = 1001:11000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                Mean     SD Naive SE Time-series SE
## (Intercept)  40.357  1.177  0.01177        0.01177
## congviecBB   -1.275  1.599  0.01599        0.01564
## congviecCC  -11.947  1.741  0.01741        0.01741
## congviecDD   -1.834  2.120  0.02120        0.02120
## congviecEE   -2.487  2.608  0.02608        0.02608
## sigma2      135.463 10.507  0.10507        0.10507
## 
## 2. Quantiles for each variable:
## 
##                2.5%     25%     50%      75%   97.5%
## (Intercept)  38.022  39.569  40.356  41.1532  42.673
## congviecBB   -4.451  -2.342  -1.262  -0.2110   1.803
## congviecCC  -15.355 -13.097 -11.944 -10.7725  -8.616
## congviecDD   -5.945  -3.286  -1.844  -0.3960   2.367
## congviecEE   -7.726  -4.192  -2.477  -0.7854   2.655
## sigma2      116.134 128.227 134.943 142.0363 157.667

 apply(posterior, 2, quantile, probs=c(0.025, 0.5, 0.975))

##       (Intercept) congviecBB congviecCC congviecDD congviecEE   sigma2
## 2.5%     38.02230  -4.451126 -15.354733  -5.944620  -7.726471 116.1341
## 50%      40.35624  -1.261680 -11.944058  -1.843784  -2.476791 134.9429
## 97.5%    42.67262   1.802598  -8.615542   2.366555   2.655258 157.6668

library(nnet)
require(brms)

## Loading required package: brms

## Loading 'brms' package (version 2.18.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').

## 
## Attaching package: 'brms'

## The following objects are masked from 'package:MCMCpack':
## 
##     ddirichlet, rdirichlet

## The following objects are masked from 'package:rstanarm':
## 
##     dirichlet, exponential, get_y, lasso, ngrps

## The following object is masked from 'package:stats':
## 
##     ar

m.2 = multinom(uwes_tol ~ congviec, data=thuan)

## # weights:  270 (220 variable)
## initial  value 1313.298559 
## iter  10 value 1100.219781
## iter  20 value 1078.549280
## iter  30 value 1073.955833
## iter  40 value 1073.587182
## iter  50 value 1072.741989
## iter  60 value 1072.397154
## iter  70 value 1072.244956
## iter  80 value 1072.209214
## final  value 1072.207948 
## converged

summary(m.2)

## Warning in sqrt(diag(vc)): NaNs produced

## Call:
## multinom(formula = uwes_tol ~ congviec, data = thuan)
## 
## Coefficients:
##      (Intercept) congviecBB congviecCC congviecDD congviecEE
## 9  -3.472838e+01  46.763203  -5.825211  -6.043696  -3.659674
## 11 -3.537491e+01  -7.244208  54.119134  -6.837061  -4.446413
## 12 -3.546169e+01 -10.118898  53.107362  -6.526137  -3.638332
## 13 -3.832267e+01  49.664292  57.354527 -12.093524  -8.142385
## 14 -3.599972e+01 -14.008224  55.947925  56.873378  -9.237040
## 15 -3.546169e+01 -10.118898  53.107362  -6.526137  -3.638332
## 16 -3.948881e+01 -18.627433  57.134191 -13.198976  71.952808
## 17 -3.988825e+01 -15.468212  58.632607 -13.732375  73.045039
## 18 -4.036045e+01  51.701951  59.104670  61.234228 -11.693535
## 19 -1.253521e-03  11.343143  19.256338  20.875420 -23.071226
## 21  6.931965e-01 -54.390345  18.051046 -43.073013 -31.416023
## 22 -4.093610e+01  52.970875  59.274749  61.809943 -11.927092
## 23 -1.531421e-04  12.440407  17.645916 -32.386215 -22.051813
## 24  6.931761e-01 -52.712321  18.338721 -42.351705 -31.100883
## 25  1.237439e-03  12.726620  18.742758 -31.134918  32.462591
## 26  1.098645e+00  10.936149  17.240038 -48.367652 -38.027886
## 27  1.791750e+00  11.747026  17.463266  20.468309  30.672481
## 28 -4.139494e+01  53.429668 -19.261468  62.268818  74.552335
## 29 -3.896763e+01  50.308759 -14.435440 -12.742043  71.431419
## 30  6.931490e-01  12.257892 -42.056804  20.873757 -30.517610
## 31  6.209797e-05  11.341563 -41.995088 -34.373197 -21.707369
## 32  6.931888e-01  10.648275  17.645457 -42.621901  31.770943
## 33  6.209797e-05  11.341563 -41.995088 -34.373197 -21.707369
## 34  6.931508e-01  11.341604 -45.511240  20.873866 -31.459173
## 35  6.209797e-05  11.341563 -41.995088 -34.373197 -21.707369
## 36  4.349610e-04  11.340950  18.338215  21.566425 -23.234772
## 37 -3.935782e+01  52.085666 -13.058612  60.231720  -8.331284
## 38  6.209797e-05  11.341563 -41.995088 -34.373197 -21.707369
## 39 -1.217638e-03  11.342940  17.646987  20.875223 -23.530778
## 40 -4.154444e+01  54.272260 -19.093941  62.418398  74.008790
## 41  6.931828e-01  11.341421 -45.797098  20.180660  31.771047
## 42 -4.081906e-04  12.728306 -35.929195  20.874425 -21.875828
## 43  1.609462e+00  11.341584 -50.542229  20.650648 -44.285799
## 44  6.931848e-01  12.440125 -38.828816  21.566954 -29.405826
## 45  2.485030e+00  11.159069  16.769957  20.180482  30.672743
## 46  1.791795e+00  10.648372 -52.490916  19.082051  32.281902
## 47  1.609472e+00  10.830835  17.134642 -51.779626  31.547908
## 48  1.945880e+00  10.782037  17.086063  18.927877  30.518272
## 49  6.934898e-01  12.034087  17.645695  20.180443 -31.363787
## 50  2.197260e+00  10.753781  16.141489 -57.990378 -51.287659
## 51  6.939792e-01  12.033872  18.049915  20.872974  31.769837
## 52  1.580051e-04  12.034631  18.338424  20.873636 -23.268218
## 53  6.932853e-01 -59.979824 -52.061091 -44.866032 -31.552047
## 54  1.791836e+00  11.852301  15.853865  19.082127  30.672534
## 
## Std. Errors:
##    (Intercept)   congviecBB   congviecCC   congviecDD   congviecEE
## 9    0.5808906 5.808906e-01          NaN          NaN 3.102231e-13
## 11   0.5195173 9.020741e-14 5.195173e-01 1.801733e-13 1.726122e-13
## 12   0.6553439 2.655777e-13 6.553439e-01          NaN          NaN
## 13   0.7016393 7.510491e-01 5.366352e-01 1.198930e-14          NaN
## 14   0.6891257 4.715820e-14 4.758497e-01 7.401265e-01 5.214743e-14
## 15   0.6553439 5.918435e-14 6.553439e-01          NaN          NaN
## 16   0.7475079          NaN 7.750824e-01          NaN 7.832667e-01
## 17   0.6586100          NaN 5.293971e-01 4.580452e-15 5.939063e-01
## 18   0.7757964 8.249053e-01 5.814657e-01 8.230234e-01          NaN
## 19   1.4146895 1.382596e+00 1.074849e+00 1.367523e+00 7.674115e-14
## 21   1.2247622          NaN 9.080545e-01 1.133402e-13          NaN
## 22   0.7621350 6.520807e-01 6.295098e-01 8.112278e-01          NaN
## 23   1.4143002 1.135547e+00 1.374717e+00          NaN 9.690468e-15
## 24   1.2247663          NaN 8.650543e-01          NaN          NaN
## 25   1.4138087 1.098148e+00 1.128818e+00 4.924475e-15 1.337648e+00
## 26   1.1547211 9.065474e-01 9.076302e-01 8.002609e-16 1.349747e-15
## 27   1.0801472 5.474829e-01 6.290998e-01 6.639683e-01 1.042634e+00
## 28   0.7683420 6.507897e-01 2.169961e-15 8.094209e-01 6.409515e-01
## 29   0.7637924 7.948101e-01          NaN          NaN 7.774383e-01
## 30   1.2247717 8.340119e-01          NaN 9.844502e-01 5.376286e-16
## 31   1.4142241 1.384589e+00 3.837931e-16 1.151971e-15 5.198019e-16
## 32   1.2247636 1.202051e+00 9.867580e-01 1.361417e-16 1.169919e+00
## 33   1.4142241 1.384589e+00 7.196753e-16          NaN          NaN
## 34   1.2247713 9.884332e-01 4.299312e-16 9.844241e-01 1.512853e-16
## 35   1.4142241 1.384589e+00 5.506926e-17 7.870138e-17 6.257877e-17
## 36   1.4140923 1.382198e+00 1.195140e+00 1.189991e+00 4.398138e-24
## 37   0.7178863 5.589500e-01 1.122972e-29 7.418113e-01          NaN
## 38   1.4142241 1.384589e+00 8.315916e-26 4.668707e-24 1.655982e-23
## 39   1.4146768 1.382640e+00 1.373767e+00 1.367564e+00 3.302165e-24
## 40   0.7768155 5.648318e-01 8.259971e-34 8.166138e-01 8.026633e-01
## 41   1.2247648 9.874266e-01 3.364777e-27 1.191639e+00 1.170003e+00
## 42   1.4143904 1.099964e+00 3.834772e-23 1.368600e+00 1.529825e-23
## 43   1.0954664 6.444660e-01 6.533012e-29 6.871269e-01 9.848591e-33
## 44   1.2247644 8.150446e-01 3.467834e-24 8.617715e-01 1.272161e-26
## 45   1.0408500 4.604425e-01 5.660106e-01 5.392943e-01 7.690402e-01
## 46   1.0801438 7.131699e-01 1.113615e-29 1.058647e+00 6.378737e-01
## 47   1.0954655 7.347734e-01 7.407265e-01 4.489246e-31 8.308042e-01
## 48   1.0690697 6.386979e-01 6.480510e-01 1.048302e+00 1.033046e+00
## 49   1.2247022 8.609172e-01 9.867920e-01 1.191620e+00 1.871982e-27
## 50   1.0541130 5.766155e-01 7.879919e-01 1.541303e-33 1.573299e-35
## 51   1.2246023 8.595597e-01 9.050277e-01 9.823468e-01 1.169055e+00
## 52   1.4141902 1.200608e+00 1.195259e+00 1.367182e+00 4.242023e-24
## 53   1.2247440 6.748217e-31 6.012110e-30 1.930339e-28 1.421933e-27
## 54   1.0801406 5.378667e-01 1.061967e+00 1.058313e+00 1.042572e+00
## 
## Residual Deviance: 2144.416 
## AIC: 2584.416

exp(coefficients(m.2))

##     (Intercept)   congviecBB   congviecCC   congviecDD   congviecEE
## 9  8.272884e-16 2.037047e+20 2.952182e-03 2.372772e-03 2.574089e-02
## 11 4.333820e-16 7.142998e-04 3.188902e+23 1.073253e-03 1.172053e-02
## 12 3.973601e-16 4.031054e-05 1.159403e+23 1.464653e-03 2.629616e-02
## 13 2.273400e-17 3.706191e+21 8.105032e+24 5.595636e-06 2.909426e-04
## 14 2.320172e-16 8.247183e-07 1.985523e+24 5.009496e+24 9.736541e-05
## 15 3.973601e-16 4.031054e-05 1.159403e+23 1.464653e-03 2.629616e-02
## 16 7.083138e-18 8.132208e-09 6.502251e+24 1.852496e-06 1.772995e+31
## 17 4.750674e-18 1.915319e-07 2.909495e+25 1.086689e-06 5.285149e+31
## 18 2.962648e-18 2.843624e+22 4.664789e+25 3.923623e+26 8.347612e-06
## 19 9.987473e-01 8.438488e+04 2.306330e+08 1.164339e+09 9.556387e-11
## 21 2.000099e+00 2.390965e-24 6.909866e+07 1.966202e-19 2.270881e-14
## 22 1.666017e-18 1.011484e+23 5.529632e+25 6.977777e+26 6.608910e-06
## 23 9.998469e-01 2.528135e+05 4.608121e+07 8.606875e-15 2.648617e-10
## 24 2.000058e+00 1.280354e-23 9.213087e+07 4.044717e-19 3.112132e-14
## 25 1.001238e+00 3.365896e+05 1.379991e+08 3.007993e-14 1.254076e+14
## 26 3.000097e+00 5.617059e+04 3.070812e+07 9.867239e-22 3.052803e-17
## 27 5.999946e+00 1.263771e+05 3.838839e+07 7.749498e+08 2.093579e+13
## 28 1.052941e-18 1.600332e+23 4.313708e-09 1.104089e+27 2.385986e+32
## 29 1.192813e-17 7.060189e+21 5.379822e-07 2.925506e-06 1.052619e+31
## 30 2.000004e+00 2.106371e+05 5.432027e-19 1.162405e+09 5.576612e-14
## 31 1.000062e+00 8.425158e+04 5.777836e-19 1.180077e-15 3.737741e-10
## 32 2.000083e+00 4.211987e+04 4.606007e+07 3.087049e-19 6.279793e+13
## 33 1.000062e+00 8.425158e+04 5.777836e-19 1.180077e-15 3.737741e-10
## 34 2.000007e+00 8.425506e+04 1.716799e-20 1.162532e+09 2.174977e-14
## 35 1.000062e+00 8.425158e+04 5.777836e-19 1.180077e-15 3.737741e-10
## 36 1.000435e+00 8.419998e+04 9.208430e+07 2.323695e+09 8.114592e-11
## 37 8.074505e-18 4.173662e+22 2.131654e-06 1.439805e+26 2.408625e-04
## 38 1.000062e+00 8.425158e+04 5.777836e-19 1.180077e-15 3.737741e-10
## 39 9.987831e-01 8.436775e+04 4.613058e+07 1.164110e+09 6.035498e-11
## 40 9.067291e-19 3.716575e+23 5.100429e-09 1.282230e+27 1.385508e+32
## 41 2.000071e+00 8.423961e+04 1.289951e-20 5.812318e+08 6.280445e+13
## 42 9.995919e-01 3.371578e+05 2.489711e-16 1.163182e+09 3.158266e-10
## 43 5.000120e+00 8.425335e+04 1.121473e-22 9.299563e+08 5.846846e-20
## 44 2.000075e+00 2.527421e+05 1.370438e-17 2.324925e+09 1.695166e-13
## 45 1.200148e+01 7.019759e+04 1.919109e+07 5.811280e+08 2.094127e+13
## 46 6.000212e+00 4.212397e+04 1.597662e-23 1.937445e+08 1.046771e+14
## 47 5.000171e+00 5.055588e+04 2.763634e+07 3.253823e-23 5.024380e+13
## 48 6.999791e+00 4.814812e+04 2.632588e+07 1.660629e+08 1.794390e+13
## 49 2.000685e+00 1.683983e+05 4.607105e+07 5.811057e+08 2.392656e-14
## 50 9.000323e+00 4.680669e+04 1.023669e+07 6.532791e-26 5.321726e-23
## 51 2.001665e+00 1.683620e+05 6.902053e+07 1.161495e+09 6.272848e+13
## 52 1.000158e+00 1.684899e+05 9.210356e+07 1.162265e+09 7.847678e-11
## 53 2.000276e+00 8.934976e-27 2.455589e-23 3.272879e-20 1.982076e-14
## 54 6.000457e+00 1.404070e+05 7.677961e+06 1.937593e+08 2.093689e+13

table1(~ uwes_tol | congviec, data=thuan)

	AA (N=91)	BB (N=112)	CC (N=78)	DD (N=41)	EE (N=23)	Overall (N=345)
uwes_tol
Mean (SD)	41.2 (10.4)	39.1 (11.0)	28.2 (14.1)	38.5 (9.98)	37.7 (11.8)	37.0 (12.5)
Median [Min, Max]	45.0 [5.00, 54.0]	42.5 [9.00, 54.0]	24.0 [11.0, 54.0]	43.0 [14.0, 54.0]	45.0 [16.0, 54.0]	42.0 [5.00, 54.0]

library(bayesplot)

## This is bayesplot version 1.9.0

## - Online documentation and vignettes at mc-stan.org/bayesplot

## - bayesplot theme set to bayesplot::theme_default()

##    * Does _not_ affect other ggplot2 plots

##    * See ?bayesplot_theme_set for details on theme setting

library(ggplot2)

p = ggplot(thuan, aes(x = uwes_tol, group=congviec))
p + geom_density(aes(fill = congviec), alpha=0.5) + theme_light(base_size = 12)

library(rstan)

## Loading required package: StanHeaders

## 
## rstan version 2.26.13 (Stan version 2.26.1)

## 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)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)

## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file

## 
## Attaching package: 'rstan'

## The following object is masked from 'package:coda':
## 
##     traceplot

## The following object is masked from 'package:tidyr':
## 
##     extract

library(shinystan)

## Loading required package: shiny

## 
## This is shinystan version 2.6.0

library(rstanarm)
library(psych)

## 
## Attaching package: 'psych'

## The following object is masked from 'package:rstan':
## 
##     lookup

## The following object is masked from 'package:brms':
## 
##     cs

## The following object is masked from 'package:rstanarm':
## 
##     logit

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

library(brms)

set.seed(123)

# Có 

prior0 = get_prior(uwes_tol ~ congviec, family = gaussian, data=thuan)
bayesfit1 = brm(data=thuan, uwes1 ~ congviec, prior0, 
                family = gaussian, chains=2, iter = 10000)

## Compiling Stan program...

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 3.8e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.38 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.199 seconds (Warm-up)
## Chain 1:                0.223 seconds (Sampling)
## Chain 1:                0.422 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 8e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.182 seconds (Warm-up)
## Chain 2:                0.253 seconds (Sampling)
## Chain 2:                0.435 seconds (Total)
## Chain 2:

summary(bayesfit1)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: uwes1 ~ congviec 
##    Data: thuan (Number of observations: 345) 
##   Draws: 2 chains, each with iter = 10000; warmup = 5000; thin = 1;
##          total post-warmup draws = 10000
## 
## Population-Level Effects: 
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      4.58      0.15     4.30     4.87 1.00     9019     7266
## congviecBB    -0.41      0.19    -0.79    -0.03 1.00     9355     8024
## congviecCC    -1.27      0.21    -1.69    -0.86 1.00     9166     7895
## congviecDD    -0.24      0.26    -0.74     0.27 1.00    10420     7759
## congviecEE    -0.54      0.33    -1.18     0.10 1.00    11398     7820
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.38      0.05     1.28     1.49 1.00    11476     6804
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

bayesfit1$fit

## Inference for Stan model: anon_model.
## 2 chains, each with iter=10000; warmup=5000; thin=1; 
## post-warmup draws per chain=5000, total post-warmup draws=10000.
## 
##                 mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff
## b_Intercept     4.58    0.00 0.15    4.30    4.48    4.58    4.68    4.87  8996
## b_congviecBB   -0.41    0.00 0.19   -0.79   -0.54   -0.41   -0.28   -0.03  9337
## b_congviecCC   -1.27    0.00 0.21   -1.69   -1.42   -1.27   -1.13   -0.86  9128
## b_congviecDD   -0.24    0.00 0.26   -0.74   -0.42   -0.24   -0.07    0.27 10384
## b_congviecEE   -0.54    0.00 0.33   -1.18   -0.76   -0.54   -0.32    0.10 11359
## sigma           1.38    0.00 0.05    1.28    1.34    1.38    1.41    1.49 11442
## lprior         -9.23    0.00 0.01   -9.24   -9.23   -9.23   -9.22   -9.21 13973
## lp__         -608.26    0.03 1.76 -612.61 -609.18 -607.93 -606.96 -605.86  4621
##              Rhat
## b_Intercept     1
## b_congviecBB    1
## b_congviecCC    1
## b_congviecDD    1
## b_congviecEE    1
## sigma           1
## lprior          1
## lp__            1
## 
## Samples were drawn using NUTS(diag_e) at Wed Nov  2 08:57:50 2022.
## 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).

 # Không 

prior1 = get_prior(uwes_tol ~ congviec - 1, family = gaussian, data=thuan)

bayesfit2 = brm(data = thuan, uwes_tol ~ congviec-1, prior1, 
                family = gaussian, chains = 2, iter = 10000)

## Compiling Stan program...

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 3.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.32 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.222 seconds (Warm-up)
## Chain 1:                0.201 seconds (Sampling)
## Chain 1:                0.423 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.219 seconds (Warm-up)
## Chain 2:                0.197 seconds (Sampling)
## Chain 2:                0.416 seconds (Total)
## Chain 2:

  summary(bayesfit2)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: uwes_tol ~ congviec - 1 
##    Data: thuan (Number of observations: 345) 
##   Draws: 2 chains, each with iter = 10000; warmup = 5000; thin = 1;
##          total post-warmup draws = 10000
## 
## Population-Level Effects: 
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## congviecAA    41.23      1.22    38.81    43.66 1.00    13599     7319
## congviecBB    39.06      1.07    36.98    41.15 1.00    14395     8253
## congviecCC    28.19      1.31    25.61    30.78 1.00    14234     6953
## congviecDD    38.45      1.82    34.92    41.94 1.00    14029     7577
## congviecEE    37.68      2.42    32.99    42.34 1.00    12631     8321
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    11.62      0.44    10.79    12.54 1.00    13772     7629
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

  pairs(bayesfit2)

  bayesfit2$fit

## Inference for Stan model: anon_model.
## 2 chains, each with iter=10000; warmup=5000; thin=1; 
## post-warmup draws per chain=5000, total post-warmup draws=10000.
## 
##                  mean se_mean   sd     2.5%      25%      50%      75%    97.5%
## b_congviecAA    41.23    0.01 1.22    38.81    40.42    41.22    42.02    43.66
## b_congviecBB    39.06    0.01 1.07    36.98    38.32    39.06    39.80    41.15
## b_congviecCC    28.19    0.01 1.31    25.61    27.30    28.18    29.08    30.78
## b_congviecDD    38.45    0.02 1.82    34.92    37.23    38.45    39.70    41.94
## b_congviecEE    37.68    0.02 2.42    32.99    36.00    37.68    39.34    42.34
## sigma           11.62    0.00 0.44    10.79    11.31    11.60    11.91    12.54
## lprior          -3.34    0.00 0.04    -3.41    -3.36    -3.33    -3.31    -3.27
## lp__         -1336.30    0.02 1.69 -1340.42 -1337.21 -1335.99 -1335.05 -1333.94
##              n_eff Rhat
## b_congviecAA 13549    1
## b_congviecBB 14403    1
## b_congviecCC 14317    1
## b_congviecDD 14012    1
## b_congviecEE 12608    1
## sigma        13584    1
## lprior       13517    1
## lp__          5060    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Nov  2 08:58:44 2022.
## 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(bayesfit2, ignore_prior = T, theme = ggplot2::theme())

launch_shinystan(bayesfit2, rstudio = getOption("shinystan.rstudio"))

## 
## Launching ShinyStan interface... for large models this  may take some time.

## 
## Listening on http://127.0.0.1:7837

marginal_effects(bayesfit2, probs = c(0.05, 0.95))

## Warning: Method 'marginal_effects' is deprecated. Please use
## 'conditional_effects' instead.

## Warning: Argument 'probs' is deprecated. Please use 'prob' instead.

marginal_effects(bayesfit2, probs=c(0.05,0.95), conditions=congviec)

## Warning: Method 'marginal_effects' is deprecated. Please use
## 'conditional_effects' instead.

## Warning: Argument 'probs' is deprecated. Please use 'prob' instead.

## Warning: The following variables in 'conditions' are not part of the model:
## 'conditions'

 # Kiểm tra giả thuyết
hypothesis(bayesfit2, "congviecAA > congviecCC")

## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (congviecAA)-(con... > 0    13.04      1.78    10.11    15.91        Inf
##   Post.Prob Star
## 1         1    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(bayesfit2, "congviecBB > congviecCC")

## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (congviecBB)-(con... > 0    10.87       1.7     8.07    13.67        Inf
##   Post.Prob Star
## 1         1    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(bayesfit2, "congviecAA > congviecEE")

## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (congviecAA)-(con... > 0     3.55      2.73    -0.97     8.12       9.24
##   Post.Prob Star
## 1       0.9     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

20221102 Bài tập số 3

Lê Minh Thuận

2022-11-02

R Markdown

KẾT THÚC