1 / 33

Univariate Model Fitting

Univariate Model Fitting. Sarah Medland QIMR – openMx workshop Brisbane 16/08/10. Univariate Twin Models. = using the twin design to assess the aetiology of one trait ( univariate ) Path Diagrams Basic ACE R Script. 1. Path Diagrams for Univariate Models. 1.0. 1.0. 1. 1. 1.

kineta
Download Presentation

Univariate Model Fitting

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Univariate Model Fitting Sarah Medland QIMR – openMx workshop Brisbane 16/08/10

  2. Univariate Twin Models = using the twin design to assess the aetiology of one trait (univariate) • Path Diagrams • Basic ACE R Script

  3. 1. Path Diagrams for Univariate Models

  4. 1.0 1.0 1 1 1 1 1 1 A C E E C A c c e a a e Twin 1 Trait Twin 2 Trait Basic Twin Model - MZ

  5. 0.5 1.0 1 1 1 1 1 1 A C E E C A c c e a a e Twin 1 Trait Twin 2 Trait Basic Twin Model - DZ

  6. rMZ = 1.0; rDZ = 0.5 rMZ = 1.0; rDZ = 1.0 1 1 1 1 1 1 A C E E C A c c e a a e Twin 1 Trait Twin 2 Trait Basic Twin Model

  7. The Variance • Since the variance of a variable is the covariance of the variable with itself, the expected variance will be the sum of all paths from the variable to itself, which follow Wright’s rules

  8. Variance for Twin 1 - A 1 1 a*1*a = a2 1 A C E c e a Twin 1 Trait

  9. Variance for Twin 1 - C 1 1 a*1*a = a2 c*1*c = c2 1 A C E c e a Twin 1 Trait

  10. Variance for Twin 1 - E a*1*a = a2 c*1*c = c2 e*1*e = e2 1 1 1 A C E c e a Twin 1 Trait

  11. 1 1 1 A C E c e a Twin 1 Trait Total Variance for Twin 1 a*1*a = a2 c*1*c = c2 e*1*e = e2 Total Variance =a2+c2+e2

  12. Covariance - MZ 1.0 a2+c2 1.0 1 1 1 1 1 1 A C E E C A c c e a a e Twin 1 Trait Twin 2 Trait

  13. Covariance - DZ 0.5 0.5a2+c2 1.0 1 1 1 1 1 1 A C E E C A c c e a a e Twin 1 Trait Twin 2 Trait

  14. Variance-Covariance Matrices

  15. Variance-Covariance Matrices

  16. Variance-Covariance Matrices

  17. Because, e2 refers to environmental influences UNIQUE to each twin. Therefore, this cannot explain why there is similarity between twins. Why isn’t e2 included in the covariance? Why is a2 only .5 for DZs but not MZs? Because DZ twins share on average half of their genes, whereas MZs share all of their genes.

  18. 2. Basic openMx ACE Script

  19. Overview • OpenMx script • Running the script • Describing the output

  20. 1 1 1 A C E c e a Twin 1 Trait ACE model Specify the matrices you need To build the model # Fit ACE Model with RawData and Matrices Input # ----------------------------------------------------------------------- univACEModel <- mxModel("univACE", mxModel("ACE", # Matrices a, c, and e to store a, c, and e path coefficients mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=10, label="a11", name="a" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=10, label="c11", name="c" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=10, label="e11", name="e" ), # Matrices A, C, and E compute variance components mxAlgebra( expression=a %*% t(a), name="A" ), mxAlgebra( expression=c %*% t(c), name="C" ), mxAlgebra( expression=e %*% t(e), name="E" ), # Algebra to compute total variances and standard deviations (diagonal only) mxAlgebra( expression=A+C+E, name="V" ), Do some algebra to get the variances

  21. 1 1 1 A C E c e a Twin 1 Trait Start values? a2 = additive genetic variance (A) c2 = Shared E variance (C) e2 = Non-shared E variance (E) Sum is modelled to be expected Total Variance Start Values for a, c, e:  (Total Variance / 3)

  22. 1 A a = * 1 V V V Twin 1 Trait inv = a/SD * SD = a 1/SD 1/SD Standardize parameter estimates # Algebra to compute total variances and standard deviations (diagonal only) mxAlgebra( expression=A+C+E, name="V" ), mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"), mxAlgebra( expression=solve(sqrt(I*V)), name="iSD"), # Algebra to compute standardized path estimares and variance components mxAlgebra( expression=a%*%iSD, name="sta"), mxAlgebra( expression=c%*%iSD, name="stc"), mxAlgebra( expression=e%*%iSD, name="ste"), mxAlgebra( expression=A/V, name="h2"), mxAlgebra( expression=C/V, name="c2"), mxAlgebra( expression=E/V, name="e2"), • The regression coefficient a is standardized by: (a * SD(A)) / SD(Trait) • where SD(Trait) is the standard deviation of the dependent variable, and SD(A) is the standard deviation of the predictor, the latent factor ‘A’ (=1)

  23. 1 A sta Twin 1 Trait 2 = / A V “h2” “sta” Standardize parameter estimates # Algebra to compute total variances and standard deviations (diagonal only) mxAlgebra( expression=A+C+E, name="V" ), mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"), mxAlgebra( expression=solve(sqrt(I*V)), name="iSD"), # Algebra to compute standardized path estimares and variance components mxAlgebra( expression=a%*%iSD, name="sta"), mxAlgebra( expression=c%*%iSD, name="stc"), mxAlgebra( expression=e%*%iSD, name="ste"), mxAlgebra( expression=A/V, name="h2"), mxAlgebra( expression=C/V, name="c2"), mxAlgebra( expression=E/V, name="e2"), The heritability ‘h2’ is the proportion of the total variance due to A (additive genetic effects; = A / V. Note: this will be “sta” squared. The standardized variance components for C and E are: C / V; E / V N

  24. Standardising data V=A+C+ E A/V= 73/233 = .31 V= [73] + [90] + [70] C/V = 90/233 = .39 V= [233] E/V= 70/233 = .30 a=8.55 c=9.49 e = 8.35 SD =sqrt(V)= 15.26 sta=8.55/ 15.26 =0.56 squared = .31 stc=9.49/ 15.26 =0.62squared =.39 ste=8.35/ 15.26 = 0.55 squared = .30

  25. 1/ 0.5 1.0 1 1 1 1 1 1 A C E E C A c c e a a e Twin 1 Trait Twin 2 Trait # Algebra for expected variance/covariance matrix in MZ mxAlgebra( expression= rbind ( cbind(A+C+E , A+C), cbind(A+C , A+C+E)), name="expCovMZ" ), # Algebra for expected variance/covariance matrix in DZ mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C), cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )

  26. mxModel("MZ", mxData( observed=mzData, type="raw" ), mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", dimnames=selVars ) ), mxModel("DZ", mxData( observed=dzData, type="raw" ), mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars ) ), mxAlgebra( expression=MZ.objective + DZ.objective, name="m2ACEsumll" ), mxAlgebraObjective("m2ACEsumll"), mxCI(c('ACE.A', 'ACE.C', 'ACE.E')) ) univACEFit <- mxRun(univACEModel, intervals=T) univACESumm <- summary(univACEFit) univACESumm

  27. A E E A Twin 1 Trait Twin 2 Trait Sub-models # Fit AE model # ----------------------------------------------------------------------- univAEModel <- mxModel(univACEFit, name="univAE", mxModel(univACEFit$ACE, mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="c11", name="c" ) ) ) univAEFit <- mxRun(univAEModel) univAESumm <- summary(univAEFit) univAESumm You can fit sub-models to test the significance of your parameters -you simply drop the parameter and see if the model fit is significantly worse than full model

  28. E E C E E C Twin 1 Trait Twin 2 Trait Twin 1 Trait Twin 2 Trait Sub-models # Fit CE model # ----------------------------------------------------------------------- univCEModel <- mxModel(univACEFit, name="univCE", mxModel(univACEFit$ACE, mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="a11", name="a" ) ) ) univCEFit <- mxRun(univCEModel) univCESumm <- summary(univCEFit) univCESumm # Fit E model # ----------------------------------------------------------------------- univEModel <- mxModel(univAEFit, name="univE", mxModel(univAEFit$ACE, mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="a11", name="a" ) ) ) univEFit <- mxRun(univEModel) univESumm <- summary(univEFit) univESumm The E parameter can never not be dropped because it includes measurement error

  29. OpenMx Output

  30. univACESumm free parameters: name matrix row col Estimate 1 a11 ACE.a 1 1 8.546504 2 c11 ACE.c 1 1 9.488454 3 e11 ACE.e 1 1 8.352197 4 mean ACE.Mean 1 1 94.147803 observed statistics: 2198 estimated parameters: 4 degrees of freedom: 2194 -2 log likelihood: 17639.91 saturated -2 log likelihood: NA number of observations: 1110 chi-square: NA p: NA AIC (Mx): 13251.91 BIC (Mx): 1127.663 adjusted BIC: RMSEA: NA

  31. tableFitStatisticsmodels compared to saturated model Name ep -2LL df AIC diffLLdiffdf p M1 : univTwinSat 10 17637.98 2188 13261.98 - - - M2 : univACE 4 17639.91 2194 13251.91 1.93 6 0.93 M3 : univAE 3 17670.38 2195 13280.38 32.4 7 0 M4 : univCE 3 17665.27 2195 13275.27 27.3 7 0 M5 : univE 2 18213.28 2196 13821.28 575.3 8 0 Smaller -2LL means better fit. -2LL of sub-model is always higher (worse fit). The question is: is it significantly worse. Chi-sq test: dif in -2LL is chi-square distributed. Evaluate sig of chi-sq test. A non-sig p-value means that the model is consistent with the data.

  32. Nested.fitmodels compared to ACE model Name ep -2LL df diffLL diffdf p univACE 4 17639.91 2194 NA NA NA univAE 3 17670.38 2195 32.47 1 0 univCE 3 17665.27 2195 25.36 1 0 univE 2 18213.28 2196 573.37 2 0 Smaller -2LL means better fit. -2LL of sub-model is always higher (worse fit). The question is: is it significantly worse. Chi-sq test: dif in -2LL is chi-square distributed. Evaluate sig of chi-sq test. Critical Chi-sq value for 1 DF = 3.84 A non-sig p-value means that the dropped parameter(s) are non-significant.

  33. Estimates ACE model > univACEFit$ACE.h2 [,1] [1,] 0.3137131 > univACEFit$ACE.c2 [,1] [1,] 0.3866762 > univACEFit$ACE.e2 [,1] [1,] 0.2996108

More Related