主座標分析

• 対象間の類似度を表す対称行列を固有値・固有ベクトルに分解する
• 参考サイトは
  • こちらと
  • こちら
    • 以下のソースは、こちら（群馬　青木先生）からの改変(ほとんど変わってないですけれど)
• Rで実施
  • 類似度行列は、距離行列の値を負にしたものでもよい
    • 距離行列の条件を満足する行列の簡易作成法はこちら

princoRY <- function(s,thres=0,noPlotAxis=2,outfile,...)#類似度行列、有効固有値下限、プロットの軸数、出力ファイル名、plot関数の引数
{
        n <- nrow(s)# 行列のサイズ
        m <- colSums(s)/n #列別の平均。対称行列なので、行別の平均でもある
        h <- s+sum(s)/n^2-outer(m, m, "+") 
         #２重中心化
         # 第i,j要素から、第i行,第j列の平均を引き、(引きすぎた分である)、全体平均を足し戻している
         # 第２項 sum(s)/n^2は、全体の平均
         # 第３項 outer(m,m,"+") は、第i行平均と第j列平均の和
        res <- eigen(h) # 固有値・固有ベクトルを求める
        res$values <- res$values[res$values > thres]  # 固有値のうち、与えた値を下限として、解の個数を決める
        ax <- length(res$values) #解の個数
        res$vectors <- t(t(res$vectors[,1:ax])*sqrt(res$values)) # 解(軸)に対してサンプルの座標を与える
        colnames(res$vectors) <- names(res$values) <- paste("Axis", 1:ax) # 軸に名前を与える
        rownames(res$vectors) <- paste("Object", 1:n) # サンプルに名前を与える
        write.table(res$vectors,outfile,sep="\t") #出力して、
        table<-read.table(outfile,sep="\t",header=TRUE) # そのファイルを読み込んでいる
        plot(table[1:noPlotAxis],...) # 第１軸から、指定軸数のplotをさせる
        return(table) #tableの形で、返り値を返すので、値参照、プロットのしなおしはそれを用いて適宜行う
}

• • 使い方(関数内での表示プロットと、返り値 "answer"から軸を指定しなおしてプロットした図とを掲載)

>distance<-scan("file.txt")
>n<-sqrt(length(distance)) # 行数・列数
>matdistance<-matrix(distance,n,n) #距離行列
>s<-(-1)*matdistance #類似度行列への変換
>thres <- 0 # 有効な固有値の下限閾値
>noPlotAxis<-4
>outfile<-"hoge.txt"
>answer<-princoRY(s,thres,noPlotAxis,outfile)
> plot(answer[8:10])

• 返り値を増やして

princoRY2 <- function(s,thres=0,noPlotAxis=2,outfile,...)#類似度行列、有効固有値下限、プロットの軸数、出力ファイル名、plot関数の引数
{
        n <- nrow(s)# 行列のサイズ
        m <- colSums(s)/n #列別の平均。対称行列なので、行別の平均でもある
        h <- s+sum(s)/n^2-outer(m, m, "+") 
         #２重中心化
         # 第i,j要素から、第i行,第j列の平均を引き、(引きすぎた分である)、全体平均を足し戻している
         # 第２項 sum(s)/n^2は、全体の平均
         # 第３項 outer(m,m,"+") は、第i行平均と第j列平均の和
        res <- eigen(h) # 固有値・固有ベクトルを求める
        res$values <- res$values[res$values > thres]  # 固有値のうち、与えた値を下限として、解の個数を決める
        ax <- length(res$values) #解の個数
        res$vectors <- t(t(res$vectors[,1:ax])*sqrt(res$values)) # 解(軸)に対してサンプルの座標を与える
        colnames(res$vectors) <- names(res$values) <- paste("Axis", 1:ax) # 軸に名前を与える
        rownames(res$vectors) <- paste("Object", 1:n) # サンプルに名前を与える
        write.table(res$vectors,outfile,sep="\t") #出力して、
        table<-read.table(outfile,sep="\t",header=TRUE) # そのファイルを読み込んでいる
        plot(table[1:noPlotAxis],...) # 第１軸から、指定軸数のplotをさせる
        ret<-list(NULL)
        ret[[1]]<-res
        ret[[2]]<-table
        return(ret) #tableの形で、返り値を返すので、値参照、プロットのしなおしはそれを用いて適宜行う
}