#
#ランキングコンジョイント分析用、Ｒプログラム
#
# データの読み込み
#
dlist <- read.table("conj1215.data", header=TRUE)
# 読み込んだデータリストからデータを切り分ける
attach(dlist, warn.conflicts=FALSE)  
#データの名前の一覧
cat("\n＜読み込んだデータの一覧＞\n")
print(names(dlist))
#
#買わない（４）が選ばれた場合、それを無視して一つ順位を詰める
#
rv1 <- (rnk1 !=4)*rnk1+
	(rnk1==4)*((rnk2 ==1)*1+(rnk2 ==2)*2+(rnk2 ==3)*3)
rv2 <- (rnk1 != 4 & rnk2 != 4)*rnk2
rv3 <- (rnk1 != 4 & rnk2 != 4 & rnk1 != 0 & rnk2 != 0)*((rnk1 != 1 & rnk2 != 1)*1+
				(rnk1 != 2 & rnk2 != 2)*2+
				(rnk1 != 3 & rnk2 != 3)*3)
   

#print(rv1)
#print(rv2)
#print(rv3)

#
# 対数尤度関数 ll の構成
#
ll <- function(p){ 	#関数の始まり
	v1 <- exp(p[1]*dis1/100+ p[2]*rec1+p[3]*eco1/100+p[4]*pri1+p[5]*var1)
	v2 <- exp(p[1]*dis2/100+ p[2]*rec2+p[3]*eco2/100+p[4]*pri2+p[5]*var2)
	v3 <- exp(p[1]*dis3/100+ p[2]*rec3+p[3]*eco3/100+p[4]*pri3+p[5]*var3)
#
#選択順位の順にvを並び替える
#
	rk1 <- (rv1 ==1)*v1+(rv1 ==2)*v2+(rv1 ==3)*v3
	rk0 <- (rv1 != 1)*v1+(rv1 != 2)*v2+(rv1 != 3)*v3
	rk2 <- (rv2 ==1)*v1+(rv2 ==2)*v2+(rv2 ==3)*v3
	rk3 <- (rv3 ==1)*v1+(rv3 ==2)*v2+(rv3 ==3)*v3
#
#１位と２位の選択確率を構成
#
	vv1 <- (rk1 !=  0)*(log(rk1/(v1+v2+v3)+(rk1 == 0)*1))
	vv2 <- (rk2 != 0)*(log(rk2/rk0+(rk2 == 0)*1))
#
#対数尤度を負の数にする（nlmが最小化するため）
#
	-sum(vv1 + vv2)
} 	#関数の終わり
#
#初期値を与えて、最大対数尤度の計算、ヘシアンを出力する
#
out <- nlm(ll, p =c(0.1,0.1,0.1,0.1,0.1), hessian = TRUE)
#
#出力準備
#
alabel <- c("dis","rec","eco","pri","var")
cat("\n＜ランキングコンジョイント　結果の概要＞ \n\n")
print(out)  #計算概要出力
#出力を切り分ける
attach(out, warn.conflicts=FALSE)
sterr = sqrt(diag(solve(hessian))) 	# 標準偏差の計算
tval = estimate/sterr 			#ｔ値の計算
#出力形式のデータフレームの作成
degreefree <- 960-2
LstOut <- data.frame(Parameter=alabel, Estimated=estimate,St_Error=sterr,
	T_value=tval,P_Value=pt(abs(tval),degreefree,lower.tail = FALSE))
cat("＜パラメータの推計結果＞ \n\n")
print(LstOut)		#推計結果の出力
estimate2 <- estimate/abs(estimate[5])
LstOut2 <- data.frame(Parameter=alabel, Yen=estimate2)
cat("\n＜効用の価格\表\示：円＞ \n\n")
print(LstOut2)
cat("\n※　警告は無視してかまいません。 \n\n")