My notepad

Diff動作解析チュートリアル（DIFFマーカーセットで測定したc3dデータを使って動作分析を行う）

2026-10-07T19:32:31.098+09:00

DIFF 動作解析（Web 版）は、三次元動作分析装置（c3d）と床反力計のデータから、DIFF モデルに基づいて関節角度・身体重心・関節モーメント・関節パワーを計算する Web アプリです。この記事では使い方を説明します。

アプリを開く

目次：はじめに／準備するデータ／解析の手順／ビューアの使い方／出力ファイル／動作ごとの解析／主な設定／困ったときは

はじめに

データは外部に送信されません。 計算はすべてお使いのブラウザの中で行われます。c3d ファイルやその内容がサーバーに送られることはありません。

推奨ブラウザ：Google Chrome、Microsoft Edge の最新版（パソコン）。
初回の読み込み：計算環境（約 30 MB）を読み込むため、数十秒かかります。インターネット接続が必要です。2 回目以降はブラウザに保存されるので速くなります。
解析の設定は、このブラウザに保存されます（c3d などのデータは保存されません）。

起動直後の画面。右上に「準備ができました」と表示されたら使えます。

準備するデータ

データ	内容
動作データ（必須）	解析する試行の c3d（マーカー座標と床反力計のアナログデータ）。
静止立位データ（推奨）	DIFF 直立肢位（両足を肩幅に開き、つま先開き角 0°）で計測した c3d。体重の推定、足部座標系の補正、足関節角度の 0° の基準などに使います。
被験者の計測値	性別、身長、膝関節左右径・足関節左右径（ノギスなどで計測）。体重は静止立位の床反力から推定できます。

マーカーセットは DIFF（検証データはViconで計測）を前提としています：左右の肩峰（SHD）、股関節（HIP）、膝関節（KNE）、外果（ANK）、第 5 中足骨頭（MP）、上腕骨外側上顆（ELB）、尺骨茎状突起（WRA）。名前が異なる場合は「DIFF モデルの設定」で変更できます。

膝関節・足関節左右径を入力しないと、その関節点はマーカー位置（DIFF 付録 A1 の第 1 近似）になり、とくに前額面の角度・モーメントの精度が下がります。

解析の手順

左側の設定欄

動作データの「ファイルを選択」で、解析する c3d を選びます。
静止立位データに、同じ被験者・同じ日の静止立位の c3d を選びます。
被験者情報で、性別・身長・膝関節左右径・足関節左右径を入力します。体重は 0 のままで静止立位から推定されます。
解析条件の「動作」で、歩行・立ち上がり・しゃがみ込みなどを選びます。動作に応じてイベントや相が検出されます。
必要に応じてフィルターや DIFF モデルの設定を変えます（通常は初期値のままで構いません）。
一番下の「解析する」を押します。数秒〜十数秒で結果が表示されます。

解析結果。上から、注意事項、ダウンロードボタン、相と所要時間（立ち上がりの場合）、床反力計の判定結果、ビューアの順に並びます。

結果で確認すること

注意事項（黄色）：左右径が未入力、足が床反力計の外で接地した、など。計算しなかった区間の理由もここに出ます。
床反力計の判定結果：どの床反力計にどの足（または椅子など）が乗ったかを自動で判定しています。違っている場合は、左の「床反力計の割り当て」で手動に変えて、もう一度「解析する」を押します。

ビューアの使い方

結果の下にあるビューアで、スティックピクチャーと各指標のグラフを同期させて確認できます。ビューアは「ビューア（.html）をダウンロード」で保存すると、このアプリがなくても開けます。

ビューア。左に矢状面・前額面、右に 3D のスティックピクチャーとグラフ。

再生と描画区間

▶ で再生、Space キーで再生・停止、← → で 1 フレームずつ送ります。速度は 0.1〜1 倍。

描画区間：開始・終了時刻を入力するか、「現在→開始」「現在→終了」で決めます。「相・周期から選ぶ」で、立ち上がりの各相や歩行周期をそのまま区間にできます。スティックピクチャー・グラフ・床面投影図がこの区間に限定されます。
時間軸の上には、相（帯）とイベント（離殿・踵接地など）が表示されます。
「被験者情報」を押すと、身長・体重などが表示されます。

3D スティックピクチャー

操作	方法
回転	左ドラッグ
移動	右ドラッグ、または `Shift`＋左ドラッグ
拡大・縮小	マウスホイール、または「＋」「−」
視点の切り替え	「視点」から斜め・右側方・左側方・前方・後方・上方
元に戻す	「リセット」またはダブルクリック

マーカーにポインタを合わせると、名前と座標が表示されます。「マーカー名」「床反力ベクトル（右・左）」「床反力ベクトル（左右合成・破線）」「重心の軌跡」などはチェックボックスで切り替えます。左右合成の床反力ベクトルは、左右の足の床反力を合成したもので、左右合成の作用点から白の破線で描きます（矢状面・前額面・3D）。

床面投影図

上から見た身体重心と床反力作用点（右・左・左右合成）の軌跡です。ホイールや「＋」「−」で拡大・縮小、ドラッグで移動、「全体」で元に戻ります。点をクリックすると、その時刻へ移動します。

時系列グラフ

青が右、橙が左、緑が全身・左右合成です。グラフをクリック・ドラッグするとその時刻へ移動し、ポインタを合わせると値が表示されます。表示する指標は「表示する指標を選ぶ」で選びます。

出力ファイル

ファイル	内容
Excel（.xlsx）	関節角度、体節絶対角度、重心、床反力、各床反力計、関節モーメント、関節パワー、関節点座標、マーカー座標、イベント、解析情報（設定の記録）。動作によって「相・所要時間」「歩行パラメータ」「歩行周期_右／左」も入ります。
縦長 CSV（R 向け）	1 行 1 値の tidy 形式（trial, frame, time_s, category, variable, side, unit, value）。
歩行周期 CSV	歩行で、1 歩行周期を 0〜100 % に正規化した値（各周期）。
ビューア（.html）	スティックピクチャーとグラフ。単独で開けるので、共同研究者への共有にも使えます。

R での読み込み例：

library(tidyverse)
d <- read_csv("sample02_DIFF_long.csv")
d |>
  filter(category == "関節モーメント", unit == "N·m/kg",
         str_detect(variable, "^膝関節モーメント (伸展|屈曲)")) |>
  ggplot(aes(time_s, value, colour = side)) + geom_line()

動作ごとの解析

動作	検出されるイベント・相
立ち上がり	動作開始（体幹前傾の開始）、離殿、足関節最大背屈、股関節伸展終了。第 1 相（屈曲相）・第 2 相（重心移動相）・第 3 相（伸展相）の所要時間。
しゃがみ込み	動き始め（重心の下降開始）、着臀、座位到達。下降相・着座相の所要時間。
歩行	踵接地・つま先離地（床反力計の接地が近ければそちらを採用）。歩行パラメータ（立脚期 %、ケイデンス、重複歩距離、歩行速度、歩幅、歩隔）。グラフの横軸を「歩行周期 %」に切り替えられます。

離殿・着臀は、椅子の下の床反力計で判定します。椅子の下に床反力計がない場合は、股関節点の高さの変化から推定します。

主な設定

設定	説明（初期値）
カットオフ周波数	バターワース低域通過フィルター（6 Hz）。初期設定では床反力にもマーカーと同じ周波数を使い、すべての値をこの周波数で求めます。
床反力のカットオフ（用途別）	「床反力にもマーカーと同じカットオフ周波数を使う」のチェックを外すと、「床反力のカットオフ周波数」（関節モーメント・パワーの計算に使用）と「床反力・作用点の出力用」を別々に指定できます（周波数の違いによる影響の検証用）。
床反力の閾値	これ未満を無負荷とする鉛直力（20 N）。
関節点の近似法	DIFF 付録 A1。股関節 3B（左右マーカー間距離の 18 % 内挿）、膝・足関節 3A（左右径の 1/2 内挿）。
静止立位の足関節角度を 0° とする	マーカーの貼り方の影響を除くため、静止立位の背屈・底屈角を差し引きます（オン）。
関節モーメントの符号	初期値は股関節伸展＋・外転＋、膝関節伸展＋、足関節底屈＋。「DIFF 解説書の表記」を選ぶと膝関節のみ屈曲＋になります。
関節モーメントの座標系	遠位体節座標系（初期値）、近位体節座標系、実験室座標系。

計算方法の詳細（座標系、角度の定義、慣性係数、逆動力学など）は、リポジトリの計算仕様書 docs/methods.md にまとめています。論文の方法の記載に使ってください。

困ったときは

「準備しています…」から進まない

初回は計算環境の読み込みに数十秒かかります。数分たっても進まない場合は、インターネット接続を確認してください。所属機関のネットワークで外部サイト（cdn.jsdelivr.net）が遮断されていると読み込めません。

「マーカー ○○ が見つかりません」と表示される

c3d のマーカー名が DIFF の名前と違います。表示された「c3d にあるマーカー」を参考に、「DIFF モデルの設定」のマーカー名を合わせてください（左右の接頭辞 R/L を除いた部分を入力します）。

「体重が不明です」と表示される

体重を入力するか、床反力計に乗った静止立位データを指定してください。

関節モーメントが一部の区間で空欄になる

足が床反力計の外で接地している区間、両足が同じ床反力計に乗っている区間、椅子などが床反力計に乗っている区間（股関節のみ）は、正しく計算できないため空欄にしています。注意事項に理由が表示されます。

床反力計の判定が違う

左の「床反力計の割り当て」で、床反力計ごとに右足・左足・その他（椅子など）・無視を選び、もう一度「解析する」を押してください。

Plot Configurator User Manual

2026-01-12T23:36:00.003+09:00

Plot Configurator User Manual

