#hierarchical_linear.R
# Basic inference.
library('rjags')
## Loading required package: coda
## Linked to JAGS 3.4.0
## Loaded modules: basemod,bugs
setwd("~/new/Bayes/Beyes/rjags/rjags1")
df <- read.csv("hierarchical_linear.csv")
head(df)
## g t x y a b
## 1 1 1 2.016819 8.905680 9.373546 -13.81636
## 2 1 2 6.607978 49.112353 9.373546 -13.81636
## 3 1 3 2.059746 2.705064 9.373546 -13.81636
## 4 1 4 3.841037 24.402757 9.373546 -13.81636
## 5 1 5 7.176185 60.663906 9.373546 -13.81636
## 6 1 6 7.774452 63.593174 9.373546 -13.81636
tail(df)
## g t x y a b
## 995 10 95 3.96388928 21.567688 9.849356 -14.39401
## 996 10 96 2.86152202 14.403575 9.849356 -14.39401
## 997 10 97 4.94417340 33.229822 9.849356 -14.39401
## 998 10 98 2.81181986 12.735653 9.849356 -14.39401
## 999 10 99 8.29269655 66.201074 9.849356 -14.39401
## 1000 10 100 0.07259715 -8.116238 9.849356 -14.39401
#library(GGally)
#pm = ggpairs(data=df,
# columns=1:4,
# upper = list(continuous = "density"),
# lower = list(combo = "facetdensity"),
# title="tips data",
# colour = "y")
#print(pm)
hierarchical_linear_code <- '
model
{
for (i in 1:N)
{
y[i] ~ dnorm(mu[i], tau)
mu[i] <- a[g[i]] * x[i] + b[g[i]]
}
for (j in 1:K)
{
a[j] ~ dnorm(mu.a, tau.a)
b[j] ~ dnorm(mu.b, tau.b)
}
mu.a ~ dnorm(0, 0.0001)
mu.b ~ dnorm(0, 0.0001)
tau <- pow(sigma, -2)
sigma ~ dunif(0, 1000)
tau.a <- pow(sigma.a, -2)
tau.b <- pow(sigma.b, -2)
sigma.a ~ dunif(0, 1000)
sigma.b ~ dunif(0, 1000)
}
'
writeLines(hierarchical_linear_code,con="hierarchical_linear.txt")
jags <- jags.model("hierarchical_linear.txt",
data = list('x' = with(df, x),
'y' = with(df, y),
'g' = with(df, g),
'N' = nrow(df),
'K' = with(df, max(g))),
n.chains = 4,
n.adapt = 5000)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 5035
##
## Initializing model
mcmc.samples <- coda.samples(jags,
c('mu.a',
'mu.b',
'tau.a',
'tau.b'),
5000)
#png(file.path('graphs', 'hierarchical', 'linear_plot1.png'))
plot(mcmc.samples)
dev.off()
## null device
## 1
summary(mcmc.samples)
##
## Iterations = 5001:10000
## Thinning interval = 1
## Number of chains = 4
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## mu.a 10.1130 0.4859 0.003436 0.003436
## mu.b -14.2447 0.3941 0.002787 0.004927
## tau.a 0.5627 0.2824 0.001997 0.002784
## tau.b 4.1039 63.5996 0.449717 1.825838
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## mu.a 9.1316 9.8176 10.1153 10.4102 11.078
## mu.b -15.0360 -14.4880 -14.2413 -14.0003 -13.458
## tau.a 0.1573 0.3535 0.5148 0.7199 1.231
## tau.b 0.2553 0.6577 1.0522 1.7144 5.904
mcmc.samples <- coda.samples(jags,
c('a',
'b',
'mu.a',
'mu.b',
'tau.a',
'tau.b'),
5000)
summary(mcmc.samples)
##
## Iterations = 10001:15000
## Thinning interval = 1
## Number of chains = 4
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## a[1] 9.4407 0.09928 0.0007020 0.001660
## a[2] 9.5103 0.09897 0.0006998 0.001882
## a[3] 11.7610 0.09405 0.0006650 0.001650
## a[4] 10.3071 0.09441 0.0006676 0.001717
## a[5] 12.1167 0.09663 0.0006833 0.001531
## a[6] 10.1247 0.09461 0.0006690 0.001562
## a[7] 8.4287 0.08850 0.0006258 0.001435
## a[8] 8.5173 0.11329 0.0008010 0.002556
## a[9] 11.1798 0.09744 0.0006890 0.001697
## a[10] 9.7636 0.09360 0.0006619 0.001588
## b[1] -14.0771 0.52555 0.0037162 0.008856
## b[2] -13.8221 0.57325 0.0040535 0.011048
## b[3] -13.4580 0.53516 0.0037841 0.009770
## b[4] -14.7887 0.54454 0.0038505 0.009697
## b[5] -14.7644 0.50337 0.0035594 0.008260
## b[6] -13.6158 0.51169 0.0036182 0.008546
## b[7] -13.8072 0.50124 0.0035443 0.008058
## b[8] -16.0060 0.63912 0.0045193 0.015441
## b[9] -13.9727 0.55048 0.0038925 0.009770
## b[10] -14.1686 0.52352 0.0037018 0.009162
## mu.a 10.1148 0.48907 0.0034583 0.003554
## mu.b -14.2491 0.40602 0.0028710 0.005182
## tau.a 0.5673 0.28472 0.0020132 0.002811
## tau.b 5.4675 79.98656 0.5655904 5.093179
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## a[1] 9.2426 9.3751 9.441 9.5075 9.633
## a[2] 9.3120 9.4450 9.512 9.5773 9.702
## a[3] 11.5748 11.6978 11.761 11.8249 11.944
## a[4] 10.1290 10.2427 10.306 10.3696 10.497
## a[5] 11.9293 12.0512 12.115 12.1811 12.309
## a[6] 9.9368 10.0620 10.126 10.1888 10.308
## a[7] 8.2554 8.3685 8.429 8.4887 8.600
## a[8] 8.2932 8.4413 8.517 8.5934 8.737
## a[9] 10.9864 11.1139 11.181 11.2462 11.369
## a[10] 9.5789 9.7013 9.764 9.8265 9.947
## b[1] -15.0929 -14.4332 -14.082 -13.7297 -13.039
## b[2] -14.9178 -14.2139 -13.834 -13.4405 -12.668
## b[3] -14.4923 -13.8254 -13.469 -13.1003 -12.387
## b[4] -15.8936 -15.1465 -14.768 -14.4130 -13.779
## b[5] -15.7797 -15.0996 -14.753 -14.4221 -13.801
## b[6] -14.5880 -13.9707 -13.624 -13.2753 -12.592
## b[7] -14.7557 -14.1524 -13.811 -13.4724 -12.811
## b[8] -17.2505 -16.4362 -16.010 -15.5779 -14.722
## b[9] -15.0299 -14.3470 -13.978 -13.6076 -12.879
## b[10] -15.1953 -14.5170 -14.168 -13.8214 -13.138
## mu.a 9.1439 9.8179 10.113 10.4067 11.095
## mu.b -15.0583 -14.5022 -14.248 -13.9995 -13.431
## tau.a 0.1535 0.3580 0.521 0.7258 1.240
## tau.b 0.2357 0.6308 1.027 1.6879 5.714
# Compare with true parameter values and OLS solution(s).
library('plyr')
ddply(df, 'g', function (df) {with(df, unique(a))})
## g V1
## 1 1 9.373546
## 2 2 9.681932
## 3 3 11.654145
## 4 4 10.509369
## 5 5 12.001719
## 6 6 9.988719
## 7 7 8.531118
## 8 8 8.589625
## 9 9 11.277952
## 10 10 9.849356
ddply(df, 'g', function (df) {with(df, unique(b))})
## g V1
## 1 1 -13.81636
## 2 2 -14.92936
## 3 3 -12.48779
## 4 4 -16.09788
## 5 5 -13.77830
## 6 6 -13.38032
## 7 7 -13.84731
## 8 8 -15.83453
## 9 9 -14.29116
## 10 10 -14.39401
lm(y ~ x, data = df)
##
## Call:
## lm(formula = y ~ x, data = df)
##
## Coefficients:
## (Intercept) x
## -14.17 10.10
lm(y ~ x : factor(g) + factor(g) - 1, data = df)
##
## Call:
## lm(formula = y ~ x:factor(g) + factor(g) - 1, data = df)
##
## Coefficients:
## factor(g)1 factor(g)2 factor(g)3 factor(g)4 factor(g)5
## -13.979 -13.555 -13.182 -15.009 -15.025
## factor(g)6 factor(g)7 factor(g)8 factor(g)9 factor(g)10
## -13.384 -13.603 -16.676 -13.866 -14.121
## x:factor(g)1 x:factor(g)2 x:factor(g)3 x:factor(g)4 x:factor(g)5
## 9.424 9.469 11.723 10.339 12.161
## x:factor(g)6 x:factor(g)7 x:factor(g)8 x:factor(g)9 x:factor(g)10
## 10.090 8.397 8.618 11.165 9.755
#hierarchical_logit.R
# Basic inference.
library('rjags')
setwd("~/new/Bayes/Beyes/rjags/rjags1")
df <- read.csv("hierarchical_logit.csv")
head(df)
## g t x y a b
## 1 1 1 2.016819 1 -0.1264538 -1.816357
## 2 1 2 9.446753 0 -0.1264538 -1.816357
## 3 1 3 6.291140 0 -0.1264538 -1.816357
## 4 1 4 2.059746 0 -0.1264538 -1.816357
## 5 1 5 6.870228 0 -0.1264538 -1.816357
## 6 1 6 7.698414 0 -0.1264538 -1.816357
tail(df)
## g t x y a b
## 1245 25 45 9.998635 0 0.2444284 -3.288627
## 1246 25 46 4.784666 0 0.2444284 -3.288627
## 1247 25 47 2.520198 0 0.2444284 -3.288627
## 1248 25 48 6.721320 0 0.2444284 -3.288627
## 1249 25 49 3.424355 0 0.2444284 -3.288627
## 1250 25 50 7.373781 0 0.2444284 -3.288627
hierarchical_logit_code <- '
model
{
for (i in 1:N)
{
y[i] ~ dbern(p[i])
logit(p[i]) <- a[g[i]] * x[i] + b[g[i]]
}
for (j in 1:K)
{
a[j] ~ dnorm(mu.a, tau.a)
b[j] ~ dnorm(mu.b, tau.b)
}
mu.a ~ dnorm(0, 0.0001)
mu.b ~ dnorm(0, 0.0001)
tau.a <- pow(sigma.a, -2)
tau.b <- pow(sigma.b, -2)
sigma.a ~ dunif(0, 1000)
sigma.b ~ dunif(0, 1000)
}
'
writeLines(hierarchical_logit_code,con="hierarchical_logit.txt")
jags <- jags.model("hierarchical_logit.txt",
data = list('x' = with(df, x),
'y' = with(df, y),
'g' = with(df, g),
'N' = nrow(df),
'K' = with(df, max(g))),
n.chains = 4,
n.adapt = 5000)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 7563
##
## Initializing model
mcmc.samples <- coda.samples(jags,
c('mu.a',
'mu.b',
'tau.a',
'tau.b'),
5000)
#png(file.path('graphs', 'hierarchical', 'logit_plot1.png'))
plot(mcmc.samples)
dev.off()
## null device
## 1
summary(mcmc.samples)
##
## Iterations = 5001:10000
## Thinning interval = 1
## Number of chains = 4
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## mu.a 0.6818 0.1871 0.001323 0.002128
## mu.b -1.9448 0.2828 0.002000 0.005771
## tau.a 1.6417 0.6804 0.004811 0.011963
## tau.b 3.4320 40.1071 0.283600 1.913829
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## mu.a 0.3267 0.5589 0.6768 0.801 1.063
## mu.b -2.5241 -2.1261 -1.9368 -1.756 -1.404
## tau.a 0.6529 1.1558 1.5309 2.004 3.274
## tau.b 0.3636 0.7221 1.0834 1.701 6.226
mcmc.samples <- coda.samples(jags,
c('a',
'b',
'mu.a',
'mu.b',
'tau.a',
'tau.b'),
5000)
summary(mcmc.samples)
##
## Iterations = 10001:15000
## Thinning interval = 1
## Number of chains = 4
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## a[1] -0.14381 0.1501 0.0010614 0.002779
## a[2] 0.83459 0.1827 0.0012916 0.004420
## a[3] 0.01724 0.1184 0.0008375 0.002323
## a[4] 1.14601 0.3123 0.0022085 0.004962
## a[5] 0.53345 0.1310 0.0009265 0.002596
## a[6] 1.09677 0.2829 0.0020001 0.005032
## a[7] -0.08632 0.1761 0.0012449 0.003080
## a[8] 0.40281 0.1153 0.0008154 0.002455
## a[9] 1.21680 0.3096 0.0021894 0.005003
## a[10] 1.81583 0.6639 0.0046945 0.011983
## a[11] -0.60732 0.3876 0.0027410 0.005244
## a[12] 1.35901 0.3248 0.0022970 0.005715
## a[13] 1.19855 0.2819 0.0019932 0.006006
## a[14] -0.75668 0.4903 0.0034668 0.006749
## a[15] 0.96851 0.2584 0.0018273 0.004061
## a[16] -0.34975 0.1648 0.0011656 0.002829
## a[17] 1.18834 0.2945 0.0020822 0.005162
## a[18] 1.76373 0.5990 0.0042356 0.009861
## a[19] 0.54969 0.1244 0.0008794 0.002806
## a[20] 1.35120 0.3413 0.0024137 0.005638
## a[21] 1.27697 0.3397 0.0024020 0.005252
## a[22] 0.18560 0.1058 0.0007483 0.002276
## a[23] 1.30135 0.3179 0.0022478 0.005840
## a[24] 0.60825 0.1380 0.0009757 0.003308
## a[25] 0.12199 0.1216 0.0008601 0.003182
## b[1] -1.68885 0.7185 0.0050809 0.013387
## b[2] -2.95649 0.7510 0.0053101 0.019012
## b[3] -1.35511 0.5621 0.0039745 0.011243
## b[4] -1.71104 0.6763 0.0047823 0.010943
## b[5] -1.72429 0.5649 0.0039947 0.011502
## b[6] -1.97501 0.6966 0.0049259 0.012396
## b[7] -2.67374 0.8255 0.0058373 0.015930
## b[8] -1.55792 0.5641 0.0039887 0.011966
## b[9] -1.93342 0.6184 0.0043729 0.009900
## b[10] -0.18626 0.8870 0.0062721 0.020208
## b[11] -2.17942 0.8447 0.0059732 0.011187
## b[12] -2.55388 0.7340 0.0051902 0.013448
## b[13] -2.77799 0.7091 0.0050138 0.016113
## b[14] -2.20453 0.8110 0.0057344 0.010437
## b[15] -1.78871 0.5847 0.0041341 0.009358
## b[16] -0.69828 0.6030 0.0042638 0.012360
## b[17] -2.03642 0.6899 0.0048782 0.012086
## b[18] -0.70604 0.8050 0.0056925 0.015722
## b[19] -2.40422 0.6328 0.0044748 0.014471
## b[20] -1.98299 0.6505 0.0046000 0.011253
## b[21] -1.81937 0.6608 0.0046727 0.010215
## b[22] -2.20194 0.6045 0.0042745 0.013233
## b[23] -2.07055 0.6981 0.0049364 0.012653
## b[24] -2.76698 0.7005 0.0049536 0.017459
## b[25] -2.71015 0.8134 0.0057516 0.022351
## mu.a 0.67815 0.1822 0.0012884 0.002119
## mu.b -1.94695 0.2826 0.0019984 0.005688
## tau.a 1.68026 0.6859 0.0048501 0.012053
## tau.b 1.89620 8.9320 0.0631587 0.409531
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## a[1] -0.46319 -0.23716 -0.13643 -0.04268 0.1368179
## a[2] 0.51294 0.70446 0.82173 0.94929 1.2282817
## a[3] -0.22548 -0.05935 0.02095 0.09734 0.2425429
## a[4] 0.60409 0.92647 1.12137 1.34062 1.8264670
## a[5] 0.28736 0.44413 0.52925 0.61529 0.8082004
## a[6] 0.60958 0.89660 1.07459 1.27127 1.7128746
## a[7] -0.45377 -0.19680 -0.07967 0.03103 0.2396713
## a[8] 0.18600 0.32292 0.39973 0.47971 0.6389534
## a[9] 0.67591 0.99740 1.19242 1.41262 1.8818299
## a[10] 0.67375 1.34325 1.76635 2.23169 3.2387731
## a[11] -1.52125 -0.82441 -0.55531 -0.33388 -0.0004804
## a[12] 0.79335 1.13075 1.33425 1.55842 2.0729175
## a[13] 0.72277 0.99734 1.17297 1.37109 1.8336626
## a[14] -1.95255 -1.01576 -0.67116 -0.40792 -0.0416342
## a[15] 0.53537 0.78487 0.94451 1.12715 1.5392791
## a[16] -0.70935 -0.45315 -0.33867 -0.23259 -0.0637120
## a[17] 0.67939 0.97829 1.16780 1.37544 1.8289126
## a[18] 0.73841 1.34102 1.71909 2.13370 3.0769842
## a[19] 0.32325 0.46334 0.54363 0.63087 0.8102618
## a[20] 0.74734 1.11411 1.33156 1.56041 2.0846720
## a[21] 0.68898 1.03525 1.25050 1.48758 2.0268495
## a[22] -0.02358 0.11479 0.18530 0.25497 0.3940425
## a[23] 0.74548 1.08037 1.27781 1.49612 1.9873561
## a[24] 0.36138 0.51216 0.60165 0.69479 0.8978831
## a[25] -0.10736 0.03966 0.11752 0.20119 0.3729608
## b[1] -3.14764 -2.14691 -1.69252 -1.22106 -0.2723121
## b[2] -4.58291 -3.42166 -2.90271 -2.42720 -1.6576505
## b[3] -2.46912 -1.72936 -1.35853 -0.98064 -0.2462790
## b[4] -3.05232 -2.14875 -1.70879 -1.27423 -0.3709368
## b[5] -2.87455 -2.08912 -1.71730 -1.34560 -0.6409497
## b[6] -3.41437 -2.41350 -1.94755 -1.51153 -0.6607409
## b[7] -4.50795 -3.16794 -2.59510 -2.09723 -1.2736467
## b[8] -2.66790 -1.93262 -1.56475 -1.18033 -0.4433711
## b[9] -3.19331 -2.33335 -1.91739 -1.51689 -0.7621501
## b[10] -1.78172 -0.80525 -0.23016 0.38359 1.6886480
## b[11] -4.00512 -2.68522 -2.12770 -1.62620 -0.6078677
## b[12] -4.12714 -3.01722 -2.50945 -2.04895 -1.2250629
## b[13] -4.33650 -3.21380 -2.71740 -2.28298 -1.5588040
## b[14] -3.94891 -2.70862 -2.15045 -1.65627 -0.7244231
## b[15] -2.99014 -2.16326 -1.77587 -1.38877 -0.6647169
## b[16] -1.86623 -1.10602 -0.69734 -0.29236 0.4748786
## b[17] -3.48179 -2.47634 -2.00912 -1.57167 -0.7531840
## b[18] -2.16956 -1.27149 -0.74898 -0.18545 0.9851871
## b[19] -3.71873 -2.81518 -2.37525 -1.97226 -1.2446812
## b[20] -3.30331 -2.40405 -1.96146 -1.55202 -0.7196094
## b[21] -3.16383 -2.24451 -1.80484 -1.38388 -0.5407680
## b[22] -3.46406 -2.58551 -2.17975 -1.78922 -1.0627083
## b[23] -3.54983 -2.50783 -2.03603 -1.60374 -0.7670339
## b[24] -4.24596 -3.20893 -2.72994 -2.27895 -1.5251974
## b[25] -4.48496 -3.21800 -2.63818 -2.13712 -1.2920037
## mu.a 0.32822 0.55754 0.67351 0.79406 1.0491863
## mu.b -2.52540 -2.12850 -1.93860 -1.76154 -1.4045579
## tau.a 0.67655 1.18950 1.56476 2.05051 3.3225406
## tau.b 0.36697 0.71926 1.06477 1.63973 4.8089136
# Compare with true parameter values and OLS solution(s).
library('plyr')
ddply(df, 'g', function (df) {with(df, unique(a))})
## g V1
## 1 1 -0.12645381
## 2 2 0.84111969
## 3 3 -0.15458464
## 4 4 1.50002880
## 5 5 0.48660048
## 6 6 1.45101281
## 7 7 0.04586309
## 8 8 0.33068167
## 9 9 1.51775423
## 10 10 1.98021396
## 11 11 -0.23732753
## 12 12 1.30570229
## 13 13 1.18960022
## 14 14 -0.55518502
## 15 15 1.00156839
## 16 16 -0.25984566
## 17 17 1.51144018
## 18 18 1.60364102
## 19 19 0.55838350
## 20 20 1.28141698
## 21 21 0.97048945
## 22 22 0.16230884
## 23 23 1.29517402
## 24 24 0.51551524
## 25 25 0.24442837
ddply(df, 'g', function (df) {with(df, unique(b))})
## g V1
## 1 1 -1.8163567
## 2 2 -3.1293631
## 3 3 -0.2327127
## 4 4 -2.6212667
## 5 5 -1.4898916
## 6 6 -2.3892372
## 7 7 -2.6557819
## 8 8 -1.3877818
## 9 9 -2.2511646
## 10 10 -0.9165701
## 11 11 -1.7093334
## 12 12 -2.4676009
## 13 13 -3.1007406
## 14 14 -2.7331119
## 15 15 -2.3547813
## 16 16 -0.8510409
## 17 17 -2.5637166
## 18 18 -0.2479644
## 19 19 -3.3549144
## 20 20 -1.3235473
## 21 21 -1.7217814
## 22 22 -2.0091490
## 23 23 -2.4882025
## 24 24 -2.4423177
## 25 25 -3.2886271
glm(y ~ x : factor(g) + factor(g) - 1,
data = df,
family = binomial(link = 'logit'))
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
##
## Call: glm(formula = y ~ x:factor(g) + factor(g) - 1, family = binomial(link = "logit"),
## data = df)
##
## Coefficients:
## factor(g)1 factor(g)2 factor(g)3 factor(g)4 factor(g)5
## -1.097e+00 -4.139e+00 -9.665e-01 -1.492e+00 -1.482e+00
## factor(g)6 factor(g)7 factor(g)8 factor(g)9 factor(g)10
## -2.012e+00 -3.679e+00 -1.235e+00 -1.983e+00 1.341e+00
## factor(g)11 factor(g)12 factor(g)13 factor(g)14 factor(g)15
## -1.430e+00 -3.868e+00 -3.908e+00 3.041e+02 -1.595e+00
## factor(g)16 factor(g)17 factor(g)18 factor(g)19 factor(g)20
## 5.006e-02 -2.253e+00 -6.930e-02 -2.590e+00 -2.275e+00
## factor(g)21 factor(g)22 factor(g)23 factor(g)24 factor(g)25
## -1.831e+00 -2.180e+00 -2.428e+00 -3.289e+00 -3.608e+00
## x:factor(g)1 x:factor(g)2 x:factor(g)3 x:factor(g)4 x:factor(g)5
## -2.264e-01 1.075e+00 -4.446e-02 1.031e+00 4.750e-01
## x:factor(g)6 x:factor(g)7 x:factor(g)8 x:factor(g)9 x:factor(g)10
## 1.080e+00 9.159e-02 3.375e-01 1.219e+00 1.391e+00
## x:factor(g)11 x:factor(g)12 x:factor(g)13 x:factor(g)14 x:factor(g)15
## -9.635e-01 1.833e+00 1.575e+00 -3.996e+03 8.658e-01
## x:factor(g)16 x:factor(g)17 x:factor(g)18 x:factor(g)19 x:factor(g)20
## -4.835e-01 1.243e+00 1.768e+00 5.710e-01 1.484e+00
## x:factor(g)21 x:factor(g)22 x:factor(g)23 x:factor(g)24 x:factor(g)25
## 1.273e+00 1.849e-01 1.426e+00 6.859e-01 2.429e-01
##
## Degrees of Freedom: 1250 Total (i.e. Null); 1200 Residual
## Null Deviance: 1733
## Residual Deviance: 700.4 AIC: 800.4