qt1 <- function(file)
{
data <- read.table(file, header=TRUE)
ncat <- dim(data)[2]-1
mx <- apply(data[1:ncat], 2, max)
start <- c(0, cumsum(mx)[1:(ncat-1)])
nobe <- sum(mx)
# ダミー変数を使ったデータ行列に変換
x <- t(apply(data, 1,
function(obs)
{
zeros <- rep(0, nobe)
zeros[start+obs[1:ncat]] <- 1
c(zeros[-start-1], obs[ncat+1])
}
))
# 重回帰分析として解く
a <- cov(x)
ndim <- nobe-ncat
B <- solve(a[1:ndim,1:ndim], a[ndim+1,1:ndim])
m <- apply(x, 2, mean)
const <- m[ndim+1]-sum(B*m[1:ndim])
# 数量化 I 類としての解に変換
ncase <- dim(data)[1]
s <- apply(x, 2, sum)
name <- coef <- NULL
en <- 0
for (i in 1:ncat) {
st <- en+1
en <- st+mx[i]-2
target <- st:en
temp.mean <- sum(s[target]*B[target])/ncase
const <- const+temp.mean
coef <- c(coef, -temp.mean, B[target]-temp.mean)
name <- c(name, paste("Cat", i, 1:mx[i], sep="-"))
}
coef <- c(coef,const)
names(coef) <- c(name, "Constant")
res <- as.matrix(coef)
colnames(res) <- c("Coefficient")
return(res)
}
# ==========================
# テストデータ qt1.data
d1 d2 d3 y
1 2 2 6837
3 2 2 7397
1 2 2 7195
1 1 1 6710
2 3 2 6670
1 3 3 6279
2 2 2 6601
1 1 1 4929
3 2 2 5471
1 1 2 6164
1 1 1 5095
1 1 1 4766
2 1 1 6525
1 1 1 5087
1 3 2 6060
# 実行結果
> qt1("qt1.data")
Coefficient
Cat-1-1 -67.3037
Cat-1-2 450.0296
Cat-1-3 -338.5259
Cat-2-1 -168.8741
Cat-2-2 372.3481
Cat-2-3 -226.5407
Cat-3-1 -450.4444
Cat-3-2 281.1111
Cat-3-3 453.7778
Constant 6119.0667