Hatena::ブログ(Diary)

arupaka-_-arupakaの日記

2018-09-13

多変量回帰と相関

多変量解析では、

(0)単独変数の調整

 (0-1)カテゴリ化すると楽 cut関数

(0-2)少なすぎる変数の除去

 

(1)説明変数間の相関を調べる

  →とりあえず、主要っぽい変数を10個くらいとりだすとわかりやすい

http://kusanagi.hatenablog.jp/entry/2014/07/23/220951

(1-1)ライブラリ

     library(psych)
     pairs.panels(dat)

cor.plot(tmp2)

→各変数がどういう情報をもっているか把握

(1-2)因子分析や主成分分析等で整理

 (1-3) 交互作用の調査

(1-4) 相関係数、ファイ、相関比を使う

 (1−5)3変数の依存

i tmp2<-d2
     d<-data.frame(y=tmp2$y,x1=as.numeric(paste(tmp2$MENSEKI)),x2=as.numeric(paste(tmp2$TIKUNENSUU)),x3=as.numeric(paste(tmp2$TINRYOU)))
     
     d2<-d%>% dplyr::filter( x1<=300 & x2<=100) %>% dplyr::group_by(x1=20*round(x1/20),x2=2*round(x2/2)) %>% dplyr:: summarise(z=length(y[!is.na(y)]),y=median(y,na.rm=T)) %>% dplyr::arrange(x2,x1)
     
     d3<-data.frame(d2[-3]) %>% tidyr::spread(key=x1,value=y)
     x_name<-d3[,1]
     d3[d3>=200]<-200
     #d3[d3<=10]<-10
     #image(x=1:length(d3[[1]]), y=as.numeric(names(d3[-1])),z=as.matrix(d3[,-1]),col=terrain.colors(100))
     #image.plot(x=1:length(d3[[1]]), y=as.numeric(names(d3[-1])),z=as.matrix(d3[,-1]),col=terrain.colors(100),zlim=c(0,200))
     image.plot(x=1:length(d3[[1]]), y=as.numeric(names(d3[-1])),z=as.matrix(d3[,-1]),col=terrain.colors(100),zlim=c(0,200),xlab="Ages",ylab="Area (m^2)")
      mtext(at=1:length(d3[[1]]),side=1,text=x_name)
     text(x=as.numeric(factor(d2[[2]])),y=d2[[1]],labels=d2[[3]])
     #contour(x=as.numeric(rownames(d3[,1])),y=as.numeric(colnames(d3[1,])[-1]),as.matrix(d3),nlevel=10,add=T)
     contour(x=1:length(d3[[1]]),y=as.numeric(names(d3[-1])),z=as.matrix(d3[,-1]),nlevel=10,add=T)
     ###

pairs.panels(tmp[,1:10])

fa.parallel(tmp[,1:10])

fa(r=tmp2[-1,-1],nfactors=7)

cv.test(y, v1)

interaction.plot

(2)説明変数と被説明変数の関係

(2-1) とりあえず、すべてカテゴリ変数化すると非線形までいける

(2-2) glm を使う anova 等で説明力チェック

(2-3) stan 等でベイズ的に調整

(例1):隣同士は同じあたりb(t+1)=b(t)+e(t))

(例2): 事前分布→データ数少ない場合は平均でおきかえ

(2-4)カテゴリ変数は一つ0にしなきゃいけない。

(2-5)ラッソやグループラッソ。

(2-6)ダミー変数化にはdumiesパッケージが便利

2017-10-09

google trend データの収集

(1)python を導入 参考:https://qiita.com/yuta_h3/items/2988c4d0811bf8c344c0

git clone git://github.com/yyuu/pyenv.git ~/.pyenv

~/.bash_profile に以下を書きこむ

export PYENV_ROOT="${HOME}/.pyenv"

if [ -d "${PYENV_ROOT}" ]; then

export PATH=${PYENV_ROOT}/bin:$PATH

eval "$(pyenv init -)"

fi

設定の読み込み

source ~/.bash_profile

インストールできる一覧を確認

pyenv install --list

バージョンを指定してインストール

pyenv install 3.6.3