Welcome to Plot Configurator!
This application(https://plot-app-seven.vercel.app/) is a tool that allows you to easily create high-quality charts (Line Charts, Bar Charts, Scatter Plots) suitable for papers and presentations simply by inputting data.

Introduction
Quick Start
Detailed Usage
Tips

Introduction

In this app, you can instantly render charts just by copying and pasting data from Excel or spreadsheets. You can finley tune colors, fonts, layouts, and save the result as an image.

Quick Start

Enter Data: Paste your data into the "Data Input" area on the left.
Select Chart: Choose from "Line", or "Bar", "Scatter".
Adjust: Use the panels below to adjust colors and axis ranges.
Save: Click the "Download" button below the chart to save the image.

Detailed Usage

1. Data Entry

Paste your input data into the text box on the left. It supports copy & paste from Excel or Google Spreadsheets.
The data format differs slightly depending on the type of chart you want to create.

Common Rules:

The first line is recognized as the "Header (Item Names)" (e.g., x, y, group).
Data supports tab-delimited (copy-paste from Excel) or comma-delimited (CSV) formats.

A. For Line / Bar Chart

Columns are required in the following order:
X Value | Y Value | Error (Optional) | Group (Optional)

Input Example:

x	y	error	group
1	10	0.5	Control
2	20	1.0	Control
1	15	0.8	Treatment

B. For Scatter Plot

Columns are required in the following order:
X Value | Y Value | X Error (Optional) | Y Error (Optional) | Group (Optional)

Input Example:

x	y	x_error	y_error	group
1	10	0.5	0.2	A
2	15	0.8	0.3	A
1	14	0.5	0.2	B

Note: If the data format is incorrect, the chart may not appear. In that case, please refer to the format of the sample data.

2. Selecting Chart Type

You can change the chart type using the selector at the top left of the screen.

Line: Suitable for viewing time-series data or trends.
Bar: Suitable for comparing categories.
Scatter: Suitable for visualizing the relationship between two variables. You can set both X and Y errors.

3. Customizing Appearance

You can adjust the design of the chart in the "Appearance" section.

Aspect Ratio: Choose the aspect ratio from 16:9, 4:3, or Custom (fit to browser width).
Layout: Choose whether to orient the chart Horizontally or Vertically.
Margins: Adjust the top, bottom, left, and right margins in pixels.
Color Palette: Batch change the chart colors from presets (Default, Nature, Vibrant, etc.).
Font: Select Japanese fonts (Hiragino/Meiryo, etc.) or Serif English fonts. Font size can also be changed.
Grid and Legend: Toggle the display of grid lines and the legend. The legend can be dragged and dropped to your preferred position.

4. Axis Configuration

In the "Axes Config" section, you can configure detailed settings for the X and Y axes.

Show: Whether to display the axis itself.
Title: The name of the axis (e.g., "Time (s)", "Concentration (%)") and its position (Center, Left/Top, Right/Bottom).
Min / Max: Manually fix the axis range. Leave blank to automatically adjust based on the data (Auto).
Tick Step: Specify the interval of the tick marks.

5. Saving Images

You can download the created chart from the bar below the preview screen.

File Name: Enter a file name for the image to be saved.
Save: Click the button to download the image in PNG format, etc.

Tips

Reordering Groups: If your data has multiple groups (e.g., Control, Treatment), you can change the display order or customize individual colors and names in the "Series Settings" at the bottom of the settings panel.
Jitter: If data points overlap in a Line chart making it hard to see, turning on "Jitter" will shift the points slightly to prevent overlap.
Language Switching: You can switch between English and Japanese using the button at the top left of the screen.

Created by: SK

Plot Configuratorの使い方

2026-01-11T21:35:00.005+09:00

Plot Configurator 取扱説明書

Plot Configuratorへようこそ！
このアプリケーション(https://plot-app-seven.vercel.app/)は、データを入力するだけで、論文やプレゼンテーションに使用できる高品質なグラフ（折れ線グラフ、棒グラフ、散布図）を簡単に作成できるツールです。

はじめに

このアプリでは、Excelやスプレッドシートからデータをコピー＆ペーストするだけで、即座にグラフを描画できます。色、フォント、レイアウトなどを細かく調整し、画像として保存することができます。

クイックスタート

データを入力: 左側の「データ入力」エリアにデータを貼り付けます。
グラフを選択: 「折れ線グラフ」「棒グラフ」「散布図」から種類を選びます。
調整する: 下のパネルで色や軸の範囲を調整します。
保存する: グラフ下の「ダウンロード」ボタンで画像を保存します。

詳細な使い方

1. データの入力

左側のテキストボックスに入力データを貼り付けます。ExcelやGoogleスプレッドシートからのコピー＆ペーストに対応しています。
データの形式は、作成したいグラフの種類によって少し異なります。

共通のルール:

1行目は「ヘッダー（項目名）」として認識されます（x, y, groupなど）。
データはタブ区切り（Excelからのコピペ）またはカンマ区切り（CSV）に対応しています。

A. 折れ線グラフ・棒グラフ (Line / Bar) の場合

以下の順序で列が必要です：
Xの値 | Yの値 | 誤差(省略可) | グループ(省略可)

入力例:

x	y	error	group
1	10	0.5	Control
2	20	1.0	Control
1	15	0.8	Treatment

B. 散布図 (Scatter Plot) の場合

以下の順序で列が必要です：
Xの値 | Yの値 | X方向の誤差(省略可) | Y方向の誤差(省略可) | グループ(省略可)

入力例:

x	y	x_error	y_error	group
1	10	0.5	0.2	A
2	15	0.8	0.3	A
1	14	0.5	0.2	B

ポイント: データ形式が正しくない場合、グラフが表示されないことがあります。その場合は、サンプルデータの形式を参考にしてください。

2. グラフの種類の選択

画面左上のセレクターでグラフの種類を変更できます。

Line (折れ線グラフ): 時系列データや変化を見るのに適しています。
Bar (棒グラフ): カテゴリごとの比較に適しています。
Scatter (散布図): 2つの変数の関係を見るのに適しています。X誤差とY誤差の両方を設定できます。

3. 見た目のカスタマイズ

「外観設定 (Appearance)」セクションでグラフのデザインを調整できます。

アスペクト比: グラフの縦横比を 16:9、4:3、あるいは Custom（ブラウザの幅に合わせる）から選べます。
レイアウト: グラフを横向き（Horizontal）にするか縦向き（Vertical）にするか選べます。
余白 (Margins): 上下左右の余白をピクセル単位で調整できます。
カラーパレット: グラフの色使いをプリセット（Default, Nature, Vibrantなど）から一括変更できます。
フォント: 日本語フォント（Hiragino/Meiryoなど）やSerif英語フォントなどを選択できます。文字サイズも変更可能です。
グリッドと凡例: グリッド線や凡例（Legend）の表示・非表示を切り替えられます。凡例はドラッグ＆ドロップで好きな位置に動かすことができます。

4. 軸の設定

「軸の設定 (Axes Config)」セクションで、X軸・Y軸ごとの詳細設定が可能です。

表示: 軸そのものを表示するかどうか。
タイトル: 軸の名前（例：「時間 (s)」「濃度 (%)」）と、その位置（中央、左/上、右/下）。
最小値 / 最大値: 軸の範囲を手動で固定できます。空白にするとデータに合わせて自動調整されます（Auto）。
目盛りの間隔 (Tick Step): 目盛りの刻み幅を指定できます。

5. 画像の保存

プレビュー画面の下にあるバーから、作成したグラフをダウンロードできます。

ファイル名: 保存する画像のファイル名を入力します。
保存: ボタンをクリックすると、PNG形式などで画像がダウンロードされます。

ヒント

グループの並び替え: データに複数のグループ（例：Control, Treatment）がある場合、設定パネルの下部にある「シリーズ設定」で、表示順序を入れ替えたり、個別の色や名前を変更したりできます。
ジッター (Jitter): 折れ線グラフでデータ点が重なって見にくい場合、「Jitter (少しずらす)」をオンにすると、点が重ならないように少しずらして表示されます。
言語切り替え: 画面左上のボタンで、英語(English)と日本語(Japanese)を切り替えることができます。

作成: SK

重回帰分析（モデル選択を含む）

2024-03-16T00:06:00.001+09:00

ある患者群の主観的幸福感（QOL；VASで取得）、栄養状態（MNA）、日中の身体活動（Mets）、痛み（VASで取得）を想定したサンプルデータです。ここからQOLが、栄養状態、身体活動、痛みで説明できるかを検証します。

QOL	mna	mets	pain
83	24	2.4	36
84	26	2.2	40
中略	中略	中略	中略
79	24	2.2	40
93	26	2.5	55

コードと結果は以下の通りです。

> model1<-lm(qol ~ mna + mets + pain, data = data1)
> model1 |> summary()

Call:
lm(formula = qol ~ mna + mets + pain, data = data1)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2957 -1.5591  0.0094  1.5762  5.7335 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 30.95741    4.67349   6.624 1.44e-08 ***
mna          0.96400    0.12507   7.708 2.34e-10 ***
mets         6.53627    1.83639   3.559 0.000766 ***
pain         0.34312    0.03929   8.734 4.87e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.757 on 56 degrees of freedom
Multiple R-squared:  0.7574,	Adjusted R-squared:  0.7444 
F-statistic: 58.28 on 3 and 56 DF,  p-value: < 2.2e-16

重回帰のモデルにmodel1という名前を付けてコードを記載しました。lm()関数を使うことで、単回帰だけでなく重回帰のモデルを作成することができます。今回はQOLが目的変数で、MNA、Mets、painが説明変数です。目的変数と説明変数は ~ で結び、説明変数は+で加えていき、モデル式をつくります。今回はqol ~ mna + mets + painとなります。その後、 , で区切ってデータの名前を指定します。そして、model1をパイプ演算子でsummary()関数とつなぎます。

　結果です。Coefficients:を確認します。Estimateは回帰係数、Std. Errorは標準誤差、t valueがt値で、Prがp値です。ただし、InterceptのEstimateは回帰係数ではなく、切片の値です。ここから以下のようなモデル式が導かれます。

qol = 0.97329×mna + 16.90631×mets + 0.37702×pain + 5.35034

　また、5%を有意水準と定めると、すべての回帰係数は有意です。ここから、回帰係数は０ではないと結論づけられ、すべての説明変数は、QOLに正の影響を及ぼしていると言えそうです（もし係数が負の値なら負の影響を及ぼすことになります）。注意が必要ですが、回帰係数だけみると約16.9のmetsの影響が抜き出て大きいようにも見えます。Metsの値は、0.1刻みで変化する値で、１変化することは考えにくいと言えます。一方で、VASで測定している痛みの測定値は、最大で100となるため1の変化はあまり大きいと言えません。このように、説明変数の値の意味が異なるため、単純な比較はできないことに注意しましょう。どの回帰係数がどれだけの影響があるのか、を考えるときには、その説明変数の値が持っている意味をドメイン知識に照らし、判断していくのが大切です。

　また、モデル式のデータに対する当てはまりをチェックするための決定係数は、Multiple R-squaredで、Adjusted R-squaredは調整済みの決定係数になります。調整済み決定係数は0.7444と非常に高い値となっています。

　多重共線性を確認します。各説明変数の相関をチェックすることも行われますが、ここではVIFを確認します。carパッケージのvif()関数が必要になります。carパッケージは、rstatixパッケージをインストールした際に同時にインストールされているパッケージです。ですので改めてインストールする必要はありません。コードと結果は以下のようになります。

> model1 |> vif()
     mna     mets     pain 
1.094704 1.097706 1.003791

library()関数でcarパッケージを読み込みます。それからmodel1とvif()関数つなぐ、もしくはvif()関数の()内にモデル名を指定します。結果は10よりも小さくすべての説明変数を投入しても問題なさそうです。

　回帰モデルのモデル選択について解説します。複数のモデルから予測性能の高いモデルを選択するために赤池情報量基準を用います。model1以外のモデルを作成します。mna、mets、painを一つずつ減らした説明変数が2つだけのモデルを作りました。コードは以下の通りで、全部で4つになります。

> model2<-lm(qol ~ mna + mets, data = data1)
> model3<-lm(qol ~ mna + pain, data = data1)
> model4<-lm(qol ~ mets + pain, data = data1)
> 
> model1 |> AIC()
[1] 297.8261
> model2 |> AIC()
[1] 347.4006
> model3 |> AIC()
[1] 308.0626
> model4 |> AIC()
[1] 339.2155

作成したモデルに対し、それぞれAIC()関数を実行します。()内にモデル名を入力しても構いません。model1のAICが最小となったため、model1を採用するのが良い、という結果になりました。

おまけです。各説明変数間での目的変数に対する影響を分析したい場合、目的変数・説明変数のすべてを平均０，標準偏差１に標準化し、重回帰分析を行うことで、比較ができます。Rではscale()関数で簡単に標準化したデータに変換できます。標準化したデータに対し重回帰分析を行うことで、説明変数間の影響の大きさを比較できるようになります。

単回帰分析

2024-03-15T23:52:00.001+09:00

回帰分析にはrstatixパッケージを使いません。10m歩行の方向速度(gait)を目的変数、1日活動量（1日の歩数; steps）を説明変数にした単回帰分析を行います。

gait	steps
1	7938
0.96	6172
1.03	9653
中略	中略
0.96	8783
1.48	8861
1.42	9292

回帰分析のコードは以下の通りです。

> model1<-lm(gait ~ steps, data = data1)
> model1 |> summary()

Call:
lm(formula = gait ~ steps, data = data1)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.44208 -0.23944  0.00866  0.19000  0.47315 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.659e-01  2.306e-01   1.153 0.253635    
steps       1.137e-04  2.815e-05   4.041 0.000159 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2543 on 58 degrees of freedom
Multiple R-squared:  0.2197,	Adjusted R-squared:  0.2062 
F-statistic: 16.33 on 1 and 58 DF,  p-value: 0.0001589

model1という名前を付けてコードを記載しました。lm()関数を使い、モデルの構造と対象データを指定します。目的変数と説明変数は ~ で結びます。すなわちgait ~ stepsとなります。そして、, で区切ってデータの名前を指定します。

model1の結果をsummary()関数で出力します。

model1をパイプ演算子でsummary()関数とつなぎます。summary()関数を用いることで、モデルの詳細を表示できます。summary(model1)としても同様にモデルの詳細を表示できます。

結果です。Coefficients:を確認します。Estimateは回帰係数、Std. Errorは標準誤差、t valueがt値で、Prがp値です。ただし、InterceptのEstimateは回帰係数ではなく、切片の値です。これらの結果から以下のモデル式が導かれます。

gait(10m歩行速度) = 0.0001137 × steps（1日の歩数） + 0.2659

ここから、1歩歩数がふえると0.0001137m/秒歩行速度が速くなることが理解できます。つまり、1000歩増えると、0.137m/秒速くなることわけです。10m歩行の平均的な速度が1.4m/秒であることを考慮するとそれなりに影響する、ということが言えそうです。また、有意水準を5%と設定した場合、回帰係数は有意ということになります。これは、回帰係数=0という、帰無仮説は棄却されたことを意味します。調整済み決定係数はおよそ0.2でした。そしてf統計量に対して、p値は有意水準未満となっていました。モデルとしても有意と言えそうです。

対応のないt検定（rstatixパッケージ）

2024-03-15T22:23:00.005+09:00

rstatixパッケージを用いて解析を行います。

以下にサンプルデータを示します。このデータは1列目にグループ、2列目に値をとっています。valuesは握力とします。この2つのグループの握力の平均値に有意な差が有るのかを対応のないt検定（Welch）で検証してみます。

groups	values
a	32.1
a	29.5
a	34.8
中略	中略
b	35.2
b	36.4
b	35.1
中略	中略

上記のデータの名前はdata1としています。

rstatixパッケージを読み込んでt検定のコードを入力していきます。

> library(rstatix)

 次のパッケージを付け加えます: ‘rstatix’ 

 以下のオブジェクトは ‘package:stats’ からマスクされています:

    filter

> t_test(data = data1, values ~ groups, paired = FALSE)
# A tibble: 1 × 8
  .y.    group1 group2    n1    n2 statistic    df         p
* <chr>  <chr>  <chr>  <int> <int>     <dbl> <dbl>     <dbl>
1 values a      b         20    20     -5.35  25.9 0.0000133

このコードでは、t_test関数を( )内の内容にしたがって実行しています。( )内では、データの指定と対応の有り無しの指定を行います。最初にデータの指定です。ここでは、data1です。valuesとgroupsを ~ で結んでおり、valuesの値についてgroupsで分けて分析するという意味になります。対応の有無は、paired=で指定します。TRUEとすれば対応あり、FALSEとすれば対応なし、となります。なお、t_test関数の対応のないt検定では、指定しない限り等分散を仮定しないWelchのt検定が実行されます。もし、等分散を仮定したStudentのt検定を実行したい場合は、var.equal = TRUEを( )内に加えると実行できます。

n1とn2は各グループのサンプルサイズ、statisticはt値、dfが自由度、pがp値になります。有意水準を5%と設定した場合、帰無仮説は棄却されるため、2つのグループに有意な差があると結論づけることができます。

順位和検定と符号付順位和検定（rsatixパッケージ）

2024-03-15T22:22:00.005+09:00

rstatixパッケージのt_test関数をwilcox_testに変えるだけで、( )内の書き方を変える必要はありません。

順位和検定

サンプルデータです。このデータは1列目にグループ、2列目に値をとっています。valuesは握力とします。

groups	values
a	32.1
a	29.5
a	34.8
中略	中略
b	35.2
b	36.4
b	35.1

コードと結果は以下の通りです。

> wilcox_test(data = data1, values ~ groups, paired = FALSE)
# A tibble: 1 × 7
  .y.    group1 group2    n1    n2 statistic           p
* <chr>  <chr>  <chr>  <int> <int>     <dbl>       <dbl>
1 values a      b         20    20        34 0.000000925

p値から、有意水準を5%と設定した場合、帰無仮説が棄却され、2つのグループに有意な差があると結論づけることができます。

符号付順位和検定

対応のあるデータで符号順位和検定を実行します。データはここでもt_test関数（対応のあるt検定のコード）をwilcox_testに変更するだけで実行可能です。

サンプルデータです。同一の被験者で筋トレ前と筋トレ後で握力を測定した結果として載せています。groups列はpre(筋トレ前)とpost（筋トレ後）で、subjectsは被験者、valuesは握力の値となっています。対応のある検定ではsubjectsのデータは使わないが、対応ありということをわかりやすくするために載せました。

groups	subjects	values
pre	s1	32.3
pre	s2	29.5
pre	s3	34.8
中略	中略	中略
post	s1	34.6
post	s2	34.9
post	s3	36.2
中略	中略	中略

コードと結果は以下の通りです。

> wilcox_test(data = data1, values ~ groups, paired = TRUE)
# A tibble: 1 × 7
  .y.    group1 group2    n1    n2 statistic          p
* <chr>  <chr>  <chr>  <int> <int>     <dbl>      <dbl>
1 values post   pre       20    20       209 0.00000381

p値から、有意水準を5%と設定した場合、帰無仮説は棄却されるため、筋トレ前後に有意な差があると結論づけることができます。

相関分析（rstatixパッケージ）

2024-03-15T22:21:00.005+09:00

握力とジャンプのデータを想定したサンプルデータで、相関分析を行います。

grip	jump
30.1	56.7
27.5	48.3
32.8	57.1
中略	中略
32.2	58.6
24.4	53.6
30.7	55.1

コードと結果は以下の通りです。

> cor_test(data1, grip, jump, method = "pearson")
# A tibble: 1 × 8
  var1  var2    cor statistic       p conf.low conf.high method 
  <chr> <chr> <dbl>     <dbl>   <dbl>    <dbl>     <dbl> <chr>  
1 grip  jump   0.65      3.60 0.00205    0.287     0.847 Pearson

コードは、最初に解析対象であるdata1を指定します。grip列とjump列の相関を観たいので、gripとjumpを , で区切って記述します。相関分析では、順番はどちらが先でも結果は変わらないので、gripとjumpが逆でも構いません。method=で”pearson”を指定するとpearsonの積率相関係数を計算します。

結果です。相関係数はcorとして出力されます。今回は0.65ですので中等度の相関がありそうです。そして無相関検定の結果として、pが0.00205と出力されており、有意水準を5％と設定した場合、帰無仮説は棄却されました。ということで、この2つの変数の間には中等度の相関がある、と結論づけることができます。なお、conf.lowとconf.highは95%信頼区間であり、0.287-0.847ということになり、0を跨いでおらず、無相関検定と一致した結果となっています。

ノンパラメトリックな方法としてスピアマンを指定すると以下のような結果となります。

> cor_test(data1, grip, jump, method = "spearman")
# A tibble: 1 × 6
  var1  var2    cor statistic      p method  
  <chr> <chr> <dbl>     <dbl>  <dbl> <chr>   
1 grip  jump   0.65       472 0.0027 Spearman

今回用いたデータでは、pearson、spearmanの方法どちらでも相関係数はあまり変わりがない、結果となっています。

一元配置分散分析（rstatixパッケージ）

2024-03-15T22:20:00.001+09:00

分散分析

3つのグループに対しbmiを測定したデータについて、グループ間の平均値に下がるのかを分散分析で調べます。読み込んだデータは以下の通りです。

groups	bmi
a	18.3
a	17.89
a	17.55
中略	中略
b	20.23
b	21.03
b	19.1
中略	中略
c	22.87
c	21.96
c	22.37

コードと結果は以下の通りです。

> anova_test(data1, bmi ~ groups)
ANOVA Table (type II tests)

  Effect DFn DFd       F        p p<.05  ges
1 groups   2  57 121.351 2.86e-21     * 0.81

このコードでは、anova_test関数を( )内の内容にしたがって実行しています。( )内では、データの指定と解析対象の目的変数と説明変数の指定を行います。最初にデータの指定です。ここでは、data1です。次に、bmiとgroupsを ~ で結びます。bmiについてgroupsの影響を分析するという意味になります。

結果ですが、DFnとDFdは自由度で、FはF値、pはｐ値です。gesは効果量です。ｐ値は極めて小さい値となっており、0.05を有意水準と設定した場合、主効果がある、と判断できます。

多重比較

ここからさらに水準間での差を確認したい場合、多重比較を行ってきます。2群間の比較になるので、基本的にはt検定を行い、p値を補正していきます。rstatixパッケージでは、多重比較補正はholm法がデフォルトです。なおholm法はRのベースパッケージでもデフォルトになっています。

コードと結果は以下の通りです。

> t_test(data1, bmi ~ groups)
# A tibble: 3 × 10
  .y.   group1 group2    n1    n2 statistic    df        p    p.adj p.adj.signif
* <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>    <dbl>    <dbl> <chr>       
1 bmi   a      b         20    20     -9.30  31.4 1.55e-10 3.10e-10 ****        
2 bmi   a      c         20    20    -16.4   30.4 1.18e-16 3.54e-16 ****        
3 bmi   b      c         20    20     -6.29  37.9 2.29e- 7 2.29e- 7 ****

コードは、対応のないt検定と全く変わりません。多重比較の指定をしなくても、t検定の関数を知っているだけで、シンプルにできてしまうのがこの関数のすばらしいところです。水準が３（グループが３つ）の場合、組み合わせは３つ（a vs b、 a vs c、 b vs c）になりますが、特にそのことをコードで指定しなくても、結果として組み合わせを検出し、ｔ検定を行います。

　結果です。group1とgroup2は比較した2群を、n1とn2を、statisticはt値、dfは、pはp値、p.adjはholm補正後のp値です。この結果から、すべての群間に有意差があることが理解できます。

割合の検定（rstatixパッケージ）

2024-03-15T22:19:00.001+09:00

二項検定

50名に幸せだと思うか、思わないかを聴取し、幸せと思う人の割合が多いのかを二項検定で検討してみます。50名中、幸せだと思う人が35人で、思わない人が15人だったとします。コードと結果は以下のようになります。

最初にrstatixパッケージを読み込みます。

> library(rstatix)
> binom_test(35, 50, p = 0.5)
# A tibble: 1 × 6
      n estimate conf.low conf.high       p p.signif
* <dbl>    <dbl>    <dbl>     <dbl>   <dbl> <chr>   
1    50      0.7    0.554     0.821 0.00660 **

( )内には、幸せだと思う人の数（35人）、その次に全体の人数（50人）を入力します。また、p= は帰無仮説の比率になります。今回は割合が拮抗している（5割=0.5）というのが帰無仮説になるので0.5となります。帰無仮説が棄却されると、割合は拮抗しておらず偏っている、ということになります。

35人と15人と入力する場合、c(35, 15)をbinom_testに入れる必要があり、以下のように記述します。

> library(rstatix)

> binom_test(c(35, 15), p = 0.5)
# A tibble: 1 × 6
      n estimate conf.low conf.high       p p.signif
* <dbl>    <dbl>    <dbl>     <dbl>   <dbl> <chr>   
1    50      0.7    0.554     0.821 0.00660 **

上と下のコードは同じことですので、当然結果も同じになります。

カイ二乗検定（独立性の検定）

独立性の検定では、分析に先立って分割表を作成する必要があります。

例えば，自分は幸福だと思う人と思わない人が高齢者と若年者で変わるのかを検討したかったとします。幸福だと思う高齢者は58人で若年者は45人、幸福だと思わない高齢者は25人で若年者は38人だとします。このデータを作成するには、以下のコードを入力します。

> data1<-matrix(c(58, 45, 25, 38), ncol = 2, byrow = 2)
> data1
     [,1] [,2]
[1,]   58   45
[2,]   25   38

データの名前をdata1としています。matrixで2行2列のデータにします。( )内では、4つの数値の指定と行列数の指定を行います。4つの数値はc(数値, 数値, … )でくくります。Rでは数値が並んだデータをベクトルと呼びますが、こういったベクトルを作るときはc( )の( )の中に数値を , で区切って入力していきます。ncol=とbyrow= でそれぞれ列と行の数を指定できます。今回はいずれも2となります。data1を実行してみると2行2列の分割表が出力されます。

結果とコードは以下の通りです。

> chisq_test(data1)
# A tibble: 1 × 6
      n statistic      p    df method          p.signif
* <dbl>     <dbl>  <dbl> <int> <chr>           <chr>   
1   166      3.68 0.0549     1 Chi-square test ns

chisq.test()関数の中にデータ名を指定するだけです。

結果です。nはサンプルサイズ（全体で166人）、statisticはカイ二乗値、pはp値、dfは自由度になります。P値は0.05より大きいので帰無仮説は棄却できませんでした。

おまけ

以下のように、自分は幸福だと思うか、思わないか、を高齢者、中高年、若年者で調べ、どの世代に差が有るかを調べたいとします。

	older	middle	younger
yes	58	45	23
no	25	38	71

エクセル等の表計算ソフトで表を作り、コピーしてデータを読み込みます。以下がデータです。データの読み込みにはread.table()関数を用います。

> data1<-read.table("clipboard", header = TRUE)
> data1
    older middle younger
yes    58     45      23
no     25     38      71

read.delim関数と間違えないようにしましょう。

rstatixパッケージを使って、カイ二乗検定と、多重比較は行っていきます。コードと結果は以下の通りです。多重比較にはpairwise_prop_test関数を用います。prop_testは母比率の差の検定です。多重比較補正はholm補正です。

> chisq_test(data1)
# A tibble: 1 × 6
      n statistic             p    df method          p.signif
* <int>     <dbl>         <dbl> <int> <chr>           <chr>   
1   260      38.0 0.00000000557     2 Chi-square test ****    
> pairwise_prop_test(data1)
# A tibble: 3 × 5
  group1 group2              p        p.adj p.adj.signif
* <chr>  <chr>           <dbl>        <dbl> <chr>       
1 older  middle  0.0549        0.0549       ns          
2 older  younger 0.00000000362 0.0000000109 ****        
3 middle younger 0.0000939     0.000188     ***

多重比較は2×2に分割してカイ二乗検定を繰り返す方法もあるとは思います。

> pairwise_chisq_gof_test(data1)
# A tibble: 3 × 8
      n group1 group2  statistic             p    df        p.adj p.adj.signif
* <int> <chr>  <chr>       <dbl>         <dbl> <int>        <dbl> <chr>       
1   166 older  middle       3.68 0.0549            1 0.0549       ns          
2   177 older  younger     34.8  0.00000000362     1 0.0000000109 ****        
3   177 middle younger     15.3  0.0000939         1 0.000188     ***

二元配置分散分析（rstatixパッケージ）

2024-03-15T22:17:00.002+09:00

炎症状態（C反応性たんぱく；CRP）について、グループ（3グループ）や栄養状態（問題の無い群と低い群）の影響があるのかを調べます。読み込んだデータは以下の通りです。

groups	nutrition	crp
a	low	0.845
a	low	0.562
a	low	0.638
中略	中略	中略
b	low	0.469
b	low	0.286
b	low	0.393
中略	中略	中略
c	hight	0.276
c	hight	0.092
c	hight	0.216

コードと結果は以下の通りとなります。

> anova_test(data1, crp ~ groups * nutrition)
ANOVA Table (type II tests)

            Effect DFn DFd      F       p p<.05   ges
1           groups   2  54 40.368 1.9e-11     * 0.599
2        nutrition   1  54 27.590 2.6e-06     * 0.338
3 groups:nutrition   2  54  6.670 3.0e-03     * 0.198

入力したコードは、一元配置分散分析と同じanova_test関数で、( )内で二元配置の構造を指定します。最初にデータの指定です。ここでは、data1です。次に、crpとgroupsまたはnutrition を~ で結びます。ここで、2つの要因は*で結びます。crpについてgroupsのnutritionの影響を分析するという意味になります。

出力結果ですが、自由度、F値、ｐ値、効果量に関する情報が提示さされていることに変わりはありません。ただし、それぞれ3行に渡って結果が出力されています。groupsはグループの違いによる効果、nutritionは栄養状態の違いによる効果、groups:nutritionは交互作用の結果です。有意水準を0.05とすると、グループや栄養状態による主効果は確認されましたが、同時に交互作用が有意となっています。こうなった場合、要因の水準毎に効果が異なっている可能性があるので、単純主効果を調べる必要が出てきます。そこで、今回は、栄養状態に問題の無い群と低い群で分けて、グループ毎に一元配置分散分析を行っていくことにします。

ここで、パイプ演算子とgroup_by()関数が活躍します。データをgroup_by()で行うのです。実際には以下のようなコードになります。

> data1 |> group_by(nutrition) |> anova_test(crp ~ groups)
# A tibble: 2 × 8
  nutrition Effect   DFn   DFd     F        p `p<.05`   ges
* <chr>     <chr>  <dbl> <dbl> <dbl>    <dbl> <chr>   <dbl>
1 hight     groups     2    27  5.20 1.2 e- 2 *       0.278
2 low       groups     2    27 65.1  4.65e-11 *       0.828

data1について、group_by()で栄養状態に問題の無い群と、低い群で分けることを指定するために、data1とgroup_by()をパイプ演算子で結びます。そししてgroup_by()関数の()内に、nutritionを指定します。ここまででデータの切り分け方を指定した上で、パイプ演算子でanova_test()関数をつないで分散分析を行います。これにより、指定したデータを切り分けながらの分析が可能となります。

結果です。1行目は、栄養状態が問題の無い群での単純主効果の結果、2行目は栄養状態が低い群での単純主効果の結果です。いずれも有意な主効果となっています。主効果が確認されたので、さらに、多重比較へ進みます。

多重比較のコードと結果は以下の通りとなります。なお、ここではt検定に対し、栄養状態に問題の無い群と低い群ごとにholm補正をした結果が出力されます。

> data1 |> group_by(nutrition) |> t_test(crp ~ groups)
# A tibble: 6 × 11
  nutrition .y.   group1 group2    n1    n2 statistic    df            p        p.adj p.adj.signif
* <chr>     <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>        <dbl>        <dbl> <chr>       
1 hight     crp   a      b         10    10      2.00  12.1 0.069        0.137        ns          
2 hight     crp   a      c         10    10      2.72  11.1 0.02         0.059        ns          
3 hight     crp   b      c         10    10      1.30  17.4 0.209        0.209        ns          
4 low       crp   a      b         10    10      7.14  15.2 0.00000311   0.00000622   ****        
5 low       crp   a      c         10    10     10.3   17.0 0.0000000105 0.0000000315 ****        
6 low       crp   b      c         10    10      4.55  17.2 0.000276     0.000276     ***

入力したコードの前半は、パイプ演算子とgroup_by()関数を使用しており、単純主効果の分散分析と同じです。そして、t_test関数でt検定を行っていきます。すでにgroup_by()関数をつかって、栄養状態の水準毎に切ることを指定してありますので、t_test関数では、crpとgroupsを~ で結ぶだけです。これで、栄養状態の水準毎にcrpについてgroupsの多重比較をt検定で行う、ということになります。

結果です。1列目のnutritionは栄養状態に問題の無い群と低い群に群ごとの結果を行ごとわかるように示してあります。group1とgroup2は比較を行った群が記載されています。n1とn2は各群のサンプルサイズ、statisticはｔ値、dfは自由度、pは多重比較補正前のｐ値、p.adj pはholm法による多重比較補正後のp値を示しています。この結果をみると、栄養状態に問題の無い群ではどのグループ間でも有意差は認められませんでした。一方で、栄養状態の低い群では、全てのグループ間に有意差がありました。

次は、a、b、cの3つのグループ毎にデータを区切って、栄養状態の水準で比較を行います。ただし栄養状態の水準は2で、問題の無い群と低い群の比較になるため、分散分析を使わずにt検定を行っていきます。

コードと結果は以下の通りです。

> data1 |> group_by(groups) |> t_test(crp ~ nutrition)
# A tibble: 3 × 9
  groups .y.   group1 group2    n1    n2 statistic    df        p
* <chr>  <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>    <dbl>
1 a      crp   hight  low       10    10    -4.00   13.7 0.00138 
2 b      crp   hight  low       10    10    -4.06   17.2 0.000804
3 c      crp   hight  low       10    10    -0.857  17.4 0.403

コードは、前述の栄養状態の水準毎にグループ間の多重比較を行った時のコードと同じで、group_by関数とt_test関数内のnutritionとgroupsを入れ替えただけです。

結果です。ｃグループでは、栄養状態による有意差は確認できませんでしたが、aとcグループでは有意差が確認できました。

対応のあるt検定（rstatixパッケージ）

2024-03-15T22:17:00.001+09:00

以下のデータは、同一の被験者で筋トレ前と筋トレ後で握力を測定した結果です。groups列はpre(筋トレ前)とpost（筋トレ後）で、subjectsは被験者、valuesは握力の値となっています。対応のある検定ではsubjectsのデータは使わないが、対応ありということをわかりやすくするために載せました。

groups	subjects	values
pre	s1	32.3
pre	s2	29.5
pre	s3	34.8
中略	中略	中略
post	s1	34.6
post	s2	34.9
post	s3	36.2
中略	中略	中略

rstatixパッケージを読み込んでから、t検定を行います。

> library(rstatix)
> t_test(data = data1, values ~ groups, paired = TRUE)
# A tibble: 1 × 8
  .y.    group1 group2    n1    n2 statistic    df         p
* <chr>  <chr>  <chr>  <int> <int>     <dbl> <dbl>     <dbl>
1 values post   pre       20    20      5.76    19 0.0000148

対応なしのｔ検定と基本的には変わりません。t_test関数の( )内では、データの指定と対応の有り無しの指定を行います。このデータの指定の仕方は対応無しの時と同じです。valuesとgroupsを ~ で結んでおり、valuesの値についてgroupsで分けて分析するという意味になります。対応の有無は、paired=で指定し、TRUEとします。実行すると、以下の結果が出力されます。

結果です。見方は対応の無いt検定と変わりません。有意水準を5%と設定した場合、帰無仮説は棄却されるため、筋トレ前後で平均値に有意な差があると結論づけることができます。

一元配置反復測定分散分析（rstatixパッケージ）

2024-03-15T22:17:00.000+09:00

同一の被験者が3条件下で血圧を測定し、測定環境の違いが血圧に及ぼす影響を調べたかったとします。データは以下の通りです。このようなデータは、対応のあるデータとなり、反復測定分散分析が使われます。

conditions	subjects	bp
a	s1	102
a	s2	102
a	s3	110
中略	中略	中略
b	s1	112
b	s2	105
b	s3	106
中略	中略	中略
c	s18	120
c	s19	117
c	s20	120

分散分析のコードと結果は以下の通りです。

> anova_test(data1, bp ~ conditions + Error(subjects|conditions))
ANOVA Table (type III tests)

$ANOVA
      Effect DFn DFd      F        p p<.05   ges
1 conditions   2  38 49.099 2.93e-11     * 0.592

$`Mauchly's Test for Sphericity`
      Effect     W     p p<.05
1 conditions 0.984 0.862      

$`Sphericity Corrections`
      Effect   GGe      DF[GG]    p[GG] p[GG]<.05   HFe      DF[HF]    p[HF] p[HF]<.05
1 conditions 0.984 1.97, 37.39 4.11e-11         * 1.097 2.19, 41.67 2.93e-11         *

anova_test関数を使い、目的変数（アウトカム指標、ここでは血圧）と説明変数（要因、ここでは測定環境）を ~ で結ぶところまでは対応のないデータの分散分析と一緒です。対応があるデータということを記述するために、Error(subjects | conditions)を追加します。Error()の中は、| で区切られています。測定環境は、同一の被験者で測定した、ということを記載している、と理解しておきましょう。

$ANOVAは通常の反復測定のための分散分析（typeⅢ）で、$`Mauchly's Test for Sphericity`は球面性の検定であり、分散分析ではありません。$`Sphericity Corrections`は自由度を補正した分散分析の結果で、Greenhouse-GeisserとHuynh-Feldtの方法についての結果が示されています。球面性の仮定が成立しないような場合、自由度補正法を使います。Greenhouse-Geisserの方が結果は保守的になるようです。

自由度補正法のGGeはGreenhouse-Geisserの統計量ε、DF[GG]は自由度、p[GG]はp値、HFeはHuynh-Feldtの統計量ε、DF[HF]は自由度、p[HF]はp値です。今回の結果は通常のANOVA、自由度を補正したANOVA、どちらでも主効果が確認されできる結果となっています。

主効果が確認された場合、多重比較を行うことがあります。多重比較のコードと結果は以下の通りです。

> data1 |> t_test(bp ~ conditions, paired = TRUE)
# A tibble: 3 × 10
  .y.   group1 group2    n1    n2 statistic    df            p        p.adj p.adj.signif
* <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>        <dbl>        <dbl> <chr>       
1 bp    a      b         20    20     -4.24    19 0.000441     0.000441     ***         
2 bp    a      c         20    20     -9.45    19 0.0000000129 0.0000000387 ****        
3 bp    b      c         20    20     -5.95    19 0.00000989   0.0000198    ****

パイプ演算子を使ったコードで表示しています（このくらいのコードではパイプ演算子を使わなくてもそれほど、解読しにくいわけではありません）。対応のあるt検定のときの入力規則と全く変わりません。水準が3つ以上になり、複数の組み合わせになったので、自動で多重比較補正の結果も返してくれます。

結果です。多重比較はt検定のホルム補正法です。補正されたp値を確認すると、すべての組み合わせで有差があることが理解できます。

Rで散布図を描く（ggplotなし）

2023-01-29T21:32:00.004+09:00

以下のデータを視覚化する。これは、相関分析で使ったデータと同じ。今回のデータは計画行列ではない。架空のデータで、a列が握力、ｂ列が垂直飛びの結果だとする。

a	b
30.1	56.7
27.5	48.3
32.8	57.1
29.6	52.2
30.6	53.4
32.1	57.7
23.3	49.2
29.4	56.9
28.2	56.6
28.8	53.5
25.7	54.2
30.2	51
29.3	52.8
28	53.9
25.8	51.3
30.9	54.9
32.5	57.4
32.2	58.6
24.4	53.6
30.7	55.1

　このデータのa列をx軸とb列をy軸に割り当てて散布図を描いてみる。

　なお、学会発表や論文のために作図する場合は、より自由度が高く、美しい作図が可能なggplotパッケージを用いた方が良い。

　data1としてデータを読みこんでいる。

　コードは以下の通り。plotが関数の名前で () 内に必要な情報を入力していく。最初にx軸(data1$a)、そのあとにy軸(data1$b)を指定する。

plot(data1$a, data1$b)

　x軸とy軸のラベルを任意に変更する場合にはxlab=、ylab=で指定してやればよい。

plot(data1$a, data1$b, xlab = "a", ylab = "b")

　ここからは、ちょっと難しいので、初級っぽくないので流しても良い。

　複数の散布図を重ねることもできる。以下のように2つのグループの握力と垂直飛びのデータがあったとして、それを散布図で重ねて表示したいとする。

a	b	group
30.1	56.7	group_1
27.5	48.3	group_1
32.8	57.1	group_1
29.6	52.2	group_1
30.6	53.4	group_1
32.1	57.7	group_1
23.3	49.2	group_1
29.4	56.9	group_1
28.2	56.6	group_1
28.8	53.5	group_1
25.7	54.2	group_1
30.2	51	group_1
29.3	52.8	group_1
28	53.9	group_1
25.8	51.3	group_1
30.9	54.9	group_1
32.5	57.4	group_1
32.2	58.6	group_1
24.4	53.6	group_1
30.7	55.1	group_1
37.2	62.2	group_2
31	49.7	group_2
35.8	58.4	group_2
30	54.3	group_2
32.2	53	group_2
32.3	59.7	group_2
23.9	52.5	group_2
33.2	58.5	group_2
31.4	60.1	group_2
29.3	53.3	group_2
28	53.7	group_2
34.1	50.8	group_2
29.6	56.4	group_2
30.6	56.4	group_2
28.5	54.3	group_2
34.6	57.3	group_2
35.9	56.6	group_2
33.8	60.4	group_2
26	53.6	group_2
31.3	55.6	group_2

　コードは以下の通り。上のデータをdata2として読み込んでいる。ドットに色を付けて区別させたいので、色がわかりやすいドットの形に変更する。pch=で変更する。19番で塗りつぶしのドットを出力できる。また、col=で色分けを指定する。factor()として色分けしたい列の名前を入れる。

plot(data2$a, data2$b, xlab = "a", ylab = "b", pch = 19, col = factor(data2$group))

ただし、これだとどの色がどのグループかわからない。そこで、plot関数を実行した後に、legends関数を実行する。コードと図は以下の通り。 () 内は最初に凡例を表示する場所、凡例として表示するデータ、凡例に表示するドットの形、凡例の色、を指定していく。表示場所は左上にしたいので"topleft"とする。凡例を表示するためにはlegend=levels(factor())と入力し、factorの () 内に列名data2$groupを指定する。pch=は19番としてドットの形を指定する。col=factor()で色を指定する列data2$groupを指定する（ただし、factorをlevelsで囲む）。

plot(data2$a, data2$b, xlab = "a", ylab = "b", pch = 19, col = factor(data2$group))
legend("topleft", legend = levels(factor(data2$group)), pch = 19, col = factor(levels(factor(data2$group))))

Rで独立性の検定（χ二乗検定）

2023-01-22T17:40:00.000+09:00

例えば，自分は幸福だと思う人と思わない人が高齢者と若年者で変わるのかを検討したかったとする。そんなときはχ二乗検定が役に立つ。幸福だと思う高齢者は58人で若年者は45人、幸福だと思わない高齢者は25人で若年者は38人だったとする。

まず、データの読み込み。matrixという関数で2行2列のデータにする。ncol=とbyrow= でそれぞれ列と行の数を指定する。

data1<-matrix(c(58, 45, 25, 38), ncol = 2, byrow = 2)
data1 #以下が読み込んだデータ
     [,1] [,2]
[1,]   58   45
[2,]   25   38

次に、χ二乗検定を実行する。chisq.testで実行できる。 () 内に作成したデータを記述する。p-valueから有意性検定の結果を読み取ることができる。

chisq.test(data1)

	Pearson's Chi-squared test with Yates' continuity correction

data:  data1
X-squared = 3.6838, df = 1, p-value = 0.05494

Rで二項検定

2023-01-20T11:57:00.000+09:00

2値のカテゴリカルな値の比率に差があるかを調べたい時に二項検定を用いる。例えば、ある集団に対し、自分は幸せだと思う人と思わない人の数に差があるかを調べたい時に用いることができる。

幸せだと思う人と思わない人の比率が1:1（全体の中の好きな人の比率が0.5）になっている時には、どちらも拮抗していてどちらが多いとは言えない（帰無仮説）、ということになる。この帰無仮説を棄却できるかを検証するのがこの検定の特徴になる。

全体が50名で、幸せだと思う人が35人で、思わない人が15人だったとすると、コードは以下のように記述する。()内の35は、幸せだと思う人の数で50は全体の数となる。また、p= は帰無仮説の比率となっている。今回は比率は拮抗しているというのが帰無仮説だから0.5となる。

binom.test(35, 50, p = 0.5)

結果は、以下の通り。p-valueは0.0066なので結構小さい。95 percent confidence intervalは、幸せだと思う人の全体に対する比率の信頼区間をあらわしている。

	Exact binomial test

data:  35 and 50
number of successes = 35, number of trials = 50, p-value = 0.0066
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.5539177 0.8213822
sample estimates:
probability of success 
                   0.7

Rで箱ひげ図と1次元散布図を描く（ggplotなし）

2023-01-19T17:21:00.003+09:00

以下の計画行列のデータを視覚化する。aとbの2つのグループの握力のデータとして表している。groups列はaとbのグループが、valuse列には握力の値がそれぞれ記述されている。

このデータをaとbに分けて箱ひげ図や散布図を描いてみる。なお、学会発表や論文のために作図する場合は、より自由度が高く、美しい作図が可能なggplotパッケージを用いた方が良い。

groups	values
a	32.1
a	29.5
a	34.8
a	31.6
a	32.6
a	34.1
a	25.3
a	31.4
a	30.2
a	30.8
a	27.7
a	32.2
a	31.3
a	30.0
a	27.8
a	32.9
a	34.5
a	34.2
a	26.4
a	32.7
b	34.6
b	34.9
b	36.2
b	34.3
b	35.7
b	35.8
b	32.3
b	32.8
b	35.3
b	33.8
b	33.1
b	33.3
b	34.7
b	35.0
b	35.9
b	33.5
b	34.4
b	35.2
b	36.4
b	35.1

箱ひげ図

箱ひげ図は、以下のコードになる。boxplotが関数の名前で () 内に必要な情報を入力していく。valuesとgroupsを ~ で結んでいるが、これはvaluesの値についてgroupsで分けて箱ひげ図にして、という意味。ylab=とxlab=はy軸とx軸の名前を指定するという意味で、特に入力しなくても図は出力される。なお、y軸の最大値と最小値を指定したい場合は、 , で区切ってylim=c(0, 36)を付け加えればよい。この意味は0が最小で36が最大という意味。数字は任意で、、

boxplot(data1$values ~ data1$groups, ylab = "Grip_strength", xlab = "Groups")

図は以下の通りに出力される。

1次元散布図

以下のコードとなる。stripchartが関数の名前で、 () 内に必要な情報を入力していく。データも軸の名前もboxplotと同じ記述方法になっている。ただし、デフォルトでは、x軸にvaluesが記載される。

stripchart(data1$values ~ data1$groups)

そこで、以下のようにvertical=を加え、TRUEにすると、valuesがy軸の図になる。

stripchart(data1$values ~ data1$groups, vertical = TRUE)

ドットが重なって見にくい場合は、method= を加えて、jitterを指定すると以下のように、ドットを横に散らしてくれる。

stripchart(data1$values ~ data1$groups, vertical = TRUE, method = "jitter")

さらに、各軸に名前を付けるときには以下のようにする。xlab=とylab=で指定すればよい。

stripchart(data1$values ~ data1$groups, vertical = TRUE, method = "jitter", xlab = "group", ylab = "Grip_strength")

図は以下のようになる。

箱ひげ図と散布図を重ねる

箱ひげ図と散布図を重ねて表示するとさらにデータの見通しが立ちやすくなるかもしれない。その場合には、まずは箱ひげ図のコードを実行し、図を出力してから、散布図を重ねる処理を行う。散布図のコードにadd=TRUEを加えればよい。

boxplot(data1$values ~ data1$groups, xlab = "group", ylab = "Grip_strength")

stripchart(data1$values ~ data1$groups, vertical = TRUE, method = "jitter", xlab = "group", ylab = "Grip_strength", add = TRUE)

図は以下の通りとなる。

Rで符号順位和検定（対応のある2群の差のノンパラメトリック検定）

2023-01-18T15:20:00.004+09:00

読み込んだデータから対応のある2群間の差の検定を行ってみる。今回は、ノンパラメトリックな方法で、符号順位和検定を使う。

読み込んだデータは表の通り。読み込み方はここ。

このデータは、計画行列（デザイン行列）となっている。同一の被験者で筋トレ前と筋トレ後で握力を測定したことを想定した仮想のデータとなっている。groups列はpre(前)とpost（後）で、subjectsは被験者、valuesは握力の値となっている。今回の検定ではsubjectsのデータは使わないが、対応ありということでわかりやすくするため載せている。

groups	subjects	values
pre	s1	32.1
pre	s2	29.5
pre	s3	34.8
pre	s4	31.6
pre	s5	32.6
pre	s6	34.1
pre	s7	25.3
pre	s8	31.4
pre	s9	30.2
pre	s10	30.8
pre	s11	27.7
pre	s12	32.2
pre	s13	31.3
pre	s14	30
pre	s15	27.8
pre	s16	32.9
pre	s17	34.5
pre	s18	34.2
pre	s19	26.4
pre	s20	32.7
post	s1	34.6
post	s2	34.9
post	s3	36.2
post	s4	34.3
post	s5	35.7
post	s6	35.8
post	s7	32.3
post	s8	32.8
post	s9	35.3
post	s10	33.8
post	s11	33.1
post	s12	33.3
post	s13	34.7
post	s14	35
post	s15	35.9
post	s16	33.5
post	s17	34.4
post	s18	35.2
post	s19	36.4
post	s20	35.1

でもって、以下のコードを実行する

なお、このデータの名前はdata1という名前になっている。~の前にvaluesを~の後にgroupsを指定する。なお、data1の後に$をつけてからvalusesまたはgroupsを指定しているが、これはdata1の中にあるvaluesまたはgroupsという意味になる。

paired=は、対応があるかないかを指定する。TRUEにすると対応ありを指定できる。

 wilcox.test(data1$values ~ data1$groups, paired = TRUE)

そうすると以下の結果が返ってくる。

	Wilcoxon signed rank exact test

data:  data1$values by data1$groups
V = 1, p-value = 3.815e-06
alternative hypothesis: true location shift is not equal to 0

p-valueがp値になる。3.815e-05は指数表示で3.815*10のマイナス5乗ということで、すごく小さい値になっている。

なお、以下のようなデータの並びにして分析をする場合も記述しておく。

これは1列目被験者、2列目にpreの握力、3列目にpostの握力を記述している。1列目のsubjectsはわかりやすくするために記述しているが、無くても良い。

subjects	pre	post
s1	32.1	34.6
s2	29.5	34.9
s3	34.8	36.2
s4	31.6	34.3
s5	32.6	35.7
s6	34.1	35.8
s7	25.3	32.3
s8	31.4	32.8
s9	30.2	35.3
s10	30.8	33.8
s11	27.7	33.1
s12	32.2	33.3
s13	31.3	34.7
s14	30	35
s15	27.8	35.9
s16	32.9	33.5
s17	34.5	34.4
s18	34.2	35.2
s19	26.4	36.4
s20	32.7	35.1

このような場合のt検定のコードは以下の通り。~を , にして、2列目にpre、3列目にpostを配置している。結果は同じになる。

wilcox.test(data1$pre, data1$post, paired = TRUE)

Rで相関分析（pearsonとspearman）

2023-01-17T23:59:00.005+09:00

Rで相関分析を行う。データは以下の通り。架空のデータで、a列が握力、ｂ列が垂直飛びの結果だとする。

a	b
30.1	56.7
27.5	48.3
32.8	57.1
29.6	52.2
30.6	53.4
32.1	57.7
23.3	49.2
29.4	56.9
28.2	56.6
28.8	53.5
25.7	54.2
30.2	51
29.3	52.8
28	53.9
25.8	51.3
30.9	54.9
32.5	57.4
32.2	58.6
24.4	53.6
30.7	55.1

Ｒのコードは以下の通り。データの名前はdata1。a列とb列のデータの相関を観たいので、data$aとdata$bを , で区切って記述する。相関分析では、順番はどちらが先でも結果は変わらない。method=でpearsonを指定するとpearsonの積率相関係数を計算してくれる。

cor.test(data1$a,data1$b, method = "pearson")　

結果は以下の通り。p-valueは無相関検定のp値で、相関係数はcorで、約0.65となっている。95 percent confidence intervalは、相関係数の95%信頼区間。

	Pearson's product-moment correlation

data:  data1$a and data1$b
t = 3.6002, df = 18, p-value = 0.002046
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2865442 0.8470224
sample estimates:
      cor 
0.6470215

ノンパラメトリックな方法としてspearmanを使いたい場合、method=をspearmanとすればよい。なお、結果は以下の通りとなる。

	Spearman's rank correlation rho

data:  data1$a and data1$b
S = 472, p-value = 0.002699
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.6451128

Rで順位和検定（対応のない2群の差のノンパラメトリック検定）

2023-01-17T16:51:00.002+09:00

読み込んだデータから対応のない2群間の差の検定を行ってみる。今回は、ノンパラメトリックな方法で、順位和検定を使う。

読み込んだデータは表の通りで、ｔ検定の時と一緒。読み込み方はここ。

データは、計画行列（デザイン行列）になっている。このデータは勝手に作ったデータで、2グループの握力を表しているものとする。groups列は2つのグループで、aとbのグループに分けられている。valuesは分析対象の値で、ここでは握力となっている。1行に1人分のデータがあることになり、各列から所属のグループと握力がわかることになる。

groups	values
a	32.1
a	29.5
a	34.8
a	31.6
a	32.6
a	34.1
a	25.3
a	31.4
a	30.2
a	30.8
a	27.7
a	32.2
a	31.3
a	30.0
a	27.8
a	32.9
a	34.5
a	34.2
a	26.4
a	32.7
b	34.6
b	34.9
b	36.2
b	34.3
b	35.7
b	35.8
b	32.3
b	32.8
b	35.3
b	33.8
b	33.1
b	33.3
b	34.7
b	35.0
b	35.9
b	33.5
b	34.4
b	35.2
b	36.4
b	35.1

順位和検定を実行するために以下のコードを実行する。

paired=は、対応があるかないかを指定する。FALSEにすると対応なしを指定できる。

wilcox.test(data1$values ~ data1$groups, paired = FALSE)

そうすると以下の結果が返ってくる。

	Wilcoxon rank sum exact test

data:  data1$values by data1$groups
W = 34, p-value = 9.249e-07
alternative hypothesis: true location shift is not equal to 0

p-valueがp値になる。9.249e-07は9.249*10のマイナス5乗ということで、すごく小さい値になっている。95%confidence intervalは差の信頼区間を表している。

なお、以下のようなデータの並びにして分析をする場合も記述しておく。

これは1列目にaグループの握力、2列目にbグループの握力を記述している。

group_a	group_b
32.1	34.6
29.5	34.9
34.8	36.2
31.6	34.3
32.6	35.7
34.1	35.8
25.3	32.3
31.4	32.8
30.2	35.3
30.8	33.8
27.7	33.1
32.2	33.3
31.3	34.7
30.0	35.0
27.8	35.9
32.9	33.5
34.5	34.4
34.2	35.2
26.4	36.4
32.7	35.1

このような場合のt検定のコードは以下の通り。~を , にし、group_aとgroup_bを配置している。結果は同じになる。

wilcox.test(data1$group_a, data1$group_b, paired = FALSE)

Rでヒストグラムを描く（ggplotなし）

2023-01-16T14:33:00.006+09:00

分析対象のデータの特徴をつかむためヒストグラムを描くことがある。Rのベース関数で描画してみる。

データは表のとおり。データは計画行列で記述されている。

groups	values
a	32.1
a	29.5
a	34.8
a	31.6
a	32.6
a	34.1
a	25.3
a	31.4
a	30.2
a	30.8
a	27.7
a	32.2
a	31.3
a	30.0
a	27.8
a	32.9
a	34.5
a	34.2
a	26.4
a	32.7
b	34.6
b	34.9
b	36.2
b	34.3
b	35.7
b	35.8
b	32.3
b	32.8
b	35.3
b	33.8
b	33.1
b	33.3
b	34.7
b	35.0
b	35.9
b	33.5
b	34.4
b	35.2
b	36.4
b	35.1

ここから、グループbのデータを取り出して、ヒストグラムを描く。

bのデータの抽出についてはここを参照。

コードは以下の通り。histと記述されているのがヒストグラムを出力する関数で、()の中に必要な情報を記述する。data1[21:40, 2]はデータのことになる。2列目の21～40行目のデータを指定していれている。 ()内は、違う情報を記載するたびに , で区切る必要がある。データ以外の上方は指定しなくてもヒストグラムは出力される。よりきれいに出力したい場合には色々指定する必要がある。breaks=は図の最大・最小・間隔を記載する。main=となっているのは、図のタイトル。ylabとかxlabとかはy軸やx軸につける名前のこと。

hist(data1[21:40,2], breaks = seq(32, 38, 1), main = "HIstogram", ylab = "n", xlab = "group_b")

下のような図が描画される。

Rで任意の行や列を抽出する（tidyverseなし）

2023-01-16T13:42:00.006+09:00

以下の表から、任意のデータを取り出す。

こんな計画行列のデータがあったとする。データの詳細はここを参照。ここからgroups列のaだけを抽出したい。なお、今回はデータの名前をdata1としている。

groups	values
a	32.1
a	29.5
a	34.8
a	31.6
a	32.6
a	34.1
a	25.3
a	31.4
a	30.2
a	30.8
a	27.7
a	32.2
a	31.3
a	30.0
a	27.8
a	32.9
a	34.5
a	34.2
a	26.4
a	32.7
b	34.6
b	34.9
b	36.2
b	34.3
b	35.7
b	35.8
b	32.3
b	32.8
b	35.3
b	33.8
b	33.1
b	33.3
b	34.7
b	35.0
b	35.9
b	33.5
b	34.4
b	35.2
b	36.4
b	35.1

任意の行を取り出す。

コードは以下の通りとなる。データの名前はdata1としている。データの中の抜き出したい行数を指定する。1から20行目までがグループaのデータになるため、1:20として、抜き出したい行数を指定する。[]の中の , より左は行、右は列を指定することになる。今回は任意の行を抜き出したいので、 , より左側で行数を指定する。

data1[1:20, ]

出てくる結果は以下の通り。

   groups values
1       a   32.1
2       a   29.5
3       a   34.8
4       a   31.6
5       a   32.6
6       a   34.1
7       a   25.3
8       a   31.4
9       a   30.2
10      a   30.8
11      a   27.7
12      a   32.2
13      a   31.3
14      a   30.0
15      a   27.8
16      a   32.9
17      a   34.5
18      a   34.2
19      a   26.4
20      a   32.7

bの行を抽出したかったら、[]の中を21:40に指定すればよい。

行数を指定するのではなく、行の文字（aとかb）を指定して抽出する方法もある。コードは以下の通りだが、可読性があまり高くない。

今回はtidyverseパッケージを使用していないが、tidyverseの方がコードの可読性が高い。

data1[data1$groups == "a", ]

任意の列を取り出す。

コードは以下の通りとなる。データの名前はdata1としている。データの中の抜き出したい列数を指定する。[]の中の , より右に列を指定することになる。もし、1から2列目までを指定したかったら1:2とする。2列目だけを抽出したければ、2とする。
1から2列目までを抽出するコードと結果は以下の通りとなる。

 data1[ ,1:2] #コード
   groups values # ここから結果
1       a   32.1
2       a   29.5
3       a   34.8
4       a   31.6
5       a   32.6
6       a   34.1
7       a   25.3
8       a   31.4
9       a   30.2
10      a   30.8
11      a   27.7
12      a   32.2
13      a   31.3
14      a   30.0
15      a   27.8
16      a   32.9
17      a   34.5
18      a   34.2
19      a   26.4
20      a   32.7
21      b   34.6
22      b   34.9
23      b   36.2
24      b   34.3
25      b   35.7
26      b   35.8
27      b   32.3
28      b   32.8
29      b   35.3
30      b   33.8
31      b   33.1
32      b   33.3
33      b   34.7
34      b   35.0
35      b   35.9
36      b   33.5
37      b   34.4
38      b   35.2
39      b   36.4
40      b   35.1

2列目だけを抽出するコードと結果は以下の通り。出てくる結果は以下の通り。

data1[ ,2] #ここはコードで下は結果
 [1] 32.1 29.5 34.8 31.6 32.6 34.1 25.3 31.4 30.2 30.8 27.7 32.2 31.3 30.0 27.8 32.9 34.5 34.2 26.4 32.7 34.6 34.9 36.2
[24] 34.3 35.7 35.8 32.3 32.8 35.3 33.8 33.1 33.3 34.7 35.0 35.9 33.5 34.4 35.2 36.4 35.1

任意の行と列を取り出す。

こんどは行と列を組み合わせるだけ。たとえば、bグループの2列目だけを取り出す。bグループは21から40行目までになるので21:40と指定する。2列目は2と指定すればよい。コードと結果は以下の通り。

data1[21:40, 2] #ここはコードで下は結果
 [1] 34.6 34.9 36.2 34.3 35.7 35.8 32.3 32.8 35.3 33.8 33.1 33.3 34.7 35.0 35.9 33.5 34.4 35.2 36.4 35.1

Rで対応のあるt検定

2023-01-16T11:33:00.012+09:00

読み込んだデータから対応のあるｔ検定を行ってみる。

読み込んだデータは表の通り。読み込み方はここ。

このデータは、計画行列（デザイン行列）といって、列ごとに得られた情報をまとめて配置したデータの配置方法となっている。回帰分析では計画行列が使われる。今回はt検定だが、計画行列を使ってみることにする。同一の被験者で筋トレ前と筋トレ後で握力を測定したことを想定している。groups列はpre(前)とpost（後）で、subjectsは被験者、valuesは握力の値となっている。今回の検定ではsubjectsのデータは使わないが、対応ありということで載せている。

groups	subjects	values
pre	s1	32.1
pre	s2	29.5
pre	s3	34.8
pre	s4	31.6
pre	s5	32.6
pre	s6	34.1
pre	s7	25.3
pre	s8	31.4
pre	s9	30.2
pre	s10	30.8
pre	s11	27.7
pre	s12	32.2
pre	s13	31.3
pre	s14	30
pre	s15	27.8
pre	s16	32.9
pre	s17	34.5
pre	s18	34.2
pre	s19	26.4
pre	s20	32.7
post	s1	34.6
post	s2	34.9
post	s3	36.2
post	s4	34.3
post	s5	35.7
post	s6	35.8
post	s7	32.3
post	s8	32.8
post	s9	35.3
post	s10	33.8
post	s11	33.1
post	s12	33.3
post	s13	34.7
post	s14	35
post	s15	35.9
post	s16	33.5
post	s17	34.4
post	s18	35.2
post	s19	36.4
post	s20	35.1

でもって、以下のコードを実行する

paired=は、対応があるかないかを指定する。TRUEにすると対応ありを指定できる。

t.test(data1$values ~ data1$groups, paired = TRUE)

そうすると以下の結果が返ってくる。

	Paired t-test

data:  data1$values by data1$groups
t = 5.8586, df = 19, p-value = 1.214e-05
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 2.256032 4.763968
sample estimates:
mean difference 
           3.51

p-valueがp値になる。1.214e-05は指数表示で1.214*10のマイナス5乗ということで、すごく小さい値になっている。95%confidence intervalは差の信頼区間を表している。

なお、以下のようなデータの並びにして分析をする場合も記述しておく。

subjects	pre	post
s1	32.1	34.6
s2	29.5	34.9
s3	34.8	36.2
s4	31.6	34.3
s5	32.6	35.7
s6	34.1	35.8
s7	25.3	32.3
s8	31.4	32.8
s9	30.2	35.3
s10	30.8	33.8
s11	27.7	33.1
s12	32.2	33.3
s13	31.3	34.7
s14	30	35
s15	27.8	35.9
s16	32.9	33.5
s17	34.5	34.4
s18	34.2	35.2
s19	26.4	36.4
s20	32.7	35.1

このような場合のt検定のコードは以下の通り。~を , にして、2列目にpre、3列目にpostを配置している。結果は同じになる。

t.test(data1$pre, data1$post, paired = TRUE)

Rで対応のないt検定

2023-01-11T18:39:00.015+09:00

読み込んだデータから対応のないｔ検定を行ってみる。正しくはにはWelchのt検定になる。等分散を仮定しないt検定になる。

読み込んだデータは表の通り。読み込み方はここ。

このデータは、計画行列（デザイン行列）といって、列ごとに得られた情報をまとめて配置したデータの配置方法となっている。回帰分析では計画行列が使われる。今回はt検定だが、計画行列を使ってみることにする。このデータは勝手に作ったデータで、2グループの握力を表しているものとする。groups列は2つのグループで、aとbのグループに分けられている。valuesは分析対象の値で、ここでは握力となっている。1行に1人分のデータがあることになり、各列から所属のグループと握力がわかることになる。

groups	values
a	32.1
a	29.5
a	34.8
a	31.6
a	32.6
a	34.1
a	25.3
a	31.4
a	30.2
a	30.8
a	27.7
a	32.2
a	31.3
a	30
a	27.8
a	32.9
a	34.5
a	34.2
a	26.4
a	32.7
b	34.6
b	34.9
b	36.2
b	34.3
b	35.7
b	35.8
b	32.3
b	32.8
b	35.3
b	33.8
b	33.1
b	33.3
b	34.7
b	35
b	35.9
b	33.5
b	34.4
b	35.2
b	36.4
b	35.1

でもって、以下のコードを実行する。

paired=は、対応があるかないかを指定する。FALSEにすると対応なしを指定できる。

t.test(data1$values ~ data1$groups, paired = FALSE)

そうすると以下の結果が返ってくる。

	Welch Two Sample t-test

data:  data1$values by data1$groups
t = -5.3544, df = 25.931, p-value = 1.333e-05
alternative hypothesis: true difference in means between group a and group b is not equal to 0
95 percent confidence interval:
 -4.857642 -2.162358
sample estimates:
mean in group a mean in group b 
         31.105          34.615

p-valueがp値になる。1.333e-05は指数表示で1.333*10のマイナス5乗ということで、すごく小さい値になっている。95%confidence intervalは差の信頼区間を表している。

なお、以下のようなデータの並びにして分析をする場合も記述しておく。

これは1列目にaグループの握力、2列目にbグループの握力を記述している。

group_a	group_b
32.1	34.6
29.5	34.9
34.8	36.2
31.6	34.3
32.6	35.7
34.1	35.8
25.3	32.3
31.4	32.8
30.2	35.3
30.8	33.8
27.7	33.1
32.2	33.3
31.3	34.7
30	35
27.8	35.9
32.9	33.5
34.5	34.4
34.2	35.2
26.4	36.4
32.7	35.1

このような場合のt検定のコードは以下の通り。~を , にし、group_aとgroup_bを配置している。結果は同じになる。

t.test(data1$group_a, data1$group_b, paired = FALSE)

エクセルで作成したデータをコピーしてRに読み込む

2023-01-11T17:38:00.009+09:00

エクセルで作ってあるデータをRに読み込みたい。Windowsでの操作は以下の通り。

取り込みたい部分のデータを選択してコピーする。

Rのコンソールで以下のコードを実行する。

  # data1という名前でエクセルからコピーした値を読み込む。
  data1<-read.delim("clipboard", header = TRUE) 
  
  # 読み込んだデータを確認するためにdata1を入力して実行する。
  data1 
     groups subjects values
  1       a       s1   32.3
  2       a       s2   29.5
  3       a       s3   34.8
  4       a       s4   31.6
  5       a       s5   32.6
  6       a       s6   34.6
  7       a       s7   25.3
  8       a       s8   31.4
  9       a       s9   30.2
  10      a      s10   30.8
  11      a      s11   27.8
  12      a      s12   32.2
  13      a      s13   31.3
  14      a      s14   30.0
  15      a      s15   27.8
  16      a      s16   32.9
  17      a      s17   34.5
  18      a      s18   34.2
  19      a      s19   26.4
  20      a      s20   32.7
  21      b       s1   34.6
  22      b       s2   34.9
  23      b       s3   36.2
  24      b       s4   34.3
  25      b       s5   35.9
  26      b       s6   35.8
  27      b       s7   32.3
  28      b       s8   32.8
  29      b       s9   35.2
  30      b      s10   33.8
  31      b      s11   33.1
  32      b      s12   33.3
  33      b      s13   34.7
  34      b      s14   35.0
  35      b      s15   35.9
  36      b      s16   33.5
  37      b      s17   34.4
  38      b      s18   35.2
  39      b      s19   36.4
  40      b      s20   35.1

  # 読み込んだデータを確認するため頭の部分だけを表示する。
  data1 |> head() #
      groups subjects values
  1      a       s1   32.3
  2      a       s2   29.5
  3      a       s3   34.8
  4      a       s4   31.6
  5      a       s5   32.6
  6      a       s6   34.6

My notepad

Diff動作解析チュートリアル（DIFFマーカーセットで測定したc3dデータを使って動作分析を行う）

はじめに

準備するデータ

解析の手順

結果で確認すること

ビューアの使い方

再生と描画区間

3D スティックピクチャー

床面投影図

時系列グラフ

出力ファイル

動作ごとの解析

主な設定

困ったときは

「準備しています…」から進まない

「マーカー ○○ が見つかりません」と表示される

「体重が不明です」と表示される

関節モーメントが一部の区間で空欄になる

床反力計の判定が違う

Plot Configurator User Manual

Plot Configurator User Manual

Table of Contents

Introduction

Quick Start

Detailed Usage

1. Data Entry

A. For Line / Bar Chart

B. For Scatter Plot

2. Selecting Chart Type

3. Customizing Appearance

4. Axis Configuration

5. Saving Images

Tips

Plot Configuratorの使い方

Plot Configurator 取扱説明書

目次

はじめに

クイックスタート

詳細な使い方

1. データの入力

A. 折れ線グラフ・棒グラフ (Line / Bar) の場合

B. 散布図 (Scatter Plot) の場合

2. グラフの種類の選択

3. 見た目のカスタマイズ

4. 軸の設定

5. 画像の保存

ヒント

重回帰分析（モデル選択を含む）

単回帰分析

対応のないt検定（rstatixパッケージ）

順位和検定と符号付順位和検定（rsatixパッケージ）

相関分析（rstatixパッケージ）

一元配置分散分析（rstatixパッケージ）

割合の検定（rstatixパッケージ）

二元配置分散分析（rstatixパッケージ）

対応のあるt検定（rstatixパッケージ）

一元配置反復測定分散分析（rstatixパッケージ）

Rで散布図を描く（ggplotなし）

Rで独立性の検定（χ二乗検定）

Rで二項検定

Rで箱ひげ図と1次元散布図を描く（ggplotなし）

Rで符号順位和検定（対応のある2群の差のノンパラメトリック検定）

Rで相関分析（pearsonとspearman）

Rで順位和検定（対応のない2群の差のノンパラメトリック検定）

Rでヒストグラムを描く（ggplotなし）

Rで任意の行や列を抽出する（tidyverseなし）

Rで対応のあるt検定

Rで対応のないt検定

エクセルで作成したデータをコピーしてRに読み込む