理学療法士がまとめたノート

統計学備忘録(R言語のメモ)

since2016 ときどきTEXのメモ

グラフ

ヒストグラム

library(ggplot2) data <- c(rnorm(n=100, mean=5, sd=1),rnorm(n=100, mean=3, sd=3)) cat <- c(rep("x",100), rep("y",100)) dat1 <- data.frame(cat, data) colnames(dat1) <- c("name","value") #乱数なのでデータは変わります 単純なヒストグラム HG1= …

対応のあるデータのグラフ

id <- 1:5 x <- c(rnorm(5, 20, 2)) y <- c(rnorm(5, 25, 1)) data <- data.frame(id, x, y) id x y 1 1 22.83957 23.13902 2 2 15.05851 23.15993 3 3 19.49754 24.70116 4 4 17.43654 26.60801 5 5 17.01846 25.09940 x <- c(1,2) t_data <- t(data[,2:3])…

回帰直線と交互作用

パッケージ無しで描く方法 dat <- data.frame(y=c(1,2,3,4), x=c(4,5,8,12), a=c(1,1,0,0)) dat$a <- as.factor(dat$a) plot(dat$x, dat$y) fit <- summary(lm(y~x, data=dat)) lines(range(dat$x), fit$coef[1]+ fit$coef[2]*range(dat$x)) pchAB <- ifelse…

混合モデル(マルチレベルモデル)の基礎

投稿日:2021.5.2 忘れないように基本的な部分のみ貼っておきますが、名称だけでも統一していただきたいです・・・ 線形混合モデル () 線形混合効果モデル () 階層線形モデル () マルチレベルモデル ()・・・etc サンプルはパッケージlme4の「cbpp」を使用し…

stripchartの横にエラーバー

