TSEM user specified model

[Izel Erdogan]

Hi, 

First of all, thank you so much for your work!

My question is:

Im running TSEM at the moment and wanted to ask if someone know more about this error: 

This is the model: 

model_3factor <-"F1=~ BMI + FIns

F2=~ adhd_total_symp + BMI + FGlu + HbA1c

F3=~ ocs_ocd + HbA1c

F1~~F2

F1~~F3

F2~~F3

F1 ~ Gene

F2 ~ Gene

F3 ~ Gene

This is the command:

insulin_factors <- userGWAS(covstruc=LDSCoutput_TSEM,SNPs=fusion_stats, estimation = "DWLS", model = model3, printwarn = TRUE, sub=c("F1~Gene","F2~Gene"), toler = FALSE, SNPSE = FALSE, parallel = TRUE,  GC="standard", smooth_check=TRUE,TWAS=TRUE)

This is the error :

      Error in crossprod(wls.diff, wls.v) : non-conformable arguments

which I remember about not specifying ~SNP effects in the model. 

To see if it works, I also ran commonfactorGWAS function but then I get a row error that I cannot debug where it comes from (because none of the constructs have 32 cols or rows):

Error in `$<-.data.frame`(`*tmp*`, "free", value = c(0, 0, 0, 0, 0, 0,  :

  replacement has 24 rows, data has 32

Calls: commonfactorGWAS -> .commonfactorGWAS_main -> $<- -> $<-.data.frame

Thank you so much in advance!!

Best regards,

Izel 

[agro...@gmail.com]

Hi Izel, 

Is there any chance you have some variables in your LDSCoutput that you are not including in your model? 

I'm noticing also that your model is labeled model_3factor but the code just calls model=model3. Neither of these would explain why you get the error for commonfactorGWAS. Can you paste the code you ran for that part? 

Best, 

  Andrew

[Izel Erdogan]

Hi Andrew,

So sorry, that was a mistake.

This is the right model:

model3 <-"F1=~ adhd_total_symp + BMI + HOMA_IR + FIns
  F2=~ adhd_total_symp + ocs_ocd + HbA1c
  F3=~ FGlu

  F1 ~ Gene
  F2 ~ Gene
  F3 ~ Gene

  adhd_total_symp ~~ a*adhd_total_symp
  a > .001
  HbA1c ~~b*HbA1c
  b > .001
  FIns ~~c*FIns
  c > .001
"

Variables are the same, I am getting this error if I try on first 100:

fusion_stats100=fusion_stats[1:100,]
insulin_factors <- userGWAS(covstruc=LDSCoutput_TSEM,SNPs=fusion_stats100, estimation = "DWLS", model = model3, printwarn = TRUE, sub=c("F1~Gene","F2~Gene"), toler = FALSE, SNPSE = FALSE, parallel = FALSE,  GC="standard", smooth_check=TRUE,TWAS=TRUE)

[1] "Please note that an update was made to userGWAS on 11/21/19 so that it combines addSNPs and userGWAS."
found:  adhd_total_symp ocs_ocd BMI HOMA_IR FIns FGlu HbA1c
expected:  adhd_total_symp BMI HOMA_IR FIns ocs_ocd HbA1c FGlu Gene
Error in .rearrange(k = ncol(S_LD), fit = ReorderModelnoSNP, names = colnames(S_LD)) :
  object 'ReorderModelnoSNP' not found

which expects''Gene' in LDSC I guess?  or is writing   F1 ~ Gene in the model unnecessary?

Thank you so much!

Best,

Izel

[agro...@gmail.com]

Hi Izel, 

Thanks for raising this! You are specifying the effects of gene expression on the factors correctly. It turns out this was a bug with a recent update to the userGWAS function that will fix the measurement model (factor correlations and factor loadings) prior to estimating the SNP or Gene effects on the factors. As an aside, the reason for this update is that it really speeds up the model estimation and it improves interpretability in so far as you know that the same exact model is being applied across all SNPs or genes. 

The bug was that this update was coded to pick-up on estimate SNP effects in the model, but was a big oversight on my part that it should also be looking for Gene effects (as you correctly specified) in the model when TWAS=TRUE. I just put in an update that should fix this error if you could redownload the package and try your code after restarting R. The update also prints about this new fix_measurement update, which I plan to also have documented on the github wiki in the next week or so. 

Unrelated to the bug you were getting, one thing I noticed about your model is that you have a single-item latent for F3. Since that part of the model is only defined by the one variable, you only have that variables variance (the SNP-based heritability) as the degree of freedom. Your model estimates both the factor variance and the the residual variance of FGlu, so the df is -1. To avoid this I would recommend adding the line below to the model that fixes the residual variance to 0 so that F3 can be interpreted as entirely defined by that indicator. 

-Andrew

[Izel Erdogan]

Hi Andrew,

Thank you (and the team) so much for fixing it, everything runs smoothly now.

I indeed corrected the residual variance, thanks for pointing it out.

I have another question on a different topic so I will open a new conversation for that.

Wish you a nice week!

Best,

Izel