トップ > サポート > ご利用での Q&A > データハンドリング、解析技術編

データハンドリング、解析技術編

は Windows 版固有、 は UNIX 版固有、 は Linux 版固有、 は S-PLUS 2000 版固有の問題です。

質問 一覧

S-PLUSでの数値の有効桁数はどうなっているのか。表示桁を変えられるか。

S-PLUSは内部でデータを倍精度（8byte）で保持していますので、比較的高い精度で計算ができます。倍精度浮動小数点数(double)は、一般に符号部1ビット、指数部11ビット、仮数部52ビットで構成される場合が多く、-10306〜10306の範囲の数を表現可能で、計算精度は15桁程度です。（CPUによって異なる場合があります）

表示桁数を変更するには、関数optionsを用いてdigitsを変更します。デフォルトでは7桁表示です。

Example :

> pi
[1] 3.141593
> options(digits=20)
> pi
[1] 3.14159265358979311

量的なデータを、"H","M","L"のようなカテゴリカルデータに変換するには。

S-PLUSではfactorというデータクラスがあるので、関数cutとfactorを用いてこれに変換します。任意の水準で分類したり、各水準に好きな名称をつけることも可能です。いったんfactorクラスに変更された後は、関数tableを用いて単純集計することや、crosstabsでクロス集計と検定を行うことも可能です。

> x <- factor(cut(lottery.payoff,3),label=c("L","M","H"))
# 宝くじの賞金を高、中、低の3つに分類する
> table(x)
L M H
177 72 5

関数interpはどのように不規則な点を格子状の点に補間するのか。

オプション引数ncp=0のとき（デフォルト時）は隣接する3角形による線形補間を行います。この方法は以下のようになります。

N 個のデータ (x_1, y_1, z_1), ... (x_N, y_N, z_N) があるとします。この時，座標 (x, y) の（線形）補間値 z を求めたいとします。

  1. N 個のデータの x-y 座標から作れるすべての三角形を考えます。このうち，座標 (x, y) を内点として持つ，最小の三角形を構成している 3 つのデータを取り出します。
  2. 次に，この 3 点を通る平面の式を求めます（平面の方程式は、z = ax + by + cと書けますので，3 点から一つの平面が一意的に決まります）。
  3. 2. の平面の式に座標 (x, y) を代入して補間値 z を求めます。

これを求めたい (x, y) について繰り返し行ないます。ただし，1. で (x, y) を内点に持つ三角形がない場合には，その点の補間値は，NAとなります。

また、2≦ncp<観測値の数で、ncpが指定されているときは、近傍のncp個の点を使って求めた偏微分による3次補間を行います。またこの時、extrap=TRUEとすれば、補外を行うことが出来ます。

（参考）三角形の内点であるかどうかを決定するアルゴリズムは以下のようになります。

xy-平面が、データの (x, y) 座標を頂点とするいくつかの三角形で分割できたとします。ただし、三角形の頂点は、反時計回りの順に保管しているものとします。

例えば、(x1, y1),(x2, y2),(x3, y3) を三角形の頂点の座標（この順で反時計まわりに並んでいるものとします）。補間を行ないたい点の座標 (xi, yi) がこの三角形の内点であるかどうかは、以下のアルゴリズムで決定します。

  1. xi が c(x1, x2, x3) の最小値、最大値の間にある。
  2. 同様に、yi が c(y1, y2, y3) の最小値、最大値の間にある。
  3. f(u1,v1,u2,v2,u3,v3) = (u1 - u3) * (v2 - v3) - (v1 - v3) * (u2 - u3)

とし、

f(x1, y1, x2, y2, xi, yi) > 0
f(x2, y2, x3, y3, xi, yi) > 0
f(x3, y3, z1, y1, zi, yi) > 0

の 3 つの式が全て成り立つ。

1、2、3 がすべて成り立つ場合、(xi, yi) は、この三角形の内点であることがわかります（3.は、三角形の頂点を反時計回りに並べたことによります）。また、3. で等号が成り立つ時には、(xi,yi) が辺上にあることになります。

どうすればS-PLUSで数量化理論1類ができるか。

S-PLUS における数量化理論は既存の多変量解析のための関数により、実現可能です。この関数では質的データ（因子、カテゴリカルデータ）を量的データと同様に扱うことができます。ダミー変数を作成する必要はありません。クロス集計の必要もありません。

例えば次のようなデータ suryo があったとします。

年収    性別    年代    職業     最終学歴
300    男    20代    会社員     高校
...

目的変数（外的規準）を年収、説明変数（アイテム）を性別、年代、職業、最終学歴とします。説明変数のいずれもが質的な変数です。数量化1類の場合、線形モデルをあてはめる関数 lm を用いて次のように指定します。

>options(contrasts=c("contr.treatment","contr.poly")) # 注意!
>ichirui <- lm(年収-mean(suryo$年収) ‾ -1 + 性別 + 年代 + 職業 + 最終学歴, data = suryo)

注意! S-PLUS では質的な変数を数値に変換する規則がいくつかあります。デフォルトは Helmert 対比（Helmert contrast、「S と統計モデル」P34をご覧ください）という規則を使います。これを数量化理論の参考書に出てくるような 1,0 のダミー変数に変換するためにはこの options 関数を使います。

