ATEとしてハザード比を求める必要があったので、IPWを使った。 ダミーデータで備忘録を残しておく。 library(survminer) library(survival) library(cobalt) library(tidyverse) # データの準備 diabetic_mod <- diabetic %>% filter(eye == "left") %>% mutate(flag_highrisk = if_else(risk >= 10, 1, 0)) # KM曲線の描画 surv <- survfit(Surv(time, status) ~ trt, data = diabetic_mod) ggsurvplot(surv) # Cox回帰 fit <- coxph(Surv(time, status) ~ trt, data = diabetic_mod) summa

写真から顔を検出したら集中線を描くのが紳士のたしなみ。 以下のアルゴリズムを参考にRでタカヤナギ＝サンに集中線を描きます。 http://stamefusa.hateblo.jp/entry/20120115/1326636719 顔領域の検出については以前の記事を参考のこと。 ひとまずこんな感じになります。 まだ集中線とは言い難いのでここからいい感じのパラメータを探したい。 library("imager") library("httr") MS_FACE_KEY <- "あなたのスクリプトキー" ### 顔領域の検出 url_image <- "https://raw.githubusercontent.com/dichika/ojisan/master/takayanagi_dj.jpg" # タカヤナギ=サン data_image <- load.image(url_image)

プロバイダわからん。 とりあえず下りの速度が遅すぎなければいい。 ということで口コミサイトの結果から比較する。 以下書きなぐりコード。 library(rvest) library(dplyr) library(tidyr) library(stringi) # 以下よりソースをダウンロード # http://netspeed.studio-radish.com/cgi-bin/netspeed/openresult/openresult.cgi dat <- read_html("~/downloads/result.html") numbers <- dat %>% html_nodes(xpath="//*[contains(@class,'graphmarker')]") %>% html_text() tags <- dat %>% html_nodes(xpath="//nob

いつも忘れるのでメモ。 round_date()は気をきかせすぎな気もするのでfloor_date()をとりあえず使っていこう。 library("lubridate") > floor_date(as.Date("2017-03-22"),"month") [1] "2017-03-01" > round_date(as.Date("2017-03-22"),"month") [1] "2017-04-01"

下記サイトにあるようなbeeswarm（蜂群図）を描きたい。 https://www.yodosha.co.jp/jikkenigaku/statistics_pitfall/pitfall_5.html そんな時はbeeswarmパッケージを使えば良い。 http://www.cbs.dtu.dk/~eklund/beeswarm/ だがggplot2でやりたい。 position_jitterで調整してやることでそれっぽいものはできる。 widthで横方向、heightで縦方向に散らす程度を設定できる。 geom_dotplotで解決する。 詳しくは以下のリンクを見られたし。（id:kaz_yosさんありがとうございます！） http://rpubs.com/kaz_yos/1406 ggplot(iris, aes(x=Species, y=Sepal.Length))+geom_p

ggplot2でラベルを日本語にしてたりするとPDFで出力した時に文字化けする。 解決策としてはfamilyをJapan1に指定しておく。 こんな感じ。 pdf("hoge.pdf", family=Japan1) 当たり前のことすぎるのかググッても情報を見つけられなかったので忘れないうちにメモしておく。 なお、他の書体（たとえばゴシック）とかも使えるので詳しくはヘルプのpdfFontsを参照のこと。 ※この件、twitterでなにげなくつぶやいたら親切に教えてくださった方がいて感激した。

