チラシの裏の落書き日記

統計とか,研究とか,日常についての備忘録的なもの。

よく必要になるけど,よく忘れるR(とExcel)の小技メモ

文字(数字)の桁数を揃えたい。
参考URL:
http://d.hatena.ne.jp/teramonagi/20100902/1283429021

formatC("hoge", widh = 5, flag="0")

flagは '0 + - #' だけが指定可能とのこと。


不要な空白を消したい。

# いらない空白を全て削除する
# 全角半角の指定に注意。
gsub(" ", "", x, fixed = TRUE)
# 文字列の前後の空白を削除する
stringr::str_trim(x, side="both")

全角半角を変換したい。
Nipponパッケージを利用するといいらしいが,上手く動作しない場合がある。
そんなときは,Excelを使う。
参考URL:全角を半角に、半角を全角に変換する エクセルの関数技

文字を半角にするにはASC関数を使い,半角文字を全角にするにはJIS関数を使う。

不要な全角半角スペースを削除するには,=SUBSTITUTE(SUBSTITUTE(A1," ","")," ","")という関数を利用するとよい。
(trimは文字列前後の空白を削除するための関数。)


Rで図を出力する場合のあれこれ。
参考URL
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/53.html
R のグラフで日本語を使う - 知に至る病
https://oku.edu.mie-u.ac.jp/~okumura/stat/first.html
Rのグラフで軸の目盛りの刻み幅を変更する方法 - Rプログラミングの小ネタ

# 文字のフォントを指定する。
# これを指定しておけば,とりあえず文字化けしない。
pdf(file="hoge.pdf",  family="Japan1GothicBBB")

# これはMac向き
par(family="HiraKakuProN-W3")  # ヒラギノ角ゴシックProN W3

# 空のプロットを作成する。
 plot.new()

# tclで軸を内側に向ける。
# 正の値が内側で,枠の大きさに対する割合で指定で指定する。0.3ぐらいがけっこうちょうどいい。
plot(0,1,tcl = 0.3)

# メモリを独自に振り直したいときに使う。
#メモリを消す
curve(dnorm, xlim=c(-3, 3), xaxt="n") 
# atに書き込みたい目盛りのベクトルを指定する
axis(side=1, at=-3:3) 
# sideは下から時計回りに1,2,3,4になっている。

# ラベルを各軸に並行して描く
las = 0
# ラベルをすべて水平に描く
las = 1
# ラベルを軸に対して垂直に描く
las = 2
# ラベルをすべて垂直に描く
las = 3

Mplusでターゲット行列を指定した回転を行う方法

Analysisコマンドで,rotation = targetを指定し,
modelコマンド内で目標とする値を~で指定する。
また,各因子の最後の(*t)などとつけて,各因子負荷が1つの因子負荷行列で表すことができることを示さないといけない。
Mplusのマニュアルは,versin 8のp681などに記載されている。正直,いろいろと回転がありすぎて何が何だかわからない。
本文をtargetで検索するのがはやい。
これを使えば,任意の値を指定したターゲット行列に向けた回転を行ってくれる。
(つまり,プロクラステス(procrustes)回転ができる。)
確認的因子分析よりもゆるく,探索的因子分析よりも仮説を反映した分析ができるので,もう少し利用されてもいい方法のように感じる。


使用したデータ
https://www.statmodel.com/usersguide/chap4/ex4.1a.dat

TITLE:
DATA:
FILE IS ex4.1a.dat;
VARIABLE:
NAMES ARE y1-y12;
ANALYSIS: 
ESTIMATOR IS ML;
 ROTATION=TARGET;
 MODEL:
f1 BY y1-y12 y1~0 y2~0 y3~0 y4~0 y5~1 y6~1(*t);
f2 BY y1-y12 y7~0 y8~0 y9~0 y10~0 y11~1 y12~1(*t);
  OUTPUT:
MODINDICES CINTERVAL;
  

参考PDF
http://www.restore.ac.uk/appliedpsychometrics/Tutorial.materials/Introduction.to.Mplus/Introduction%20to%20Mplus%20session%202.pdf

