#### 第5章 ####
#### 5.2節 ####
# 線形回帰モデル
fp <- read.table("foreignpurchase.csv", sep=",", header=T)
fp.lm <- lm(internet~age, data=fp)
summary(fp.lm)

# マルチレベル・線形回帰モデル　最尤推定法
library(lme4)

# 固定効果：年齢，ランダム効果：定数項
fp.lmer1 <- lmer(internet~age+(1|nation), data=fp)
summary(fp.lmer1)
# 固定効果のみを表示
fixef(fp.lmer1) 
# ランダム効果のみを表示
ranef(fp.lmer1)

# 固定効果：定数項，ランダム効果：年齢
fp.lmer2 <- lmer(internet~(0+age|nation), data=fp)
summary(fp.lmer2)
# 固定効果のみを表示
fixef(fp.lmer2) 
# ランダム効果のみを表示
ranef(fp.lmer2)

# 分散分析
anova(fp.lmer1,fp.lmer2)

# ランダム効果モデル
fp.lmer3 <- lmer(internet~gend+(age|nation), data=fp)
summary(fp.lmer3)
# 固定効果のみを表示
fixef(fp.lmer3) 
# ランダム効果のみを表示
ranef(fp.lmer3)

# マルチレベル・線形回帰モデル　ベイズ推定法
# lme4のmcmcsamp()を使う例
# library(coda)
# set.seed(66)
fp.mcmc1 <- mcmcsamp(fp.lmer1,n=10000)
summary(fp.mcmc1,burnin=1000)
str(fp.mcmc1, burnin=1000)

# densityplot(fp.mcmc1)
densityplot(fp.mcmc1,burnin=1000)
xyplot(fp.mcmc1)

fp.mcmc2 <- mcmcsamp(fp.lmer2, n=50000)
summary(fp.mcmc2)

fp.mcmc3 <- mcmcsamp(fp.lmer3, n=50000)
summary(fp.mcmc3)

# baysmパッケージを使う例
library(bayesm)
# mcmcの設定
R=10000
keep=1
# 事前分布を設定
reg=levels(factor(fp$area))
nreg=length(reg)
nobs=nrow(fp)
nvar=3 	#2説明変数＋定数項
# 変数を設定
regdata=NULL
for (j in 1:nreg) {
        y=fp$internet[fp$area==reg[j]]
        iota=c(rep(1,length(y)))
        X=cbind(iota, fp$age[fp$area==reg[j]], fp$gend[fp$area==reg[j]])
        regdata[[j]]=list(y=y,X=X)}

Z=matrix(c(rep(1,nreg)),ncol=1)
Data1=list(regdata=regdata,Z=Z)
Mcmc1=list(R=R,keep=1)
set.seed(66)
out1=rhierLinearModel(Data=Data1,Mcmc=Mcmc1)

summary(out1$Deltadraw)
summary(out1$Vbetadraw)
summary(t(out1$betadraw[1,,]))

# 稼働検査期間を除く場合
summary(t(out1$betadraw[1,,]), burnin=1000)

#### 5.2.2 二項ロジットモデル ####
# library(bayesm)
#選択肢属性
ChoiAttr<-read.table("blogitChoiceAttr.csv",sep=",",header=T)
#個人（グループ属性）
IndAttr<-read.table("blogitIndAttr.csv",sep=",",header=T)
#MCMC設定
R=10000
keep=1
#個人（グループ数）
reg=levels(factor(ChoiAttr$id))
nreg=length(reg)
#個人（グループ）毎のサンプル数
nobs=(nrow(ChoiAttr)/nreg)
#nvar=6 #説明変数＋定数項
nz=ncol(IndAttr) #個人属性数

lgtdata=NULL
for (j in 1:nreg) {
	y=ChoiAttr$choice[ChoiAttr$id==reg[j]]
	iota=c(rep(1,length(y)))
	X=as.matrix(ChoiAttr[ChoiAttr[,1]==reg[j],c(4:8)])
	X=cbind(iota,X)
	lgtdata[[j]]=list(y=y,X=X)
}
Z=NULL
Z=as.matrix(cbind(c(rep(1,nreg)),IndAttr))
Data2=list(lgtdata=lgtdata,Z=Z)
Mcmc2=list(R=R,keep=1,sbeta=0.2)

