Nothing
inst/doc/grmsem.R
In grmsem: Genetic-Relationship-Matrix Structural Equation Modelling (GRMSEM)
## ----eval = FALSE-------------------------------------------------------------
# > G <- grm.input("large.gcta.grm.gz")
# > G[1:3,1:3] #Relationships among the first three individuals
# [,1] [,2] [,3]
# [1,] 0.99354762 0.02328514 0.01644197
# [2,] 0.02328514 0.99406837 0.01021175
# [3,] 0.01644197 0.01021175 1.02751472
## ----eval = FALSE-------------------------------------------------------------
# > load("ph.large.RData")
# > ph.large[1:2,]
# Y1 Y2 Y3 Y4
# [1,] -0.7640819 -0.6016908 -0.3981901 -0.3169821
# [2,] -0.5099606 0.6671311 -1.3119328 -0.5601261
## ----path-cholesky, echo=FALSE, out.width="65%", fig.align='center', fig.cap="Quad-variate Cholesky decomposition model. Observed phenotypes (P) are represented by squares and factors by circles. Single-headed arrows (paths) define relationships between variables. The variance of latent variables is constrained to unit variance; this is omitted from the path diagram to improve clarity."----
knitr::include_graphics("model_quadvariate_path_diagram_Cholesky.png")
## ----path-indep, echo = FALSE, out.width="65%", fig.align='center', fig.cap="Quad-variate Independent Pathway model. Observed phenotypes (P) are represented by squares and factors by circles. Single-headed arrows (paths) define relationships between variables. The variance of latent variables is constrained to unit variance; this is omitted from the path diagram to improve clarity."----
knitr::include_graphics("model_quadvariate_path_diagram_IP.png")
## ----path-ipc, echo = FALSE, out.width="65%", fig.align='center', fig.cap="Quad-variate IPC model. Observed phenotypes (P) are represented by squares and factors by circles. Single-headed arrows (paths) define relationships between variables. The variance of latent variables is constrained to unit variance; this is omitted from the path diagram to improve clarity."----
knitr::include_graphics("model_quadvariate_path_diagram_IPC.png")
## ----eval = FALSE-------------------------------------------------------------
# #If files were downloaded externally, please load here:
# > load("G.small.RData")
# > load("ph.small.RData")
# > fit <- grmsem.fit(ph.small,
# G.small,
# model="Cholesky",
# LogL=TRUE,
# estSE=TRUE)
# > fit
# $model.in
# part label value freepar
# 1 a a11 0.6971323 1
# 2 a a21 -0.2357364 1
# 3 a a31 0.3048627 1
# 4 a a22 0.2067668 1
# 5 a a32 -0.3823634 1
# 6 a a33 -0.1557505 1
# 7 e e11 -0.4958634 1
# 8 e e21 -0.5924593 1
# 9 e e31 0.1700410 1
# 10 e e22 0.3570643 1
# 11 e e32 -0.3220280 1
# 12 e e33 0.5024512 1
#
# $model.fit$LL
# [1] -54.39143
#
# $model.out
# label estimates gradient se Z p
# 1 a11 -0.5618594 -4.490283e-05 0.18730897 -2.999640 2.702993e-03
# 2 a21 -0.5172780 6.160157e-05 0.25150203 -2.056755 3.970982e-02
# 3 a31 -0.6022167 -4.596719e-05 0.20040106 -3.005057 2.655308e-03
# 4 a22 -0.7791023 -1.154262e-05 0.13696258 -5.688432 1.282113e-08
# 5 a32 -0.3313097 6.118543e-05 0.14755775 -2.245289 2.474961e-02
# 6 a33 -0.4201447 7.389654e-05 0.07000431 -6.001697 1.952654e-09
# 7 e11 -0.8266200 2.460492e-04 0.11621420 -7.112900 1.136299e-12
# 8 e21 -0.2344285 -1.605268e-04 0.12464049 -1.880838 5.999398e-02
# 9 e31 -0.4124524 6.209974e-05 0.12534761 -3.290469 1.000207e-03
# 10 e22 0.3862229 -2.636815e-05 0.12634091 3.056990 2.235717e-03
# 11 e32 0.3410509 8.582604e-05 0.10625220 3.209824 1.328162e-03
# 12 e33 -0.2511355 1.911040e-04 0.05125420 -4.899803 9.593276e-07
#
# $k
# [1] 3
#
# $n
# [1] 300
#
# $n.obs
# [1] 100 100 100
#
# $n.ind
# [1] 100
#
# $model
# [1] "Cholesky"
## ----eval = FALSE-------------------------------------------------------------
# > var.fit <- grmsem.var(fit)
# > print(var.fit)
# $VA
# 1 2 3
# Y1 0.3156860 0.2906375 0.3383611
# Y2 0.2906375 0.8745770 0.5696376
# Y3 0.3383611 0.5696376 0.6489526
#
# $VA.se
# 1 2 3
# Y1 0.2104826 0.1832943 0.1850112
# Y2 0.1832943 0.2583423 0.2220882
# Y3 0.1850112 0.2220882 0.2227517
#
# $VE
# 1 2 3
# Y1 0.6833006 0.1937833 0.3409414
# Y2 0.1937833 0.2041249 0.2284123
# Y3 0.3409414 0.2284123 0.3495017
#
# $VE.se
# 1 2 3
# Y1 0.1921300 0.1180125 0.1387263
# Y2 0.1180125 0.1269780 0.1222845
# Y3 0.1387263 0.1222845 0.1373590
#
# $VP
# 1 2 3
# Y1 0.9989866 0.4844208 0.6793025
# Y2 0.4844208 1.0787019 0.7980499
# Y3 0.6793025 0.7980499 0.9984543
#
# $RG
# 1 2 3
# Y1 1.0000000 0.5531272 0.7475604
# Y2 0.5531272 1.0000000 0.7561242
# Y3 0.7475604 0.7561242 1.0000000
#
# $RG.se
# 1 2 3
# Y1 0.0000000 0.22348027 1.575650e-01
# Y2 0.2234803 0.00000000 8.697161e-02
# Y3 0.1575650 0.08697161 4.449797e-17
#
# $RE
# 1 2 3
# Y1 1.0000000 0.5188747 0.6976686
# Y2 0.5188747 1.0000000 0.8551589
# Y3 0.6976686 0.8551589 1.0000000
#
# $RE.se
# 1 2 3
# Y1 5.160947e-17 2.051003e-01 1.173883e-01
# Y2 2.051003e-01 5.610664e-17 8.730948e-02
# Y3 1.173883e-01 8.730948e-02 2.359273e-17
#
## ----eval = FALSE-------------------------------------------------------------
# > print(st.fit)
# $stand.model.out
# label stand.estimates stand.se Z p
# 1 a11 -0.5621443 0.17186763 -3.270798 1.072444e-03
# 2 a21 -0.4980504 0.22959865 -2.169222 3.006584e-02
# 3 a31 -0.6026826 0.17316262 -3.480443 5.005852e-04
# 4 a22 -0.7501425 0.12682263 -5.914895 3.320880e-09
# 5 a32 -0.3315661 0.14661954 -2.261404 2.373422e-02
# 6 a33 -0.4204698 0.07358378 -5.714164 1.102447e-08
# 7 e11 -0.8270392 0.11681963 -7.079625 1.445454e-12
# 8 e21 -0.2257147 0.12967922 -1.740562 8.176046e-02
# 9 e31 -0.4127715 0.13122602 -3.145500 1.658030e-03
# 10 e22 0.3718667 0.14164900 2.625269 8.658054e-03
# 11 e32 0.3413148 0.11390419 2.996507 2.730917e-03
# 12 e33 -0.2513298 0.05457774 -4.604987 4.124919e-06
#
## ----eval=FALSE---------------------------------------------------------------
# > st.var.fit <- grmsem.var(st.fit)
# > print(st.var.fit)
# $VA
# 1 2 3
# Y1 0.3160062 0.2799762 0.3387946
# Y2 0.2799762 0.8107680 0.5488882
# Y3 0.3387946 0.5488882 0.6499573
#
# $VA.se
# 1 2 3
# Y1 0.1932288 0.1619891 0.1613852
# Y2 0.1619891 0.1379212 0.1565957
# Y3 0.1613852 0.1565957 0.1544661
#
# $VE
# 1 2 3
# Y1 0.6839938 0.1866749 0.3413782
# Y2 0.1866749 0.1892320 0.2200922
# Y3 0.3413782 0.2200922 0.3500427
#
# $VE.se
# 1 2 3
# Y1 0.1932288 0.1218645 0.1436571
# Y2 0.1218645 0.1379212 0.1334911
# Y3 0.1436571 0.1334911 0.1544661
#
# $VP
# 1 2 3
# Y1 1.0000000 0.4666511 0.6801728
# Y2 0.4666511 1.0000000 0.7689803
# Y3 0.6801728 0.7689803 1.0000000
#
# $RG
# 1 2 3
# Y1 1.0000000 0.5531272 0.7475604
# Y2 0.5531272 1.0000000 0.7561242
# Y3 0.7475604 0.7561242 1.0000000
#
# $RG.se
# 1 2 3
# Y1 0.0000000 0.23210789 1.574432e-01
# Y2 0.2321079 0.00000000 8.945803e-02
# Y3 0.1574432 0.08945803 1.633888e-17
#
# $RE
# 1 2 3
# Y1 1.0000000 0.5188747 0.6976686
# Y2 0.5188747 1.0000000 0.8551589
# Y3 0.6976686 0.8551589 1.0000000
#
# $RE.se
# 1 2 3
# Y1 0.0000000 2.130183e-01 1.172975e-01
# Y2 0.2130183 6.290479e-17 8.880531e-02
# Y3 0.1172975 8.880531e-02 5.827606e-17
## ----eval=FALSE---------------------------------------------------------------
# > fit.DS <- grmsem.fit(ph.small,
# G.small,
# model="DS",
# LogL=TRUE,
# estSE=TRUE)
# [1] 300
# [1] 100 100 100
# [1] 1 101 201
# [1] 100 200 300
# [1] "Parameters for a direct symmetric model (variance components):"
# [1] "Vector of genetic variances:"
# [1] "A11" "A21" "A31" "A22" "A32" "A33"
# [1] "Vector of residual factor loadings:"
# [1] "E11" "E21" "E31" "E22" "E32" "E33"
#
# > fit.DS
# $model.in
# part label value freepar
# 1 a A11 1.0000000 1
# 2 a A21 0.4602739 1
# 3 a A31 0.6738565 1
# 4 a A22 1.0000000 1
# 5 a A32 0.7527887 1
# 6 a A33 1.0000000 1
# 7 e E11 1.0000000 1
# 8 e E21 0.4602739 1
# 9 e E31 0.6738565 1
# 10 e E22 1.0000000 1
# 11 e E32 0.7527887 1
# 12 e E33 1.0000000 1
#
# $model.fit$LL
# [1] -54.39145
#
# $model.fit$calls
# function gradient
# 122 26
#
# $model.fit$convergence
# [1] 0
#
# $model.fit$message
# NULL
#
#
# $model.out
# label estimates gradient se Z p
# 1 A11 0.3161061 0.0009860302 0.2106008 1.500972 0.1333627052
# 2 A21 0.2912542 0.0039066047 0.1834675 1.587497 0.1124001236
# 3 A31 0.3389444 -0.0036802953 0.1852176 1.829979 0.0672530362
# 4 A22 0.8754990 0.0025629175 0.2584953 3.386905 0.0007068577
# 5 A32 0.5705239 -0.0051792201 0.2222947 2.566520 0.0102724710
# 6 A33 0.6497819 0.0048300247 0.2229978 2.913849 0.0035700307
# 7 E11 0.6828573 0.0009099445 0.1921358 3.554035 0.0003793684
# 8 E21 0.1932362 0.0081933028 0.1180255 1.637241 0.1015800561
# 9 E31 0.3403784 -0.0046169601 0.1387849 2.452561 0.0141843201
# 10 E22 0.2034536 0.0015010666 0.1268743 1.603584 0.1088057273
# 11 E32 0.2277079 0.0414495710 0.1222574 1.862529 0.0625285683
# 12 E33 0.3487737 0.0049090361 0.1373909 2.538551 0.0111312504
#
# $k
# [1] 3
#
# $n
# [1] 300
#
# $n.obs
# [1] 100 100 100
#
# $n.ind
# [1] 100
#
# $model
# [1] "DS"
## ----eval=FALSE---------------------------------------------------------------
# > var.fit.ds <- grmsem.var(fit.DS)
# > var.fit.ds
# $VA
# 1 2 3
# Y1 0.3161061 0.2912542 0.3389444
# Y2 0.2912542 0.8754990 0.5705239
# Y3 0.3389444 0.5705239 0.6497819
#
# $VA.se
# 1 2 3
# Y1 0.2106008 0.1834675 0.1852176
# Y2 0.1834675 0.2584953 0.2222947
# Y3 0.1852176 0.2222947 0.2229978
#
# $VE
# 1 2 3
# Y1 0.6828573 0.1932362 0.3403784
# Y2 0.1932362 0.2034536 0.2277079
# Y3 0.3403784 0.2277079 0.3487737
#
# $VE.se
# 1 2 3
# Y1 0.1921358 0.1180255 0.1387849
# Y2 0.1180255 0.1268743 0.1222574
# Y3 0.1387849 0.1222574 0.1373909
#
# $VP
# 1 2 3
# Y1 0.9989633 0.4844903 0.6793228
# Y2 0.4844903 1.0789526 0.7982318
# Y3 0.6793228 0.7982318 0.9985556
#
# $RG
# 1 2 3
# Y1 1.0000000 0.5536406 0.7478738
# Y2 0.5536406 1.0000000 0.7564186
# Y3 0.7478738 0.7564186 1.0000000
#
# $RG.se
# 1 2 3
# Y1 0.0000000 0.22323218 0.15738360
# Y2 0.2232322 0.00000000 0.08684325
# Y3 0.1573836 0.08684325 0.00000000
#
# $RE
# 1 2 3
# Y1 1.0000000 0.5184307 0.6974693
# Y2 0.5184307 1.0000000 0.8548176
# Y3 0.6974693 0.8548176 1.0000000
#
# $RE.se
# 1 2 3
# Y1 0.0000000 0.20548054 0.11753303
# Y2 0.2054805 0.00000000 0.08761989
# Y3 0.1175330 0.08761989 0.00000000
## ----eval = FALSE-------------------------------------------------------------
# #Do not run!
# #Please downloaded externally
# #load("G.large.RData")
# #load("ph.large.RData")
# #fit <- grmsem.fit(ph.large,
# # G.large,model="Cholesky",
# # LogL=TRUE,
# # estSE=TRUE)
#
# > load("fit.large.RData")
# > fit.large
# $model.in
# part label value freepar
# 1 a a11 0.01686961 1
# 2 a a21 0.75022252 1
# 3 a a31 -0.10352849 1
# 4 a a41 -0.81445295 1
# 5 a a22 -0.34964589 1
# 6 a a32 0.70013464 1
# 7 a a42 0.30439698 1
# 8 a a33 0.87503509 1
# 9 a a43 -0.37434617 1
# 10 a a44 0.86376688 1
# 11 e e11 -0.28008894 1
# 12 e e21 -0.92449828 1
# 13 e e31 -0.72312602 1
# 14 e e41 -0.70924504 1
# 15 e e22 0.51178943 1
# 16 e e32 0.66207502 1
# 17 e e42 -0.64331456 1
# 18 e e33 0.16870994 1
# 19 e e43 0.30222504 1
# 20 e e44 0.21833386 1
#
# $model.fit$LL
# [1] -4556.647
#
# $model.fit$calls
# function gradient
# 488 118
#
# $model.fit$convergence
# [1] 0
#
# $model.fit$message
# NULL
#
# $model.out
# label estimates gradient se Z p
# 1 a11 -0.590283919 7.196924e-03 0.018936941 -31.1710286 2.632163e-213
# 2 a21 -0.406228554 4.760323e-04 0.022331488 -18.1908415 6.099640e-74
# 3 a31 -0.412282704 -1.306064e-03 0.021811928 -18.9017085 1.104185e-79
# 4 a41 -0.251181508 4.515862e-03 0.025355870 -9.9062469 3.910610e-23
# 5 a22 -0.611871764 7.519419e-03 0.013081543 -46.7736703 0.000000e+00
# 6 a32 -0.412494347 5.020765e-05 0.015153798 -27.2205248 3.712770e-163
# 7 a42 -0.406496327 2.858898e-03 0.020248170 -20.0757069 1.203583e-89
# 8 a33 0.443344860 -2.806537e-04 0.013549030 32.7215213 7.719589e-235
# 9 a43 -0.015668913 1.437494e-03 0.024698450 -0.6344088 5.258141e-01
# 10 a44 0.655953940 -5.747793e-03 0.014789225 44.3535029 0.000000e+00
# 11 e11 -0.794244501 2.325624e-02 0.012513836 -63.4693048 0.000000e+00
# 12 e21 -0.290249172 -8.449568e-04 0.013042715 -22.2537381 1.037729e-109
# 13 e31 -0.445243952 -5.908758e-03 0.012154714 -36.6313812 9.056960e-294
# 14 e41 0.006719942 7.527450e-03 0.011697368 0.5744833 5.656408e-01
# 15 e22 0.521077401 -1.251593e-02 0.008841668 58.9342896 0.000000e+00
# 16 e32 -0.054613913 1.163049e-03 0.010020443 -5.4502495 5.029922e-08
# 17 e42 -0.018229268 -8.354153e-03 0.011967918 -1.5231779 1.277142e-01
# 18 e33 0.405590017 -2.845990e-02 0.007471397 54.2856995 0.000000e+00
# 19 e43 -0.022769655 7.889148e-04 0.012392180 -1.8374213 6.614772e-02
# 20 e44 -0.467128074 1.039490e-02 0.009007167 -51.8618176 0.000000e+00
#
# $k
# [1] 4
#
# $n
# [1] 20000
#
# $n.obs
# [1] 5000 5000 5000 5000
#
# $n.ind
# [1] 5000
#
# $model
# [1] "Cholesky"
## ----eval = FALSE-------------------------------------------------------------
# > var.fit <- grmsem.var(fit.large)
# > print(var.fit)
# $VA
# 1 2 3 4
# Y1 0.3484351 0.2397902 0.2433639 0.1482684
# Y2 0.2397902 0.5394087 0.4198747 0.3507607
# Y3 0.2433639 0.4198747 0.5366833 0.2642885
# Y4 0.1482684 0.3507607 0.2642885 0.6588525
#
# $VA.se
# 1 2 3 4
# Y1 0.02235634 0.01660915 0.01767332 0.01514868
# Y2 0.01660915 0.02017283 0.01672469 0.01562045
# Y3 0.01767332 0.01672469 0.02042099 0.01515375
# Y4 0.01514868 0.01562045 0.01515375 0.02078249
#
# $VE
# 1 2 3 4
# Y1 0.630824327 0.23052881 0.35363256 -0.005337277
# Y2 0.230528809 0.35576624 0.10077361 -0.011449317
# Y3 0.353632560 0.10077361 0.36572812 -0.011231587
# Y4 -0.005337277 -0.01144932 -0.01123159 0.219104559
#
# $VE.se
# 1 2 3 4
# Y1 0.019878092 0.012210347 0.013679555 0.009287867
# Y2 0.012210347 0.012054860 0.009085111 0.007093293
# Y3 0.013679555 0.009085111 0.012413525 0.007220553
# Y4 0.009287867 0.007093293 0.007220553 0.008358356
#
# $VP
# 1 2 3 4
# Y1 0.9792594 0.4703190 0.5969964 0.1429311
# Y2 0.4703190 0.8951749 0.5206483 0.3393114
# Y3 0.5969964 0.5206483 0.9024114 0.2530569
# Y4 0.1429311 0.3393114 0.2530569 0.8779571
#
# $RG
# 1 2 3 4
# Y1 1.0000000 0.553110 0.5627767 0.3094522
# Y2 0.5531100 1.000000 0.7803720 0.5883800
# Y3 0.5627767 0.780372 1.0000000 0.4444523
# Y4 0.3094522 0.588380 0.4444523 1.0000000
#
# $RG.se
# 1 2 3 4
# Y1 0.00000000 0.02488356 2.233207e-02 3.056075e-02
# Y2 0.02488356 0.00000000 1.557784e-02 1.976462e-02
# Y3 0.02233207 0.01557784 3.008489e-18 2.206251e-02
# Y4 0.03056075 0.01976462 2.206251e-02 4.557191e-18
#
# $RE
# 1 2 3 4
# Y1 1.00000000 0.48661851 0.73623906 -0.01435621
# Y2 0.48661851 1.00000000 0.27937355 -0.04100828
# Y3 0.73623906 0.27937355 1.00000000 -0.03967676
# Y4 -0.01435621 -0.04100828 -0.03967676 1.00000000
#
# $RE.se
# 1 2 3 4
# Y1 1.111452e-17 0.01770789 1.099568e-02 2.499190e-02
# Y2 1.770789e-02 0.00000000 2.164591e-02 2.549970e-02
# Y3 1.099568e-02 0.02164591 5.397777e-18 2.555311e-02
# Y4 2.499190e-02 0.02549970 2.555311e-02 7.971451e-18
#
## ----eval = FALSE-------------------------------------------------------------
# > st.fit <- grmsem.stpar(fit.large)
# > print(st.fit)
#
# $stand.model.out
# label stand.estimates stand.se Z p
# 1 a11 -0.596502229 0.016233733 -36.7446127 1.417362e-295
# 2 a21 -0.429354966 0.020456160 -20.9890304 8.261383e-98
# 3 a31 -0.434003098 0.019542251 -22.2084494 2.845907e-109
# 4 a41 -0.268071735 0.024887619 -10.7712891 4.703688e-27
# 5 a22 -0.646705353 0.011739313 -55.0888607 0.000000e+00
# 6 a32 -0.434225891 0.014761620 -29.4158701 3.441759e-190
# 7 a42 -0.433830407 0.018983259 -22.8533157 1.354554e-115
# 8 a33 0.466701710 0.013216780 35.3113009 3.939089e-273
# 9 a43 -0.016722540 0.024694677 -0.6771718 4.982969e-01
# 10 a44 0.700062328 0.012710221 55.0786896 0.000000e+00
# 11 e11 -0.802611420 0.012064939 -66.5242831 0.000000e+00
# 12 e21 -0.306772929 0.012940717 -23.7060220 3.125266e-124
# 13 e31 -0.468700852 0.012081674 -38.7943628 0.000000e+00
# 14 e41 0.007171812 0.011697423 0.6131104 5.398033e-01
# 15 e22 0.550742107 0.009617573 57.2641473 0.000000e+00
# 16 e32 -0.057491151 0.010016353 -5.7397289 9.482821e-09
# 17 e42 -0.019455060 0.011964987 -1.6259993 1.039498e-01
# 18 e33 0.426957819 0.008330528 51.2521942 0.000000e+00
# 19 e43 -0.024300758 0.012391052 -1.9611537 4.986110e-02
# 20 e44 -0.498539222 0.010093805 -49.3906111 0.000000e+00
#
# $k
# [1] 4
#
# $n
# [1] 20000
#
# $model
# [1] "Cholesky"
## ----eval=FALSE---------------------------------------------------------------
# > st.var.fit <- grmsem.var(st.fit)
# > print(st.var.fit)
# $VA
# 1 2 3 4
# Y1 0.3558149 0.2561112 0.2588838 0.1599054
# Y2 0.2561112 0.6025735 0.4671576 0.3956584
# Y3 0.2588838 0.4671576 0.5947213 0.2969199
# Y4 0.1599054 0.3956584 0.2969199 0.7504382
#
# $VA.se
# 1 2 3 4
# Y1 0.01936692 0.01501072 0.01589623 0.01481065
# Y2 0.01501072 0.01344962 0.01252746 0.01320059
# Y3 0.01589623 0.01252746 0.01377554 0.01417036
# Y4 0.01481065 0.01320059 0.01417036 0.01000981
#
# $VE
# 1 2 3 4
# Y1 0.644185091 0.24621946 0.37618466 -0.005756178
# Y2 0.246219456 0.39742650 0.11212194 -0.012914839
# Y3 0.376184656 0.11212194 0.40527870 -0.012618339
# Y4 -0.005756178 -0.01291484 -0.01261834 0.249561817
#
# $VE.se
# 1 2 3 4
# Y1 0.019366916 0.012115870 0.013465601 0.009385665
# Y2 0.012115870 0.013449618 0.009581551 0.007494475
# Y3 0.013465601 0.009581551 0.013775535 0.007600587
# Y4 0.009385665 0.007494475 0.007600587 0.010009812
#
# $VP
# 1 2 3 4
# Y1 1.0000000 0.5023307 0.6350685 0.1541492
# Y2 0.5023307 1.0000000 0.5792795 0.3827435
# Y3 0.6350685 0.5792795 1.0000000 0.2843016
# Y4 0.1541492 0.3827435 0.2843016 1.0000000
#
# $RG
# 1 2 3 4
# Y1 1.0000000 0.553110 0.5627767 0.3094522
# Y2 0.5531100 1.000000 0.7803720 0.5883800
# Y3 0.5627767 0.780372 1.0000000 0.4444523
# Y4 0.3094522 0.588380 0.4444523 1.0000000
#
# $RG.se
# 1 2 3 4
# Y1 7.209226e-18 2.354325e-02 2.121442e-02 2.863522e-02
# Y2 2.354325e-02 4.542180e-18 1.479096e-02 1.851946e-02
# Y3 2.121442e-02 1.479096e-02 5.454486e-18 2.069598e-02
# Y4 2.863522e-02 1.851946e-02 2.069598e-02 4.092585e-18
#
# $RE
# 1 2 3 4
# Y1 1.00000000 0.48661851 0.73623906 -0.01435621
# Y2 0.48661851 1.00000000 0.27937355 -0.04100828
# Y3 0.73623906 0.27937355 1.00000000 -0.03967676
# Y4 -0.01435621 -0.04100828 -0.03967676 1.00000000
#
# $RE.se
# 1 2 3 4
# Y1 0.00000000 1.675409e-02 1.044538e-02 2.341725e-02
# Y2 0.01675409 4.271060e-18 2.052641e-02 2.389309e-02
# Y3 0.01044538 2.052641e-02 6.732149e-18 2.394314e-02
# Y4 0.02341725 2.389309e-02 2.394314e-02 4.512840e-18
#
## ----eval=FALSE---------------------------------------------------------------
# > grmsem.fcoher(fit.large)
# > grmsem.fcoher(fit.large)$fcoher.model.out[,c(1,7:10)]
#
# label Vi Vi.se FCOHER FCOHER.se
# 1 a11 3.484351e-01 0.0223563431 1.0000000000 0.000000000
# 2 a21 1.650216e-01 0.0181433759 0.3059306238 0.027526691
# 3 a31 1.699770e-01 0.0179853613 0.3167175772 0.025135935
# 4 a41 6.309215e-02 0.0127378513 0.0957606594 0.018914180
# 5 a22 3.743871e-01 0.0160084532 0.6940693762 0.027526691
# 6 a32 1.701516e-01 0.0125017124 0.3170428312 0.024459178
# 7 a42 1.652393e-01 0.0164616136 0.2507985687 0.023183443
# 8 a33 1.965547e-01 0.0120137852 0.3662395917 0.021786386
# 9 a43 2.455148e-04 0.0007739957 0.0003726401 0.001174465
# 10 a44 4.302756e-01 0.0194021010 0.6530681318 0.023536859
# 11 e11 6.308243e-01 0.0198780916 NA NA
# 12 e21 8.424458e-02 0.0075712747 NA NA
# 13 e31 1.982422e-01 0.0108236256 NA NA
# 14 e41 4.515762e-05 0.0001572113 NA NA
# 15 e22 2.715217e-01 0.0092143864 NA NA
# 16 e32 2.982679e-03 0.0010945112 NA NA
# 17 e42 3.323062e-04 0.0004363328 NA NA
# 18 e33 1.645033e-01 0.0060606482 NA NA
# 19 e43 5.184572e-04 0.0005643313 NA NA
# 20 e44 2.182086e-01 0.0084150015 NA NA
#
## ----eval=FALSE---------------------------------------------------------------
# > grmsem.fcoher(fit.large)
# > grmsem.fcoher(fit.large)$fcoher.model.out[,c(1,7:8,13:14)]
#
# label Vi Vi.se FCOENV FCOENV.se
# 1 a11 3.484351e-01 0.0223563431 NA NA
# 2 a21 1.650216e-01 0.0181433759 NA NA
# 3 a31 1.699770e-01 0.0179853613 NA NA
# 4 a41 6.309215e-02 0.0127378513 NA NA
# 5 a22 3.743871e-01 0.0160084532 NA NA
# 6 a32 1.701516e-01 0.0125017124 NA NA
# 7 a42 1.652393e-01 0.0164616136 NA NA
# 8 a33 1.965547e-01 0.0120137852 NA NA
# 9 a43 2.455148e-04 0.0007739957 NA NA
# 10 a44 4.302756e-01 0.0194021010 NA NA
# 11 e11 6.308243e-01 0.0198780916 1.0000000000 0.0000000000
# 12 e21 8.424458e-02 0.0075712747 0.2367975722 0.0172339750
# 13 e31 1.982422e-01 0.0108236256 0.5420479497 0.0161908960
# 14 e41 4.515762e-05 0.0001572113 0.0002061008 0.0007175779
# 15 e22 2.715217e-01 0.0092143864 0.7632024278 0.0172339750
# 16 e32 2.982679e-03 0.0010945112 0.0081554556 0.0029995986
# 17 e42 3.323062e-04 0.0004363328 0.0015166558 0.0019944184
# 18 e33 1.645033e-01 0.0060606482 0.4497965946 0.0162294819
# 19 e43 5.184572e-04 0.0005643313 0.0023662548 0.0025782116
# 20 e44 2.182086e-01 0.0084150015 0.9959109886 0.0034506100
## ----eval = FALSE-------------------------------------------------------------
# > fit.var <- grmsem.var(fit)
# > load("ph.large.RData")
# > grmsem.biher(ph.large, fit.var)
# $VPO
# 1 2 3 4
# Y1 1.0000000 0.5061277 0.6278417 0.1792348
# Y2 0.5061277 1.0000000 0.6167281 0.4352479
# Y3 0.6278417 0.6167281 1.0000000 0.3331679
# Y4 0.1792348 0.4352479 0.3331679 1.0000000
#
# $VA
# 1 2 3 4
# Y1 0.3484351 0.2397902 0.2433639 0.1482684
# Y2 0.2397902 0.5394087 0.4198747 0.3507607
# Y3 0.2433639 0.4198747 0.5366833 0.2642885
# Y4 0.1482684 0.3507607 0.2642885 0.6588525
#
# $BIHER
# 1 2 3 4
# Y1 0.3484351 0.4737741 0.3876197 0.8272299
# Y2 0.4737741 0.5394087 0.6808100 0.8058872
# Y3 0.3876197 0.6808100 0.5366833 0.7932591
# Y4 0.8272299 0.8058872 0.7932591 0.6588525
#
# $BIHER.se
# 1 2 3 4
# Y1 0.02235634 0.03281612 0.02814933 0.08451862
# Y2 0.03281612 0.02017283 0.02711842 0.03588862
# Y3 0.02814933 0.02711842 0.02042099 0.04548383
# Y4 0.08451862 0.03588862 0.04548383 0.02078249
#
# $BIHER.Z
# 1 2 3 4
# Y1 15.585514 14.43723 13.77012 9.787546
# Y2 14.437234 26.73937 25.10508 22.455231
# Y3 13.770125 25.10508 26.28097 17.440466
# Y4 9.787546 22.45523 17.44047 31.702288
#
# $BIHER.p
# 1 2 3 4
# Y1 9.132538e-55 3.017077e-47 3.855470e-43 1.273488e-22
# Y2 3.017077e-47 1.641347e-157 4.377379e-139 1.137653e-111
# Y3 3.855470e-43 4.377379e-139 3.165351e-152 4.067453e-68
# Y4 1.273488e-22 1.137653e-111 4.067453e-68 1.444882e-220
#
## ----eval = FALSE-------------------------------------------------------------
# #Do not run!
# #Please downloaded externally
# #load("G.large.RData")
# #load("ph.large.RData")
# #> out<-grmsem.fit(ph.large,
# # G.large,
# # model="DS",
# # LogL=TRUE,
# # estSE=TRUE)
#
# > out
# $model.in
# part label value freepar
# 1 a A11 1.0000000 1
# 2 a A21 0.5061277 1
# 3 a A31 0.6278417 1
# 4 a A41 0.1792348 1
# 5 a A22 1.0000000 1
# 6 a A32 0.6167281 1
# 7 a A42 0.4352479 1
# 8 a A33 1.0000000 1
# 9 a A43 0.3331679 1
# 10 a A44 1.0000000 1
# 11 e E11 1.0000000 1
# 12 e E21 0.5061277 1
# 13 e E31 0.6278417 1
# 14 e E41 0.1792348 1
# 15 e E22 1.0000000 1
# 16 e E32 0.6167281 1
# 17 e E42 0.4352479 1
# 18 e E33 1.0000000 1
# 19 e E43 0.3331679 1
# 20 e E44 1.0000000 1
#
#
# $model.fit$LL
# [1] -4556.647
#
# $model.fit$calls
# function gradient
# 173 38
#
# $model.fit$convergence
# [1] 0
#
# $model.fit$message
# NULL
#
#
# $model.out
# label estimates gradient se Z p
# 1 A11 0.348458001 -0.052226894 0.022356069 15.5867299 8.960447e-55
# 2 A21 0.239787813 0.073256044 0.016608680 14.4374998 3.005481e-47
# 3 A31 0.243373776 0.046398011 0.017672850 13.7710544 3.806174e-43
# 4 A41 0.148252661 0.045240115 0.015148427 9.7866704 1.284564e-22
# 5 A22 0.539403719 -0.038929159 0.020172359 26.7397447 1.625068e-157
# 6 A32 0.419867512 0.053261468 0.016724145 25.1054698 4.334317e-139
# 7 A42 0.350754104 -0.026684946 0.015620099 22.4553060 1.135724e-111
# 8 A33 0.536686987 -0.032016186 0.020420456 26.2818311 3.094103e-152
# 9 A43 0.264276207 0.014920668 0.015153366 17.4400997 4.093589e-68
# 10 A44 0.658845787 0.011119193 0.020782130 31.7025151 1.434490e-220
# 11 E11 0.630774354 -0.070278866 0.019875702 31.7359534 4.961474e-221
# 12 E21 0.230506722 0.166057579 0.012209173 18.8797985 1.672229e-79
# 13 E31 0.353594799 0.308389144 0.013677755 25.8518155 2.321384e-147
# 14 E41 -0.005335735 -0.015130111 0.009287206 -0.5745253 5.656124e-01
# 15 E22 0.355757271 -0.040844105 0.012054375 29.5127094 1.977804e-191
# 16 E32 0.100760965 -0.097429983 0.009084347 11.0917126 1.376273e-28
# 17 E42 -0.011450535 0.016775106 0.007093098 -1.6143207 1.064579e-01
# 18 E33 0.365701992 -0.075051659 0.012412202 29.4631041 8.554205e-191
# 19 E43 -0.011230439 0.004879057 0.007220099 -1.5554412 1.198410e-01
# 20 E44 0.219102869 0.026457431 0.008358236 26.2140086 1.839860e-151
#
# $k
# [1] 4
#
# $n
# [1] 20000
#
# $n.obs
# [1] 5000 5000 5000 5000
#
# $n.ind
# [1] 5000
#
# $model
# [1] "DS"
## ----eval = FALSE-------------------------------------------------------------
# > var.out<-grmsem.var(out)
# > print(var.out)
# $VA
# 1 2 3 4
# Y1 0.3484580 0.2397878 0.2433738 0.4198675
# Y2 0.2397878 0.5394037 0.1482527 0.3507541
# Y3 0.2433738 0.4198675 0.5366870 0.2642762
# Y4 0.1482527 0.3507541 0.2642762 0.6588458
#
# $VA.se
# 1 2 3 4
# Y1 0.02235607 0.01660868 0.01767285 0.01672414
# Y2 0.01660868 0.02017236 0.01514843 0.01562010
# Y3 0.01767285 0.01672414 0.02042046 0.01515337
# Y4 0.01514843 0.01562010 0.01515337 0.02078213
#
# $VE
# 1 2 3 4
# Y1 0.630774354 0.23050672 0.353594799 0.10076097
# Y2 0.230506722 0.35575727 -0.005335735 -0.01145053
# Y3 0.353594799 0.10076097 0.365701992 -0.01123044
# Y4 -0.005335735 -0.01145053 -0.011230439 0.21910287
#
# $VE.se
# 1 2 3 4
# Y1 0.019875702 0.012209173 0.013677755 0.009084347
# Y2 0.012209173 0.012054375 0.009287206 0.007093098
# Y3 0.013677755 0.009084347 0.012412202 0.007220099
# Y4 0.009287206 0.007093098 0.007220099 0.008358236
#
# $VP
# 1 2 3 4
# Y1 0.9792324 0.4702945 0.5969686 0.5206285
# Y2 0.4702945 0.8951610 0.1429169 0.3393036
# Y3 0.5969686 0.5206285 0.9023890 0.2530458
# Y4 0.1429169 0.3393036 0.2530458 0.8779487
#
# $RG
# 1 2 3 4
# Y1 1.0000000 0.5530889 0.5627792 0.8762846
# Y2 0.5530889 1.0000000 0.2755402 0.5883746
# Y3 0.5627792 0.7803596 1.0000000 0.4444323
# Y4 0.3094107 0.5883746 0.4444323 1.0000000
#
# $RG.se
# 1 2 3 4
# Y1 0.00000000 0.02488339 0.02233109 0.03667066
# Y2 0.02488339 0.00000000 0.02690300 0.01976462
# Y3 0.02233109 0.01557817 0.00000000 0.02206248
# Y4 0.03055988 0.01976462 0.02206248 0.00000000
#
# $RE
# 1 2 3 4
# Y1 1.00000000 0.48659729 0.73621590 0.27103868
# Y2 0.48659729 1.00000000 -0.01479291 -0.04101331
# Y3 0.73621590 0.27935199 1.00000000 -0.03967428
# Y4 -0.01435269 -0.04101331 -0.03967428 1.00000000
#
# $RE.se
# 1 2 3 4
# Y1 4.413292e-18 0.01770803 1.099634e-02 0.02340894
# Y2 1.770803e-02 0.00000000 2.576743e-02 0.02549942
# Y3 1.099634e-02 0.02164575 5.512125e-18 0.02555251
# Y4 2.499120e-02 0.02549942 2.555251e-02 0.00000000
## ----eval = FALSE-------------------------------------------------------------
# gcta64 --reml-bivar 1 2 --grm-gz large.gcta --pheno large.gcta.phe
# --reml-bivar-lrt-rg 0 --thread-num 4 --out Y1_Y2 --reml-maxit 1000 >
# large.gcta.Y1Y2.log`
## ----eval = FALSE-------------------------------------------------------------
# Summary result of REML analysis (shortened output):
# Source Variance SE
# V(G)/Vp_tr1 0.361387 0.020470
# V(G)/Vp_tr2 0.635426 0.015053
# rG 0.557367 0.026039
## ----eval = FALSE-------------------------------------------------------------
# load("ph.large.RData")
# load("G.large.RData")
# ph.biv<-ph.large[,c(1,2)]
# fit <- grmsem.fit(ph.biv, G.large, LogL = TRUE, estSE = TRUE)
# print(fit)
# st.fit <- grmsem.stpar(fit)
# print(st.fit)
# st.var.fit <- grmsem.var(st.fit)
# print(st.var.fit)
## ----eval = FALSE-------------------------------------------------------------
# $VA
# 1 2
# Y1 0.3613708 0.2669526
# Y2 0.2669526 0.6349136
#
# $VA.se
# 1 2
# Y1 0.01943807 0.01555813
# Y2 0.01555813 0.01355490
#
# $RG
# 1 2
# Y1 1.0000000 0.5573145
# Y2 0.5573145 1.0000000
#
# $RG.se
# 1 2
# Y1 7.179876e-18 0.0239077
# Y2 2.390770e-02 0.0000000
#
#
Try the grmsem package in your browser
Any scripts or data that you put into this service are public.
grmsem documentation built on Jan. 29, 2021, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----eval = FALSE-------------------------------------------------------------
# > G <- grm.input("large.gcta.grm.gz")
# > G[1:3,1:3] #Relationships among the first three individuals
# [,1] [,2] [,3]
# [1,] 0.99354762 0.02328514 0.01644197
# [2,] 0.02328514 0.99406837 0.01021175
# [3,] 0.01644197 0.01021175 1.02751472
## ----eval = FALSE-------------------------------------------------------------
# > load("ph.large.RData")
# > ph.large[1:2,]
# Y1 Y2 Y3 Y4
# [1,] -0.7640819 -0.6016908 -0.3981901 -0.3169821
# [2,] -0.5099606 0.6671311 -1.3119328 -0.5601261
## ----path-cholesky, echo=FALSE, out.width="65%", fig.align='center', fig.cap="Quad-variate Cholesky decomposition model. Observed phenotypes (P) are represented by squares and factors by circles. Single-headed arrows (paths) define relationships between variables. The variance of latent variables is constrained to unit variance; this is omitted from the path diagram to improve clarity."----
knitr::include_graphics("model_quadvariate_path_diagram_Cholesky.png")
## ----path-indep, echo = FALSE, out.width="65%", fig.align='center', fig.cap="Quad-variate Independent Pathway model. Observed phenotypes (P) are represented by squares and factors by circles. Single-headed arrows (paths) define relationships between variables. The variance of latent variables is constrained to unit variance; this is omitted from the path diagram to improve clarity."----
knitr::include_graphics("model_quadvariate_path_diagram_IP.png")
## ----path-ipc, echo = FALSE, out.width="65%", fig.align='center', fig.cap="Quad-variate IPC model. Observed phenotypes (P) are represented by squares and factors by circles. Single-headed arrows (paths) define relationships between variables. The variance of latent variables is constrained to unit variance; this is omitted from the path diagram to improve clarity."----
knitr::include_graphics("model_quadvariate_path_diagram_IPC.png")
## ----eval = FALSE-------------------------------------------------------------
# #If files were downloaded externally, please load here:
# > load("G.small.RData")
# > load("ph.small.RData")
# > fit <- grmsem.fit(ph.small,
# G.small,
# model="Cholesky",
# LogL=TRUE,
# estSE=TRUE)
# > fit
# $model.in
# part label value freepar
# 1 a a11 0.6971323 1
# 2 a a21 -0.2357364 1
# 3 a a31 0.3048627 1
# 4 a a22 0.2067668 1
# 5 a a32 -0.3823634 1
# 6 a a33 -0.1557505 1
# 7 e e11 -0.4958634 1
# 8 e e21 -0.5924593 1
# 9 e e31 0.1700410 1
# 10 e e22 0.3570643 1
# 11 e e32 -0.3220280 1
# 12 e e33 0.5024512 1
#
# $model.fit$LL
# [1] -54.39143
#
# $model.out
# label estimates gradient se Z p
# 1 a11 -0.5618594 -4.490283e-05 0.18730897 -2.999640 2.702993e-03
# 2 a21 -0.5172780 6.160157e-05 0.25150203 -2.056755 3.970982e-02
# 3 a31 -0.6022167 -4.596719e-05 0.20040106 -3.005057 2.655308e-03
# 4 a22 -0.7791023 -1.154262e-05 0.13696258 -5.688432 1.282113e-08
# 5 a32 -0.3313097 6.118543e-05 0.14755775 -2.245289 2.474961e-02
# 6 a33 -0.4201447 7.389654e-05 0.07000431 -6.001697 1.952654e-09
# 7 e11 -0.8266200 2.460492e-04 0.11621420 -7.112900 1.136299e-12
# 8 e21 -0.2344285 -1.605268e-04 0.12464049 -1.880838 5.999398e-02
# 9 e31 -0.4124524 6.209974e-05 0.12534761 -3.290469 1.000207e-03
# 10 e22 0.3862229 -2.636815e-05 0.12634091 3.056990 2.235717e-03
# 11 e32 0.3410509 8.582604e-05 0.10625220 3.209824 1.328162e-03
# 12 e33 -0.2511355 1.911040e-04 0.05125420 -4.899803 9.593276e-07
#
# $k
# [1] 3
#
# $n
# [1] 300
#
# $n.obs
# [1] 100 100 100
#
# $n.ind
# [1] 100
#
# $model
# [1] "Cholesky"
## ----eval = FALSE-------------------------------------------------------------
# > var.fit <- grmsem.var(fit)
# > print(var.fit)
# $VA
# 1 2 3
# Y1 0.3156860 0.2906375 0.3383611
# Y2 0.2906375 0.8745770 0.5696376
# Y3 0.3383611 0.5696376 0.6489526
#
# $VA.se
# 1 2 3
# Y1 0.2104826 0.1832943 0.1850112
# Y2 0.1832943 0.2583423 0.2220882
# Y3 0.1850112 0.2220882 0.2227517
#
# $VE
# 1 2 3
# Y1 0.6833006 0.1937833 0.3409414
# Y2 0.1937833 0.2041249 0.2284123
# Y3 0.3409414 0.2284123 0.3495017
#
# $VE.se
# 1 2 3
# Y1 0.1921300 0.1180125 0.1387263
# Y2 0.1180125 0.1269780 0.1222845
# Y3 0.1387263 0.1222845 0.1373590
#
# $VP
# 1 2 3
# Y1 0.9989866 0.4844208 0.6793025
# Y2 0.4844208 1.0787019 0.7980499
# Y3 0.6793025 0.7980499 0.9984543
#
# $RG
# 1 2 3
# Y1 1.0000000 0.5531272 0.7475604
# Y2 0.5531272 1.0000000 0.7561242
# Y3 0.7475604 0.7561242 1.0000000
#
# $RG.se
# 1 2 3
# Y1 0.0000000 0.22348027 1.575650e-01
# Y2 0.2234803 0.00000000 8.697161e-02
# Y3 0.1575650 0.08697161 4.449797e-17
#
# $RE
# 1 2 3
# Y1 1.0000000 0.5188747 0.6976686
# Y2 0.5188747 1.0000000 0.8551589
# Y3 0.6976686 0.8551589 1.0000000
#
# $RE.se
# 1 2 3
# Y1 5.160947e-17 2.051003e-01 1.173883e-01
# Y2 2.051003e-01 5.610664e-17 8.730948e-02
# Y3 1.173883e-01 8.730948e-02 2.359273e-17
#
## ----eval = FALSE-------------------------------------------------------------
# > print(st.fit)
# $stand.model.out
# label stand.estimates stand.se Z p
# 1 a11 -0.5621443 0.17186763 -3.270798 1.072444e-03
# 2 a21 -0.4980504 0.22959865 -2.169222 3.006584e-02
# 3 a31 -0.6026826 0.17316262 -3.480443 5.005852e-04
# 4 a22 -0.7501425 0.12682263 -5.914895 3.320880e-09
# 5 a32 -0.3315661 0.14661954 -2.261404 2.373422e-02
# 6 a33 -0.4204698 0.07358378 -5.714164 1.102447e-08
# 7 e11 -0.8270392 0.11681963 -7.079625 1.445454e-12
# 8 e21 -0.2257147 0.12967922 -1.740562 8.176046e-02
# 9 e31 -0.4127715 0.13122602 -3.145500 1.658030e-03
# 10 e22 0.3718667 0.14164900 2.625269 8.658054e-03
# 11 e32 0.3413148 0.11390419 2.996507 2.730917e-03
# 12 e33 -0.2513298 0.05457774 -4.604987 4.124919e-06
#
## ----eval=FALSE---------------------------------------------------------------
# > st.var.fit <- grmsem.var(st.fit)
# > print(st.var.fit)
# $VA
# 1 2 3
# Y1 0.3160062 0.2799762 0.3387946
# Y2 0.2799762 0.8107680 0.5488882
# Y3 0.3387946 0.5488882 0.6499573
#
# $VA.se
# 1 2 3
# Y1 0.1932288 0.1619891 0.1613852
# Y2 0.1619891 0.1379212 0.1565957
# Y3 0.1613852 0.1565957 0.1544661
#
# $VE
# 1 2 3
# Y1 0.6839938 0.1866749 0.3413782
# Y2 0.1866749 0.1892320 0.2200922
# Y3 0.3413782 0.2200922 0.3500427
#
# $VE.se
# 1 2 3
# Y1 0.1932288 0.1218645 0.1436571
# Y2 0.1218645 0.1379212 0.1334911
# Y3 0.1436571 0.1334911 0.1544661
#
# $VP
# 1 2 3
# Y1 1.0000000 0.4666511 0.6801728
# Y2 0.4666511 1.0000000 0.7689803
# Y3 0.6801728 0.7689803 1.0000000
#
# $RG
# 1 2 3
# Y1 1.0000000 0.5531272 0.7475604
# Y2 0.5531272 1.0000000 0.7561242
# Y3 0.7475604 0.7561242 1.0000000
#
# $RG.se
# 1 2 3
# Y1 0.0000000 0.23210789 1.574432e-01
# Y2 0.2321079 0.00000000 8.945803e-02
# Y3 0.1574432 0.08945803 1.633888e-17
#
# $RE
# 1 2 3
# Y1 1.0000000 0.5188747 0.6976686
# Y2 0.5188747 1.0000000 0.8551589
# Y3 0.6976686 0.8551589 1.0000000
#
# $RE.se
# 1 2 3
# Y1 0.0000000 2.130183e-01 1.172975e-01
# Y2 0.2130183 6.290479e-17 8.880531e-02
# Y3 0.1172975 8.880531e-02 5.827606e-17
## ----eval=FALSE---------------------------------------------------------------
# > fit.DS <- grmsem.fit(ph.small,
# G.small,
# model="DS",
# LogL=TRUE,
# estSE=TRUE)
# [1] 300
# [1] 100 100 100
# [1] 1 101 201
# [1] 100 200 300
# [1] "Parameters for a direct symmetric model (variance components):"
# [1] "Vector of genetic variances:"
# [1] "A11" "A21" "A31" "A22" "A32" "A33"
# [1] "Vector of residual factor loadings:"
# [1] "E11" "E21" "E31" "E22" "E32" "E33"
#
# > fit.DS
# $model.in
# part label value freepar
# 1 a A11 1.0000000 1
# 2 a A21 0.4602739 1
# 3 a A31 0.6738565 1
# 4 a A22 1.0000000 1
# 5 a A32 0.7527887 1
# 6 a A33 1.0000000 1
# 7 e E11 1.0000000 1
# 8 e E21 0.4602739 1
# 9 e E31 0.6738565 1
# 10 e E22 1.0000000 1
# 11 e E32 0.7527887 1
# 12 e E33 1.0000000 1
#
# $model.fit$LL
# [1] -54.39145
#
# $model.fit$calls
# function gradient
# 122 26
#
# $model.fit$convergence
# [1] 0
#
# $model.fit$message
# NULL
#
#
# $model.out
# label estimates gradient se Z p
# 1 A11 0.3161061 0.0009860302 0.2106008 1.500972 0.1333627052
# 2 A21 0.2912542 0.0039066047 0.1834675 1.587497 0.1124001236
# 3 A31 0.3389444 -0.0036802953 0.1852176 1.829979 0.0672530362
# 4 A22 0.8754990 0.0025629175 0.2584953 3.386905 0.0007068577
# 5 A32 0.5705239 -0.0051792201 0.2222947 2.566520 0.0102724710
# 6 A33 0.6497819 0.0048300247 0.2229978 2.913849 0.0035700307
# 7 E11 0.6828573 0.0009099445 0.1921358 3.554035 0.0003793684
# 8 E21 0.1932362 0.0081933028 0.1180255 1.637241 0.1015800561
# 9 E31 0.3403784 -0.0046169601 0.1387849 2.452561 0.0141843201
# 10 E22 0.2034536 0.0015010666 0.1268743 1.603584 0.1088057273
# 11 E32 0.2277079 0.0414495710 0.1222574 1.862529 0.0625285683
# 12 E33 0.3487737 0.0049090361 0.1373909 2.538551 0.0111312504
#
# $k
# [1] 3
#
# $n
# [1] 300
#
# $n.obs
# [1] 100 100 100
#
# $n.ind
# [1] 100
#
# $model
# [1] "DS"
## ----eval=FALSE---------------------------------------------------------------
# > var.fit.ds <- grmsem.var(fit.DS)
# > var.fit.ds
# $VA
# 1 2 3
# Y1 0.3161061 0.2912542 0.3389444
# Y2 0.2912542 0.8754990 0.5705239
# Y3 0.3389444 0.5705239 0.6497819
#
# $VA.se
# 1 2 3
# Y1 0.2106008 0.1834675 0.1852176
# Y2 0.1834675 0.2584953 0.2222947
# Y3 0.1852176 0.2222947 0.2229978
#
# $VE
# 1 2 3
# Y1 0.6828573 0.1932362 0.3403784
# Y2 0.1932362 0.2034536 0.2277079
# Y3 0.3403784 0.2277079 0.3487737
#
# $VE.se
# 1 2 3
# Y1 0.1921358 0.1180255 0.1387849
# Y2 0.1180255 0.1268743 0.1222574
# Y3 0.1387849 0.1222574 0.1373909
#
# $VP
# 1 2 3
# Y1 0.9989633 0.4844903 0.6793228
# Y2 0.4844903 1.0789526 0.7982318
# Y3 0.6793228 0.7982318 0.9985556
#
# $RG
# 1 2 3
# Y1 1.0000000 0.5536406 0.7478738
# Y2 0.5536406 1.0000000 0.7564186
# Y3 0.7478738 0.7564186 1.0000000
#
# $RG.se
# 1 2 3
# Y1 0.0000000 0.22323218 0.15738360
# Y2 0.2232322 0.00000000 0.08684325
# Y3 0.1573836 0.08684325 0.00000000
#
# $RE
# 1 2 3
# Y1 1.0000000 0.5184307 0.6974693
# Y2 0.5184307 1.0000000 0.8548176
# Y3 0.6974693 0.8548176 1.0000000
#
# $RE.se
# 1 2 3
# Y1 0.0000000 0.20548054 0.11753303
# Y2 0.2054805 0.00000000 0.08761989
# Y3 0.1175330 0.08761989 0.00000000
## ----eval = FALSE-------------------------------------------------------------
# #Do not run!
# #Please downloaded externally
# #load("G.large.RData")
# #load("ph.large.RData")
# #fit <- grmsem.fit(ph.large,
# # G.large,model="Cholesky",
# # LogL=TRUE,
# # estSE=TRUE)
#
# > load("fit.large.RData")
# > fit.large
# $model.in
# part label value freepar
# 1 a a11 0.01686961 1
# 2 a a21 0.75022252 1
# 3 a a31 -0.10352849 1
# 4 a a41 -0.81445295 1
# 5 a a22 -0.34964589 1
# 6 a a32 0.70013464 1
# 7 a a42 0.30439698 1
# 8 a a33 0.87503509 1
# 9 a a43 -0.37434617 1
# 10 a a44 0.86376688 1
# 11 e e11 -0.28008894 1
# 12 e e21 -0.92449828 1
# 13 e e31 -0.72312602 1
# 14 e e41 -0.70924504 1
# 15 e e22 0.51178943 1
# 16 e e32 0.66207502 1
# 17 e e42 -0.64331456 1
# 18 e e33 0.16870994 1
# 19 e e43 0.30222504 1
# 20 e e44 0.21833386 1
#
# $model.fit$LL
# [1] -4556.647
#
# $model.fit$calls
# function gradient
# 488 118
#
# $model.fit$convergence
# [1] 0
#
# $model.fit$message
# NULL
#
# $model.out
# label estimates gradient se Z p
# 1 a11 -0.590283919 7.196924e-03 0.018936941 -31.1710286 2.632163e-213
# 2 a21 -0.406228554 4.760323e-04 0.022331488 -18.1908415 6.099640e-74
# 3 a31 -0.412282704 -1.306064e-03 0.021811928 -18.9017085 1.104185e-79
# 4 a41 -0.251181508 4.515862e-03 0.025355870 -9.9062469 3.910610e-23
# 5 a22 -0.611871764 7.519419e-03 0.013081543 -46.7736703 0.000000e+00
# 6 a32 -0.412494347 5.020765e-05 0.015153798 -27.2205248 3.712770e-163
# 7 a42 -0.406496327 2.858898e-03 0.020248170 -20.0757069 1.203583e-89
# 8 a33 0.443344860 -2.806537e-04 0.013549030 32.7215213 7.719589e-235
# 9 a43 -0.015668913 1.437494e-03 0.024698450 -0.6344088 5.258141e-01
# 10 a44 0.655953940 -5.747793e-03 0.014789225 44.3535029 0.000000e+00
# 11 e11 -0.794244501 2.325624e-02 0.012513836 -63.4693048 0.000000e+00
# 12 e21 -0.290249172 -8.449568e-04 0.013042715 -22.2537381 1.037729e-109
# 13 e31 -0.445243952 -5.908758e-03 0.012154714 -36.6313812 9.056960e-294
# 14 e41 0.006719942 7.527450e-03 0.011697368 0.5744833 5.656408e-01
# 15 e22 0.521077401 -1.251593e-02 0.008841668 58.9342896 0.000000e+00
# 16 e32 -0.054613913 1.163049e-03 0.010020443 -5.4502495 5.029922e-08
# 17 e42 -0.018229268 -8.354153e-03 0.011967918 -1.5231779 1.277142e-01
# 18 e33 0.405590017 -2.845990e-02 0.007471397 54.2856995 0.000000e+00
# 19 e43 -0.022769655 7.889148e-04 0.012392180 -1.8374213 6.614772e-02
# 20 e44 -0.467128074 1.039490e-02 0.009007167 -51.8618176 0.000000e+00
#
# $k
# [1] 4
#
# $n
# [1] 20000
#
# $n.obs
# [1] 5000 5000 5000 5000
#
# $n.ind
# [1] 5000
#
# $model
# [1] "Cholesky"
## ----eval = FALSE-------------------------------------------------------------
# > var.fit <- grmsem.var(fit.large)
# > print(var.fit)
# $VA
# 1 2 3 4
# Y1 0.3484351 0.2397902 0.2433639 0.1482684
# Y2 0.2397902 0.5394087 0.4198747 0.3507607
# Y3 0.2433639 0.4198747 0.5366833 0.2642885
# Y4 0.1482684 0.3507607 0.2642885 0.6588525
#
# $VA.se
# 1 2 3 4
# Y1 0.02235634 0.01660915 0.01767332 0.01514868
# Y2 0.01660915 0.02017283 0.01672469 0.01562045
# Y3 0.01767332 0.01672469 0.02042099 0.01515375
# Y4 0.01514868 0.01562045 0.01515375 0.02078249
#
# $VE
# 1 2 3 4
# Y1 0.630824327 0.23052881 0.35363256 -0.005337277
# Y2 0.230528809 0.35576624 0.10077361 -0.011449317
# Y3 0.353632560 0.10077361 0.36572812 -0.011231587
# Y4 -0.005337277 -0.01144932 -0.01123159 0.219104559
#
# $VE.se
# 1 2 3 4
# Y1 0.019878092 0.012210347 0.013679555 0.009287867
# Y2 0.012210347 0.012054860 0.009085111 0.007093293
# Y3 0.013679555 0.009085111 0.012413525 0.007220553
# Y4 0.009287867 0.007093293 0.007220553 0.008358356
#
# $VP
# 1 2 3 4
# Y1 0.9792594 0.4703190 0.5969964 0.1429311
# Y2 0.4703190 0.8951749 0.5206483 0.3393114
# Y3 0.5969964 0.5206483 0.9024114 0.2530569
# Y4 0.1429311 0.3393114 0.2530569 0.8779571
#
# $RG
# 1 2 3 4
# Y1 1.0000000 0.553110 0.5627767 0.3094522
# Y2 0.5531100 1.000000 0.7803720 0.5883800
# Y3 0.5627767 0.780372 1.0000000 0.4444523
# Y4 0.3094522 0.588380 0.4444523 1.0000000
#
# $RG.se
# 1 2 3 4
# Y1 0.00000000 0.02488356 2.233207e-02 3.056075e-02
# Y2 0.02488356 0.00000000 1.557784e-02 1.976462e-02
# Y3 0.02233207 0.01557784 3.008489e-18 2.206251e-02
# Y4 0.03056075 0.01976462 2.206251e-02 4.557191e-18
#
# $RE
# 1 2 3 4
# Y1 1.00000000 0.48661851 0.73623906 -0.01435621
# Y2 0.48661851 1.00000000 0.27937355 -0.04100828
# Y3 0.73623906 0.27937355 1.00000000 -0.03967676
# Y4 -0.01435621 -0.04100828 -0.03967676 1.00000000
#
# $RE.se
# 1 2 3 4
# Y1 1.111452e-17 0.01770789 1.099568e-02 2.499190e-02
# Y2 1.770789e-02 0.00000000 2.164591e-02 2.549970e-02
# Y3 1.099568e-02 0.02164591 5.397777e-18 2.555311e-02
# Y4 2.499190e-02 0.02549970 2.555311e-02 7.971451e-18
#
## ----eval = FALSE-------------------------------------------------------------
# > st.fit <- grmsem.stpar(fit.large)
# > print(st.fit)
#
# $stand.model.out
# label stand.estimates stand.se Z p
# 1 a11 -0.596502229 0.016233733 -36.7446127 1.417362e-295
# 2 a21 -0.429354966 0.020456160 -20.9890304 8.261383e-98
# 3 a31 -0.434003098 0.019542251 -22.2084494 2.845907e-109
# 4 a41 -0.268071735 0.024887619 -10.7712891 4.703688e-27
# 5 a22 -0.646705353 0.011739313 -55.0888607 0.000000e+00
# 6 a32 -0.434225891 0.014761620 -29.4158701 3.441759e-190
# 7 a42 -0.433830407 0.018983259 -22.8533157 1.354554e-115
# 8 a33 0.466701710 0.013216780 35.3113009 3.939089e-273
# 9 a43 -0.016722540 0.024694677 -0.6771718 4.982969e-01
# 10 a44 0.700062328 0.012710221 55.0786896 0.000000e+00
# 11 e11 -0.802611420 0.012064939 -66.5242831 0.000000e+00
# 12 e21 -0.306772929 0.012940717 -23.7060220 3.125266e-124
# 13 e31 -0.468700852 0.012081674 -38.7943628 0.000000e+00
# 14 e41 0.007171812 0.011697423 0.6131104 5.398033e-01
# 15 e22 0.550742107 0.009617573 57.2641473 0.000000e+00
# 16 e32 -0.057491151 0.010016353 -5.7397289 9.482821e-09
# 17 e42 -0.019455060 0.011964987 -1.6259993 1.039498e-01
# 18 e33 0.426957819 0.008330528 51.2521942 0.000000e+00
# 19 e43 -0.024300758 0.012391052 -1.9611537 4.986110e-02
# 20 e44 -0.498539222 0.010093805 -49.3906111 0.000000e+00
#
# $k
# [1] 4
#
# $n
# [1] 20000
#
# $model
# [1] "Cholesky"
## ----eval=FALSE---------------------------------------------------------------
# > st.var.fit <- grmsem.var(st.fit)
# > print(st.var.fit)
# $VA
# 1 2 3 4
# Y1 0.3558149 0.2561112 0.2588838 0.1599054
# Y2 0.2561112 0.6025735 0.4671576 0.3956584
# Y3 0.2588838 0.4671576 0.5947213 0.2969199
# Y4 0.1599054 0.3956584 0.2969199 0.7504382
#
# $VA.se
# 1 2 3 4
# Y1 0.01936692 0.01501072 0.01589623 0.01481065
# Y2 0.01501072 0.01344962 0.01252746 0.01320059
# Y3 0.01589623 0.01252746 0.01377554 0.01417036
# Y4 0.01481065 0.01320059 0.01417036 0.01000981
#
# $VE
# 1 2 3 4
# Y1 0.644185091 0.24621946 0.37618466 -0.005756178
# Y2 0.246219456 0.39742650 0.11212194 -0.012914839
# Y3 0.376184656 0.11212194 0.40527870 -0.012618339
# Y4 -0.005756178 -0.01291484 -0.01261834 0.249561817
#
# $VE.se
# 1 2 3 4
# Y1 0.019366916 0.012115870 0.013465601 0.009385665
# Y2 0.012115870 0.013449618 0.009581551 0.007494475
# Y3 0.013465601 0.009581551 0.013775535 0.007600587
# Y4 0.009385665 0.007494475 0.007600587 0.010009812
#
# $VP
# 1 2 3 4
# Y1 1.0000000 0.5023307 0.6350685 0.1541492
# Y2 0.5023307 1.0000000 0.5792795 0.3827435
# Y3 0.6350685 0.5792795 1.0000000 0.2843016
# Y4 0.1541492 0.3827435 0.2843016 1.0000000
#
# $RG
# 1 2 3 4
# Y1 1.0000000 0.553110 0.5627767 0.3094522
# Y2 0.5531100 1.000000 0.7803720 0.5883800
# Y3 0.5627767 0.780372 1.0000000 0.4444523
# Y4 0.3094522 0.588380 0.4444523 1.0000000
#
# $RG.se
# 1 2 3 4
# Y1 7.209226e-18 2.354325e-02 2.121442e-02 2.863522e-02
# Y2 2.354325e-02 4.542180e-18 1.479096e-02 1.851946e-02
# Y3 2.121442e-02 1.479096e-02 5.454486e-18 2.069598e-02
# Y4 2.863522e-02 1.851946e-02 2.069598e-02 4.092585e-18
#
# $RE
# 1 2 3 4
# Y1 1.00000000 0.48661851 0.73623906 -0.01435621
# Y2 0.48661851 1.00000000 0.27937355 -0.04100828
# Y3 0.73623906 0.27937355 1.00000000 -0.03967676
# Y4 -0.01435621 -0.04100828 -0.03967676 1.00000000
#
# $RE.se
# 1 2 3 4
# Y1 0.00000000 1.675409e-02 1.044538e-02 2.341725e-02
# Y2 0.01675409 4.271060e-18 2.052641e-02 2.389309e-02
# Y3 0.01044538 2.052641e-02 6.732149e-18 2.394314e-02
# Y4 0.02341725 2.389309e-02 2.394314e-02 4.512840e-18
#
## ----eval=FALSE---------------------------------------------------------------
# > grmsem.fcoher(fit.large)
# > grmsem.fcoher(fit.large)$fcoher.model.out[,c(1,7:10)]
#
# label Vi Vi.se FCOHER FCOHER.se
# 1 a11 3.484351e-01 0.0223563431 1.0000000000 0.000000000
# 2 a21 1.650216e-01 0.0181433759 0.3059306238 0.027526691
# 3 a31 1.699770e-01 0.0179853613 0.3167175772 0.025135935
# 4 a41 6.309215e-02 0.0127378513 0.0957606594 0.018914180
# 5 a22 3.743871e-01 0.0160084532 0.6940693762 0.027526691
# 6 a32 1.701516e-01 0.0125017124 0.3170428312 0.024459178
# 7 a42 1.652393e-01 0.0164616136 0.2507985687 0.023183443
# 8 a33 1.965547e-01 0.0120137852 0.3662395917 0.021786386
# 9 a43 2.455148e-04 0.0007739957 0.0003726401 0.001174465
# 10 a44 4.302756e-01 0.0194021010 0.6530681318 0.023536859
# 11 e11 6.308243e-01 0.0198780916 NA NA
# 12 e21 8.424458e-02 0.0075712747 NA NA
# 13 e31 1.982422e-01 0.0108236256 NA NA
# 14 e41 4.515762e-05 0.0001572113 NA NA
# 15 e22 2.715217e-01 0.0092143864 NA NA
# 16 e32 2.982679e-03 0.0010945112 NA NA
# 17 e42 3.323062e-04 0.0004363328 NA NA
# 18 e33 1.645033e-01 0.0060606482 NA NA
# 19 e43 5.184572e-04 0.0005643313 NA NA
# 20 e44 2.182086e-01 0.0084150015 NA NA
#
## ----eval=FALSE---------------------------------------------------------------
# > grmsem.fcoher(fit.large)
# > grmsem.fcoher(fit.large)$fcoher.model.out[,c(1,7:8,13:14)]
#
# label Vi Vi.se FCOENV FCOENV.se
# 1 a11 3.484351e-01 0.0223563431 NA NA
# 2 a21 1.650216e-01 0.0181433759 NA NA
# 3 a31 1.699770e-01 0.0179853613 NA NA
# 4 a41 6.309215e-02 0.0127378513 NA NA
# 5 a22 3.743871e-01 0.0160084532 NA NA
# 6 a32 1.701516e-01 0.0125017124 NA NA
# 7 a42 1.652393e-01 0.0164616136 NA NA
# 8 a33 1.965547e-01 0.0120137852 NA NA
# 9 a43 2.455148e-04 0.0007739957 NA NA
# 10 a44 4.302756e-01 0.0194021010 NA NA
# 11 e11 6.308243e-01 0.0198780916 1.0000000000 0.0000000000
# 12 e21 8.424458e-02 0.0075712747 0.2367975722 0.0172339750
# 13 e31 1.982422e-01 0.0108236256 0.5420479497 0.0161908960
# 14 e41 4.515762e-05 0.0001572113 0.0002061008 0.0007175779
# 15 e22 2.715217e-01 0.0092143864 0.7632024278 0.0172339750
# 16 e32 2.982679e-03 0.0010945112 0.0081554556 0.0029995986
# 17 e42 3.323062e-04 0.0004363328 0.0015166558 0.0019944184
# 18 e33 1.645033e-01 0.0060606482 0.4497965946 0.0162294819
# 19 e43 5.184572e-04 0.0005643313 0.0023662548 0.0025782116
# 20 e44 2.182086e-01 0.0084150015 0.9959109886 0.0034506100
## ----eval = FALSE-------------------------------------------------------------
# > fit.var <- grmsem.var(fit)
# > load("ph.large.RData")
# > grmsem.biher(ph.large, fit.var)
# $VPO
# 1 2 3 4
# Y1 1.0000000 0.5061277 0.6278417 0.1792348
# Y2 0.5061277 1.0000000 0.6167281 0.4352479
# Y3 0.6278417 0.6167281 1.0000000 0.3331679
# Y4 0.1792348 0.4352479 0.3331679 1.0000000
#
# $VA
# 1 2 3 4
# Y1 0.3484351 0.2397902 0.2433639 0.1482684
# Y2 0.2397902 0.5394087 0.4198747 0.3507607
# Y3 0.2433639 0.4198747 0.5366833 0.2642885
# Y4 0.1482684 0.3507607 0.2642885 0.6588525
#
# $BIHER
# 1 2 3 4
# Y1 0.3484351 0.4737741 0.3876197 0.8272299
# Y2 0.4737741 0.5394087 0.6808100 0.8058872
# Y3 0.3876197 0.6808100 0.5366833 0.7932591
# Y4 0.8272299 0.8058872 0.7932591 0.6588525
#
# $BIHER.se
# 1 2 3 4
# Y1 0.02235634 0.03281612 0.02814933 0.08451862
# Y2 0.03281612 0.02017283 0.02711842 0.03588862
# Y3 0.02814933 0.02711842 0.02042099 0.04548383
# Y4 0.08451862 0.03588862 0.04548383 0.02078249
#
# $BIHER.Z
# 1 2 3 4
# Y1 15.585514 14.43723 13.77012 9.787546
# Y2 14.437234 26.73937 25.10508 22.455231
# Y3 13.770125 25.10508 26.28097 17.440466
# Y4 9.787546 22.45523 17.44047 31.702288
#
# $BIHER.p
# 1 2 3 4
# Y1 9.132538e-55 3.017077e-47 3.855470e-43 1.273488e-22
# Y2 3.017077e-47 1.641347e-157 4.377379e-139 1.137653e-111
# Y3 3.855470e-43 4.377379e-139 3.165351e-152 4.067453e-68
# Y4 1.273488e-22 1.137653e-111 4.067453e-68 1.444882e-220
#
## ----eval = FALSE-------------------------------------------------------------
# #Do not run!
# #Please downloaded externally
# #load("G.large.RData")
# #load("ph.large.RData")
# #> out<-grmsem.fit(ph.large,
# # G.large,
# # model="DS",
# # LogL=TRUE,
# # estSE=TRUE)
#
# > out
# $model.in
# part label value freepar
# 1 a A11 1.0000000 1
# 2 a A21 0.5061277 1
# 3 a A31 0.6278417 1
# 4 a A41 0.1792348 1
# 5 a A22 1.0000000 1
# 6 a A32 0.6167281 1
# 7 a A42 0.4352479 1
# 8 a A33 1.0000000 1
# 9 a A43 0.3331679 1
# 10 a A44 1.0000000 1
# 11 e E11 1.0000000 1
# 12 e E21 0.5061277 1
# 13 e E31 0.6278417 1
# 14 e E41 0.1792348 1
# 15 e E22 1.0000000 1
# 16 e E32 0.6167281 1
# 17 e E42 0.4352479 1
# 18 e E33 1.0000000 1
# 19 e E43 0.3331679 1
# 20 e E44 1.0000000 1
#
#
# $model.fit$LL
# [1] -4556.647
#
# $model.fit$calls
# function gradient
# 173 38
#
# $model.fit$convergence
# [1] 0
#
# $model.fit$message
# NULL
#
#
# $model.out
# label estimates gradient se Z p
# 1 A11 0.348458001 -0.052226894 0.022356069 15.5867299 8.960447e-55
# 2 A21 0.239787813 0.073256044 0.016608680 14.4374998 3.005481e-47
# 3 A31 0.243373776 0.046398011 0.017672850 13.7710544 3.806174e-43
# 4 A41 0.148252661 0.045240115 0.015148427 9.7866704 1.284564e-22
# 5 A22 0.539403719 -0.038929159 0.020172359 26.7397447 1.625068e-157
# 6 A32 0.419867512 0.053261468 0.016724145 25.1054698 4.334317e-139
# 7 A42 0.350754104 -0.026684946 0.015620099 22.4553060 1.135724e-111
# 8 A33 0.536686987 -0.032016186 0.020420456 26.2818311 3.094103e-152
# 9 A43 0.264276207 0.014920668 0.015153366 17.4400997 4.093589e-68
# 10 A44 0.658845787 0.011119193 0.020782130 31.7025151 1.434490e-220
# 11 E11 0.630774354 -0.070278866 0.019875702 31.7359534 4.961474e-221
# 12 E21 0.230506722 0.166057579 0.012209173 18.8797985 1.672229e-79
# 13 E31 0.353594799 0.308389144 0.013677755 25.8518155 2.321384e-147
# 14 E41 -0.005335735 -0.015130111 0.009287206 -0.5745253 5.656124e-01
# 15 E22 0.355757271 -0.040844105 0.012054375 29.5127094 1.977804e-191
# 16 E32 0.100760965 -0.097429983 0.009084347 11.0917126 1.376273e-28
# 17 E42 -0.011450535 0.016775106 0.007093098 -1.6143207 1.064579e-01
# 18 E33 0.365701992 -0.075051659 0.012412202 29.4631041 8.554205e-191
# 19 E43 -0.011230439 0.004879057 0.007220099 -1.5554412 1.198410e-01
# 20 E44 0.219102869 0.026457431 0.008358236 26.2140086 1.839860e-151
#
# $k
# [1] 4
#
# $n
# [1] 20000
#
# $n.obs
# [1] 5000 5000 5000 5000
#
# $n.ind
# [1] 5000
#
# $model
# [1] "DS"
## ----eval = FALSE-------------------------------------------------------------
# > var.out<-grmsem.var(out)
# > print(var.out)
# $VA
# 1 2 3 4
# Y1 0.3484580 0.2397878 0.2433738 0.4198675
# Y2 0.2397878 0.5394037 0.1482527 0.3507541
# Y3 0.2433738 0.4198675 0.5366870 0.2642762
# Y4 0.1482527 0.3507541 0.2642762 0.6588458
#
# $VA.se
# 1 2 3 4
# Y1 0.02235607 0.01660868 0.01767285 0.01672414
# Y2 0.01660868 0.02017236 0.01514843 0.01562010
# Y3 0.01767285 0.01672414 0.02042046 0.01515337
# Y4 0.01514843 0.01562010 0.01515337 0.02078213
#
# $VE
# 1 2 3 4
# Y1 0.630774354 0.23050672 0.353594799 0.10076097
# Y2 0.230506722 0.35575727 -0.005335735 -0.01145053
# Y3 0.353594799 0.10076097 0.365701992 -0.01123044
# Y4 -0.005335735 -0.01145053 -0.011230439 0.21910287
#
# $VE.se
# 1 2 3 4
# Y1 0.019875702 0.012209173 0.013677755 0.009084347
# Y2 0.012209173 0.012054375 0.009287206 0.007093098
# Y3 0.013677755 0.009084347 0.012412202 0.007220099
# Y4 0.009287206 0.007093098 0.007220099 0.008358236
#
# $VP
# 1 2 3 4
# Y1 0.9792324 0.4702945 0.5969686 0.5206285
# Y2 0.4702945 0.8951610 0.1429169 0.3393036
# Y3 0.5969686 0.5206285 0.9023890 0.2530458
# Y4 0.1429169 0.3393036 0.2530458 0.8779487
#
# $RG
# 1 2 3 4
# Y1 1.0000000 0.5530889 0.5627792 0.8762846
# Y2 0.5530889 1.0000000 0.2755402 0.5883746
# Y3 0.5627792 0.7803596 1.0000000 0.4444323
# Y4 0.3094107 0.5883746 0.4444323 1.0000000
#
# $RG.se
# 1 2 3 4
# Y1 0.00000000 0.02488339 0.02233109 0.03667066
# Y2 0.02488339 0.00000000 0.02690300 0.01976462
# Y3 0.02233109 0.01557817 0.00000000 0.02206248
# Y4 0.03055988 0.01976462 0.02206248 0.00000000
#
# $RE
# 1 2 3 4
# Y1 1.00000000 0.48659729 0.73621590 0.27103868
# Y2 0.48659729 1.00000000 -0.01479291 -0.04101331
# Y3 0.73621590 0.27935199 1.00000000 -0.03967428
# Y4 -0.01435269 -0.04101331 -0.03967428 1.00000000
#
# $RE.se
# 1 2 3 4
# Y1 4.413292e-18 0.01770803 1.099634e-02 0.02340894
# Y2 1.770803e-02 0.00000000 2.576743e-02 0.02549942
# Y3 1.099634e-02 0.02164575 5.512125e-18 0.02555251
# Y4 2.499120e-02 0.02549942 2.555251e-02 0.00000000
## ----eval = FALSE-------------------------------------------------------------
# gcta64 --reml-bivar 1 2 --grm-gz large.gcta --pheno large.gcta.phe
# --reml-bivar-lrt-rg 0 --thread-num 4 --out Y1_Y2 --reml-maxit 1000 >
# large.gcta.Y1Y2.log`
## ----eval = FALSE-------------------------------------------------------------
# Summary result of REML analysis (shortened output):
# Source Variance SE
# V(G)/Vp_tr1 0.361387 0.020470
# V(G)/Vp_tr2 0.635426 0.015053
# rG 0.557367 0.026039
## ----eval = FALSE-------------------------------------------------------------
# load("ph.large.RData")
# load("G.large.RData")
# ph.biv<-ph.large[,c(1,2)]
# fit <- grmsem.fit(ph.biv, G.large, LogL = TRUE, estSE = TRUE)
# print(fit)
# st.fit <- grmsem.stpar(fit)
# print(st.fit)
# st.var.fit <- grmsem.var(st.fit)
# print(st.var.fit)
## ----eval = FALSE-------------------------------------------------------------
# $VA
# 1 2
# Y1 0.3613708 0.2669526
# Y2 0.2669526 0.6349136
#
# $VA.se
# 1 2
# Y1 0.01943807 0.01555813
# Y2 0.01555813 0.01355490
#
# $RG
# 1 2
# Y1 1.0000000 0.5573145
# Y2 0.5573145 1.0000000
#
# $RG.se
# 1 2
# Y1 7.179876e-18 0.0239077
# Y2 2.390770e-02 0.0000000
#
#
Try the grmsem package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.