こんな話がある。 兼子毅（2015）『ゼロから始めるRー四則演算から多変量解析まで』日科技連は、サンプルデータの公開をしていないのかな？？ AKB48の年齢と身長のデータとか、手で全部打つの辛い。。。一発で変換できなさそうな名前も多いし。。。 pic.twitter.com/5NApBIjd5p— Yuichiro Kobayashi (@langstat) 2016年10月17日 こういうのを見たらGoogle Cloud VisionでとりあえずOCRである。 CROWD_VISION_KEY <- "YOUR KEY" library("httr") f <- "https://pbs.twimg.com/media/Cu9cVUlUkAAiHPL.jpg" tmp <- tempfile() download.file(f, tmp) img <- readBin(tmp, "ra

最近readrパッケージが1.0になった。 https://blog.rstudio.org/2016/08/05/readr-1-0-0/ この中で「実験的だが」という断りつきでread_csv_chunked関数が最後に紹介されている。 メモリに乗らないようなデータを分割して読み込むことができる。 baseのread.**関数群も読み込み行数および読み込み開始位置は指定できたので、今までできなかったというわけではないが、１回の読み込みデータ数を指定すればあとはfor文を書かずともよしなにやってくれるので楽ではある。 しかも毎回の読み込み時に呼び出されるコールバック関数も指定できる。 以下は毎回20ずつデータを読み込み、そのたびにstr関数を呼ぶ例。 library("readr") read_csv_chunked(file=readr_example("mtcars.csv"), c

元ネタはこちら http://timelyportfolio.blogspot.jp/2013/07/all-my-roads-lead-back-to-financepimco.html google analyticsでも「ユーザーフロー」に使われているsankey-diagramという可視化手法がある。 http://en.wikipedia.org/wiki/Sankey_diagram で、d3.js関連の可視化をがしがし実装してくれているtimelyportfolioの人がそれも実装していた。 早速試してみる。 library(rCharts) # 事前に↓をzipで落として展開しておく # https://github.com/timelyportfolio/rCharts_d3_sankey # 展開した場所に移動しておく setwd("/hoge/rCharts_d3_s

ひとりアドベントカレンダー18日目。 Advanced RにもあるようにRは遅延評価を導入しているため、思ってたのんと違うわこれという動作がたまにある。 例えばlapply関数の挙動。 > add <- function(x) { + function(y) x + y + } > adders <- lapply(1:10, add) > adders[[1]](10) [1] 20 # ここは11になって欲しい。 > adders[[10]](10) [1] 20 これを防ぐにはforce関数を用いるなどして強制評価する必要があった。 add <- function(x) { force(x) function(y) x + y } adders2 <- lapply(1:10, add) adders2[[1]](10) [1] 11 adders2[[10]](10) [1] 20

ひとりアドベントカレンダー15日目。 Rはよしなにデータを処理してくれるので、非常に使い勝手が良い。 その一方、よしなに処理しすぎてこちらが泣きをみることがある。 たとえば、Rで使わない人はいないであろう［関数。 data.frameクラスに適用される［.data.frame関数は通常データフレームを返す。 > iris[1:2,1:2] Sepal.Length Sepal.Width 1 5.1 3.5 2 4.9 3.0 しかし、返すデータフレームが1列の時はベクトルを返してくる。 > iris[1:2,1] [1] 5.1 4.9 もう慣れすぎて個人的にはこの挙動に対してなんら疑問を持たないのだが、これは一貫性が無い。 一貫性をもたせる、つまりデータフレームで返すためにはdrop=FALSEとする。 > iris[1:2,1,drop=FALSE] Sepal.Length 1 5

分を抽出して10分単位でfloorして、1時間単位でfloorしたデータと足し合わせる。 たまにしかやらないので忘れてしまう。 というか以下の記事でも言及しているようにxtsパッケージとか使えばもっとスマートに書けるんじゃないのか。 http://notchained.hatenablog.com/entry/2015/09/22/072820 library("dplyr") library("lubridate") res <- df %>% group_by(m=floor_date(time, unit = "hour") + minutes(floor(minute(time) / 10) * 10)) %>% summarise(count=sum(hoge)) あと、ミリ秒が入っていたときの書式もよく忘れる。 %Sの代わりに%OSを使えばよい。 as.POSIXctを使えばこ

（rbind_allはdeprecatedとなっており、代わりに現在はbind_rowsが推奨されています） データフレームを縦に結合していく時、皆さんdo.call("rbind", list(data1, data2, data3))みたいなことやっていると思います。 しかしdplyrにC++で書き直されたrbind_allというのがあるのでそちらを使うとめちゃ速くなります。 enjoy! 追記 data.tableに適用する時はdata.table::rbindlistを使えとのことです。 https://github.com/hadley/dplyr/issues/240 追記2 コメント欄に情報いただきました。（musimasamiさんありがとうございます） factorが入っているとうまくいかないこともあるそうなのでcharacterに変換しておくと良いとのことです。 libr

私、集計はplyrでやろうぜええええええええええとか言ってる割に、データマージはmerge()でやっておりました。 もうそろそろ冬休み終わるし沢山の集計作業が待っているが俺は早く帰りたい。 ということでdplyrの*_join()のうち、inner_join()を試してみました。 ３倍速い。 これからはinner_join使います。 なお、本当はplyrのjoinとも比較したかったんだけど、なぜかRが落ちるのでやめました。 追記 ヘルプを読むと、outputのrenameとか変数名が異なるものをbyに指定した結合とかはできないとのことなので注意。 Unlike merge, preserves the order of x no matter what join type is used. If needed, rows from y will be added to the bottom

lubridate（ラブリーデイトって読むのか？）パッケージ使うと、日付の変換やらなんやら直感的に操作できそうなので少し触ってみた。 http://rpubs.com/dichika/lubridate 英語であれば充実した解説が以下に。 http://www.r-statistics.com/2012/03/do-more-with-dates-and-times-in-r-with-lubridate-1-1-0/

俺の睡眠を可視化したい。 で、Sleep CycleというiPhoneアプリを使った。 これは睡眠時の体動をiPhoneの加速度センサで測定するというもの。 Sleep Cycle alarm clock Northcube AB 価格： 100円 posted with sticky on 2015.4.1 ただ残念な事にこれはあくまでアプリ内でしか可視化できない。 runkeeperとも連携しているのでそこからデータを抽出できるかと思ったが何時間寝たとか要約データのみが渡されており俺が欲しいのはそんなんじゃない。 ということでiTunes接続時のバックアップデータを利用してデータを抽出することにした。 バックアップデータから欲しいデータを抽出する 一回抽出するだけだったら以下のアプリ(iPhone Backup Extracter)を動かしてファイル抽出しておしまいで良い。 http:

開発者の@overlastさんからインストール時にユーザ辞書を生成する方法を教えていただいた。 @dichika こんばんは。とても参考になる記事を書いて頂きありがとうございます！！ 記事の内容を踏まえ、ユーザ辞書作成機能を付け、擬似コードを書き換えました。お時間があるときに動きそうな雰囲気かご確認をお願いいたします ^^ https://t.co/khgMUf2k5Y— Toshinori Sato (@overlast) 2015, 3月 27 これでRMeCabからユーザ辞書を指定するだけで良くなった。素敵！！！ まずここにあるようにneologdインストール時にユーザ辞書を生成する。 https://github.com/neologd/mecab-ipadic-neologd/wiki/ProgrammingLanguage.ja なお、ここで書かれているYYYYMMDDは本日の

追記 （３月２８日追記）neologd側でインストール時にユーザ辞書を生成すれば対応できるので別途記事を書いた。そちらを参照のこと。 http://d.hatena.ne.jp/dichika/20150328/p1 経緯 MeCabの新しい辞書が公開された。 https://github.com/neologd/mecab-ipadic-neologd/blob/master/README.ja.md 開発者の方が今週末のTokyoRで発表してくださるらしいが俺は行けない。 https://atnd.org/events/63990 行けないのでRMeCabから使ってみてこの気持ちを昇華させる。 使い方 Mac、Linuxなら難しくない。 インストールは公式の手順に従えば良い。 https://github.com/neologd/mecab-ipadic-neologd/blob/ma

webAPIを叩くパッケージはたくさんあるし、自分でもいくつか作っている。 しかし、毎日使うわけではないし、使いたくなった時には使い方を忘れている。 ということで、統一した記法でデータを取得できるパッケージを作った。 名前は蛇口(jaguchi)。 devtools::install_github("dichika/jaguchi") library("jaguchi") # speakerdeckに上がったスライドの情報を取得 jaguchi("speakerdeck", url="https://speakerdeck.com/dichika/tokyor-extra-number-1") # ニコニコ動画にクエリを投げて情報を取得 jaguchi("niconicos", query="艦これ", size=100, type="tag") 位置づけ 分析フローを下記のように分けた時

GETで持って来れるので別に難しい話ではない。 が、また使いたくなったときに調べるのが面倒なのでパッケージにした。 install_github("dichika/hatenab") library(hatenab) # もっとも古い日付と最新の日付からURLを作成 dates <- seq(as.Date("2010-02-10"),as.Date("2015-01-22"), by=1) dates <- format(dates, "%Y%m%d") urls <- sprintf("http://d.hatena.ne.jp/dichika/%s/p1", dates) # はてブAPIを叩く result <- getBookmarkCount(urls) result <- result[order(result$count, decreasing=TRUE),] head(r