授業で使用するいいデータがなかったので，作ってしまいました。
作るのにはpsychパッケージを使います。
library(psych)
# 乱数の種
set.seed(17)
# モデルの係数を指定
fx <- matrix(c(0.9,0.8,0.7,rep(0,9),0.7,0.6,0.7,rep(0,9),0.8,0.7,0.6),ncol=3)
fy <- c(0.8,0.7,0.9)
Phi <- matrix(c(1,0.4,0.3,0.2,0.4,1,0.2,0.6,0.3,0.2,1,0.3,0.2,0.6,0.3,1),ncol=4)
# サンプルデータを作成
sample <- sim.structure(fx,Phi,fy)## Loading required package: MASSsample$model## Vx1 Vx2 Vx3 Vx4 Vx5 Vx6 Vx7 Vx8 Vx9 Vy1 Vy2
## Vx1 1.000 0.720 0.630 0.252 0.216 0.252 0.216 0.189 0.162 0.144 0.126
## Vx2 0.720 1.000 0.560 0.224 0.192 0.224 0.192 0.168 0.144 0.128 0.112
## Vx3 0.630 0.560 1.000 0.196 0.168 0.196 0.168 0.147 0.126 0.112 0.098
## Vx4 0.252 0.224 0.196 1.000 0.420 0.490 0.112 0.098 0.084 0.336 0.294
## Vx5 0.216 0.192 0.168 0.420 1.000 0.420 0.096 0.084 0.072 0.288 0.252
## Vx6 0.252 0.224 0.196 0.490 0.420 1.000 0.112 0.098 0.084 0.336 0.294
## Vx7 0.216 0.192 0.168 0.112 0.096 0.112 1.000 0.560 0.480 0.192 0.168
## Vx8 0.189 0.168 0.147 0.098 0.084 0.098 0.560 1.000 0.420 0.168 0.147
## Vx9 0.162 0.144 0.126 0.084 0.072 0.084 0.480 0.420 1.000 0.144 0.126
## Vy1 0.144 0.128 0.112 0.336 0.288 0.336 0.192 0.168 0.144 1.000 0.560
## Vy2 0.126 0.112 0.098 0.294 0.252 0.294 0.168 0.147 0.126 0.560 1.000
## Vy3 0.162 0.144 0.126 0.378 0.324 0.378 0.216 0.189 0.162 0.720 0.630
## Vy3
## Vx1 0.162
## Vx2 0.144
## Vx3 0.126
## Vx4 0.378
## Vx5 0.324
## Vx6 0.378
## Vx7 0.216
## Vx8 0.189
## Vx9 0.162
## Vy1 0.720
## Vy2 0.630
## Vy3 1.000# モデル図
structure.diagram(fx,Phi,fy,cut=0)さて，ここからデータセットを作ります。自作関数を準備(この辺の話は http://kosugitti.net/archives/112 を参照)。
# 任意の相関行列をもつデータセットを作る関数
corMatVal <- function(size,n.vals,cor.mat,corvec=FALSE){
Rmat <- matrix(rnorm(n=size*n.vals),nrow=size)
#if given cor.mat is just vector,make correlational matrix
if(corvec==TRUE){
temp.mat <- diag(n.vals)
temp.mat[(upper.tri)(temp.mat,diag=FALSE)] <- cor.mat
temp.cor <- (t(temp.mat)+temp.mat)
diag(temp.cor) <- 1
}else{
temp.cor <- cor.mat
}
Upper <- chol(temp.cor)
result <- as.data.frame(Rmat %*% Upper)
return(result)
}この関数は，サンプルサイズ，変数の数，行列を渡してやるとデータセットを作るというもの。使い方は以下の通り。
sampledata <- corMatVal(500,12,sample$model)
summary(sampledata)## Vx1 Vx2 Vx3
## Min. :-3.556136 Min. :-3.267431 Min. :-2.88423
## 1st Qu.:-0.714150 1st Qu.:-0.588818 1st Qu.:-0.78146
## Median :-0.000591 Median : 0.009451 Median :-0.09036
## Mean :-0.029510 Mean : 0.023383 Mean :-0.04252
## 3rd Qu.: 0.673071 3rd Qu.: 0.659493 3rd Qu.: 0.65921
## Max. : 3.330153 Max. : 4.258519 Max. : 2.98409
## Vx4 Vx5 Vx6
## Min. :-2.77290 Min. :-3.02830 Min. :-3.172601
## 1st Qu.:-0.62570 1st Qu.:-0.53395 1st Qu.:-0.575488
## Median : 0.02396 Median : 0.05308 Median : 0.001294
## Mean : 0.02385 Mean : 0.04205 Mean : 0.012784
## 3rd Qu.: 0.62457 3rd Qu.: 0.62317 3rd Qu.: 0.684504
## Max. : 2.98966 Max. : 2.78232 Max. : 2.418963
## Vx7 Vx8 Vx9
## Min. :-2.8903 Min. :-3.2256 Min. :-2.41139
## 1st Qu.:-0.7673 1st Qu.:-0.7430 1st Qu.:-0.73895
## Median :-0.1329 Median :-0.1395 Median :-0.08882
## Mean :-0.1276 Mean :-0.1055 Mean :-0.10181
## 3rd Qu.: 0.5445 3rd Qu.: 0.5705 3rd Qu.: 0.60653
## Max. : 2.3493 Max. : 2.7647 Max. : 2.89852
## Vy1 Vy2 Vy3
## Min. :-3.843148 Min. :-3.13183 Min. :-2.85006
## 1st Qu.:-0.670198 1st Qu.:-0.60068 1st Qu.:-0.64280
## Median :-0.047808 Median : 0.06205 Median : 0.03828
## Mean :-0.000494 Mean : 0.06736 Mean : 0.02538
## 3rd Qu.: 0.671363 3rd Qu.: 0.75519 3rd Qu.: 0.70597
## Max. : 2.968449 Max. : 3.82201 Max. : 4.18128print(cor(sampledata),digit=2)## Vx1 Vx2 Vx3 Vx4 Vx5 Vx6 Vx7 Vx8 Vx9 Vy1 Vy2
## Vx1 1.00 0.731 0.634 0.203 0.191 2.0e-01 1.7e-01 0.217 0.106 0.101 0.16
## Vx2 0.73 1.000 0.597 0.234 0.161 2.1e-01 1.4e-01 0.194 0.118 0.099 0.13
## Vx3 0.63 0.597 1.000 0.217 0.125 2.0e-01 1.4e-01 0.204 0.114 0.099 0.10
## Vx4 0.20 0.234 0.217 1.000 0.321 4.5e-01 7.4e-02 0.075 0.043 0.293 0.25
## Vx5 0.19 0.161 0.125 0.321 1.000 4.0e-01 6.2e-02 0.054 0.036 0.267 0.19
## Vx6 0.20 0.209 0.195 0.454 0.400 1.0e+00 -1.1e-06 0.029 -0.072 0.284 0.25
## Vx7 0.17 0.140 0.140 0.074 0.062 -1.1e-06 1.0e+00 0.520 0.524 0.198 0.20
## Vx8 0.22 0.194 0.204 0.075 0.054 2.9e-02 5.2e-01 1.000 0.417 0.192 0.22
## Vx9 0.11 0.118 0.114 0.043 0.036 -7.2e-02 5.2e-01 0.417 1.000 0.154 0.15
## Vy1 0.10 0.099 0.099 0.293 0.267 2.8e-01 2.0e-01 0.192 0.154 1.000 0.53
## Vy2 0.16 0.127 0.104 0.253 0.194 2.5e-01 2.0e-01 0.216 0.150 0.530 1.00
## Vy3 0.12 0.094 0.080 0.373 0.258 3.7e-01 2.2e-01 0.196 0.149 0.687 0.65
## Vy3
## Vx1 0.120
## Vx2 0.094
## Vx3 0.080
## Vx4 0.373
## Vx5 0.258
## Vx6 0.371
## Vx7 0.219
## Vx8 0.196
## Vx9 0.149
## Vy1 0.687
## Vy2 0.654
## Vy3 1.000これで完成でもいいんですが，カバーストーリーとして基礎体力が野球のパフォーマンスにどういう影響があるか，というものを考えているので，適当な変数名をつけたり，単位をそれっぽいものにしたりします。 ここで丸め誤差などが生じてくるので，理想的なモデルが再現できるとは限りません。
colnames(sampledata)<-c("Run50M","Run500M","Run1500M","Kensui","Entoh","Udetate","Hanpuku","Kiritsu","Habatobi","Hit","Homerun","Choda")
sampledata$Run50M <- round(sampledata$Run50M+8,3)
sampledata$Run500M <- round(sampledata$Run500M*10,3)+120
sampledata$Run1500M <- round(sampledata$Run1500M*10,3)+240
sampledata$Kensui <- round(sampledata$Kensui*5+30)
sampledata$Entoh <- round(sampledata$Entoh*5+30)
sampledata$Udetate <- round(sampledata$Udetate*10+34)
sampledata$Hanpuku <- round(sampledata$Hanpuku*10+30)
sampledata$Kiritsu <- round(sampledata$Kiritsu+5,3)
sampledata$Habatobi <- round(sampledata$Habatobi+10,2)
sampledata$Hit <- round(sampledata$Hit*10+100)
sampledata$Homerun <- round(sampledata$Homerun*2+20)
sampledata$Choda <- round(sampledata$Choda*10+35)
summary(sampledata)## Run50M Run500M Run1500M Kensui
## Min. : 4.444 Min. : 87.33 Min. :211.2 Min. :16.00
## 1st Qu.: 7.286 1st Qu.:114.11 1st Qu.:232.2 1st Qu.:27.00
## Median : 7.999 Median :120.09 Median :239.1 Median :30.00
## Mean : 7.970 Mean :120.23 Mean :239.6 Mean :30.12
## 3rd Qu.: 8.673 3rd Qu.:126.59 3rd Qu.:246.6 3rd Qu.:33.00
## Max. :11.330 Max. :162.59 Max. :269.8 Max. :45.00
## Entoh Udetate Hanpuku Kiritsu
## Min. :15.00 Min. : 2.00 Min. : 1.00 Min. :1.774
## 1st Qu.:27.00 1st Qu.:28.00 1st Qu.:22.00 1st Qu.:4.257
## Median :30.00 Median :34.00 Median :29.00 Median :4.861
## Mean :30.23 Mean :34.12 Mean :28.72 Mean :4.894
## 3rd Qu.:33.00 3rd Qu.:41.00 3rd Qu.:35.00 3rd Qu.:5.571
## Max. :44.00 Max. :58.00 Max. :53.00 Max. :7.765
## Habatobi Hit Homerun Choda
## Min. : 7.590 Min. : 62.00 Min. :14.00 Min. : 6.00
## 1st Qu.: 9.260 1st Qu.: 93.00 1st Qu.:19.00 1st Qu.:29.00
## Median : 9.910 Median :100.00 Median :20.00 Median :35.00
## Mean : 9.898 Mean : 99.98 Mean :20.14 Mean :35.25
## 3rd Qu.:10.605 3rd Qu.:107.00 3rd Qu.:22.00 3rd Qu.:42.00
## Max. :12.900 Max. :130.00 Max. :28.00 Max. :77.00一応，丸まったデータがどういう推定値を出すか，見てみます。SEMにはlavaanパッケージを使います。
library(lavaan)## This is lavaan 0.5-17
## lavaan is BETA software! Please report any bugs.model <- '
Ashi =~ Run50M + Run500M + Run1500M
Ude =~ Kensui + Entoh + Udetate
Shunbin =~ Hanpuku+Kiritsu+Habatobi
Performance =~ Hit + Homerun + Choda
Performance ~ Ashi + Ude + Shunbin
'
fit <- sem(model,data=sampledata)
summary(fit,standardized=TRUE)## lavaan (0.5-17) converged normally after 136 iterations
##
## Number of observations 500
##
## Estimator ML
## Minimum Function Test Statistic 44.970
## Degrees of freedom 48
## P-value (Chi-square) 0.598
##
## Parameter estimates:
##
## Information Expected
## Standard Errors Standard
##
## Estimate Std.err Z-value P(>|z|) Std.lv Std.all
## Latent variables:
## Ashi =~
## Run50M 1.000 0.928 0.876
## Run500M 8.839 0.457 19.330 0.000 8.205 0.832
## Run1500M 8.067 0.471 17.120 0.000 7.489 0.723
## Ude =~
## Kensui 1.000 3.161 0.646
## Entoh 0.769 0.087 8.814 0.000 2.430 0.521
## Udetate 2.179 0.211 10.348 0.000 6.889 0.716
## Shunbin =~
## Hanpuku 1.000 7.312 0.789
## Kiritsu 0.086 0.008 11.281 0.000 0.631 0.662
## Habatobi 0.087 0.008 11.182 0.000 0.636 0.648
## Performance =~
## Hit 1.000 7.259 0.746
## Homerun 0.188 0.012 15.180 0.000 1.363 0.697
## Choda 1.265 0.073 17.437 0.000 9.181 0.920
##
## Regressions:
## Performance ~
## Ashi -1.215 0.433 -2.804 0.005 -0.155 -0.155
## Ude 1.407 0.174 8.071 0.000 0.613 0.613
## Shunbin 0.313 0.056 5.617 0.000 0.315 0.315
##
## Covariances:
## Ashi ~~
## Ude 1.071 0.190 5.645 0.000 0.365 0.365
## Shunbin 1.762 0.389 4.526 0.000 0.260 0.260
## Ude ~~
## Shunbin 1.396 1.442 0.969 0.333 0.060 0.060
##
## Variances:
## Run50M 0.261 0.037 0.261 0.232
## Run500M 29.861 3.189 29.861 0.307
## Run1500M 51.250 3.926 51.250 0.477
## Kensui 13.984 1.223 13.984 0.583
## Entoh 15.847 1.165 15.847 0.729
## Udetate 45.118 4.789 45.118 0.487
## Hanpuku 32.370 4.410 32.370 0.377
## Kiritsu 0.510 0.044 0.510 0.561
## Habatobi 0.558 0.047 0.558 0.580
## Hit 42.043 3.427 42.043 0.444
## Homerun 1.966 0.147 1.966 0.514
## Choda 15.341 3.440 15.341 0.154
## Ashi 0.862 0.076 1.000 1.000
## Ude 9.994 1.494 1.000 1.000
## Shunbin 53.464 6.367 1.000 1.000
## Performance 30.163 3.925 0.572 0.572これでオッケー。
実は，こんな面倒なことをしなくても，sim.structure関数で仮想データを作ることも可能なのです。
# 5件法のデータセットをN=500で作る(itemsオプションをTに)
sample <- sim.structure(fx,Phi,fy,n=500,items=T,cat=5)
sample## Call: sim.structure(fx = fx, Phi = Phi, fy = fy, n = 500, items = T,
## cat = 5)
##
## $model (Population correlation matrix)
## Vx1 Vx2 Vx3 Vx4 Vx5 Vx6 Vx7 Vx8 Vx9 Vy1 Vy2 Vy3
## Vx1 1.00 0.72 0.630 0.252 0.216 0.252 0.216 0.189 0.162 0.14 0.126 0.16
## Vx2 0.72 1.00 0.560 0.224 0.192 0.224 0.192 0.168 0.144 0.13 0.112 0.14
## Vx3 0.63 0.56 1.000 0.196 0.168 0.196 0.168 0.147 0.126 0.11 0.098 0.13
## Vx4 0.25 0.22 0.196 1.000 0.420 0.490 0.112 0.098 0.084 0.34 0.294 0.38
## Vx5 0.22 0.19 0.168 0.420 1.000 0.420 0.096 0.084 0.072 0.29 0.252 0.32
## Vx6 0.25 0.22 0.196 0.490 0.420 1.000 0.112 0.098 0.084 0.34 0.294 0.38
## Vx7 0.22 0.19 0.168 0.112 0.096 0.112 1.000 0.560 0.480 0.19 0.168 0.22
## Vx8 0.19 0.17 0.147 0.098 0.084 0.098 0.560 1.000 0.420 0.17 0.147 0.19
## Vx9 0.16 0.14 0.126 0.084 0.072 0.084 0.480 0.420 1.000 0.14 0.126 0.16
## Vy1 0.14 0.13 0.112 0.336 0.288 0.336 0.192 0.168 0.144 1.00 0.560 0.72
## Vy2 0.13 0.11 0.098 0.294 0.252 0.294 0.168 0.147 0.126 0.56 1.000 0.63
## Vy3 0.16 0.14 0.126 0.378 0.324 0.378 0.216 0.189 0.162 0.72 0.630 1.00
##
## $reliability (population reliability)
## Vx1 Vx2 Vx3 Vx4 Vx5 Vx6 Vx7 Vx8 Vx9 Vy1 Vy2 Vy3
## 0.81 0.64 0.49 0.49 0.36 0.49 0.64 0.49 0.36 0.64 0.49 0.81
##
## $r (Sample correlation matrix for sample size = 500 )
## Vx1 Vx2 Vx3 Vx4 Vx5 Vx6 Vx7 Vx8 Vx9
## Vx1 1.000 0.37352 0.181 0.048 0.086 0.126 0.103 0.115 0.095
## Vx2 0.374 1.00000 0.285 0.086 0.124 0.137 0.018 0.089 0.112
## Vx3 0.181 0.28517 1.000 0.102 0.123 0.246 0.029 0.088 0.073
## Vx4 0.048 0.08556 0.102 1.000 0.183 0.131 0.049 -0.036 -0.058
## Vx5 0.086 0.12417 0.123 0.183 1.000 0.268 0.057 0.014 -0.108
## Vx6 0.126 0.13742 0.246 0.131 0.268 1.000 0.064 -0.015 -0.018
## Vx7 0.103 0.01794 0.029 0.049 0.057 0.064 1.000 0.237 0.135
## Vx8 0.115 0.08917 0.088 -0.036 0.014 -0.015 0.237 1.000 0.226
## Vx9 0.095 0.11240 0.073 -0.058 -0.108 -0.018 0.135 0.226 1.000
## Vy1 0.035 -0.00079 -0.023 0.125 0.130 0.099 0.106 0.071 -0.017
## Vy2 0.061 0.09508 0.065 0.113 0.060 0.142 0.060 0.116 -0.032
## Vy3 0.023 0.07774 0.090 0.136 0.105 0.123 -0.030 0.076 0.028
## Vy1 Vy2 Vy3
## Vx1 0.03491 0.061 0.023
## Vx2 -0.00079 0.095 0.078
## Vx3 -0.02332 0.065 0.090
## Vx4 0.12480 0.113 0.136
## Vx5 0.12989 0.060 0.105
## Vx6 0.09867 0.142 0.123
## Vx7 0.10578 0.060 -0.030
## Vx8 0.07095 0.116 0.076
## Vx9 -0.01669 -0.032 0.028
## Vy1 1.00000 0.205 0.183
## Vy2 0.20529 1.000 0.263
## Vy3 0.18276 0.263 1.000summary(sample$observed)## Vx1 Vx2 Vx3 Vx4
## Min. :0.000 Min. :0.000 Min. :0.000 Min. :0.000
## 1st Qu.:4.750 1st Qu.:1.000 1st Qu.:0.000 1st Qu.:5.000
## Median :5.000 Median :2.000 Median :0.000 Median :5.000
## Mean :4.588 Mean :2.442 Mean :0.422 Mean :4.644
## 3rd Qu.:5.000 3rd Qu.:4.000 3rd Qu.:1.000 3rd Qu.:5.000
## Max. :5.000 Max. :5.000 Max. :5.000 Max. :5.000
## Vx5 Vx6 Vx7 Vx8
## Min. :0.00 Min. :0.000 Min. :0.000 Min. :0.000
## 1st Qu.:1.00 1st Qu.:0.000 1st Qu.:5.000 1st Qu.:1.000
## Median :2.00 Median :0.000 Median :5.000 Median :2.000
## Mean :2.48 Mean :0.384 Mean :4.598 Mean :2.382
## 3rd Qu.:4.00 3rd Qu.:1.000 3rd Qu.:5.000 3rd Qu.:4.000
## Max. :5.00 Max. :5.000 Max. :5.000 Max. :5.000
## Vx9 Vy1 Vy2 Vy3
## Min. :0.000 Min. :0.000 Min. :0.000 Min. :0.000
## 1st Qu.:0.000 1st Qu.:5.000 1st Qu.:1.000 1st Qu.:0.000
## Median :0.000 Median :5.000 Median :3.000 Median :0.000
## Mean :0.342 Mean :4.672 Mean :2.626 Mean :0.412
## 3rd Qu.:0.000 3rd Qu.:5.000 3rd Qu.:4.000 3rd Qu.:1.000
## Max. :5.000 Max. :5.000 Max. :5.000 Max. :5.000まあ5件法っていってるのに0が入っていたり，全部の変数が同じスケールになってしまうという点が引っかかるのであれば，上の自作関数の法を使ってください。
あとは，こうやって作ったデータセットを保存して配布すれば終わり。
write.table(sampledata,"sports.csv",sep=",",quote=F,row.names=F,col.names=T)以上です。