set.seed(66)
out2=rhierBinLogit(Data=Data2,Mcmc=Mcmc2)

summary(out2$Deltadraw,tvalues=as.vector(Delta))
summary(out2$Vbetadraw,tvalues=as.vector(Vbeta[upper.tri(Vbeta,diag=TRUE)]))

# 稼働検査期間を除いたパラメータ計算結果
summary(t(out2$betadraw[1,,]), burnin=1000)

# 図5.4
plot(out2$Deltadraw)
# 図5.5
plot(out2$betadraw)
# 図5.6
plot(out2$Vbetadraw)

#### 二項ロジットモデル分析用データの作成 ####
# P. RossiのbaysmパッケージのrhirBinLogitに記載されている
# サンプルデータ作成部分をアレンジ
set.seed(66)
nvar=5 ## number of coefficients
nlgt=30 ## number of cross-sectional units サンプル数
nobs=20 ## number of observations per unit 質問数
nz=2
Z=NULL
Z=matrix(c(runif(nlgt,min=-1,max=1),runif(nlgt,min=-1,max=1)),nrow=nlgt,ncol=nz)
Delta=matrix(c(-2,-1,0,1,2,-1,1,-.5,.5,0),nrow=nz,ncol=nvar)
iota=matrix(1,nrow=nvar,ncol=1)
Vbeta=diag(nvar)+.5*iota%*%t(iota)
## simulate data
lgtdata=NULL
for (i in 1:nlgt)
{ beta=t(Delta)%*%Z[i,]+as.vector(t(chol(Vbeta))%*%rnorm(nvar))
X=matrix(runif(nobs*nvar),nrow=nobs,ncol=nvar)
prob=exp(X%*%beta)/(1+exp(X%*%beta))
unif=runif(nobs,0,1)
y=ifelse(unif<prob,1,0)
lgtdata[[i]]=list(y=y,X=X,beta=beta)
}

yy=NULL
xx=NULL
for(i in 1:nlgt){
yy=rbind(yy,lgtdata[[i]]$y)
xx=rbind(xx,lgtdata[[i]]$X)
}

ChoiceAttr=NULL
ChoiceAttr=matrix(0,nrow=nobs*nlgt,ncol=nvar+3)
for(i in 1:nlgt){ChoiceAttr[(((i-1)*nobs+1):(i*nobs)),1]=i}
ChoiceAttr[,2]=rep(1:nobs)
ChoiceAttr[,3]=yy
ChoiceAttr[,4:(nvar+3)]=xx
colnames(ChoiceAttr)=c("id","qid","choice", "var1","var2","var3","var4","var5")

#Z2=NULL
#Z2=matrix(c(runif(nreg,min=-1,max=1),runif(nreg,min=-1,max=1)),nrow=nobs,ncol=nz)
IndAttr=Z
colnames(IndAttr)=c("ind1","ind2")

write.table(ChoiceAttr,"blogitChoiceAttr.csv",sep=",")
write.table(IndAttr,"blogitIndAttr.csv",sep=",")

ChoiAttr=NULL
IndAttr=NULL

#### 5.2.3 多項ロジットモデル ####
# library(bayesm)
# 選択肢属性
mnlChoiAttr<-read.table("mnlChoiceAttr.csv",sep=",",header=T)
# 個人（グループ属性）
mnlIndAttr<-read.table("mnlIndAttr.csv",sep=",",header=T)
# MCMC設定
R=10000
keep=5
# 個人数
reg=levels(factor(mnlChoiAttr$id))
nreg=length(reg)
# 個人毎の質問数
nobs=(nrow(mnlChoiAttr)/nreg)
p=3  # 選択肢数
na=1 # 選択肢共通変数の数
nz=ncol(mnlIndAttr) #個人属性数

