成績評価に使うので、以下の課題に回答してください。

・提出締切日: 11/30
・メイルで  kubo@ees.hokudai.ac.jp に送信してください
・メイルの件名は kubostat としてください

不明な点は久保までおたずねください。

課題  1
あなたの研究に必要とされるデータ解析はどういったものでしょうか。
できるだけわかりやすく説明してみてください。
今まで使ったことのあるデータ解析の手法やソフトウェアについて
簡単に説明してみてください。

また授業でわかりにくかった点，とりあげてほしかった統計モデルや
手法などがあれば，教えてください．


課題 2
まず，今回の例題データ data をよみこんでください (下記のどちらかの方法で):
方法 1: 授業 web page (b) から data.RData をダウンローどして load("data.RData")
方法 2: 以下のようにして R で直接ダウンロード
load(url("http://hosho.ees.hokudai.ac.jp/~kubo/stat/2017/Fig/distribution/data.RData"))


課題 2
R で data をよみこんだあと，次の R コードを実行してください．
各行で何をやっているのか，簡単に説明してください．
説明するときに R の出力 (output) を使ってもいいです．

sample.mean <- mean(data)
print(sample.mean)
rpois(50, lambda = sample.mean)
mean(rpois(50, lambda = sample.mean))
replicate(10, mean(rpois(50, lambda = sample.mean)))
hist(replicate(3000, mean(rpois(50, lambda = sample.mean))))

上の [課題 2-1] で出力されたヒストグラムは何をあらわしている
のか，授業で説明した内容と関連づけて説明してみてください．
このようなシミュレイションが何の役にたつのでしょうか?


課題 3
次の 3 行の R の関数の挙動や最終行で得られたアウトプットについて，
説明してください．

d <- read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/2017/Fig/poisson/data3a.csv")
d$y.new <- rpois(nrow(d), lambda = exp(1 + 0.1 * d$x + c(0, 0.8)[d$f]))
summary(glm(y.new ~ x + f, data = d, family = poisson))

上の [課題 3-1] で得られた推定の結果を使って，モデル
による定量的な予測とデータを表示する図を使う R コードを書いて
みてください．plot() や lines() などを使ってください．


課題 4
次の 5 行の R コードは自動的に AIC 最小のモデルを見つける関数
stepAIC() の使用例です．4 行目で途中経過が表示されますが，そこ
で何が表示されているのか説明してみてください．

d <- read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/2017/Fig/poisson/data3a.csv")
fit <- glm(y ~ x + f, data = d, family = poisson)
library(MASS) # これは MASS package のよみこみ
fit.best <- stepAIC(fit)
summary(fit.best)


課題 5
次の R コードを実行すると，ロジスティック回帰の例題データ
に対して AIC 最良のモデルが選択されて，その結果が fit.best に格納
されます．

d <- read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/2017/Fig/binomial/data4a.csv")
fit <- glm(cbind(y, N - y) ~ x + f, data = d, family = binomial)
library(MASS)
fit.best <- stepAIC(fit)

また，R では logistic 関数を

logistic <- function(z) { 1 / (1 + exp(-z)) }

このように自分で定義できます．さて，

plot(d$x, d$y, pch = c(21, 19)[d$f])

この図の上に f = C の場合・f = T の場合のモデルの
予測 (生存種子数の期待値) を曲線として図示する
にはどのようなコードを書けばいいでしょうか．
lines() 関数や上で定義した logistic() 関数などを使っ
て R コードを書き，それに簡単な説明もつけてくだ
さい．


課題 6
下の R コードを実行してみてください．

# [Code 6-1]
logistic <- function(z) 1 / (1 + exp(-z))
n <- 200
n.seed <- 8
par(mfrow = c(2, 2)) # 作図画面を 2x2 に分割している
for (s in c(0.2, 0.5, 1.0, 2.0)) {
   # binomial
   prob.binom <- dbinom(0:n.seed, size = n.seed, prob = logistic(0))
   plot(
     0:n.seed, n * prob.binom,
     type = "b", # both, lines and points
     ylim = c(0, max(n * prob.binom) * 1.2),
     main = paste("s =", s),
     col = "blue"
  )
  # binomial + Gaussian
  r <- rnorm(n, mean = 0, sd = s)
  y <- rbinom(n, size = n.seed, prob = logistic(r))
  table.y <- table(y)
  print(table.y)
  px <- as.numeric(names(table.y))
  points(px, table.y, col = "red", pch = 19)
}

出力された図が何をあらわしているのか，簡単に説明してください．


課題 7

# [R-code 7]
# read data and plot
d1 <- read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/2017/Ees/g/fig/nested/d1.csv")

# install and load packages
install.packages("lme4") # if you have not installed lme4 package yet.
library(lme4) # load lme4

# plot data
plot(d1$pot, d1$y, col = rep(c("blue", "red"), each = 5))

# fit GLMM (i.e., hierarchical Bayesian model)
fit1 <- glmer(
  y ~ f + (1 | id) + (1 | pot), 
  family = poisson,
  data = d1
)
print(summary(fit1))