そのディレクトリ以下でバージョンをデフォルトに。

pyenv local 3.6.3

(2)pytrendの導入

pytrends でできる。

pip install pytrends

あとは、git の exampleをみるとOK

*日本語をやるときは、ロケールに注意。プリントできなくてとまる。

*URLエンコードは必要なし。

#--*-- encode: utf-8 --*--

from pytrends.request import TrendReq

import sys

# Login to Google. Only need to run this once, the rest of requests will use the same session.

pytrend = TrendReq()


pytrend.build_payload(kw_list=[u'駅', 'bagel'])

# Interest Over Time

interest_over_time_df = pytrend.interest_over_time()

print(interest_over_time_df)

~

2017-05-21

Rでスペクトル解析

Rでスペクトル解析。

正弦波以外の周期をもった波が倍音が出力されることに注意

例1)正弦波 30日周期

x<-sin((1:3000)/30*2*pi)

plot(x)

spec.pgram(x,c(10,10))


例2)30日周期のイベント(倍音がでる)

x<-rep(c(rep(0,30)+3*dnorm(1:30,mean=15,sd=2)),100)

plot(x)

spec.pgram(x,c(10,10))

https://jp.mathworks.com/help/signal/examples/practical-introduction-to-frequency-domain-analysis.html

https://jp.mathworks.com/help/signal/examples/analyzing-harmonic-distortion.html

2017-03-16

Rでパラメータ制約付き非線形最適化とパラメータ制約付の最小二乗法とロジスティック回帰(カテゴリカルデータ、質的データ)

  • Rでパラメータ制約付非線形最適化

Rsolnpパッケージをもちいる。

Rで非線形制約条件付非線形最適化を行う方法 by Rsolnp

http://d.hatena.ne.jp/teramonagi/20091217/1261048574を参考に

  • とりあえず最小二乗法でためす.

library(Rsolnp)

#答えと問題のデータを作成

#y=5*x1-8*x2-3*x3

a1<-5

a2<--8

a3<--3

x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

x3<-rnorm(1000,mean=0,sd=1)


y<-a1*x1+a2*x2+a3*x3+1*rnorm(1000,mean=0,sd=1)


###最小二乗法の目的関数

eq1 <- function(x)

{

a1<-x[1]

a2<-x[2]

a3<-x[3]

return(sum*1^2))

}

#係数a1+a2=0の非線形制約

eq2 <- function(x)

{

return((x[1]+x[2])^1)

# return(x_[1]^4)

}

#初期パラメータ

x0 <- c(1,-1,1)

#最適化 a1+a2=0の制約で最適化. eqB=0 は制約式=0という意味

solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

#solution <- solnp(x0,fun = objectiveFunc)

#結果出力

print(solution)

  • ロジスティック回帰(制約なし)

a1<-5

a2<--8

a3<--3

x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

x3<-rnorm(1000,mean=0,sd=1)


y<-a1*x1+a2*x2+a3*x3+1*rnorm(1000,mean=0,sd=1)


prob<-1/(1+exp(-y))

z<-rbinom(length(prob),prob=prob,size=1)

#対数尤度

log_li<-function(a){

y<-x1*a[1]+x2*a[2]+x3*a[3]+a[4]

p<-1/(1+exp(-y))

-sum(z*log(p)+(1-z)*log(1-p))

}


#決定係数のスタート値

x0 <- c(1,1,1,0)

#最適化

#solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

solution <- solnp(x0,fun = log_li)

#結果出力

print(solution$par)

出力

print(solution$par)

[1] 4.50399000 -6.65304329 -2.54410502 -0.01038315


一応、標準のglmと比較

d<-data.frame(z,x1,x2,x3)

glm(z~.,d,family=binomial(link="logit"))

#Call: glm(formula = z ~ ., family = binomial(link = "logit"), data = d)

#Coefficients:

#(Intercept) x1 x2 x3

# -0.01038 4.50399 -6.65304 -2.54410

#

#Degrees of Freedom: 999 Total (i.e. Null); 996 Residual

#Null Deviance: 1386

#Residual Deviance: 294 AIC: 302

#>

