Rを使って負の二項分布で回帰する Negative binomial regression

裾の長いカウントデータを扱いたい（外れ値が有る）ときに、負の二項分布を仮定したモデリングを考える。

MASS packageにはいっていたので、覚書として残しておく。

R: Fit a Negative Binomial Generalized Linear Model

以下で実行できる。

library(MASS)
quine.nb1 <- glm.nb(Days ~ Sex+Age, data = quine)
confint(quine.nb1)

仮想通貨の勢力図をRを使って可視化する

今日は少し脱線して、仮想通貨情勢をRで調べてみる。

この記事を見ながら、自分で実装してみた。

今回、使うのは3つのパッケージ。

元記事にあったとおり、coinmarketcaprで仮想通貨のデータをとってきて、treemapで可視化。

米ドルではなく日本円で計算してみるため、quantmod packageを使って、ドル円データを引っ張ってくる。

良い感じに出来た（参考までに、図は、ブログ更新時点のもの）。

library(coinmarketcapr)
market_today <- get_marketcap_ticker_all()

#install.packages("quantmod")
library(quantmod)
from <- c("USD")
to <- c("JPY")
usdjpy<-getQuote(paste0(from, to, "=X"))
usdjpy$Last

library(treemap)
df1 <- na.omit(market_today[,c('id','market_cap_usd')])
df1$market_cap_usd <- as.numeric(df1$market_cap_usd)*as.numeric(usdjpy$Last)
colnames(df1) <- c("id", "market_cap_jpy")
df1$formatted_market_cap <-  paste0(df1$id,'\n','¥',format(df1$market_cap_jpy,big.mark = ',',scientific = F, trim = T))
treemap(df1, index = 'formatted_market_cap', vSize = 'market_cap_jpy', title = 'Cryptocurrency Market Cap', fontsize.labels=c(12, 8), palette='RdYlGn')

CRANから消えた"Archived R package"をインストールしたい

質的研究に目覚めた朝。

論文を読んで、package "concord"をインストール！

しようにもCRANの表舞台から削除されている。
そんな時、アーカイブされているpackageをインストールするのが、以下の方法。

*こちらを参考にしました。

# Download package tarball from CRAN archive （以下のような感じで格納）

url <- "https://cran.r-project.org/src/contrib/Archive/concord/concord_1.4-9.tar.gz"
pkgFile <- "concord_1.4-9.tar.gz"
download.file(url = url, destfile = pkgFile)

# Install dependencies （何か、dependenciesなパッケージが有れば入れる）

install.packages(c("xxxx"))

# Install package
install.packages(pkgs=pkgFile, type="source", repos=NULL)

# Delete package
unlink(pkgFile)

これで、消えたpackage問題は解決されます！！

論文を書く際の統計的表現について

論文執筆の際、p値はどう記載するんですか？

有効数字は？

有意差なしはNSだけでよいですか？

0.05未満は<0.05でよいのですか？

・・・

このように聞かれたり、議論になったりすることがあるので、纏めておきます。

あくまで、一例ですが、Annals of Internal Medicine (IF 17)の投稿規定は以下のようになっています。

p値は0.001より小さい場合は、<0.001、0.001-0.20は有効数字3桁、0.20より大きい場合で2桁、と書かれています。

一流誌の投稿規定は、それ自身が勉強になることが多いですね。

以下はAnnals of Internal Medicineからの引用です。

Reporting guideline

##Percentages

Report percentages to one decimal place (i.e., xx.x%) when sample size is =200. 

To avoid the appearance of a level of precision that is not present with small samples, do not use decimal places (i.e., xx%, not xx.xx%) when sample size is < 200.

##Standard deviations

Use “mean (SD)” rather than “mean ± SD” notation. The ± symbol is ambiguous and can represent standard deviation or standard error.

##Standard errors

Report confidence intervals, rather than standard errors, when possible.

##P values

For P values between 0.001 and 0.20, please report the value to the nearest thousandth. For P values greater than 0.20, please report the value to the nearest hundredth. For P values less than 0.001, report as “P<0.001.”

##“Trend”

Use the word trend when describing a test for trend or dose-response.

Avoid the term trend when referring to P values near but not below 0.05. In such instances, simply report a difference and the confidence interval of the difference (if appropriate) with or without the P value.

##Statistical software

Specify in the statistical analysis section the statistical software—version, manufacturer, and the specific functions, procedures, or programs—used for analyses.

##Cox models

When reporting findings from Cox proportional hazards models:

Do not describe hazard ratios as relative risks.
Do report how the assumption of proportional hazards was tested, and what the test showed.

##Descriptive tables

In tables that simply describe the characteristics of 2 or more groups (e.g., Table 1 of a clinical trial):