Wordの小技(大文字小文字の変換,テキストボックス内の文字列の検索)

Wordで便利な小技を見つけた。
どちらも,「高度な検索と置換」を使う。
特に,ワイルドカードと,検索対象を指定するのがポイント。
「[0-9]」で検索を行うと,大文字の数字が選択される。その上で,「A->a」で「半角」を選択すればOK。

テキストボックス内の文字列の選択では,ワイルドカードの「*」を使って,検索対象を「メイン文書」ではなく「メイン文書のテキストボックス」にして検索をする。その後,検索ウィンドを閉じて「F9」をつかって,フィールドの更新をすればよい。

Wordで数字を大文字から小文字に変換する。
office-qa.com


Wordで複数のテキストボックスの文字列を選択する。
detail.chiebukuro.yahoo.co.jp


そして,フィールドについて,何処かでまとめて勉強したほうがいい気がする。
あとで,まとめて読んでみること。

pasonyu.com

pasonyu.com

フィールドについて

フィールドについて:Word(ワード)2010基本講座

平均値の信頼区間のシミュレーションプロット

Rで,信頼区間を書いていいプロットするコードのメモ。
layout関数で,上手くプット画面を分割して,さも1枚のグラフかのように見せかけるのがポイント。
けっこう可愛い感じのグラフを描画できる。

# 
# 縦向きの,信頼区間の図。
# ポイントは,レイアウトとparの使い方。
# 
layout(matrix(c(1,2),ncol=1), heights=c(1/5,4/5))
par(mar=c(0,3,1,1))
mu <- 0
sigma <- 1
n <- 5
size <- 100

curve(dnorm(x,mu,sigma/sqrt(n)),from=-3, to = 3,axes=F, xlab = "",ylab="",lty =1, col =  rgb(1,0,1,alpha = 0.7),lwd = 2)
arrows(mu, 0, mu, dnorm(mu,mu,sigma/sqrt(n)), lwd=2, angle = 0,lty =1, col =  rgb(1,0,1,alpha = 0.7))

x <- matrix(rnorm(size*n,mean = mu,sd = sigma), ncol=n)
x_mean <- rowMeans(x)
upp <- x_mean + 1.96*sigma/sqrt(n)
low <- x_mean - 1.96*sigma/sqrt(n)
is_sig <- ((upp < 0) + (low < 0)) !=1