だいたい同じ結果.

  • カテゴリカルデータ(質的データ)のロジスティック回帰(パラメータ制約なし,不定形)
    • 不定形のため初期値により最適解がかわる。
      • x3とx4が線形従属している。x3:男性1,女性0, x4:女性1, 男性0

a1<-5

a2<--8

a3<--3

a4<-2

x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

##Man flg

x3<-rbinom(1000,prob=0.1,size=1)

##Woman flg

x4<-1-x3

y<-a1*x1+a2*x2+a3*x3+a4*x4+1*rnorm(1000,mean=0,sd=1)


prob<-1/(1+exp(-y))

z<-rbinom(length(prob),prob=prob,size=1)

log_li<-function(a){

y<-x1*a[1]+x2*a[2]+x3*a[3]+x4*a[4]+a[5]

p<-1/(1+exp(-y))

-sum(z*log(p)+(1-z)*log(1-p))

}

#初期値を変えて、不定性を確認

x0a <- c(1,1,1,1,0)

x0b<-c(3,4,2,2,1.764)

#最適化

#solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

solution1 <- solnp(pars=x0a,fun = log_li)

solution2 <- solnp(pars=x0b,fun = log_li)

#結果出力

print(solution1$par)

print(solution2$par)

#glmと比較

d<-data.frame(z,x1,x2,x3,x4)

glm(z~.,d,family=binomial(link="logit"))

    • 出力結果

> #結果出力

> print(solution1$par)

[1] 4.352484439 -7.153849857 -3.101372925 1.684972724 -0.002192511

> print(solution2$par)

[1] 4.352489 -7.153857 -1.877847 2.908487 -1.225704

>

> d<-data.frame(z,x1,x2,x3,x4)

> glm(z~.,d,family=binomial(link="logit"))

Call: glm(formula = z ~ ., family = binomial(link = "logit"), data = d)

Coefficients:

(Intercept) x1 x2 x3 x4

1.683 4.352 -7.154 -4.786 NA

Degrees of Freedom: 999 Total (i.e. Null); 996 Residual

Null Deviance: 1378

Residual Deviance: 312.5 AIC: 320.5

初期値によって解が違うことがわかる.glmあx4を自動的に捨てている。


  • 例2:制約を入れてみる.まず,切片に制約をいれてみる。gslの結果にあわせる

#####Logistic regression with Categolical data and constraints####

####Example 1: Fixed thresholds

a1<-5

a2<--8

a3<--3

a4<-2

x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

##Man flg

x3<-rbinom(1000,prob=0.1,size=1)

##Woman flg

x4<-1-x3

y<-a1*x1+a2*x2+a3*x3+a4*x4+1*rnorm(1000,mean=0,sd=1)


prob<-1/(1+exp(-y))

z<-rbinom(length(prob),prob=prob,size=1)

log_li<-function(a){

y<-x1*a[1]+x2*a[2]+x3*a[3]+x4*a[4]+a[5]

p<-1/(1+exp(-y))

-sum(z*log(p)+(1-z)*log(1-p))

}

constraint <- function(x)

{

return((x[5]-1.683)^1)

# return(x_[1]^4)

}

#決定係数のスタート値

x0a <- c(1,1,1,1,0)

x0b<-c(3,4,2,2,1.764)

#最適化

#solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

solution1 <- solnp(pars=x0a,fun = log_li,eqfun = constraint,eqB =0)

solution2 <- solnp(pars=x0b,fun = log_li,eqfun = constraint,eqB =0)

#結果出力

print(solution1$par)

print(solution2$par)

d<-data.frame(z,x1,x2,x3,x4)

glm(z~.,d,family=binomial(link="logit"))

出力結果

> print(solution1$par)

[1] 4.2505029 -6.8385070 -4.4021475 -0.0818387 1.6830000

> print(solution2$par)

[1] 4.25048803 -6.83826458 -4.40227303 -0.08189154 1.68300000

> d<-data.frame(z,x1,x2,x3,x4)

> glm(z~.,d,family=binomial(link="logit"))

Call: glm(formula = z ~ ., family = binomial(link = "logit"), data = d)

Coefficients:

(Intercept) x1 x2 x3 x4

