##########################################################
#狩野裕（未定）回帰分析
#http://www.sigmath.es.osaka-u.ac.jp/~kano/
#http://www.sigmath.es.osaka-u.ac.jp/~kano/old/lecture/u_graduate/multivariate/regression.pdf
##########################################################
library(car)
library(sem)
library(wle)
library(scatterplot3d)




#--------------------------------------------------------#
#データハンドリング
#--------------------------------------------------------#
dat <- read.csv("http://blue.zero.jp/yokumura/R/reg/page.csv",header=F)
names(dat) <- c("shuppan","genko","zuhyo")	#変数名の指定
dat						#データの確認



#--------------------------------------------------------#
#Step1:相関分析
#--------------------------------------------------------#
cor(dat)					#相関行列
scatterplot.matrix(dat,ellipse=T)		#図で確認

graphics.off()					#図を閉じる



#--------------------------------------------------------#
#Step2:単回帰分析
#--------------------------------------------------------#
reg1 <- lm(shuppan~genko,data=dat)	#原稿頁数から刷上頁数を予測
summary(reg1)

##図11.3
attach(dat)
plot(shuppan~genko,xlim=c(5,15),ylim=c(3,8))
abline(reg1)
text(genko-0.1,shuppan+0.1,rownames(dat))	#連載番号を上書き
graphics.off()					#図を閉じる
detach(dat)


##図11.6
anova(reg1)	#平方和・平均平方・F値・p値
summary(reg1)	#RootMSE(Residual standard error)・決定係数・自由度修正済み決定係数


##表 11.3 SASでの残差分析の出力に近づける
reg1$model$shuppan	#従属変数
predict(reg1)		#予測値
residuals(reg1)		#残差
			#残差のSEは不明
rstandard(reg1)		#標準化残差
rstudent(reg1)		#スチューデント化残差



##残差分析の出力関数を定義
resid.reg <- function(model){		#beginFnc
round(
data.frame(
	Dep.Value=model$model[1],		#従属変数
	Predict.Value=predict(model),		#予測値
	Residual=residuals(model),		#残差
	Student.residual=rstandard(model),	#標準化残差
	Rstudent=rstudent(model)		#スチューデント化残差
),digits=3)
}					#EndFnc

resid.reg(reg1)		#テスト

##残差分析（Rで楽な方法）

op <- par(mfrow=c(2,2))
plot(reg1)
graphics.off()


##オブザベーション13を修正
dat2 <- dat
dat2[13,1] <- 6.1

reg2 <- lm(shuppan~genko,data=dat2)	#原稿頁数から刷上頁数を予測
summary(reg2)
anova(reg2)
resid.reg(reg2)


##影響診断統計量（表11.4）

cooks.distance(reg1)		#Cook'sD
hatvalues(reg1)			#てこ比
covratio(reg1)			#Cov Ratio
dffits(reg1)			#Dffits
dfbetas(reg1)			#Dfbeta


##影響診断統計量の出力関数を定義
inf.reg <- function(model){		#beginFnc
round(
data.frame(
	CooksD =cooks.distance(model),		#Cook'sD
	HatDiagH=hatvalues(model),		#てこ比
	CovRatio=covratio(model),		#Cov Ratio
	Dffits=  dffits(model),			#Dffits
	beta=dfbetas(model)				#Dfbeta
),digits=3)
}					#EndFnc

inf.reg(reg1)				#テスト



#--------------------------------------------------------#
#Step3:重回帰分析
#--------------------------------------------------------#
cor(dat2)					#相関行列
scatterplot.matrix(dat2,ellipse=T)		#図で確認
graphics.off()

reg3 <- lm(shuppan~genko+zuhyo,data=dat2)	#重回帰分析
summary(reg3)
vif(reg3)
anova(reg3)

##標準偏回帰係数は，算出されないので自分で計算する
#(1)標準偏差と推定値から計算（お勧め）
sdd <- c(0,sd(dat2)[-1])			#標準偏差ベクトル
coef(reg3)*sdd/sd(dat2)[1]			#標準偏回帰係数

#(2)元々の変数を標準化してから計算 type1
dat3 <- data.frame(scale(dat2))
coef(lm(shuppan~genko+zuhyo,data=dat3))		#標準偏回帰係数

#(3)元々の変数を標準化してから計算 type2
coef(lm(scale(shuppan)~scale(genko)+scale(zuhyo),data=dat2))	#標準偏回帰係数

#(4)よく使うので関数を定義しておく（1）を定義したもの
std.est <- function(model){
	sdd <- c(0,sd(model$model[-1]))
	std.est <- coef(model)*sdd/sd(model$model[1])[1]
	print(std.est)
}
std.est(reg3)


anova(reg2,reg3)		#モデル比較


##図11.11を確認
attach(dat2)
s2d <- scatterplot3d(x=genko,y=zuhyo,z=shuppan,type="h")
s2d$plane3d(reg3,lty.box="solid",col="red")
detach(dat2)
graphics.off()


##変数選択
kb <- c(26,24,44,39,45,45,46,53,39,40,50,40,47,40)
dat4 <- cbind(dat2,kb)
cor(dat4)
reg4 <- lm(shuppan~genko+zuhyo+kb,data=dat4)	#重回帰分析
summary(reg4)

step(reg4)			#AICステップワイズ
mle.cp(reg4)			#Cp
mle.aic(reg4)			#AIC


#--------------------------------------------------------#
#Step4:パス解析 p.98-
#--------------------------------------------------------#
dat2


##(a) Model の指定 type 1
	#    　パス	　      母数　初期値
model <- matrix(c(
	'zuhyo -> genko',       'b',  NA,
	'zuhyo -> shuppan',     'a',  NA,
	'genko -> shuppan',     'c',  NA,
	'genko <-> genko',     'e1',  NA,
	'shuppan <-> shuppan', 'e2',  NA
  ),ncol=3,byrow=T
	)
model


##(b) Model の指定 type 2
	# 　パス	　                 #母数　#初期値
model2 <- specify.model()
    zuhyo      ->     genko,      b,      NA
    zuhyo      ->   shuppan,      a,      NA
    genko      ->   shuppan,      c,      NA
    genko     <->   genko,       e1,      NA
    shuppan   <->   shuppan,     e2,      NA

model2


##(c) 推定
sem1 <- sem(model, S=cov(dat2),N=nrow(dat2), fixed.x="zuhyo")
sem2 <- sem(model2, S=cov(dat2),N=nrow(dat2), fixed.x="zuhyo")

summary(sem1)
summary(sem2)

std.coef(sem1)
std.coef(sem2)