Report averages with standard deviations, not standard errors, when data are normally distributed.
Report median (minimum, maximum) or median (25th, 75th percentile [interquartile range, or IQR] when data are not normally distributed.
Avoid reporting P values as there can be imbalance when P values are not significant (because of small sample size) and balance when P values are significant (because of large sample size).

##Tables reporting multivariable analyses

Authors sometimes present tables that compare one by one an outcome with multiple individual factors followed by a multivariable analysis that adjusts for confounding. If confounding is present, as is often the case, the one-way comparisons are simply intermediate steps that offer little useful information for the reader. In general, omit presenting these intermediate steps in the manuscript and do not focus on them in the Results or Discussion.

##Figures

When developing informative graphics, follow these simple rules of thumb:
Avoid pie charts.
Avoid simple bar plots that do not present measures of variability.
For meta-analysis forest plots, provide the raw data (numerators and denominators) in the margins.
For survival plots, provide the numbers of people at risk by group and time below the horizontal axis.

Rを使って時系列データの変化点をみつける: changepoint

時系列データを扱う時に、変化点を見つけたいことがある。

Rにchangepointというパッケージがあるので、実装してみる。

• １つのchangepointはdefaultのmethod = AMOCで良い。

• 複数のchangepointであれば、methodでPELT, SegNeighやBinSegを指定すれば良い。

#install.packages("changepoint")
library(changepoint)

# change in variance
x=c(rnorm(100,0,3),rnorm(100,0,10))
res1=cpt.var(x, )
plot(res1)
print(res1) 

# change in mean
y=c(rnorm(100,0,1),rnorm(100,7,1))
res2=cpt.mean(y)
plot(res2,cpt.col='blue')
print(res2)

# change in mean and variance
z=c(rnorm(100,0,1),rnorm(100,3,7))
res3=cpt.meanvar(z)
plot(res3,cpt.width=3)
print(res3)

# change in mean and variance, multiple changepoints
w=c(rnorm(100,0,1),rnorm(100,3,7), rnorm(100,-5,1))
res4=cpt.meanvar(w, Q=2, method = "PELT")
plot(res4)
print(res4)

計量経済分析ことはじめ~時系列データを正しく扱う~

Rによる計量経済分析を読んだ。このシリーズはとても勉強になる。

書籍のHPはこちら

キーワードと、キーとなるパッケージ・関数を列挙しておく。
時系列データを扱う時の問題点がわかりやすく解説されていて、よかった。

クロスセクションデータの回帰分析

最小二乗推定量 (OLS：ordinary least squares estimator)
決定係数 (coefficient of determination)
BLUE (the best linear unbiased estimator)

不均一分散

実世界のデータを扱う際に標準的仮定が満たされない場合がある。其の１つが不均一分散。

不均一分散の検定は以下で実行できる。

1. Breusch-Pagan 検定

2. White 検定

3. Goldfeld-Quandt 検定: データを2分割して、それらを同じモデルを回帰した場合に誤差項の分散が等しくなるかどうかを検定

bptest(), gqtest()の説明書を見ながら、以下を実行。

library(lmtest)

## generate a regressor
x <- rep(c(-10,10), 100)

## generate heteroskedastic and homoskedastic disturbances
err1 <- rnorm(100, sd=rep(c(1,2), 50))  
err2 <- rnorm(100) 

## generate a linear relationship
y1 <- 1 + x + err1
y2 <- 1 + x + err2

## perform studentized Breusch-Pagan test
bptest(y1 ~ x)
bptest(y2 ~ x)

## perform Breusch-Pagan test
bptest(y1 ~ x, studentize = F) #studentize is a logical number
bptest(y2 ~ x, studentize = F)

## perform Goldfeld-Quandt test
gqtest(y1 ~ x, point=0.5) #point is interpreated to be a breakpoint.
gqtest(y2 ~ x, point=0.5)

時系列データの回帰分析

系列相関

時系列データでは、系列相関という問題が発生する。ある時点の回帰分析の誤差項と、其の隣の誤差項が相関している状態のことを指す。時系列データは、各時点毎に観測されたものであり、観測値の順番に意味があるので、注意が必要だ。

研究者は、このことをしっかり意識する必要がある。

系列相関の検定

1. Durbin-Watson検定
2. Breusch-Godfrey検定

library(lmtest)
## perform Durbin-Watson test
dwtest(y1 ~ x)

## perform Breusch-Godfrey test for first-order serial correlation:
bgtest(y1 ~ x)
## or for fourth-order serial correlation
bgtest(y1 ~ x, order = 4)

系列相関がある場合の推定

#FGLSを用いる方法
library(nlme)
res <- gls(y1~x, corr=corARMA(p=1))
summary(res)

#Newey-Westの修正を行う
library(sandwich)
coeftest(lm(y1~x), vcov=NeweyWest)

Stationary time series analysis (定常時系列分析)

定常過程 (Stationary process): 期待値・分散・共分散が時点tに依存せず、一定である。
independent and identically distributed; IID

モデル化には、自己回帰と移動平均を考慮して考える。

1. Autoregressive model （自己回帰）: p次のAR modelはAR(p)と記述される。

2. Moving average model （移動平均）

3. Autoregressive moving average model （自己回帰移動平均）

4. Autoregressive integrated moving average model （自己回帰和分移動平均）

ARCHとGARCH

時期によって分散の大きさがことなることを考慮する場合に用いる。

1. (Generalized) autoregressive conditional heteroskedasticity model

GARCH(1,1)でモデリングしたいときは以下。

library(fGarch)
GARCH <- garchFit(formula = ~ garch(1, 1), data = data, 
         trace = FALSE, include.mean =TRUE)
summary(GARCH)

多変量時系列

複数の時系列データ（Aの株価とBの株価など）が互いに影響を及ぼしうるときに使う。

VAR vector autoregressive model

定常過程かどうか

library(tseries)
adf.test(data)
pp.test(data)

モデルの次数決定

library(vars)
VARselect(data, lag.max=5, type="const")
vars.res <- VAR(data, p=1, type="const")
summary(vars.res)
causality(vars.res, cause="A")
causality(vars.res, cause="B")

非定常時系列

単位根検定: ある時系列がランダム・ウォークに従っているかどうかの検定
見せかけの回帰: ランダム・ウォークに従っている者同士を回帰すると、見せかけの回帰（spurious regression）を来すことあり。 共和分分析: ランダム・ウォークに従っている者同士でも、共和分関係にあるものであれば、回帰にも意味がある。

library(tseries)
adf.test(data) #Augmented Dickey-Fuller
pp.test(data) #Phillips-Perron

library(urca)
ca.po(data) #Philips-Ouliaris
ca.jo(data) #Johansen

パネル

複数の同じ経済主体に関する複数時点でのデータである。

fixed effectとrandom effectの推定結果の検定: Hausman検定

library(plm)
phtest(res.fixed, res.random)

効率的なR運用を目指して

Rの基礎とプログラミング技法

Rの基礎とプログラミング技法に、Rを効率よく扱うコツについて非常にわかりやすく載っていたので、トライしてみた。

第5章：効率的なプログラミング

apply

1. applyは、 ベクトル単位で処理できる
2. lapplyはリストやデータフレームごと処理できる
3. sapplyはlapplyの簡易版（返ってくるのはlistではなくベクトル）
4. mapplyは多変量に適用できる
5. tapplyはクロス集計表を作成できる

#apply

x<-matrix(c(1:9), nr=3)
apply(x,1,mean) #第二の引数は次元、ここでは行。2だと列。

#lapply, sapply

data(anscombe)
attach(anscombe)
ans.reg<-vector(4, mode="list")

for(i in 1:4){
x<-get(paste("x", i, sep=""))
y<-get(paste("y", i, sep=""))
ans.reg[[i]] <- lm(y~x)
}

lapply(ans.reg, coef) #データフレーム、リスト、ベクトルなどの構造に、coef()関数を適用
sapply(ans.reg, coef) #sはsimplifyであり、返されるのはベクトル

lapply(ans.reg, "[", 1) #このようにも[]を関数としても使える。

#mapply

mapply(sum, 1:10, 2:11, 3) #多変量に同じ処理が出来るが、一つ目の引数が関数となる。

#tapply

data(iris)
attach(iris)
tapply(Sepal.Width, Species, mean)

data(iris)
attach(iris)
tapply(Sepal.Length, Species, range)

効率性を分析する

以下のようなツールが搭載されている。

proc.time() #Rを起動してからの経過時間を出力する
system.time() #式の実行時間を出力する

Rprof() #コレ以降の実行時間を分析開始
Rprof(NULL) #これまでで終了
summaryRprof() #summary表示

gc() #R実行中のメモリの消費状況を調べるには、gc()を使う。其の点で使われていないメモリ領域が解放される。

さらに詳しい分析は、proftools packageなどが使える。

バッチ処理（第9章：拡張）

一括処理とも呼ばれる。定期的に実行しなければならない処理や、エンジニア同士でシェアするような状況では、Rを立ち上げずとも、以下のようにOSのバックグラウンドで出来るようにすべきである。

例えば、xxxx.rというファイルをRで実行して、output.txtを出力するには、以下をターミナル（macでは）で実行すればよい。

R CMD BATCH xxxx.r output.txt