1.601 4.250 -6.838 -4.320 NA

予想通りglmと同じ結果が初期値に依存せずにえれた.

    • 例3 制約:「a3=-a4」を入れてみる(男性=-女性).

a1<-5

a2<--8

a3<--3

a4<-2

x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

##Man flg

x3<-rbinom(1000,prob=0.1,size=1)

##Woman flg

x4<-1-x3

y<-a1*x1+a2*x2+a3*x3+a4*x4+1*rnorm(1000,mean=0,sd=1)


prob<-1/(1+exp(-y))

z<-rbinom(length(prob),prob=prob,size=1)

log_li<-function(a){

y<-x1*a[1]+x2*a[2]+x3*a[3]+x4*a[4]+a[5]

p<-1/(1+exp(-y))

-sum(z*log(p)+(1-z)*log(1-p))

}

constraint <- function(a)

{

return*2

}

#決定係数のスタート値

x0a <- c(1,1,1,1,0)

x0b<-c(3,4,2,2,1.764)

#最適化

#solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

solution1 <- solnp(pars=x0a,fun = log_li,eqfun = constraint,eqB =0)

solution2 <- solnp(pars=x0b,fun = log_li,eqfun = constraint,eqB =0)

#結果出力

print(solution1$par)

print(solution2$par)

#出力 does not depends on the innitail value

#> print(solution1$par)

#[1] 4.4599154 -6.9335231 -2.5641968 2.5641968 -0.9094139

#> print(solution2$par)

#[1] 4.4599725 -6.9335964 -2.5644191 2.5644191 -0.9096364

予想通り解が初期値に依存しない.

    • 例4 平均0制約:「a3*p3+a4*p4=0」を入れてみる(ランダムサンプルして平均をとると男女効果は0になる).
      • すべてにこの制約をいれれば、切片が平均になる。

a1<-5

a2<--8

a3<--3

a4<-2

x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

##Man flg

x3<-rbinom(1000,prob=0.1,size=1)

##Woman flg

x4<-1-x3

y<-a1*x1+a2*x2+a3*x3+a4*x4+1*rnorm(1000,mean=0,sd=1)


prob<-1/(1+exp(-y))

z<-rbinom(length(prob),prob=prob,size=1)

log_li<-function(a){

y<-x1*a[1]+x2*a[2]+x3*a[3]+x4*a[4]+a[5]

p<-1/(1+exp(-y))

-sum(z*log(p)+(1-z)*log(1-p))

}

constraint <- function(a)

{

#return*3

pp3<-sum(x3)/length(x3)

pp4<-sum(x4)/length(x4)

return*4

}

#決定係数のスタート値

x0a <- c(1,1,1,1,0)

x0b<-c(3,4,2,2,1.764)

#最適化

#solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

solution1 <- solnp(pars=x0a,fun = log_li,eqfun = constraint,eqB =0)

solution2 <- solnp(pars=x0b,fun = log_li,eqfun = constraint,eqB =0)

#結果出力

print(solution1$par)

print(solution2$par)


出力

> #結果出力

> print(solution1$par)

[1] 4.3169926 -6.6144392 -3.4385721 0.3693704 1.1827882

> print(solution2$par)

[1] 4.3704808 -6.6988868 -3.4753002 0.3733157 1.2258578


  • 係数ax+ayの平均が0制約. 回帰の分散が係数の説明力より十分小さいとき,係数をy=1/(1+exp(-mu))

で変換したものがyの平均確率(p=(p1+p2+..+pN)/N)に対応する.

#####Logistic regression with Categolical data and constraints####

####Example 3: p3*a3+p4*a4=0 Mean zero constraionts

####The case of var[a1x1+a1x2] << mu0, p=1(1+exp(-mu0));p~mean(z)

a1<-5

a2<--8

a3<--3

a4<-2

a1<-0.0309115

a2<--0.07597877

a3<--0.1681380

a4<-0.137884

#a1<-0

#a2<-0

#a3<-0

#a4<-0

a5<-0.4594540


x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

##Man flg

x3<-rbinom(1000,prob=0.1,size=1)

##Woman flg

x4<-1-x3

x1<-scale(x1)