記号 ‾ の左辺が目的変数、右辺が説明変数です。数量化のため、左辺は年収の平均を引いたもの、右辺は説明変数を + でつないだものとします。右辺のマイナスは左辺で年収の平均を引いているため、定数項のないモデルにするためです。

モデルを実行したところで、ダミー変数の係数（カテゴリースコア）を参照してみます。

>coef(ichirui)

性別女     性別男     年代          職業会社員      職業教師     最終学歴
-51.18687 -45.05051 66.13636     -7.727273      -17.95455     -15

このような結果が出力されたとします。例えば、「職業会社員」という名前の係数は、名前から「職業が会社員のとき1、そうでないとき0」のダミー変数の係数であることがわかりますが、2水準、つまり年代のダミー変数は20代が 1 か、30代が 1 となっているかわかりません。

そこで次にこれをチェックします。

> ichirui$contrasts

$"性別":
[1] "性別女" "性別男"

$"年代":
        30代
20代     0
30代     1

$"職業":
         会社員    教師
医師     0          0
会社員  1          0
教師     0            1

$"最終学歴":   
        高校
高校     1
大学     0


これでこのダミー変数「年代」が30代のとき1、20代のとき0であることがわかります。
（定数項を除いたモデル式を指定しましたので、最初のカテゴリー変数はすべてのアイテムについてダミーが作られます）

これは「基準化されていない」カテゴリースコアです。必要に応じて、以下の例を参考に基準化してください。

性別を基準化

>n <- nrow(suryo)
>coef(ichirui)[1:2]-sum(coef(ichirui)[1:2]*table(suryo$性別))/n

年代を規準化

>c(0,coef(ichirui)[3])-sum(c(0,coef(ichirui)[3])*table(suryo$年代))/n

1番目のカテゴリーはすべてのアイテムについてダミーが作られていましたが、2番め以降は 1 つのアイテムが 0 ですので、これを明示的に足してください。

***プログラムに関して***

ちょっとプログラムに関する解説を付け加えます。

>c(0,coef(ichirui)[3])    # 最初の水準が 0 に設定されているので付け加えたベクトルにした。
>table(suryo$年代)     # 年代ごとの発生件数を計算。ここで最初に表示される水準が 0 に設定される。

つまり

sum(c(0,coef(ichirui)[3])*table(suryo$年代))/n

の部分で、カテゴリースコアの平均を求めているわけです。

高速フーリエ変換（FFT）やスペクトル解析はどうやってやるのか。

fft やスペクトル解析はコマンド上での操作が必要になります。

fft は関数 fft を使って行います。利用するデータはベクトルまたは配列です。配列の場合は、多次元離散フーリエ変換を行います。詳細は help(fft)もしくは「S言語2」の「fft」の項をご覧ください。

スペクトル解析は関数spectrumなどを利用して行います。デフォルトではペリオドグラムの描画を行います。詳しくは help(spectrum) またはGuide to statisticsの「Time Series Analysis」の項をご覧ください。

ペリオドグラムのスムージング方法

Example :

> fft(lynx) #lynx時系列データをfftにかける
> spectrum(lynx)
> fftdata <- scan("A:¥¥fft.dat") #テキストデータファイルfft.datから読みこむ
> fft(fftdata)
> spectrum(fftdata)

