★ R -- 数量化 I 類 ★

 283 R -- 数量化 I 類  青木繁伸  2002/01/18 (金) 12:54


283. R -- 数量化 I 類  青木繁伸  2002/01/18 (金) 12:54
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


● 「統計学関連なんでもあり」の過去ログ--- 017 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る