x2<-scale(x2)


#y<-a1*x1+a2*x2+a3*x3+a4*x4+1*rnorm(1000,mean=0,sd=1)+a5

y<-a1*x1+a2*x2+a3*x3+a4*x4+1*rnorm(1000,mean=0,sd=1)+a5


prob<-1/(1+exp(-y))

z<-rbinom(length(prob),prob=prob,size=1)

log_li<-function(a){

y<-x1*a[1]+x2*a[2]+x3*a[3]+x4*a[4]+a[5]

p<-1/(1+exp(-y))

-sum(z*log(p)+(1-z)*log(1-p))

}

constraint <- function(a)

{

#return*5

pp3<-sum(x3)/length(x3)

pp4<-sum(x4)/length(x4)

return*6

}

#決定係数のスタート値

x0a <- c(1,1,1,1,0)

x0b<-c(3,4,2,2,1.764)

#最適化

#solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

solution1 <- solnp(pars=x0a,fun = log_li,eqfun = constraint,eqB =0)

solution2 <- solnp(pars=x0b,fun = log_li,eqfun = constraint,eqB =0)

#結果出力

print(solution1$par)

print(solution2$par)

#> 1/(1+exp(-0.45194))

#[1] 0.6111004

#> mean(z)

#[1] 0.61

#>



http://documentation.statsoft.com/portals/0/formula%20guide/Logistic%20Regression%20Formula%20Guide.pdf

http://home.hiroshima-u.ac.jp/tkurita/lecture/statface/node13.html

回帰

http://ifs.nog.cc/gucchi24.hp.infoseek.co.jp/MRA2.htm

http://norimune.net/625