サンプル ID <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20) 治療 <- c("A","B","A","B","B","A","A","A","B","B","A","B","A","A","A","B","A","B","B","B") BP変化量 <- c(治療後BP - 治療前BP) dat1 <- data.frame(id, 治療, 治療前BP, 治療…

2標本のグラフ、記述統計

2標本のまとめを簡単に サンプル dat1 <- c(-0.65, -0.58, 1.98, 0.98, 1.55, -0.05, 0.10, -1.35, 1.24, -0.28) dat2 <- c(1.31, 2.39, 2.87, 2.25, 2.33, 5.32, 3.13, 1.21) まずはデータセットを作成 data <- c(dat1, dat2) group <- c(rep("A", length(d…

感度と偽陽性 ( 単回帰分析より )

X <- c(rep(0,70), rep(1,75)) #検査結果 0=陰性、1=陽性 Y <- c(rep(0, 45),rep(1, 25), rep(0, 8), rep(1, 67)) #罹患 0=無し、1=あり (b <- xtabs(~Y + X)) X Y 0 1 0 45 8 1 25 67 上記分割表より 感度 P(Y=1, X=1):67/(25+67) = 0.7282609 1-特異度 P…

エラーバー付きのグラフ ggplot2

更新日2020.3.13 使用するパッケージ install.packages("ggplot2") install.packages("ggsignif") データセットを作成します dat <- c(-0.29733004,-2.63812280,-0.90097072,1.06843016,-3.03846846,2.14097694,2.47494865,-0.02154341,0.56411223,2.8182606…

簡単なエラーバーの描き方

たぶんこれが最も簡単な書き方だと思います Aの平均 = 2 Bの平均 = 3 Aの標準偏差 = 0.4 #ここに標本標準偏差や不偏標準偏差を代入します Bの標準偏差 = 0.7 dplot <- plot(c(1,2), c(Aの平均, Bの平均), ylim = c(0,5), xlim = c(0.5, 2.5)) arrows(1:2, c(…

比較した箱ひげにアスタリスクを入れる

アスタリスク 比較したラインを箱ひげの上限から2メモリ離した位置まで伸ばす・・・ 治療 <- c("A","A","A","A","A","A","A","A","A","A","B","B","B","B","B","B","B","B","B","B") 治療前BP <- c(140,135,145,141,145,159,141,150,149,160,170,158,165,170…

層別した散布図と回帰直線

以下のようなグラフ 忘れないように記載しておきます データサンプル Group <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3) p…

独立している2群を比較するグラフ ggplot2

データの準備 dat <- c(c(rnorm(15)),c(rnorm(15,3,1.2)),c(rnorm(15,10,2)),c(rnorm(15,17,1.5))) pre_post <- c(rep("前", 30),rep("後", 30)) treat <- c(rep("A",15),rep("B",15),rep("A",15),rep("B",15)) dataf <- data.frame(dat, pre_post, treat) h…

散布図に直線を描く

x <- rnorm(20) #乱数を使用してるので毎回変わります y <- c(1:20) plot(x,y) xs <- seq(min(x), max(x), length=1000) # xの影響を含むモデル fit <- glm(y~x,family = poisson) lines(xs, exp(fit$coef[1] + xs*fit$coef[2])) # 切片のみのモデル fit.nul…

正規分布の色塗り (polygon関数 )

投稿日:2018.7.19 最終更新日:2019.5.18 polygon関数 例)y = 2*x (x,y)=(0,0), (2,0), (2,4)の△を塗りつぶす場合 x <- 0:5 y <- 2*x plot(x,y) polygon(c(0,1,2,2,1,0),c(0,0,0,4,2,0),col="gray") # c(x軸, y軸, 色) #三角形なので、のように書いても同じ…

個人的によく使うグラフ

サンプル a <- iris$Sepal.Length[51:75] b <- iris$Petal.Length[101:125] c <- a+b data <- data.frame(a,b,c) summary(data) a b c Min. :4.900 Min. :4.50 Min. :10.10 1st Qu.:5.600 1st Qu.:5.10 1st Qu.:11.10 Median :6.100 Median :5.60 Median :11…

級内相関係数 Case1

以前掲載していたのですが、以下の文献を参考に書き直しました。 主にRのプログラム (パッケージは使用してません) を忘れないようにメモしておきます。 理論は以下の文献を参照してください。 1) SHROUT, Patrick E.; FLEISS, Joseph L. Intraclass correla…

共分散分析を考えるためのグラフ

いつも勉強させていただいおります、下記サイトの共分散分析より グラフの描き方の備忘録です 我楽多頓陳館--雑学と統計学の館 以下のようなグラフになります データのサンプルは上記サイトから引用してdataに入れました # まずはこれを描いて考える・・・!…

ヒストグラムのビン調整

head(iris) Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa 変数Speciesがsetosa…

相関係数のイメージ

パッケージmvtnormを使用して相関係数0.0, 0.2, 0.5, 0.7, 0.8, 0.9のグラフを作成してみます install.packages("mvtnorm") library(mvtnorm) 共分散行列.分散を全て1に設定しているので共分散=相関係数となります. sigma00 <- matrix(c(1,0,0,1), ncol=…

2変量の正規分布をグラフでイメージ(persp)

また、ここで勉強させていただきました. http://cse.naro.affrc.go.jp/minaka/R/R-binormal.html 忘れないように要点のみ転記させていただます.まさに備忘録. 今回はRの関数perspを使用して、密度関数の数式から3Dのグラフを描いてみます 確率変数x1…

2変量の正規分布をグラフでイメージ(scatterplot3d)

ここで勉強させていただきました. http://cse.naro.affrc.go.jp/minaka/R/R-binormal.html 忘れないように要点のみ転記させていただます. 必要なパッケージをインストールします install.packages("mvtnorm") library(mvtnorm) install.packages("scatterp…

正規分布の重ね描き

text関数でグラフに文字の挿入 curve(dnorm(x, -2, 4), from=-10, to=10, ylim=c(0,0.4),ylab ="") text(-5, 0.1, "N(-2,4)") par(new=T) curve(dnorm(x,3, 1), from=-10, to=10, ylim=c(0,0.4),ylab ="") text(1.5, 0.3, "N(3,1)") par(new=T) curve(dnorm(…

信頼区間のプロット

同じサイズのデータサンプルからt分布を利用した信頼区間の作図 まずは3×4の場合(サンプルサイズ3を4回実施する) x <- matrix(NA,nrow=3,ncol=4) #3×4の空セル for (i in 1:4){ #列数分乱数を代入 x[,i] <- rnorm(3) #標準正規分布の乱数を行数分繰り返…

クラスター分析

Rのサンプルattitudeを使用して、クラスタ分析(階層的方法)を勉強します. attitudeは、管理者態度のデータです.無作為に選ばれた35名の雇用者よりアンケート.好意的な割合が数値化されています. rating全般的評価、 complaints雇用者からの苦情処理、 …

指数分布

指数分布 投稿日2017.9.4 今回はRを使って苦手な指数関数に挑戦します定義 を確率変数の確率密度関数とします は確率変数の累積分布関数とします 離散型の変数で確率密度関数の復習をします 例)正確なサイコロ 確率変数 y<-c(rep(1/6,6)) plot(y, type = "h…

残差分析

残差分析 Rのサンプル ChickWeight を使用y<-ChickWeight$weight[1:10]x<-ChickWeight$Time[1:10]summary(lm(y~x))plot(x,y,xlab = "生後日数",ylab = "体重")abline(30.327 ,7.030) #回帰直線の挿入 残差 ei をプロットすることで、回帰モデルからのズレを…

F分布の確率密度関数

F分布の確率密度関数df(x, df1, df2, ncp, log = FALSE) 自由度1、自由度5のF分布の確率密度関数のグラフを描いてみます # par()を使用する# グラフィックパラメータ:グラフを構成する点や線の形状や色、大きさや余白などを決定する# par()は余白などの土台…

並行箱ひげ図

エクセルでデータのみコピーt<-read.table("clipboard")x<-t$V1y<-t$V2boxplot(x,y,xaxt="n",main="タイトル") #x軸のラベルを消去name<-c("ラベル1","ラベル2") #x軸のラベルを定義axis(side=1,at=c(1,2),labels=name) #x軸のラベルを更新 このような小技が…

t分布からの信頼区間

自由度5のt分布より確率点:qt(0.025,5,lower.tail = F)累積分布:pt(2.570582,5,lower.tail = F)確率密度:dt(2.570582,5) 拡張期血圧を6回測定して95%信頼区間を求めるDBP<-c(86,92,88,94,89,88)vDBP<-var(DBP) t値 { mean(DBP)-X } ÷ {sqrt(vDBP/6)} 95%…

データの要約

irisx<-iris$Sepal.Length[1:50]y<-iris$Sepal.Length[51:100]xy<-data.frame(x,y)summary(xy)boxplot(x,y)t<-c(1,2)ave<-c(mean(x),mean(y)) #平均値points(t,ave,pch=16) #平均値(・)を挿入 x y Min. :4.300 Min. :4.900 1st Qu.:4.800 1st Qu.:5.600 Me…