par(mar=c(3,3,0,1))
plot(x_mean,1:size, xlim = c(-3,3),col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)), pch = 1- is_sig ,xlab="x平均")
arrows(x_mean, 1:size, upp, 1:size, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
arrows(x_mean, 1:size, low, 1:size, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
abline(v = mu , lwd = 2,lty=2,col =  rgb(1,0,1,alpha = 0.7))

出来上がる図は次のようなもの。
f:id:kazzstat:20171020204451p:plain

#
# 横向きの,信頼区間の図。
#

layout(matrix(c(2,1),nrow=1), widths =c(4/5,1/5))
par(mar=c(3,0,1,1)) # 下,左,上,右の順で,余白を設定する。上が0なので,いい感じでくっつけられる。
mu <- 0
sigma <- 1
n <- 10
size <- 100

# 右に90度,図を転置させる。
t_dnorm <- function(x, mu, sigma){
  cbind(x, dnorm(x,mu,sigma)) %*% matrix(c(0,-1,1,0),ncol=2,byrow=T)
}
plot(t_dnorm(seq(from=-3, to =3, by=0.1), mu, sigma/sqrt(n)),type="l",axes=F, xlab = "",ylab="",lty =1, col =  rgb(1,0,1,alpha = 0.7),lwd = 2)
arrows(0, mu, dnorm(mu,mu,sigma/sqrt(n)), mu, lwd=2, angle = 0,lty =3, col =  rgb(1,0,1,alpha = 0.7))


x <- matrix(rnorm(size*n,mean = mu,sd = sigma), ncol=n)
x_mean <- rowMeans(x)
upp <- x_mean + 1.96*sigma/sqrt(n)
low <- x_mean - 1.96*sigma/sqrt(n)
is_sig <- ((upp < 0) + (low < 0)) !=1

par(mar=c(3,3,1,0)) # この余白の設定で,絶妙に図をくっつけるのがポイント。
plot(x_mean, ylim = c(-3,3),col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)), pch = 1- is_sig ,xlab="x平均")
arrows(1:size,x_mean, 1:size, upp, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
arrows(1:size,x_mean, 1:size, low,lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
abline(h = mu , lwd = 2,lty=3,col =  rgb(1,0,1,alpha = 0.7))

出来上がる図は次のようなもの。
f:id:kazzstat:20171020204524p:plain

日々淡々と進むこと

最近思うこととしては、とにかく毎日進み続けることが何よりも大切。

ルーチンワークでやらなくてはならない雑務も多く、体力や気力も削られるかもしれないが、その中でも研究をし続け、論文をなんとしても書き続ける必要がある。

書き続けていれば、書くことの労力も減っていく。自分の中に書くということを続けることが大事。

 

スポーツをやるとこの辺りの感覚もわかる。

自分が練習しないと記録も伸びない。結局、自分との戦いだ。

Mplusによるモンテカルロ・シミュレーション覚書2

今回は前回紹介しなかったコマンドなどについて,紹介する。しかしあまり気合が入っていないのでグダグダな感じになっている。

SCALE OF DEPENDENT VARIABLES FOR ANALYSIS

CENSORED,CATEGORICAL,NOMINAL,COUNTオプションが分析の従属変数に利用できる。指定方法は前述したgenerateと類似していいて,「CENSORED = y1(a);」などと指定する。
この指定によって,モデルと推定での変数の扱いが決まる。
独立変数としては,2値あるいは連続値が指定できる。

OPTIONS FOR DATA ANALYSIS

CLASSES, AUXILIARY, SURVIVAL のみがanalysisコマンドで利用できる。

CLASSES

CLASSESはカテゴリカル潜在変数の名前とクラス数を指定する。
これを使うためにはTYPE =MIXTUREが必要。
Betweenレベルのカテゴリカル潜在変数を使う場合には,betweenオプションを使って,宣言をする必要がある。
使い方としては,「CLASSES = c1 (2) c2 (2) c3 (3);」などと指定する。

AUXILIARY

補助変数は分析モデル内では利用されない変数である。
TYPE=MIXUREとともに,AUXILIARYオプションは自動的に3ステップアプローチ(Asparouhov and Muthén, 2014a, b 参照)を実行するために使用される。
詳細はマニュアルp.880-881を参照。

SURVIVAL

このオプションはTIMECENSOREDオプションと同時に使われないといけない。
気が向いたら追記します。
p.881-883

VARIABLES WITH SPECIAL FUNCTIONS FOR DATA GENERATION AND ANALYSIS

WITHIN

「TYPE=TWOLEVEL:」の場合は,個人レベルの変数を「WITHIN = y1 y2 x1;」と指定する。ここで指定しない変数は,within, betweenの両方でモデル化される。
「TYPE=THREELEVEL:」の場合は,「WITHIN = y1-y3 (level3) y7-y9 (level2) y4-y6;」とする。レベル1,3, 2の順に指定することに注意。y1-y3は個人レベルで,レベル1でのみモデル化される。レベル1変数は他の変数に先立ってモデル化される必要がある。y4-y6はレベル1とレベル2でモデル化される。y7-y9はレベル1とレベル3でモデル化される。
あるいは
「WITHIN = y1-y3;
WITHIN = (level2) y4-y6;
WITHIN = (level3) y7-y9;」
などと指定する。この場合には順番は気にしない。
「TYPE=CROSSCLASSIFIED:」の場合には「WITHIN = y1-y3 (level2a) y4-y6 (level2b) y7-y9;」と指定する。

BETWEEN

クラスターレベルで測定された変数の指定。
withinオプションと同様の指定方法になる。

POPULATION, COVERAGE, AND STARTING VALUES

POPULATIONと COVERAGEは前回の記事参照。
また,「STARTING = estimates.dat;」によって初期値を指定することも出来る。

SAVING DATA AND RESULTS

「REPSAVE = 1 10-15 100;」で保存するデータを指定できる。
単一のデータを保存する場合には,「SAVE = rep1.dat;」とする。
「SAVE = rep*.dat;」とすると,レプリケーションの個別のデータを指定できる。

「RESULTS = results.sav;」で各推定結果を利用できる。
しかし,すこぶる使いにくい。outファイルに以下のような出力が得られ,パラメタの順がわかる。(PARAMETER SPECIFICATIONの部分でパラメタの順も得られるので,それを参照する必要がある。)

RESULTS SAVING INFORMATION

  Order of data

    Replication number
    Parameter estimates
     (saved in order shown in Technical 1 output)
    Standard errors
     (saved in order shown in Technical 1 output)
    Chi-square : Value
    Chi-square : Degrees of Freedom
    Chi-square : P-Value
    Wald Test : Value
    Wald Test : Degrees of Freedom
    Wald Test : P-Value
    H0 Loglikelihood
    H1 Loglikelihood
    Number of Free Parameters
    Akaike (AIC)
    Bayesian (BIC)
    Sample-Size Adjusted BIC
    RMSEA : Estimate
    SRMR

推定法にbayesを指定した場合に,全てのイテレーションのパラメタの事後分布のパラメタ値が得られる。
「BPARAMETERS = bayes.dat;」

時系列解析では,ラグを指定できる。
例:「LAGGED = y (1);」
最大のラグは2である。

Mplusのモンテカルロ・シミュレーションコマンド覚書1

「MONTECARLO:」コマンドを使った場合に一緒に利用するコマンド・オプションの一部を一覧にしておく。

GENERAL SPECIFICATIONS

「NAMES = y1-y10;」は変数名を指定するオプション。

「NOBSERVATIONS = 500;」はデータ生成と分析に利用するサンプルサイズを指定するオプション。

多母集団の場合には,「NOBSERVATIONS = 500 1000;」とするとグループ1は500人,グループ2は1000人という意味になる。
さらに,「NOBSERVATIONS = 2(100) 10(500);」とすると,サンプルサイズ100人のグループが2つ,500人のグループが10つとなる。

「NGROUPS = 3;」で3つの母集団を仮定する。NGROUPSはデフォルトでは1で,1の場合にはTYPE=MIXTUREが使用できない。
なお,グループ名はg1, g2, ...と自動で名前が振られる。
TYPE=MIXTUREとともに使う場合には,既知のクラス数を指定する。ラベル「%g#1%」がはじめの既知のクラス,「%g#2%」が2番目のクラスを表す。これらのクラスを用いて,モデルを指定する。

「NREPS =100;」はシミュレーションの反復回数を表す。デフォルトは1。

「SEED = 123;」は乱数生成をするときのシード。結果の再現をしたい場合にはきちんと指定する必要がある。デフォルトは0。

DATA GENERATION

「GENERATE」コマンドは,生成する従属変数のスケールを指定する。
その他,打ち切り(censored),2値(binary),順序カテゴリ(ordered categorical; ordinal),名義(nominal),カウント(count),time-to-event変数の生成が出来る。
例えば「GENERATE = u1-u2 (1) u3 (1 p) u4 (1 l) u5 u6 (2 p) y1 (ca 1) y2 (cbi 0) u7 (n 2) u8 (ci) t1 (s 5*1);」などのように利用する。カッコ内の文字の意味を以下に説明する。

打ち切り変数は,以上・以下,インフレーションを含む・含まない,形式で生成できる。
「y1(ca 0)」は打ち切り値0より上の値を打ち切る指定。
「y1(cb 0)」は打ち切り値0より下の値を打ち切る指定。
「y1(cai 0)」はinflationを伴って,打ち切り値0より上の値を打ち切る指定。
「y1(cbi 0)」はinflationを伴って,打ち切り値0より下の値を打ち切る指定。

カテゴリカル変数の指定は以下のように行う。
「y1(1 l )」は一文字目は数字のいち, 二文字目はローマ字小文字のエルである。これで,ロジスティックモデルによる2値のカテゴリカル変数の生成が出来る。
閾値の数はカテゴリ数マイナス1である。)
「y1(1 p )」はプロビットによる2値のカテゴリカル変数の生成が出来る。(2パラメタ正規累積モデル, two-parameter normal ogive model.)
もし,lやpの指定がない場合にはロジスティックモデルが使用される。
ちなみに,「y1(1 3pl )」や「y1(1 4pl)」とすると,3パラメタ,4パラメタロジスティックモデルを指定出来る。
順序カテゴリを生成する場合には,段階反応モデル(graded response model)と同じ比例オッズモデル(proportional odds model)が使用される。
「y1(3 gpcm)」を使うと閾値3つで一般化部分採点モデル(generalized partial credit model)に従うデータを生成出来る。
また,推定方法が重み付き最小自乗法(weighted least squares estimation)を使う場合,2値と順っ女カテゴリカル変数には,プロビットモデルのみが許容される。
「y1(3 n)」は切片数3(カテゴリ数4)の名義変数を生成出来る。最後のカテゴリがレファレンスカテゴリとして切片が0と固定される。

カウント変数としては,Poisson, zero-inflated Poisson, negative binomial , zero-inflated negative binomial, zero-truncated negative binomial, negativi binomial hurdleの6種類が利用できる。
「x(c)」あるいは,「x(p)」でPoisson,「x(ci)」「x(pi)」でzero-inflated Poissonが指定できる。
「x(nb)」でnegative binomial, 「x(nbi)」でzero-inflated negative binomial モデル,「x(nbt)」でzero-truncated negative binomial modelが指定できる。「x(nbh)」でnegativi binomial hurdleを指定できる。

Time-to-event変数を指定するためには,文字「s」とベースライハザード関数の時間間隔の長さの数を指定する。もし,文字「s」飲みが指定された場合には,観測数と同じ長さの間隔数が指定される。
例えば,「t1 (s 5*1)」は長さ1の5つの時点間隔がデータ生成に利用されることを意味する。

MODEL POPULATIONでは,打ち切り,カウント変数のinflation partは数1に続くナンバーサイン(#)変数の打ち切り・カウント変数の名前を追加することによって参照される。

CUTPOINTS

カテゴリカル独立変数を指定するために使用される。
例えば,「CUTPOINTS = x1(0) x2(1);」とすると,x1は0, x2は1をカットポイントを持つ。カットポイントより低い値は0,高い値は1となる。
多母集団分析では,記号「|」で区切って指定を行う。例えば,「CUTPOINTS = x1 (0) x2 ( 1) | x1 (1) x2 (0);」とすると,グループ1と2で異なったカットポイントを指定出来る。

GENCLASSES

「GENCLASSES」オプションはカテゴリカル潜在変数を生成するために使用される。このオプションは,TYPE = MIXUREとともに利用される。
例えば,「GENCLASSES = c1 (3) c2 (2) c3 (3);」とすると,c1, c2, c3という名前の3つのカテゴリカル潜在変数を生成する。c1は3,c2は2,c3は3つのクラスを持つことを示す。

「GENCLASSES = cb (2 b);」のように文字「b」を着けると,between-levenのカテゴリ潜在変数を指定出来る。この変数は「BETWEEN」オプションを使用するbetween-level変数として必ず指定されなければならない。

NCSIZES

「NCSIZES」は,TYPE=TWOLEVEL,TYPE=THREELEVEL,TYPE=CROSSCLASSIFIEDとともに使用されなければならない。
例えば,TYPE=TWOLEVELに対しては,「NCSIZES = 3;」はデータ生成で使用されるユニークなクラスターサイズ3を指定している。
多母集団の場合には,「NCSIZES = 3 | 2;」とすると,グループ1では3つのユニークなクラスター,グループ2では2つのユニークなクラスターを生成する。
TYPE=THREELEVELでは,学生はクラスにネストされており,クラスは学校にネストされている。例えば,「NCSIZES=3[2];」と指定される。3(レベル3,学校)と2(レベル2,クラス)はユニークなクラスサイズとしてデータ生成に使用される。
多母集団の場合には,「NCSIZES = 3 [2] | 4 [3];」とすると,グループ1(レベル3で3つ,レベル2で2つのユニークなクラスタ)とグループ2(レベル3で4つ,レベル2で3つのユニークなクラス)で異なったクラスターサイズを指定できる。

TYPE =CROSSCLASSIFIEDを指定した場合には,生徒は学校と地域にネストされる。この場合,レベル1は生徒,レベル2aは学校,レベル2bは地域となる。
ここで「NCSIZES = 3 [2];」とすると,3はレベル2a, 2はレベル2bに対応する。分析のタイプがTHREELEVELかCROSSCLASSIFIEDかによって,NCSIZESの指定の意味が異なることに注意。

CSIZES

CSIZESはクラスターとサイズを指定する。
TYPE = TWOLEVELでは,「CSIZES = 100 (10) 30 (5) 15 (1);」と指定すると,サイズ10のクラスタを100クラスター,サイズ5のクラスターを30, サイズ1のクラスターを15がデータ生成に使用される。
多母集団では,「CSIZES = 100 (10) 30 (5) 15 (1) | 80 (10) 20 (5);」とすると,グループごとにクラスサイズを別に指定できる。
TYPE=THREELEVELの場合,「CSIZES = 40 [15(2) 10(5)] 30 [6(8)] 7 [20(2)];」と指定する。このようにすると,レベル3のクラスタ数が40, 30, 7となる。さらに,レベル3が40のものに対して,レベル1が2で15このレベル2クラスターとレベル1が5でレベル2のクラスターサイズが10こ,データ生成に利用される。この場合,3200の観測が得られる。以下同様。
多母集団の場合は「CSIZES = 30 [6(8)] 7 [20(2)] | 40 [5(6)] 20 [4(2)] ;」とする。

TYPE=CROSSCLASSIFIEDの場合,「CSIZES = 40 [15(2) 10(5)] 30 [6(8)] 7 [20(2)];」とすると,レベル2bのクラスタ数は角括弧の外の数(40, 30, 7),レベル2aのクラスタ数は角括弧の中,()の外の数(15, 10, 30, 20)を示す。レベル1のサイズは()内の数となる。

HAZARDC

HAZARDCは連続時間の生存解析の打ち切りプロセスのハザードを指定するために使用される。
0が打ち切られていない,1が右側打ち切りである打ち切りインディケータを生成するために使われる。
例えば,「HAZARDC = t1 (.5);」ここで,「t1」はtime-to-event変数名,.5は打ち切りのハザードを示す。

MISSING DATA GENERATION

Missingの生成するためには,PATMISSとPATPROBSを利用する。
例えば,「PATMISS = y1 (.2) y2 (.3) y3 (.1) |
y2 (.2) y3 (.1) y4 (.3) |
y3 (.1) y4 (.3);」とすると3種類の欠測パタンを生成できる。
「PATPROBS = .4 | .3 | .3;」とすると,PATMISSで指定した3つの欠測パタンが生じる確率を指定する。

MISSINGオプションは,どの従属変数で,欠測データが生成されるかを指定する。このオプションはMODLE MISSINGコマンドと同時に使用される。
Missing dataに観測独立変数を指定することはできない。
例えば,「MISSING = y1 y2 u1;」と指定される。
MODEL MISSINGは以前の記事で説明したように,ロジスティック回帰モデルで欠測確率を指定する。

次の記事で紹介する予定?

気が向いたら,以下のものについても紹介します。

  • SCALE OF DEPENDENT VARIABLES FOR ANALYSIS
  • OPTIONS FOR DATA ANALYSIS
  • VARIABLES WITH SPECIAL FUNCTIONS FOR DATA GENERATION AND ANALYSIS
  • POPULATION, COVERAGE, AND STARTING VALUES
  • SAVING DATA AND RESULTS