(図： spectrum(lynx)の結果表示、縦軸はdb）

欠損値の入ったデータを解析するには。

S-PLUSでは、欠損値を表す NA というデータ形式が用意されています（Not Available）。

一般的に、欠損値が入ったデータをそのまま解析の入力とするとエラーになります。欠損値を含むペアを除去するなどの前処理が必要です。いくつかの関数によって、この処理方法が異なります。

  1. meanなどの関数は、引数 na.rm をTRUEにすることにより、欠損値を取り除いた上で処理します。（デフォルトはFALSEですので、結果はNAになります。）
  2. varなどの関数は、引数 na.method を"omit"にすることにより、欠損値を取り除いた上で処理します。（デフォルトは "fail" ですのでエラーになります。）
  3. lmやfactanalといった関数では、引数 na.action をna.omitにすることにより、欠損値を含むペアを取り除いた上で処理します。（デフォルトは na.fail ですのでエラーになります。）

このように、関数によって若干引数が異なります。詳しくはオンラインヘルプをご覧ください。

また、欠損値除去のオプションがない関数もあります。この場合は、関数 is.na を用いて欠損値を含むペアを除去してください。

Example :

> x <- c(1,NA,2) #NAを含むデータを作成
> x
[1] 1 NA 2
> x[x!=NA] #これでは除去できない（NAとの演算は論理演算を含め、すべて結果がNAになるため）
[1] NA NA NA
> x[!is.na(x)] #関数 is.na を使って除去する。
[2] 1 2
> x <- matrix(c(1,2,3,4,NA,5.5),ncol=2) #NAを含むデータを作成
> x <- as.data.frame(x) #データフレームに変換
> lm(x)
Error in function(object, ...): missing values not allowed: found in V2
Dumped
> lm(x, na.action=na.omit)
（出力は略）
> x[!is.na(x[,2]),] #2列目にNAを含むペアの除去
V1 V2
1 1 4.0
3 3 5.5

モデル式（Formula）とは何か。

S言語ではモデル（仮説）をソフトウエアで扱いやすい形で、しかも柔軟でなおかつ理解しやすい形に表現するために、「モデル式（Formula）」を利用します。

例えば車の燃費（Fuel）と、車の重量（Weight）と排気量（Disp）との関係を調べようとしたとき、

Fuel ‾ Weight + Disp

というモデル式は、「燃費は車の重さとエンジンの排気量で決まる」というモデル（仮説）を表します。

このモデル式を、適当な解析関数に当てはめることにより実際の処理を行います。例えば、

lm(Fuel ‾ Weight + Disp)

では、線形モデルへのあてはめ、すなわち重回帰分析を行い、

tree(Fuel ‾ Weight + Disp)

では、ツリーモデルへのあてはめを行います。

また、因子分析のように外的基準となる反応変量が存在しない場合やクロス集計を行う場合は、左辺値が空になります。

モデル式の各項は任意のS式でもかまいません。このため非常に柔軟なモデリングが可能になります。例えば

1/Fuel ‾ sqrt(Weight) + Disp

といった記述も可能になります。

この他、交互作用などの表現も簡潔にできます。詳しくは「Sと統計モデル」の第2章をご覧ください。

Example :

> attach(fuel.frame) #データフレームfuel.frameの列名（変数名）だけで参照可能にする
> lm(Fuel ‾ Weight + Disp.)
（出力略）
> attach(market.survey)
> crosstabs( ‾ pick + nonpub, na.action=na.omit) #クロス集計とカイ二乗検定

必ず原点(0,0)を通過する回帰を行いたい。

関数lmを使った回帰を行うとき、必ず原点を通る（切片が0になる）回帰を行いたい場合は、モデル式に明示的に -1 を追加します。

Note!

Windows版で、（回帰式は不要で）単に描画のみ行いたい場合は、回帰線をダブルクリックしてCurve Fitting plot のダイアログを表示し、Curve Fittingのタブで「Omit Constant」をチェックすることにより原点を通る直線となります。

Example :

> attach(air)
> lm.out1 <- lm(ozone ‾ radiation) #通常の回帰
> lm.out2 <- lm(ozone ‾ -1 + radiation) #モデル式に「-1」を追加
> plot(radiation, ozone)
> abline(lm.out1)
> abline(lm.out2,col=2)

カテゴリカルデータからいわゆる「ダミー変数行列」を作成するには。

ダミー変数行列を作成するには、関数 model.matrix を利用します。デフォルトではヘルマート対比でダミー変数化します。（通常良く使われる）0と1でのダミー変数行列に直すには、あらかじめoptions(contrasts=c("contr.treatment","contr.poly")と設定しておく必要があります。

また、デフォルトでは水準を1つ取り除いてランク落ちを防いだ適切な形でダミー変数行列化を行いますが、model.matrixの入力とするFormulaに-1を加えると全部の変数をダミー変数に変換します。詳細は「Sと統計モデル」の32ページ以降を参照してください。

Note!

S-PLUSでは一般に、こうした処理を内部で自動的に行うので、データが適切なクラスで作成されていれば、わざわざダミー化する必要はあまりありません。

Example :

> x <- factor(cut(runif(10)*3,breaks=c(0,1,2,3)),label=c("A","B","C"))
# 乱数から、3水準のカテゴリカルデータを作成する
> model.matrix(‾x)

(Intercept) x1 x2
1 1 -1 -1
2 1 1 -1
3 1 0 2
（略）

# Helmert対比で出力
> options(contrasts=c("contr.treatment","contr.poly"))
> model.matrix(‾x)

(Intercept) xB xC
1 1 0 0
2 1 1 0
3 1 0 1
（略）

# Aという水準を一つ落としてダミー変数化
> model.matrix(‾-1+x)

xA xB xC
1 1 0 0
2 0 1 0
3 0 0 1
（略）

# 全ての水準でダミー変数化

金融ポートフォリオの最適化問題は解けるか。

以下のような問題を考えてみます。n 個の銘柄の株式があり、それぞれの組み入れ比率を x1, x2, ..., xn としたときのポートフォリオのリスクは、リスク指標として分散を用いたとき、

で、これを最小にすることを考えますが、そのとき、空売りを認めない場合、制約式

を同時に満たす必要があります。

最適化問題を解くための関数として、S-PLUSには nlmin、nlminb が用意されていますが、ここで問題になるのが制約式の存在です。単純な xi >=0 という制約式は組み込めますが、これより複雑な制約式は S-PLUS 単独ではサポートされていません。

しかしながら、NTTデータ数理システムで開発しております、数理計画パッケージ（NUOPT）により、このような典型的なリスク - リターンの最適化問題はもちろん、さらに複雑な制約式の考慮、大規模問題の解法が可能です。

S-PLUS 上からこの NUOPT をご利用いただけるアドオンソフト SNUOPT もございますので、詳細に付きましては、お気軽にNTTデータ数理システムまでお問い合わせください。

スペクトル解析でペリオドグラムのスムージングはどうするのか。

spectrumなどの関数を用いる時に、オプション span を利用するとスムージングを行います。デフォルトでは span = 1 で、スムージングしていないペリオドグラムを描画します。

詳しくは help(spec.pgram) でリファレンスを参照するか、Guide to statisticsの「Time Series Analysis」の項をご覧ください。

Example :

> par(mfrow=c(2,1))
> spec.pgram(sunspots, span=1, plot=T)
> spectrum(sunspots, spans=c(19,19,19)) #スムージング

主成分分析を行う関数でprcompとprincompの2つの関数があるが、どちらを使えば良いか。

princompがお奨めです。

prcompはベル研究所の「S言語」に元々付属していた関数ですが、princompはS-PLUSとして新しく追加された関数であり、機能も増えています。欠損値処理や、モデル式での記述も可能です。

どうすれば S-PLUS で数量化理論2類ができるか。

S-PLUS における数量化2類は既存の行列演算、統計解析の関数により、実現可能です。ダミー変数を作成する必要はありません。

例えば次のようなアンケート結果のデータ nirui があったとします。

classic SF 群
1 好き 好き 1
2 好き 好き 1
3 好き どちらでもない 1
4 嫌い どちらでもない 1
5 嫌い どちらでもない 2
6 嫌い 嫌い 2
7 好き どちらでもない 2
8 好き 嫌い 3
9 好き 嫌い 3

各行（観測値）はこの質問に答えた人の回答を、列（変数）classic は「クラシック音楽が好きか嫌いか」を、SF は「SF 映画が好きか嫌いか」を、群は「それぞれの人が所属するなんらかの群」をあらわします。

このデータで、

クラシック音楽が好き a1
嫌い a2
SF 映画が好き b1
どちらでもない b2
嫌い b3

という点数を与えて、その情報だけで最大限に 1 群と 2 群を区別できるような a1,a2,b1,b2,b3 を求めることが数量化 2 類の目的です。

この点数を使って、上のデータをあらわすと

classic SF 群
1 a1 b1 1
2 a1 b1 1
3 a1 b2 1
4 a2 b2 1
5 a2 b2 2
6 a2 b3 2
7 a1 b2 2
8 a1 b3 3
9 a1 b3 3

という得点になります。それぞれの得点を z11,z12,...z39 というように z 群（群の中での順番）と表すと

z11 = a1 + b1
z12 = a1 + b1
z13 = a1 + b2
z14 = a2 + b2
z25 = a2 + b2
...
z39 = a1 + b3

となります。全データの点数の平均を Z ダッシュ、1 群の点数の平均を Z1 ダッシュ、2 群の点数の平均を Z2 ダッシュ、3 群の点数の平均を Z3 ダッシュとします。

全データの点数の変動を考えます。

ここで、1 の部分は

となります。同様に、z2,z3 についても同じことがいえますから、上の式は

と表すことができます。これを文章で表すと

全体の変動 = 群内変動 + 群間変動

と表されます。数量化 2 類の目標は（群間変動 / 群内変動）を大きくすることです。

さらに、これらの式を一般化します。

全体の変動

群間の変動

と表されます。N は全データ個数、z は総和を表します。

さて、いよいよ a1,a2,b1,b2,b3 を求めます。ここで、それぞれの変数（クラシック音楽が好き、SF 映画がすき）のうち、ひとつの係数を 0 としないと、解法できないので、ここでは a1,b1 を 0 とします。

すると得点 z は、

a2 b2 b3
0 0 0
0 0 0
0 1 0
...

という行列 D と [a2,b2,b3] という行列 a の行列積となります。次にわかりやすくするために、次のような計算をしておきます。

★計算A

Z0 として、[3,4,3]（各列の合計）b と a の行列積

Z1 として、[1,2,0]（群1での列合計）c と a の行列積

Z2 として、[2,2,1]（群2での列合計）d  と a の行列積

Z3 として、[0,0,2]（群3での列合計）e と a の行列積

すると、

となります。これらと先ほど計算した式から

全変動

群間変動

n1,n2,n3はそれぞれ、1群、2群、3群のデータ個数を表します。

目的は（群間変動 / 全変動）を最大にする a2,b2,b3 を求めることでした（a1,b1 はすでに0としています）。a2,b2,b3の取り方により、全変動を1とすることができますので、この条件の下で、群間変動を最大にします。

 ラグランジュの定理より

で、固有ベクトルがa2,b2,b3の値になります（a1=0,b1=0と仮定していました）

さて、いよいよ S-PLUS 上でコーディングしてみます。nirui というデータがすでにExcelファイルなどからインポートされているとします。まずは、行列Dを作成します。通常のプログラミング言語では、このような場合、（if classic="好き"  then classsic1 = 1 else classic = 0）などと if 文を書くと思われますが、S-PLUS には便利な機能があります。

　

>options(contrasts=c("contr.treatment","contr.poly"))…コマンド1

>matD <- model.matrix( ‾ classic + SF, data=nirui)

>matD

(Intercept)

classic

SF嫌い

SF好き

1

1

1

0

1

2

1

1

0

1

3

1

1

0

0

4

1

0

0

0

5

1

0

0

0

6

1

0

1

0

7

1

1

0

0

8

1

1

1

0

9

1

1

1

0

...

あらかじめ、コマンド1が必要になります（理由は数量化理論1類のページをご覧ください）。これで、a1,a2,b1,b2,b3 のうち、「クラシックが嫌い」と「SFはどちらでもない」を落とした行列Dをつくる準備ができました（S-PLUS では1文字の名前を付けた関数がかなりありますので、マスクを避けるため、ここでは行列の名前の頭にmat を付けるようにします）。

ところがこの行列matDには問題があります。それは定数項（Intercept）が付いていることです。これをむりやりとってしまいましょう。次のような簡単な演算です。

>matD <- matD[,-1]

続いて、A,B,C の計算をするために、いくつか必要な数字を求めておきます。

>gun <- nirui$群      # 群を表すインデックス

>vecN <- nrow(nirui) # データ件数を求める

>vecn <- table(gun)  # 各群ごとのデータ件数(n1,n2,n3)を求める

　

以下の()の中のアルファベットは計算 A におけるベクトルの名前を表します。

　

>vecb <- apply(matD,2,sum)  # 各列ごとに1の数を数える（b）

>vecc <- apply(matD[gun==1,],2,sum)

# 群が1のもので各列ごとに1の数を数える（c）

>vecd <- apply(matD[gun==2,],2,sum)

# 群が2のもので各列ごとに1の数を数える（d）

>vece <- apply(matD[gun==3,],2,sum)

# 群が3のもので各列ごとに1の数を数える（e）

>matT <- t(matD)%*%matD-vecb%*%t(vecb)/vecN

# B を計算

>matB <- vecc%*%t(vecc)/vecn[1]+vecd%*%t(vecd)/vecn[2]+vece%*%t(vece)/vecn[3]

-vecb%*%t(vecb)/vecN

                    # C を計算

>i <- eigen(solve(matT)%*%matB)

                    # D を計算

関数eigenは固有値を解くための関数です。ここで求められたiの要素のうち、vectors が固有ベクトルになります（いわゆるカテゴリースコアといわれているものです）。

さらに個々のデータの評点を求めます。

> matD%*%i$vectors[,1]

1 0.01821303

2 0.01821303

3 0.36675507

4 0.00000000

5 0.00000000

6 0.71529712

7 0.36675507

8 1.08205219

9 1.08205219

...

とそれぞれの評点が求まりました。

S-PLUS でニューラルネットワークは実行可能か。

Oxford 大学の Brian Ripley 教授、W.N. Venables 教授によるライブラリ nnet の利用により、可能です。参考書籍は

Modern Applied Statistics winth S-PLUS 
Springer
ISBN 0-387-98214-0

です。

次のコマンドにより、ライブラリ nnet を読み込み、関数 nnet が実行できます。

> library(nnet)
> result <- nnet(Kyphosis ‾ Age, data=kyphosis, size=2)

詳細は、help(nnet) で表示されるヘルプを参照してください。

ステップワイズ変数選択をしたいのだが。

S-PLUSにおけるステップワイズ変数選択の方法には2通りの方法が考えられます。

一般的に後者を指して、「ステップワイズ法」と呼ことが多いようです。それでは両者の方法を比較してみましょう。

段階的モデル選択

まず、通常の回帰分析を行います。

> stepx1 <- lm(skips ‾ ., data=solder)

続いて、段階的モデル選択を行いますが、変数の2次の交互作用項を加えるか加えないかを考えてみます。

> stepxx <- step(stepx1, ‾.^2)

このとき、step は最初に add1(stepx1, ‾ .^2) という関数を動かし、その結果、加えると一番 AIC 基準が減少するような項を探します。その一方、drop1(stepx1) として、落とすと一番 AIC が減少するような項を探します。この繰り返しステップをあらかじめ定めた限界（上ならば、2次の交互作用をすべて加える、あるいは定数項のみのモデルにする）に達するか、もうそれ以上 AIC 基準を減少できなくなるまで、繰り返します。

つまり、変数選択の基準には、AIC が使われているわけです。

RSS の値によるステップワイズ法（関数 stepwise による方法）

はじめから、回帰のかわりにステップワイズを行います。

> stepyy <- stepwise(evap.x, evap.y)

こちらは、RSS を基準にして、モデルを入れ換え、変数選択します。

S-PLUS で多項ロジスティック回帰は実行可能か。

Oxford 大学の Brian Ripley 教授、W.N. Venables 教授によるライブラリ nnet の利用により、可能です。参考書籍は

Modern Applied Statistics winth S-PLUS 
Springer
ISBN 0-387-98214-0

です。

関数multinom でご利用可能です。

> library(nnet)
> help(multinom)

で表示されるヘルプを参照してください。

ダービン=ワトソン比を出したい。

関数としてすでに組み込まれてはいませんが、次のように S 言語で簡単に算出可能です。

なんらかの回帰（標準的な回帰分析でも、AR モデルのようなものでも同じです）を行ったとします。

a <- 回帰(lm,ar など）

d.w <- sum(diff(a$resid)^2,na.rm=T)/sum(a$resid^2,na.rm=T)

で OK です。GUI で回帰などを行った場合はダイアログの Save As に入力された名前を上の a の代わりに用いて計算してください（Windows 版 S-PLUS 4.5 まではデフォルトで last.lm のような名前がついていますが、S-PLUS 2000 ではデフォルトの名前がついていません（結果は保存されません）。Save As に必ず名前を入力してください）。

ARIMA モデルを行うため arima.mle を実行すると、Cannot coerce mode list to double というエラーがでる。

分析のためのデータは Excel などからインポートされたものではありませんか？
これらは見かけ上、ベクトルのように見えますが、実際はデータ・フレームというオブジェクト・クラスで、ベクトルではありません。arima.mle など、いくつかの時系列のための関数は、ベクトル、あるいは時系列ベクトルをその引数とします。以下のように、1列だけを抜き出すか、ベクトルに型変換してください。

arima.mle(data[,1],....)

    # ベクトルとして1列だけ抽出

vec <- as.vector(unlist(data))

    # ベクトルに型変換、unlist(data) だけでも、arima.mle などは実行可能

arima.mle(vec,...)

制約つき非線形最小二乗法を行いたい。

nlregb という関数を使います。
例えば、

> set.seed(20)
　　　　# 乱数のシードを設定（解析の再現が可能です）
> x <- seq(0, 2*pi, length=100)
> noise <- rnorm(100, sd=0.5)
> y <- 4 * sin(x - pi) + noise

y=4 * sin(x - 3.14)  という式にノイズを与えました。
この、x と y から、パラメータ 4 と pi を求めるには

> my.new.fcn <- function(param, x, y)
> {
>    amp <- param[1]
>    horshft <- param[2]
>    y - amp * sin(x - horshft)
> }

という関数を定義します。最後の行が最小にしたい式になるわけです。amp = 4, horshft = 3.14 と推定されればOK です。

では、最小二乗法に移行します。このときに、amp が非負という制約条件を入れてみましょう。

> nlregb(n = 100, start = c(3, pi/2),
>     residuals = my.new.fcn,
>     lower = c(0, -Inf), x = x, y = y)

結果のうち、

$parameters:
[1] 4.112227 3.152186

が求めるパラメータとなります。

ただ、この場合の制約ですが、「非負」や「1 < param < 5」みたいなごくごく単純なものしか与えられません。ご了承ください。
さらに複雑な制約式が必要な場合は、弊社開発製品 NUOPT をご利用ください。

S-PLUS 上からこの NUOPT をご利用いただけるアドオンソフト SNUOPT もございますので、詳細に付きましては、お気軽にNTTデータ数理システムまでお問い合わせください。

コレスポンデンス分析は可能か？

ライブラリ MASS の関数 corresp により、可能です。

> library(mass)

とすることにより、関数 corresp の利用が可能になります。
コレスポンデンス分析の場合、クロス集計表を用いますが、この関数 corresp の場合、元のデータ（クロス集計する前の）でも、クロス集計済みのデータでも、用いることができます。

元のデータの場合
> corresp( ‾ gender + age, data=data)

と、 チルダの左側をブランク、右側に分析したい変数を + でつなぎます。S-PLUS が関数内部でテーブル処理を行います。

テーブルの場合

> corresp(data.table)

のように、直接テーブル処理したデータを当てはめてください。S-PLUS 上でテーブルを作る場合、関数 table を用います。上のデータの場合は、

> data.table <- table(data$gender, data$age)

とします。

どうすればS-PLUSで数量化理論3類ができるか？

統計学辞典（東洋経済）によれば、数量化3類とコレスポンデンス分析が同等ということです。「コレスポンデンス分析は可能か」のページをご覧ください。

他の統計ソフトで計算させた判別分析の結果と、S-PLUS で計算した結果が異なるようだが？

例えば、判別分析の目的変数が2群の因子（「はい」「いいえ」のどちらか）だったとすると、通常、よく使われる判別分析では、判別式が1つだけ求まります。これは、2つの群の境界を表す式です。

一方、S-PLUS で同じデータを判別分析すると、結果の判別式は2式になります。

これは、それぞれが2群の各群に対する、判別スコアになるからです。
この判別スコア式をそれぞれ、f1, f2 とすると、f1-f2=0 が境界線になり、実際、f1-f2 あるいは、f2-f1 を計算すると、よく使われる判別分析の判別境界式になります。
多群でも同様です（例えば、3群ならば、f1-f2, f1-f3, f2-f3 の3式）。

例えば、判別分析の結果を last.discrim とします（メニューから「統計 / 多変量解析 / 判別分析」を選んだときのダイアログ「Discriminant Analysis」で「Save As」を「last.discrim」とします）。

境界式は

定数項：
coef(last.discrim)[[1]][1]-coef(last.discrim)[[1]][2]
係数：
coef(last.discrim)[[2]][,1]-coef(last.discrim)[[2]][,2]

で、計算ができます。

どうすれば S-PLUS で共分散分析(Analysis of Covariance）ができるか？

共分散分析（Analysis of covariance) を特定するには、通常の 回帰の関数 lm に「入れ子」を指定するモデルを指定するだけです。具体的には

lm(post ‾ test/pre,...)

などとなります。/ の左側が効果となる因子です。

Windows の GUI から主成分分析を行ったが、そのときの主成分得点をデータとして保存したい。

主成分分析をGUI から実行する際に、下の「Save As」に適当な名前を入力、OK をクリックします。

主成分分析GUI

続いて、コマンド・ウィンドウから、次のように入力します（仮に、上の Save As で last.princomp と入力したとします）

> pr.scores <- as.data.frame(last.princomp$scores)

> names(pr.socres) <- paste("prin", 1:ncol(pr.scores), sep="")

すると、prin1, prin2, ..., という列名を持つデータが作成されます。このデータを Excel にエクスポートしたりすることが可能です。

ARIMA モデルを作成し、予測を行った。arima.forecast のヘルプの EXAMPLES に出ている通り、予測のグラフを作成したが、思ったとおりのグラフにならない。

ARIMA モデルのための関数 arima.mle では、解析できるデータが、ベクトルか時系列（ts など）です。EXAMPLES に出ている ship データは時系列（ts）です。このデータをarima.mle に当てはめ、arima.forecast で予測を行うと、日付が正しく予測（未来）のものになります。そのため、ts.plot を使うことにより、ship データに続けて、予測されたデータが表示される結果になります。

arima.mle でベクトル（時系列ではない）データを扱うと、正しく解析されますが、元のデータに時系列の属性がないため、ts.plot で予測のデータが、目盛り１から表示されてしまいます。

あらかじめ、関数 ts などで、ベクトルを時系列に直して、解析を行ってください。

コックスの比例ハザードモデルで、ベースライン関数を求めるには、また、新たなデータから、そのデータのハザード関数を求めるには？

コックスの比例ハザードモデル（Cox Proportional Hazard）を行う関数は coxph です。ここでは、サンプルデータ ovarian を取り上げます。

fit <- coxph(formula = Surv(futime, fustat) ‾ age, data = ovarian)

futime が時間、fustat が打ち切りなどのイベントです。

ここから、ベースライン関数H0(t)を求めるには

h0t <- survfit(fit)

とします。作図するには

plot(h0t)

値を得るには

summary(h0t)

とします。新しいデータを用いて、そのデータによるハザードH(t) を求めるには

summary(survfit(fit, newdata=data.frame(age=34)))

とします。上ならば、34歳の人のハザードが求まります。

因子列を含むデータから抽出を行ったが、この因子列の因子水準にもう存在しない値が含まれてしまう。これを取り除くには？

たとえば、サンプルのデータ・フレーム fuel.frame は因子列 Type を含みます。この Type が "Compact" のもの以外で新しいデーを作成します。

> fuel2 <- fuel.frame[fuel.frame$Type!="Compact",]

すると、新しいデータ fuel2 の列 Type もやはり因子ですが、因子水準を見ると、fuel.frame のものが継承されています。つまり、もはや存在しない水準 Compact が含まれています。

> levels(fuel2$Type)

[1] "Compact" "Large" "Medium" "Small" "Sporty" "Van" 

この Type から、もう存在しない水準 Compact を取り除くには、特殊な引数 drop を使って

> fuel2$Type <- fuel2$Type[drop=T]

とします。

回帰分析で "Problem in lm.fit.qr(x, y): computed fit is singular," というエラーメッセージで実行できないことがある。

たとえば、次のようなデータがあったとします。

> データ
    月のお小遣い 体重 身長 
1             30000     45 135
2             40000     50 150
3             45000     60 180

これを使って、回帰を行います（関数 lm でも、メニューバー「統計 / 回帰 / 線形」でも、結果は同じです）。

> lm(月のお小遣い‾身長+体重, data=データ)
Problem in lm.fit.qr(unset(x), unset(y)): computed fit is singular, rank 2 
Use traceback() to see the call stack

ここでは、身長=体重 * 3 という関係が説明変数の間に成り立っており、この場合、「月のお小遣い」を体重で決めてよいのか、身長できめるべきなのか、判別できない状況にあります。これを singular と呼びます。

> lm(…, singular.ok = T)

とすると、一応答えが求まります（該当の列が係数ゼロとなっていることが分かります）。

エラーメッセージの、rank 2 の意味ですが、上の回帰では説明変数が2個あり、定数項をあわせると、全体で説明変数の数が3個になります（因子（文字の情報を持っている）列があると、列の数=変数の数になりません、ご注意ください。因子水準の数-1だけ、変数ができます）。
ところが、このデータは、上のように説明変数間に線形の関係があるので、これが、2列で表せてしまう、ということを言っているわけです。

樹形モデルの結果をplotで作図して、textでラベルを表示させたが、因子データの因子が a, b, などと表示され、どれがどれに対応しているか分からない。

サンプルデータを使って、樹形モデルを作成します。

> fit <- tree(Fuel ‾ Type, data=fuel.frame)

続いて、作図を行います。

> plot(fit)
> text(fit)

すると、分岐の表示が Type:ade など表示されますが、サンプルデータ fuel.frame の因子変数 Type にはade なる因子水準は存在しません。では、この ade はどこから来たものでしょうか？これは、デフォルトの作図オプションで「説明変数が因子だった場合、分岐条件には因子水準そのものを使わず、a, b, c, ... というアルファベットに置き換える」ということをしているためです。因子水準の文字列がとても長い場合に備えての処置です。

これを因子水準そのものにするには

> plot(fit)
> text(fit, pretty=0)

とします。詳しくは、関数 text.tree のヘルプをご覧ください。

データ編集後、数値計算がおかしくなることがある。

次のようなデータを作成します。

> data <- data.frame(x=c(".", seq(0.1, 0.9, by=0.1)), y=seq(0, 0.9, by=0.1))
> data
     x   y 
 1   . 0.0
 2 0.1 0.1
 3 0.2 0.2
 4 0.3 0.3
 5 0.4 0.4
 6 0.5 0.5
 7 0.6 0.6
 8 0.7 0.7
 9 0.8 0.8
10 0.9 0.9

SAS ユーザーの方が、SAS のデータを一旦テキストになさって、S-PLUS に読み込んだときこういう状態になっていることがありますが、SAS でピリオドが欠損値を表しているためです。

S-PLUS では、欠損をシンボルNAを使って表しますから、ピリオドを NA に置き換えます。 

> data[data=="."] <- NA
> data
     x   y 
 1  NA 0.0
 2 0.1 0.1
 3 0.2 0.2
 4 0.3 0.3
 5 0.4 0.4
 6 0.5 0.5
 7 0.6 0.6
 8 0.7 0.7
 9 0.8 0.8
10 0.9 0.9

見た目は正しく置換されているように見えます。ところが 

> sum(data$y)
[1] 4.5
> # 正しい
> sum(data$x)
Problem in Summary.factor(x, na.rm = na.rm): A factor is not a numeric object 
Use traceback() to see the call stack
> # エラー

となります。原因は x 列が最初ピリオドを含んでいたため、文字が含まれるデータとして、因子データに読み込まれたためです。では、これを as.numeric を使って数値にすればよいのですが、

> as.numeric(data$x)
 [1] NA  2  3  4  5  6  7  8  9 10

と、元の数値とは似ても似つかぬ結果となります。これは、as.numeric(因子データ)とすることによって、S-PLUS は因子データを数値に変換しますが、このとき、見た目の数値ではなく、因子水準で何番目のデータか、という順位に変換するためです。このため、計算結果が予想していたものとまったく異なることがあります。

これを防ぐには

> class(data$x)
[1] "factor"
# データは因子
> as.numeric(as.character(data$x))
 [1]  NA 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

と、一旦 as.character で文字列にしてから、数値に変換します。特に、1から１０までの整数などの時に注意が必要です。因子水準はアルファベット順で作成されるため、1,2, ...., 10 ではなく、1, 10, 2, ...., 9 となっています。

多次元正規分布に従う乱数を発生させることができるか？また、2系列の乱数に相関を指定することができるか？

関数rmvnormが利用可能です。

> x<-rmvnorm(100, rho=.9) 
> cor(x[,1], x[,2])
[1] 0.8949779
# 乱数のため、結果は異なります。

引数で、サンプル個数（行数）、サンプルの平均、標準偏差、分散共分散行列、2次元正規分布の場合の系列相関を指定することが可能です。
平均、標準偏差はd次元の正規分布に従う乱数の場合はd個のベクトルを、相関については、スカラー（値1個）あるいは、サンプル個数と同じ長さの-1から1の間の数からなるベクトルの指定が可能です。

S-PLUSの解析結果のヴァリデーションはどうなっているか？

関数validateが利用可能です。

> validate("hypotest", verbose=T)

関数validateでは、既にS-PLUSに組み込まれている解析のヴァリデーション、ユーザー独自のヴァリデーションを行うことができます。使い方などは次のPDFファイルをご覧ください。

Application Developer's Guide の19章日本語訳

また、外部の機関として、米国NIST（National Institute of Standard and Technology)が、ソフトウェアの検証のためのデータなどを用意しております。

そのNIST（National Institute of Standard and Technology) が用意した統計ソフトのヴァリデーション用のデータのページはこちら。

このうち、AOV と、Linear Regression のいくつかのデータについて S-PLUSのデータとしてあらかじめ用意しました。
zipファイルのダウンロード

例えば、Analysis of Variance の SiRstv データをS-PLUS for Windows ver 6.0以上で開きます。

  1. メニューバー「統計 / 分散分析 / 固定効果」
  2. 「Dependent」から Resistance を 「Independent」から Instrument を選択、OK

とすると、分散分析が実行され、Report1 というウィンドウに結果が 表示されます。この結果がNISTのCertified Value に等しくなることで、計算結果をご確認いただけます。（結果の桁数が、NISTのものと異なりますが、S-PLUS の場合 メニューバー「オプション / 設定」の Computations にある、 Print Digits で制御可能です）

▲topへ