*1:y-(a1*x1+a2*x2+a3*x3

*2:a[3]+a[4]

*3:a[3]+a[4]

*4:a[3]*pp3+a[4]*pp4

*5:a[3]+a[4]

*6:a[3]*pp3+a[4]*pp4

2017-01-18

Rで因子分析と分散共分散分析

Rで因子分析と分散共分散分析をする。

  • 因子分析は、潜在要因をみつける方法。
    • 主成分分析と回帰分析と異なり観測変数xを説明する潜在変数をみつける方法(下のx1とx2を見つける)

y1=a*x1+a*x2

y2=a*x1+a*x2

y3=a*x1+a*x2

y4=a*x1+a*x2

    • 心理学等では、「瞬発筋力x1」や「持久力x2」や「器用さx3」など観測できない構成概念を扱うことが重要なためこの方法が発達した。
      • 観測できるのは、持久走のタイムy1、短距離走のタイムy2、ソフトボールなど記録y3、体前屈の記録y4、身長y5、体重y6など潜在要因からきまってくるもの
    • ただし、係数の決め方や変数の決め方に不定性(同じパラメータでおなじ説明力)があるので、軸の回転をおこなうことで、説明しやすい解や理解しやすい解をさぐる。
      • 各潜在要因が直交しない解がだすことができるのも因子分析の特徴
      • 自由に色々選べるので調整がききやすい
      • 代表的な自動回転法 promax回転(非直行)やvarimax回転(直行)。構造を目指す(出来る限り観測変数が一つの説明の要因で説明できるようにする回転。潜在要因と観測要因の関係が明確になるので説明しやすい)。
      • 要因の説明力をたかめると、個別データをバラバラにする野の違いもある(主成分と因子分析)
    • Rでは、fa関数(psychパッケージ)

参考:Rで因子分析 商用ソフトで実行できない因子分析のあれこれ

http://www.slideshare.net/simizu706/r-42283141

      • Step1 下処理と標準化 scale(X)関数.na.omit(X), 正規分布っぽくする.欠損値除去
      • Step2 因子数の自動推定(vss, parallel関数)

install.packages("psych")

library(psych)

#最低値がでる。BICがおすすめらしい

map <- vss(rr2, fm="minres")

#最大値がでる。BICがおすすめらしい

para <- fa.parallel(rr2)

        • ほかには主成分析→おりまがり+-1を調べる plot(eigen(cor(X))$values)
      • Step3 因子分析

#例:8要因でoblimin回転(promaxの進化版?)

v<-fa(rr3,nfactor=8,rotate = "oblimin")

#例:8要因でvarimax回転(直行)

v2<-fa(rr3,nfactor=8,rotate = "varimax")

#例:8要因でpromax回転(非直行)

v3<-fa(rr3,nfactor=8,rotate = "promax")

      • Step4: 結果の取り出し(各潜在要因のうち因子スコアが大きい観測量をとりだし)。

for(i in 1:8){

cat(labels*12[i],"\n")

s<-(v2$loadings[,i][rev(order(abs(v2$loadings[,i])))][1:10])

print(s)

#print(paste(names(s),collapse=", "))

#print(paste(sprintf("%.3f",s),collapse=", "))

}

      • RMSEA(root mean square error of approximation) 0.05 OK 0.1> No
      • 共通因子スコアと独自因子スコア

#install.pakcages("lavaan")

#install.packages("semPlot")

#install.packages("knitr")

#install.packages("data.table")

#library(semPlot)

library(semPlot)

library(lavaan)

    • (2)練習1
      • 例1:線形回帰
        • 構造方程式(xがyを説明する場合)はy~x
        • 観測方程式(潜在要因fが観測量xを説明する場合)はy=~xとかく
        • 相関なしは潜在要因f1とf2が無相関はf1~0*f2とする

#通常の回帰

x<-rnorm(1000,mean=0,sd=2)

y<-3*x+rnorm(1000,mean=0,sd=0.1)

model1<-"y~x"

X<-cbind(x,y)

getCov(X)

#データ直接入力モード

fit<-sem(model=model1,data=X,estimator="ML")

#サマリー

summary(fit)

#サマリー、精度評価指標つき

summary(fit,fit.measure=TRUE)

#解析結果のグラフ化

semPaths(fit,whatLabels="stand",optimizeLatRes=T)

#解析結果のとりだし

parameterEstimates(object=fit)

#解析結果のとりだし標準化

standardizedSolution(object=fit)

      • 例3:二潜在要因を探す(体重(g)の列だけ除いて解析[相関係数1の要因があるとcov(X)が特異になり正定値行列でなく解析不明になるため])

#体重gの列を除く(体重(kg)と相関係数1の関係)

rr2b<-(rr2[,-which(colnames(rr2)=="体重(g)")])

#分散共分散行列

r_d<-cov(rr2b)

r_d3<-r_d

#データのラベルを取得

n1<-labels(rr2b)2

#すべてのデータラベルの観測方程式の(y1=~長距離走; y1=~ 短距離走; y1~=ソフトボールなげなど)を作る。

model_test<-paste("y1=~",n1,sep="")

#すべてのデータラベルの観測方程式の(y2=~長距離走; y2=~ 短距離走; y2~=ソフトボールなげなど)を作る。

model_testb<-paste("y2=~",n1,sep="")

#潜在要因y1とy2は無相関という条件をいれる

model_testc<-'y1~0*y2'

#モデルをまとめる

model_test2<-paste(c(model_test,model_testb,model_testc),collapse="\n")

#データ数の取得

data_no<-length(rr2b[,1])

#SEMを実行(分散共分散行列モード、データ数の指定が必要)

fit<-sem(model=model_test2,sample.cov=r_d3,estimator="ML",sample.nobs=data_no)

#semPaths(fit,whatLabels="stand",optimizeLatRes=T,nCharNodes=12)

#parameterEstimates(object=fit)

#標準化回帰係数の取り出し

standardizedSolution(object=fit)

    • 実際の解析
      • とりあえず、探索的因子分析(因子分析の結果)を利用して観測方程式を作ってみる。
      • 潜在要因をさらに朝ごはんや学校成績で説明するモデル(構造方程式)

長距離走<-(持久力),(瞬発力)

短距離走<-(瞬発力)

ソフトボール投げ<-(瞬発力)

朝ごはん->(持久力)

学校成績->(持久力)

運動習慣->(持久力)

朝ごはん->(瞬発力)

学校成績->(瞬発力)

運動習慣->(瞬発力)

朝ごはん->(筋力)

学校成績->(筋力)

運動習慣->(筋力)

*1:v3$loadings