lgtdata=NULL
for (j in 1:nreg) {
	y=mnlChoiAttr$choice[mnlChoiAttr$id==reg[j]]
	Xa=as.matrix(mnlChoiAttr[mnlChoiAttr[,1]==reg[j],c(4:6)])
	X=createX(p, na=na, nd=NULL, Xa=Xa, Xd=NULL, base=1)
	lgtdata[[j]]=list(y=y,X=X)
}
Z=NULL
Z=as.matrix(mnlIndAttr)
Data3=list(p=p,lgtdata=lgtdata,Z=Z)
Prior3=list(ncomp=3)
Mcmc3=list(R=R,keep=1)

set.seed(66)
out3=rhierMnlRwMixture(Data=Data3,Mcmc=Mcmc3,Prior=Prior3)

summary(out3$Deltadraw)
summary(out3$Vbetadraw)
summary(t(out3$betadraw[1,,]), burnin=20000)

# 図5.7
plot(out3$Deltadraw)
# 図5.8
plot(out3$betadraw)

#### 多項ロジットモデル分析用データの作成 ####
set.seed(66)
p=3 # num of choice alterns
ncoef=3 # num of coefficient
nlgt=30 # num of cross sectional units サンプル数
nobs=20 # 質問数
na=1 #選択肢共通変数の数
nz=2 #個人属性の数
Z=matrix(runif(nz*nlgt),ncol=nz)
Z=t(t(Z)-apply(Z,2,mean)) # demean Z
ncomp=3 # no of mixture components
Delta=matrix(c(1,0,1,0,1,2),ncol=2)
comps=NULL
comps[[1]]=list(mu=c(0,-1,-2),rooti=diag(rep(1,3)))
comps[[2]]=list(mu=c(0,-1,-2)*2,rooti=diag(rep(1,3)))
comps[[3]]=list(mu=c(0,-1,-2)*4,rooti=diag(rep(1,3)))
pvec=c(.4,.2,.4)
simmnlwX= function(n,X,beta) {
## simulate from MNL model conditional on X matrix
k=length(beta)
Xbeta=X%*%beta
j=nrow(Xbeta)/n
Xbeta=matrix(Xbeta,byrow=TRUE,ncol=j)
Prob=exp(Xbeta)
iota=c(rep(1,j))
denom=Prob%*%iota
Prob=Prob/as.vector(denom)
y=vector("double",n)
ind=1:j
for (i in 1:n)
{yvec=rmultinom(1,1,Prob[i,]); y[i]=ind%*%yvec}
return(list(y=y,X=X,beta=beta,prob=Prob))
}
## simulate data
simlgtdata=NULL
ni=rep(nobs,nlgt)
for (i in 1:nlgt)
{ betai=Delta%*%Z[i,]+as.vector(rmixture(1,pvec,comps)$x)
Xa=matrix(runif(ni[i]*p,min=0,max=2),ncol=p)
X=createX(p,na=1,nd=NULL,Xa=Xa,Xd=NULL,base=1)
outa=simmnlwX(ni[i],X,betai)
simlgtdata[[i]]=list(y=outa$y,X=X,beta=betai)
}

yy=NULL
xx=NULL
for(i in 1:nlgt){
yy=cbind(yy,simlgtdata[[i]]$y)
xx=rbind(xx,t(matrix(simlgtdata[[i]]$X[,3],nrow=p)))
}
yy=c(yy)
na=1

mnlChoiceAttr=NULL
mnlChoiceAttr=matrix(0,nrow=nobs*nlgt,ncol=na*p+3)
for(i in 1:nlgt){mnlChoiceAttr[(((i-1)*nobs+1):(i*nobs)),1]=i}
mnlChoiceAttr[,2]=rep(1:nobs)
mnlChoiceAttr[,3]=yy
mnlChoiceAttr[,4:(na*p+3)]=xx
colnames(mnlChoiceAttr)=c("id","qid","choice", "costA","costB","costC")
summary(mnlChoiceAttr)

mnlIndAttr=Z
colnames(mnlIndAttr)=c("ind1","ind2")

write.table(mnlChoiceAttr,"C:/Program Files/R/R-2.7.2/mnlChoiceAttr.csv",sep=",")
write.table(mnlIndAttr,"C:/Program Files/R/R-2.7.2/mnlIndAttr.csv",sep=",")

