Metric 予測 by 複数 Metric：複数 Predictors の実装例

「冗長 Predictors」からの続き。

ここでは、複数 predictors による予測モデルの実装例を、本書提供のスクリプトもとに見る。

以下は、本章の SAT スコアに関するデータ分析モデル（図 18.4）。

左の図17.2「単純直線回帰」モデルとの違いは、複数 predictors のみ（Σ_j β_j x_j,i の部分）。

以下は、本モデルが使用するデータ（Jags-Ymet-XmetMulti-Mrobust-Example.R 実行後に参照可能）。

> head(myData,3)
    State Spend StuTeaRat Salary PrcntTake SATV SATM SATT
1 Alabama 4.405      17.2 31.144         8  491  538 1029
2  Alaska 8.963      17.6 47.951        47  445  489  934
3 Arizona 4.778      19.3 32.175        27  448  496  944

predictors を Spend と PrcntTake にする場合、以下のように指定。

yName = "SATT" ; xName = c("Spend","PrcntTake")

その後、以下のように Jags-Ymet-XmetMulti-Mrobust-Example.R から Jags-Ymet-XmetMulti-Mrobust.R に定義されている関数 getMCMC を呼び出す。

mcmcCoda = genMCMC( data=myData, xName=xName, yName=yName, 

    numSavedSteps=numSavedSteps, thinSteps=thinSteps,

    saveName=fileNameRoot )

以下は Jags-Ymet-XmetMulti-Mrobust.R に定義されている data ブロック。

  y = data[,yName]

  x = as.matrix(data[,xName],ncol=length(xName))

  data {

    ym <- mean(y)

    ysd <- sd(y)

    for ( i in 1:Ntotal ) { # Ntotal is the number of data rows

      zy[i] <- ( y[i] - ym ) / ysd

    }

    for ( j in 1:Nx ) {     # Nx is the number of x predictors

      xm[j]  <- mean(x[,j]) # x is a matrix, each column a predictor

      xsd[j] <-   sd(x[,j])

      for ( i in 1:Ntotal ) {

        zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]

      }

    }

  }

matrix の x, zx は、それぞれ「標準化前」「標準化後」のデータ。x は以下のようにして再現できる。

> head(as.matrix(myData[,c("Spend","PrcntTake")],ncol=length(c("Spend","PrcntTake"))))

     Spend PrcntTake

[1,] 4.405         8

[2,] 8.963        47

[3,] 4.778        27

[4,] 4.459         6

[5,] 4.992        45

[6,] 5.443        29

以下は model ブロック。

  model {

    for ( i in 1:Ntotal ) {

      zy[i] ~ dt( zbeta0 + sum( zbeta[1:Nx] * zx[i,1:Nx] ) , 1/zsigma^2 , nu )

    }

    # Priors vague on standardized scale:

    zbeta0 ~ dnorm( 0 , 1/2^2 )  

    for ( j in 1:Nx ) {

      zbeta[j] ~ dnorm( 0 , 1/2^2 )

    }

    zsigma ~ dunif( 1.0E-5 , 1.0E+1 )

    nu ~ dexp(1/30.0)

    # Transform to original scale:

    beta[1:Nx] <- ( zbeta[1:Nx] / xsd[1:Nx] )*ysd

    beta0 <- zbeta0*ysd  + ym - sum( zbeta[1:Nx] * xm[1:Nx] / xsd[1:Nx] )*ysd

    sigma <- zsigma*ysd

  }

主に、脚色した点（複数 predictors 処理）を除けば、predictors が 1 つの場合と同じ。

以下がスクリプト実行結果（図 18.5 に相当）。

これは、前回の投稿の「predictors 冗長なし」に相当する。

冗長 predictor の PropNotTake の追加は、以下のようにする。

 myData = read.csv( file="Guber1999data.csv" )
 PropNotTake = (100-myData[,"PrcntTake"])/100
 myData = cbind( myData , PropNotTake )
 yName = "SATT" ; xName = c("Spend","PrcntTake","PropNotTake")

「掛算の効果」に続く。