OpenStreetMap
OpenStreetMap Japan
Creative Commonsのライセンスにて提供されている地図。QGISからはOpen Layers Pluginから利用可能。
地図データは、OpenStreetMapにより、CC-BY-SA 2.0のライセンスで提供されています。
]]>
![[コハクチョウ]](http://ito-hi.blog.so-net.ne.jp/_images/blog/_5c3/ito-hi/P1295315.jpg)
![[ハシビロガモ]](http://ito-hi.blog.so-net.ne.jp/_images/blog/_5c3/ito-hi/P1295308.jpg)
![[ハイイロチュウヒ]](/_images/blog/_5c3/ito-hi/m_P1285297.jpg)
library(boot)
n <- 10000
m <- 5000
c <- 3
x <- vector("list", c)
for (i in 1:c) {
x[[i]] <- runif(n, 0, 1)
}
boot.fun <- function(data, index) {
mean(data[index])
}
library(parallel)
r <- vector("list", c)
for (i in 1:c) {
r[[i]] <- mcparallel(boot(x[[i]], boot.fun, m))
}
# wait
res <- mccollect(r)
print(res)
$`99644`
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = x[[i]], statistic = boot.fun, R = m)
Bootstrap Statistics :
original bias std. error
t1* 0.4946006 2.003117e-05 0.002872614
$`99645`
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = x[[i]], statistic = boot.fun, R = m)
Bootstrap Statistics :
original bias std. error
t1* 0.4990066 -2.169221e-05 0.002894828
$`99646`
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = x[[i]], statistic = boot.fun, R = m)
Bootstrap Statistics :
original bias std. error
t1* 0.4987062 2.583383e-05 0.002893329
]]>library(MASS) data(bacteria)
library(lme4)
bac.lmer <- lmer(y ~ week + trt + (1|ID), family = binomial,
data = bacteria)
print(bac.lmer)
Generalized linear mixed model fit by the Laplace approximation
Formula: y ~ week + trt + (1 | ID)
Data: bacteria
AIC BIC logLik deviance
207.8 224.7 -98.89 197.8
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 1.3146 1.1465
Number of obs: 220, groups: ID, 50
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.1440 0.5290 5.943 2.79e-09 ***
week -0.1437 0.0483 -2.975 0.00293 **
trtdrug -1.3202 0.6252 -2.111 0.03473 *
trtdrug+ -0.7956 0.6401 -1.243 0.21395
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) week trtdrg
week -0.552
trtdrug -0.628 0.071
trtdrug+ -0.596 0.039 0.489
library(glmmML)
bac.ml <- glmmML(I(y == "y") ~ week + trt, cluster = ID, family=binomial,
data = bacteria)
print(bac.ml)
Call: glmmML(formula = I(y == "y") ~ week + trt, family = binomial, data = bacteria, cluster = ID)
coef se(coef) z Pr(>|z|)
(Intercept) 3.1440 0.61916 5.078 3.82e-07
week -0.1437 0.05092 -2.822 4.77e-03
trtdrug -1.3202 0.64199 -2.056 3.97e-02
trtdrug+ -0.7955 0.65172 -1.221 2.22e-01
Scale parameter in mixing distribution: 1.147 gaussian
Std. Error: 0.3782
LR p-value for H_0: sigma = 0: 0.007011
Residual deviance: 197.8 on 215 degrees of freedom AIC: 207.8
library(MCMCglmm)
bac.mcmc <- MCMCglmm(y ~ week + trt, random = ~ID, family = "categorical",
data = bacteria,
prior = list(R = list(V = 1, nu = 0),
G = list(G1 = list(V = 1, nu = 0))),
nitt = 60000, thin = 50, burnin = 10000,
verbose = FALSE)
summary(bac.mcmc)
Iterations = 10001:59951
Thinning interval = 50
Sample size = 1000
DIC: 5.125227
G-structure: ~ID
post.mean l-95% CI u-95% CI eff.samp
ID 5145 1.448e-12 17051 10.68
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 21231 3293 39949 129.4
Location effects: y ~ week + trt
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 231.131 121.768 355.445 220.3 <0.001 ***
week -10.284 -18.973 -2.641 501.6 0.006 **
trtdrug -94.833 -187.894 -7.489 586.1 0.030 *
trtdrug+ -54.059 -146.717 40.176 1000.0 0.206
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
bac.mcmc2 <- MCMCglmm(y ~ week + trt, random = ~ID, family = "categorical",
data = bacteria,
prior = list(R = list(V = 1, nu = 1, fix = 1),
G = list(G1 = list(V = 1, nu = 0))),
nitt = 60000, thin = 50, burnin = 10000,
verbose = FALSE)
summary(bac.mcmc2)
Iterations = 10001:59951
Thinning interval = 50
Sample size = 1000
DIC: 192.108
G-structure: ~ID
post.mean l-95% CI u-95% CI eff.samp
ID 2.874 0.002028 7.083 93.78
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 1 1 1 0
Location effects: y ~ week + trt
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 3.90023 2.32113 5.55957 117.7 <0.001 ***
week -0.17786 -0.31257 -0.06714 325.2 0.002 **
trtdrug -1.62696 -3.40207 -0.06575 291.0 0.044 *
trtdrug+ -0.93034 -2.56224 0.94574 501.9 0.276
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
library(R2jags)
set.seed(31415)
model.file <- system.file(package="R2WinBUGS", "model", "schools.txt")
data(schools)
J <- nrow(schools)
y <- schools$estimate
sigma.y <- schools$sd
data <- list(J = J, y = y, sigma.y = sigma.y)
parameters <- c("theta", "mu.theta", "sigma.theta")
n.chain <- 3
inits <- vector("list", n.chain)
inits[[1]] <- list(.RNG.seed = 123,
.RNG.name = "base::Mersenne-Twister",
theta=rnorm(J, 0, 100),
mu.theta=rnorm(1, 0, 100),
sigma.theta=runif(1, 0, 100))
inits[[2]] <- list(.RNG.seed = 1234,
.RNG.name = "base::Mersenne-Twister",
theta=rnorm(J, 0, 100),
mu.theta=rnorm(1, 0, 100),
sigma.theta=runif(1, 0, 100))
inits[[3]] <- list(.RNG.seed = 12345,
.RNG.name = "base::Mersenne-Twister",
theta=rnorm(J, 0, 100),
mu.theta=rnorm(1, 0, 100),
sigma.theta=runif(1, 0, 100))
fit <- jags(data = data, inits = inits,
parameters.to.save = parameters,
model.file = model.file,
n.chain = 3, n.iter = 1000)
samples <- autojags(fit)
print(samples)
Inference for Bugs model at "/Library/Frameworks/R.framework/Versions/2.14/Resources/library/R2WinBUGS/model/schools.txt", fit using jags,
3 chains, each with 1000 iterations (first 1000 discarded)
n.sims = 3000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
mu.theta 7.759 5.128 -2.174 4.677 7.711 10.805 17.824 1.012 560
sigma.theta 6.535 6.974 0.083 2.157 4.918 8.248 26.692 1.049 68
theta[1] 11.071 8.244 -2.292 6.355 9.560 14.787 31.065 1.021 170
theta[2] 7.699 5.855 -3.918 4.047 7.649 11.079 19.663 1.006 1300
theta[3] 5.895 7.670 -11.632 1.959 7.037 10.425 19.763 1.011 3000
theta[4] 7.547 6.209 -5.447 3.840 7.702 11.176 19.743 1.005 1200
theta[5] 5.467 6.351 -9.158 1.769 6.273 9.409 16.885 1.003 2200
theta[6] 6.229 6.634 -8.146 2.383 6.993 10.211 18.600 1.005 3000
theta[7] 10.322 6.596 -1.372 6.399 9.397 13.938 25.456 1.011 300
theta[8] 8.070 7.608 -6.682 3.769 7.924 11.827 24.152 1.012 530
deviance 60.357 2.187 57.077 59.124 59.897 61.210 65.730 1.012 810
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = var(deviance)/2)
pD = 2.4 and DIC = 62.7
DIC is an estimate of expected predictive error (lower deviance is better).
]]>![[オオハクチョウ]](http://ito-hi.blog.so-net.ne.jp/_images/blog/_5c3/ito-hi/P1085249.jpg)
![[ジョウビタキ]](http://ito-hi.blog.so-net.ne.jp/_images/blog/_5c3/ito-hi/P1085253-798b2.jpg)
![[ホシハジロ]](http://ito-hi.blog.so-net.ne.jp/_images/blog/_5c3/ito-hi/P1085258-b9755.jpg)
![[キンクロハジロ]](http://ito-hi.blog.so-net.ne.jp/_images/blog/_5c3/ito-hi/P1085259.jpg)
![[牛久沼]](/_images/blog/_5c3/ito-hi/m_Q1075552.jpg)
Y ~ dbin(p, N)で、Nを推定するときのバグを修正。
CFLAGS='-arch x86_64 -O2' \ CXXFLAGS='-arch x86_64 -O2' \ FFLAGS='-arch x86_64 -O2' \ LDFLAGS='-arch x86_64 -O2' \ ./configure --with-blas='-lblas' --with-lapack='-llapack' make -j4 sudo make install]]>
![[牛久沼から筑波山]](/_images/blog/_5c3/ito-hi/m_P1035220.jpg)
![[マガモ]](http://ito-hi.blog.so-net.ne.jp/_images/blog/_5c3/ito-hi/P1035225-4e573.jpg)
##
## データは狩野・三浦(2002) p.74より
##
library(sem)
n <- 220
R <- readMoments(diag = FALSE,
names = c("X1", "X2", "X3", "X4", "X5", "X6"))
0.439
0.410 0.351
0.288 0.354 0.164
0.329 0.320 0.190 0.595
0.248 0.329 0.181 0.470 0.484
# model 1
model1 <- specifyModel()
X1 <- F1, lambda1
X2 <- F1, lambda2
X3 <- F1, lambda3
X4 <- F1, lambda4
X5 <- F1, lambda5
X6 <- F1, lambda6
F1 <-> F1, NA, 1
X1 <-> X1, e1
X2 <-> X2, e2
X3 <-> X3, e3
X4 <-> X4, e4
X5 <-> X5, e5
X6 <-> X6, e6
s1 <- sem(model1, S = R, N = n)
summary(s1)
# model 2
model2 <- specifyModel()
X1 <- F1, lambda1
X2 <- F1, lambda2
X3 <- F1, lambda3
X4 <- F2, lambda4
X5 <- F2, lambda5
X6 <- F2, lambda6
F1 <-> F1, NA, 1
F2 <-> F2, NA, 1
F1 <-> F2, phi12
X1 <-> X1, e1
X2 <-> X2, e2
X3 <-> X3, e3
X4 <-> X4, e4
X5 <-> X5, e5
X6 <-> X6, e6
s2 <- sem(model2, S = R, N = n)
summary(s2)
結果
> summary(s1) Model Chisquare = 53.245 Df = 9 Pr(>Chisq) = 2.6267e-08 Chisquare (null model) = 318.29 Df = 15 Goodness-of-fit index = 0.91393 Adjusted goodness-of-fit index = 0.79916 RMSEA index = 0.14983 90% CI: (0.11244, 0.1899) Bentler-Bonnett NFI = 0.83272 Tucker-Lewis NNFI = 0.75687 Bentler CFI = 0.85412 SRMR = 0.084726 AIC = 77.245 AICc = 54.752 BIC = 117.97 CAIC = -4.298 Normalized Residuals Min. 1st Qu. Median Mean 3rd Qu. Max. -1.2500 -0.8320 -0.0273 0.1730 0.2630 3.5200 R-square for Endogenous Variables X1 X2 X3 X4 X5 X6 0.2406 0.2832 0.1188 0.5310 0.5511 0.3913 Parameter Estimates Estimate Std Error z value Pr(>|z|) lambda1 0.49050 0.071142 6.8947 5.3996e-12 X1 <--- F1 lambda2 0.53214 0.070323 7.5672 3.8150e-14 X2 <--- F1 lambda3 0.34471 0.073521 4.6886 2.7509e-06 X3 <--- F1 lambda4 0.72866 0.066294 10.9913 4.2063e-28 X4 <--- F1 lambda5 0.74235 0.066047 11.2397 2.6039e-29 X5 <--- F1 lambda6 0.62557 0.068334 9.1545 5.4602e-20 X6 <--- F1 e1 0.75941 0.079194 9.5893 8.8693e-22 X1 <--> X1 e2 0.71683 0.076513 9.3687 7.3449e-21 X2 <--> X2 e3 0.88117 0.087286 10.0953 5.7980e-24 X3 <--> X3 e4 0.46905 0.064961 7.2205 5.1814e-13 X4 <--> X4 e5 0.44892 0.064531 6.9566 3.4848e-12 X5 <--> X5 e6 0.60866 0.070311 8.6567 4.8567e-18 X6 <--> X6 Iterations = 13
> summary(s2) Model Chisquare = 7.7844 Df = 8 Pr(>Chisq) = 0.45481 Chisquare (null model) = 318.29 Df = 15 Goodness-of-fit index = 0.98816 Adjusted goodness-of-fit index = 0.96892 RMSEA index = 0 90% CI: (NA, 0.077836) Bentler-Bonnett NFI = 0.97554 Tucker-Lewis NNFI = 1.0013 Bentler CFI = 1 SRMR = 0.031955 AIC = 33.784 AICc = 9.5513 BIC = 77.902 CAIC = -43.365 Normalized Residuals Min. 1st Qu. Median Mean 3rd Qu. Max. -1.0900 -0.2600 -0.0312 -0.0180 0.1260 1.1300 R-square for Endogenous Variables X1 X2 X3 X4 X5 X6 0.4727 0.4506 0.2842 0.5730 0.6069 0.3935 Parameter Estimates Estimate Std Error z value Pr(>|z|) lambda1 0.68756 0.075701 9.0826 1.0599e-19 X1 <--- F1 lambda2 0.67129 0.075566 8.8836 6.4756e-19 X2 <--- F1 lambda3 0.53312 0.075615 7.0504 1.7834e-12 X3 <--- F1 lambda4 0.75700 0.067113 11.2794 1.6581e-29 X4 <--- F2 lambda5 0.77901 0.066942 11.6371 2.6692e-31 X5 <--- F2 lambda6 0.62732 0.068494 9.1587 5.2508e-20 X6 <--- F2 phi12 0.59481 0.071772 8.2876 1.1558e-16 F2 <--> F1 e1 0.52726 0.082350 6.4027 1.5270e-10 X1 <--> X1 e2 0.54937 0.081687 6.7253 1.7520e-11 X2 <--> X2 e3 0.71578 0.082010 8.7279 2.5938e-18 X3 <--> X3 e4 0.42695 0.067242 6.3495 2.1607e-10 X4 <--> X4 e5 0.39314 0.067591 5.8164 6.0135e-09 X5 <--> X5 e6 0.60647 0.070501 8.6022 7.8212e-18 X6 <--> X6 Iterations = 17]]>