OSCR

Downregulated transcription in chromosomal domains of midbrain dopamine neurons linked to schizophrenia.

Code ↔ Paper

3 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 3 matches
  1. [1] § Methods › Data normalization, selection of covariates, and statistical analysis of differences in gene expression of cases and controls ↔ R/evalCriterion.R, lines 4–152 · score 0.65 · mvForwardStepwise, model selection, BIC, Covariates, gene
  2. [2] § Methods › Data normalization, selection of covariates, and statistical analysis of differences in gene expression of cases and controls ↔ vignette/seqc.Rmd, lines 115–136 · score 0.59 · mvForwardStepwise, model selection, BIC, gene
  3. [3] § Methods › remaCorr analysis of differences in Tx expression ↔ R/evalCriterion.R, lines 161–244 · score 0.51 · variancePartition, voom, limma, transformation, DREAM, linear

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 628 lines · 23 KB · no license · 2 matches

  1. # Gabriel Hoffman
  2. # April 12, 2020
  4. #' Multivariate forward stepwise regression
  5. #'
  6. #' Multivariate forward stepwise regression evluated by multivariate BIC
  7. #'
  8. #' @param exprObj matrix of expression data (g genes x n samples), or ExpressionSet, or EList returned by voom() from the limma package
  9. #' @param baseFormula specifies baseline variables for the linear (mixed) model. Must only specify covariates, since the rows of exprObj are automatically used as a response. e.g.: \code{~ a + b + (1|c)} Formulas with only fixed effects also work, and \code{lmFit()} followed by \code{contrasts.fit()} are run.
  10. #' @param data data.frame with columns corresponding to formula
  11. #' @param variables array of variable names to be considered in the regression. If variable should be considered as random effect, use '(1|A)'.
  12. #' @param criterion multivariate criterion ('AIC', 'BIC') or summing score assuming independence of reponses ('sum AIC', 'sum BIC')
  13. #' @param shrink.method Shrink covariance estimates to be positive definite. Using "var_equal" assumes all variance on the diagonal are equal. This method is the fastest because it is linear time. Using "var_unequal" allows each response to have its own variance term, however this method is quadratic time. Using "none" does not apply shrinkge, but is only valid when there are very few responses
  14. #' @param nparamsMethod "edf": effective degrees of freedom. "countLevels" count number of levels in each random effect. "lme4" number of variance compinents, as used by lme4. See description in \code{\link{nparam}}
  15. #' @param deltaCutoff stop interating of the model improvement is less than deltaCutoff. default is 5
  16. #' @param pca use PCA to transform variables
  17. #' @param verbose Default TRUE. Print messages
  18. #' @param ... additional arguements passed to logDet
  19. #'
  20. #' @return list with formula of final model, and trace of iterations during model selection
  21. #' @examples
  22. #'
  23. #' Y = with(iris, rbind(Sepal.Width, Sepal.Length))
  24. #'
  25. #' # fit forward stepwise regression starting with model: ~1.
  26. #' bestModel = mvForwardStepwise( Y, ~ 1, data=iris, variables=colnames(iris)[3:5])
  27. #'
  28. #' bestModel
  29. #'
  30. #' @import variancePartition
  31. #' @importFrom stats as.formula
  32. #' @importFrom dplyr as_tibble
  33. #' @export
  34. mvForwardStepwise = function( exprObj, baseFormula, data, variables, criterion = c( "BIC", "sum BIC", "AIC", "AICC", "CAIC", "sum AIC"), shrink.method = c( "EB", "none", "var_equal", "var_unequal"), nparamsMethod = c("edf", "countLevels", "lme4"), deltaCutoff = 5, pca = TRUE, verbose=TRUE,... ){
  36. criterion = match.arg(criterion)
  37. shrink.method = match.arg(shrink.method)
  38. nparamsMethod = match.arg(nparamsMethod)
  39. baseFormula = as.formula( baseFormula )
  41. if( ! is.data.frame(data) ){
  42. data = as.data.frame(data, stringsAsFactors=FALSE)
  43. }
  44. data = droplevels(data)
  46. # Apply PCA only once, instead of in each call to mvIC_fit
  47. if( pca ){
  48. if( is(exprObj, "matrix") ){
  49. exprObj = t(pcTransform(t(exprObj)))
  50. }else{
  51. exprObj = t(pcTransform(t(exprObj$E)))
  52. }
  53. }
  55. # score base model
  56. suppressWarnings({
  57. baseScore = mvIC_fit( exprObj, baseFormula, data, criterion=criterion, shrink.method=shrink.method, nparamsMethod=nparamsMethod, verbose=FALSE, pca=FALSE,...)
  58. })
  60. resTrace = data.frame( iter = 0,
  61. variable = as.character(baseFormula)[-1],
  62. delta = NA,
  63. score = as.numeric(baseScore),
  64. isBest = "yes",
  65. isAdded = "yes",
  66. stringsAsFactors=FALSE)
  67. resTrace = cbind(resTrace, baseScore@params)
  69. iteration = 1
  71. # run until break
  72. while(1){
  74. if( verbose ) message(paste0("Base model: ", paste0(baseFormula, collapse=' ')))
  76. # evaluate score of each potential model
  77. score = lapply( variables, function(feature){
  78. if( verbose ) message(paste("\r\tevaluating: +", feature, ' '), appendLF=FALSE)
  80. # formula of model to try
  81. form = as.formula(paste(paste0(baseFormula, collapse=' '), "+", feature))
  83. suppressWarnings({
  84. # evaluate multivariate score
  85. mvIC_fit( exprObj, form, data, criterion=criterion, shrink.method=shrink.method, nparamsMethod=nparamsMethod, verbose=FALSE, pca=FALSE, ...)
  86. })
  87. })
  89. # get index of minumum score
  90. i = which.min(unlist(score))
  92. # get difference between best and second best score
  93. delta = as.numeric(score[[i]] - baseScore)
  94. if( verbose ) message("\nBest model delta: ", format(delta, nsmall=1, digits=1))
  96. isBest = rep("", length(score))
  97. isAdded = rep("", length(score))
  98. isBest[i] = "yes"
  99. if( delta < -deltaCutoff ){
  100. isAdded[i] = "yes"
  101. }
  102. resNew = data.frame(iter = iteration,
  103. variable = variables,
  104. delta = as.numeric(unlist(score) - baseScore),
  105. score = unlist(score),
  106. isBest = isBest,
  107. isAdded = isAdded,
  108. stringsAsFactors=FALSE)
  110. # get summary stats from model fits
  111. params = lapply(score, function(x) x@params)
  112. params = do.call(rbind, params)
  114. # combine results from this iteration
  115. resNew = cbind(resNew, params)
  117. # combine with results from previous interations
  118. resTrace = rbind(resTrace, resNew)
  120. iteration = iteration + 1
  122. # evaluate of best model is better than existing model
  123. if( delta < -deltaCutoff ){
  124. if( verbose ) message("Add variable to model: ", variables[i], '\n')
  126. # set new model, baseScore and possible variable list
  127. baseFormula = as.formula(paste(paste0(baseFormula, collapse=' '), "+", variables[i]))
  128. baseScore = score[[i]]
  130. variables = variables[-i]
  132. # if there are no additional variables to try
  133. if( length(variables) == 0){
  134. break
  135. }
  137. }else{
  138. if( verbose ){
  139. message(paste0("\nFinal model:\n ", paste0(baseFormula, collapse=' ')))
  140. }
  141. break
  142. }
  143. }
  145. # remove some columns from resTrace that are constant
  146. idx = colnames(resTrace) %in% c("n", 'p', 'criterion', 'shrink.method')
  148. # return as mvIC_result object
  149. new("mvIC_result", list(formula = baseFormula,
  150. settings= resTrace[1,idx],
  151. trace = as_tibble(resTrace[,!idx]) ))
  152. }
  161. #' Evaluate multivariate BIC
  162. #'
  163. #' Evaluate multivariate BIC on matrix of response variables. Smaller is better.
  164. #'
  165. #' @param exprObj matrix of expression data (g genes x n samples), or ExpressionSet, or EList returned by voom() from the limma package
  166. #' @param formula specifies variables for the linear (mixed) model. Must only specify covariates, since the rows of exprObj are automatically used as a response. e.g.: \code{~ a + b + (1|c)} Formulas with only fixed effects also work, and \code{lmFit()} followed by \code{contrasts.fit()} are run.
  167. #' @param data data.frame with columns corresponding to formula
  168. #' @param criterion multivariate criterion ('AIC', 'BIC') or summing score assuming independence of reponses ('sum AIC', 'sum BIC')
  169. #' @param shrink.method Shrink covariance estimates to be positive definite. Using "var_equal" assumes all variance on the diagonal are equal. This method is the fastest because it is linear time. Using "var_unequal" allows each response to have its own variance term, however this method is quadratic time. Using "none" does not apply shrinkge, but is only valid when there are very few responses
  170. #' @param nparamsMethod "edf": effective degrees of freedom. "countLevels" count number of levels in each random effect. "lme4" number of variance compinents, as used by lme4. See description in \code{\link{nparam}}
  171. #' @param pca use PCA to transform variables
  172. #' @param verbose Default TRUE. Print messages
  173. #' @param ... additional arguements passed to logDet
  174. #'
  175. #' @description
  176. #' Evaluate multivariate BIC while considering correlation between response variables. For n samples, p responses and m parameters for each model, evaluate the multivariate BIC as \deqn{n * logDet(\Sigma) + log(n) * (p*m + 0.5*p*(p+1))}
  177. #' where \eqn{\Sigma} is the residual covariance matrix. This formula extends the standard univariate BIC to the multivariate case.
  178. #' For one response the standard penalty is \eqn{log(n)*m}, this just adds a \eqn{log(n)} to that value, but only the differences between two models is important. Estimating the \eqn{p x p} covariance matrix requires \eqn{0.5*p*(p+1)} parameters.
  179. #' When \eqn{p > m} the residual covariance matrix Sigma is not full rank. In this case the psudo-determinant is used instead.
  180. #'
  181. #' @references{
  182. #' \insertRef{pauler1998schwarz}{mvIC}
  183. #'
  184. #' \insertRef{bedrick1994model}{mvIC}
  185. #'
  186. #' \insertRef{wu2013weighted}{mvIC}
  187. #' }
  188. #'
  189. #' @return multivariate BIC value
  190. #' @examples
  191. #'
  192. #' # create matrix of responses
  193. #' Y = with(iris, rbind(Sepal.Width, Sepal.Length))
  194. #'
  195. #' # Evaluate model 1
  196. #' mvIC_fit( Y, ~ Species, data=iris)
  197. #'
  198. #' # Evaluate model 2
  199. #' # smaller mvIC means better model
  200. #' mvIC_fit( Y, ~ Petal.Width + Petal.Length + Species, data=iris)
  201. #'
  202. #' @import variancePartition Rdpack
  203. #'
  204. #' @export
  205. mvIC_fit = function( exprObj, formula, data, criterion = c( "BIC", "sum BIC", "AIC", "AICC", "CAIC", "sum AIC"), shrink.method = c( "EB", "none", "var_equal", "var_unequal"), nparamsMethod = c("edf", "countLevels", "lme4"), pca = TRUE, verbose=FALSE,... ){
  207. criterion = match.arg(criterion)
  208. shrink.method = match.arg(shrink.method)
  209. nparamsMethod = match.arg(nparamsMethod)
  210. formula = as.formula(formula)
  212. if( ! is.data.frame(data) ){
  213. data = as.data.frame(data, stringsAsFactors=FALSE)
  214. }
  216. # if pca
  217. if( pca ){
  218. if( is(exprObj, "matrix") ){
  219. exprObj = t(pcTransform(t(exprObj)))
  220. }else{
  221. exprObj = t(pcTransform(t(exprObj$E)))
  222. }
  223. }
  225. # fit model and compute residuals
  226. suppressWarnings({
  227. modelFit = dream( exprObj, formula, data, REML=FALSE, computeResiduals=TRUE, quiet=!verbose)
  228. })
  230. # extract residuals
  231. residMatrix = residuals( modelFit )
  233. # effect fixed
  234. if( modelFit$method == "ls"){
  235. m <- ncol(coef(modelFit)) #+ 1
  236. }else{
  238. # effective number of parameters is returned by dream
  239. m = mean(attr(modelFit, "edf"))
  240. attr(m,"nparamsMethod") = "edf"
  241. }
  243. mvIC_from_residuals( residMatrix, m, criterion=criterion, shrink.method=shrink.method,... )
  244. }
  248. #' residuals for MArrayLM
  249. #'
  250. #' residuals for MArrayLM
  251. #'
  252. #' @param object MArrayLM object from dream
  253. #' @param ... other arguments, currently ignored
  254. #'
  255. #' @return results of residuals
  256. #' @importFrom limma residuals.MArrayLM
  257. #' @rdname residuals-method
  258. #' @aliases residuals,MArrayLM-method
  259. setMethod("residuals", "MArrayLM",
  260. function( object, ...){
  261. if( is.null(object$residuals) ){
  262. # use residuals computed by limma
  263. res = residuals.MArrayLM( object, ...)
  264. }else{
  265. # use precomputed residuals
  266. res = object$residuals
  267. }
  268. res
  269. })
  271. #Matrix of residual from list of model fits
  272. #
  273. # Matrix of residual from list of model fits
  274. #
  275. # @param fitList list of model fits from \code{lm()} or \code{lmer()}
  276. #
  277. getResids = function(fitList){
  279. if(is(fitList, "list")){
  280. residMatrix = lapply(fitList, residuals)
  281. residMatrix = do.call(rbind, residMatrix)
  282. }else if(is(fitList, "mlm")){
  283. residMatrix = t(residuals( fitList))
  284. }else{
  285. residMatrix = t(residuals(fitList))
  286. }
  287. residMatrix
  288. }
  290. #' Number of parameters in model
  291. #'
  292. #' Number of parameters in model from \code{lm()} or \code{lmer()}
  293. #'
  294. #' @param object model fit by \code{lm()} or \code{lmer()}
  295. #' @param nparamsMethod "edf": effective degrees of freedom. "countLevels" count number of levels in each random effect. "lme4" number of variance compinents, as used by lme4. See description in \code{\link{nparam}}
  296. #'
  297. #' @description
  298. #' In the case of \code{lm()}, the result is the number of coefficients For a linear mixed model fit with \code{lmer()} there are 3 options. "edf": effective degrees of freedom as computed by sum of diagonal values of the hat matrix return by \code{lmer()} . "countLevels", returns the number of fixed effects + number of levels in random effects + 1 for residual variance term. This treats each level of a random effect as a parameter. "lme4", returns number of fixed effects + number of variance components. Here a random effect with 10 levels is only counted as 1 parameter. This tends to underpenalize.
  299. # , + 1 for the variance term.
  300. #' @return number of parameters
  301. #' @importFrom stats coef
  302. #' @importFrom methods is
  303. #' @importFrom stats hatvalues
  304. nparam = function( object, nparamsMethod = c("edf", "countLevels", "lme4")){
  306. nparamsMethod = match.arg(nparamsMethod)
  308. if( is(object, "list") & nparamsMethod == "edf" ){
  310. if( all(sapply(object, function(fit) is(fit, "merMod"))) ){
  311. # mean of effective degrees of freedom across all responses
  312. m = mean(sapply( object, function(fit) sum(hatvalues(fit))))
  313. attr(m, "nparamsMethod") = "edf"
  314. return(m)
  315. }
  316. }
  318. # if object is not any of these
  319. if( !is(object, "lm") & !is(object, "mlm") & !is(object, 'merMod') ){
  320. # see if element of list is valid model fit
  321. object = object[[1]]
  322. }
  324. if( is(object, "mlm") ){
  325. # must be evaluated first because if object is 'mlm', it is also 'lm'
  326. # need 'mlm' to take presidence
  327. m = nrow(coef(object)) #+ 1
  328. attr(m, "nparamsMethod") = "lm"
  329. }else if( is(object, "lm") ){
  330. m = length(coef(object)) #+ 1
  331. attr(m, "nparamsMethod") = "lm"
  332. }else if( is(object, "merMod") ){
  334. m = switch( nparamsMethod,
  335. # effective degrees of freedom
  336. # Add term for residual variance
  337. "edf" = sum(hatvalues(object)),
  339. # fixed + number of random levels
  340. "countLevels" = length(object@beta) + object@devcomp[["dims"]][['q']] + object@devcomp[["dims"]][["useSc"]],
  342. # lme4:::npar.merMod
  343. # counts each random effect as a single parameter
  344. "lme4" = length(object@beta) + length(object@theta) + object@devcomp[["dims"]][["useSc"]],
  346. "already estimated" = m
  347. )
  349. attr(m, "nparamsMethod") = nparamsMethod
  351. }else{
  352. stop("object is not a valid model fit from lm() or lmer()")
  353. }
  355. m
  356. }
  358. #' Evaluate multivariate BIC
  359. #'
  360. #' Evaluate multivariate BIC from a list of regression fits
  361. #'
  362. #' @param fitList list of model fits with \code{lm()} or \code{lmer()}. All models must have same data, response and formula.
  363. #' @param criterion multivariate criterion ('AIC', 'BIC') or summing score assuming independence of reponses ('sum AIC', 'sum BIC')
  364. #' @param shrink.method Shrink covariance estimates to be positive definite. Using "var_equal" assumes all variance on the diagonal are equal. This method is the fastest because it is linear time. Using "var_unequal" allows each response to have its own variance term, however this method is quadratic time. Using "none" does not apply shrinkge, but is only valid when there are very few responses
  365. #' @param nparamsMethod "edf": effective degrees of freedom. "countLevels" count number of levels in each random effect. "lme4" number of variance compinents, as used by lme4. See description in \code{\link{nparam}}
  366. #' @param ... additional arguements passed to logDet
  367. #'
  368. #' @description
  369. #' Evaluate multivariate BIC while considering correlation between response variables. For n samples, p responses and m parameters for each model, evaluate the multivariate BIC as \deqn{n * logDet(\Sigma) + log(n) * (p*m + 0.5*p*(p+1))}
  370. #' where \eqn{\Sigma} is the residual covariance matrix. This formula extends the standard univariate BIC to the multivariate case.
  371. #' For one response the standard penalty is \eqn{log(n)*m}, this just adds a \eqn{log(n)} to that value, but only the differences between two models is important. Estimating the \eqn{p x p} covariance matrix requires \eqn{0.5*p*(p+1)} parameters.
  372. #' When \eqn{p > m} the residual covariance matrix Sigma is not full rank. In this case the psudo-determinant is used instead.
  373. #'
  374. #' See References
  375. #'
  376. #' Pauler, DK. The Schwarz criterion and related methods for normal linear models. Biometrika (1998), 85, 1, pp. 13-27
  377. #'
  378. #' Edward J. Bedrick and Chih-Ling Tsai. Model Selection for Multivariate Regression in Small Samples. Biometrics, 50:1 1994 226-231
  379. #'
  380. #' TJ Wu, P Chen, Y Yan. The weighted average information criterion for multivariate regression model selection. Signal Processing 93.1 (2013): 49-55.
  381. #'
  382. #' @return multivariate BIC value
  383. #' @examples
  384. #' # Predict Sepal width and Length given Species
  385. #' # Evaluate model fit
  386. #' fit1 = lm( cbind(Sepal.Width, Sepal.Length) ~ Species, data=iris)
  387. #' mvIC( fit1 )
  388. #'
  389. #' # add Petal width and length
  390. #' # smaller mvIC means better model
  391. #' fit2 = lm( cbind(Sepal.Width, Sepal.Length) ~ Petal.Width + Petal.Length + Species, data=iris)
  392. #' mvIC( fit2 )
  393. #'
  394. #' @importFrom methods is
  395. #' @export
  396. mvIC = function( fitList, criterion = c( "BIC", "sum BIC", "AIC", "AICC", "CAIC", "sum AIC"), shrink.method = c( "EB", "none", "var_equal", "var_unequal"), nparamsMethod = c("edf", "countLevels", "lme4"), ...){
  398. criterion = match.arg(criterion)
  399. shrink.method = match.arg(shrink.method)
  400. nparamsMethod = match.arg(nparamsMethod)
  402. # get residuals for 'mlm', 'lm', or list of 'lm' or 'lmer'
  403. residMatrix = getResids( fitList )
  405. # get number of parameters for multiple forms of fitList
  406. m = nparam( fitList, nparamsMethod=nparamsMethod )
  408. mvIC_from_residuals( residMatrix, m, criterion=criterion, shrink.method=shrink.method,...)
  409. }
  413. #' Evaluate multivariate BIC from matrix of residuals
  414. #'
  415. #' Evaluate multivariate BIC from matrix of residuals
  416. #'
  417. #' @param residMatrix matrix of residuals where rows are features
  418. #' @param m number of parameters for each model
  419. #' @param criterion multivariate criterion ('AIC', 'BIC') or summing score assuming independence of reponses ('sum AIC', 'sum BIC')
  420. #' @param shrink.method Shrink covariance estimates to be positive definite. Using "var_equal" assumes all variance on the diagonal are equal. This method is the fastest because it is linear time. Using "var_unequal" allows each response to have its own variance term, however this method is quadratic time. Using "none" does not apply shrinkge, but is only valid when there are very few responses
  421. #' @param ... other arguments passed to logDet
  422. #'
  423. #' @importFrom methods new
  424. mvIC_from_residuals = function( residMatrix, m, criterion = c( "BIC", "sum BIC", "AIC", "AICC", "CAIC", "sum AIC"), shrink.method = c( "EB", "none", "var_equal", "var_unequal"), ... ){
  426. criterion = match.arg(criterion)
  427. shrink.method = match.arg(shrink.method)
  429. n = ncol(residMatrix) # number of samples
  430. p = nrow(residMatrix) # number of response variables
  432. if( criterion %in% c("AIC", "BIC", "AICC", "CAIC") ){
  434. if( criterion %in% c("AICC", "CAIC") & n < p){
  435. stop(paste("Criterion", criterion, "cannot be evaluated when n < p"))
  436. }
  438. # compute log determinant explicitly
  439. # slower and not defined for low rank matrices
  440. # dataTerm = n * determinant(crossprod(residMatrix), log=TRUE)$modulus[1]
  442. if( shrink.method == "EB"){
  444. # est_param = shrinkcovmat.equal_lambda( residMatrix )
  446. # responses are *rows*
  447. res = eclairs(t(residMatrix))
  448. lambda = res$lambda
  450. # b = beam::beam(t(residMatrix), verbose=FALSE)
  451. # lambda = b@alphaOpt
  452. # res = list(logLik = b@valOpt)
  454. # dataTerm = -2*res$logLik
  455. dataTerm = -2*res$logML
  456. gdf_cov = p + (1-lambda)*p*(p-1)/2
  458. }else{
  459. # Evaluate logDet based on shrink.method
  460. logDet = rlogDet( residMatrix, shrink.method,... )
  461. dataTerm = n * logDet
  463. # get effective number of parameter used to estimate covariance by shrinkage
  464. gdf_cov = attr(logDet, "param")$gdf
  465. lambda = attr(logDet, "param")$lambda
  466. }
  468. # see Yanagihara, et al. 2015
  469. # doi:10.1214/15-EJS1022
  470. # penalty = switch( criterion,
  471. # "AIC" = 2 * (p*(m-1) + gdf_cov),
  472. # "BIC" = log(n) * (p*(m-1) + gdf_cov),
  473. # "AICC" = 2 * n*(p*(m-1) + gdf_cov) / (n-(m-1) - p - 1),
  474. # "CAIC" = (1+log(n)) * (p*(m-1) + gdf_cov))
  476. penalty = switch( criterion,
  477. "AIC" = 2 * (p*m + gdf_cov),
  478. "BIC" = log(n) * (p*m + gdf_cov),
  479. "AICC" = 2 * n*(p*m + gdf_cov) / (n-m - p - 1),
  480. "CAIC" = (1+log(n)) * (p*m + gdf_cov))
  482. # retrun data term plus penalty
  483. res = dataTerm + penalty
  484. attr(res, 'params') = data.frame( n = n,
  485. p = p,
  486. m = as.numeric(m),
  487. dataTerm = dataTerm,
  488. penalty = penalty,
  489. lambda = lambda,
  490. df_cov = gdf_cov,
  491. criterion = criterion,
  492. shrink.method = shrink.method,
  493. stringsAsFactors=FALSE)
  495. if( ! is.null( attr(m, 'nparamsMethod') )){
  496. attr(res, 'nparamsMethod') = attr(m, 'nparamsMethod')
  497. }else{
  498. attr(res, 'nparamsMethod') = "lm"
  499. }
  500. }else{
  501. # Naive metric summing BIC from all models independently
  502. rss = apply(residMatrix, 1, function(x) sum(x^2))
  503. dataTerm = n*sum(log(rss/n))
  505. penalty = switch(criterion,
  506. "sum AIC" = 2 * m*p,
  507. "sum BIC" = log(n) * m*p)
  510. # retrun data term plus penalty
  511. res = dataTerm + penalty
  512. attr(res, 'params') = data.frame( n = n,
  513. p = p,
  514. m = as.numeric(m),
  515. dataTerm = dataTerm,
  516. penalty = penalty,
  517. lambda = NA,
  518. df_cov = NA,
  519. criterion = criterion,
  520. shrink.method = "none",
  521. stringsAsFactors=FALSE)
  522. attr(res, 'nparamsMethod') = "naive"
  523. }
  525. new("mvIC", as.numeric(res),
  526. nparamsMethod = attr(res, 'nparamsMethod'),
  527. params = attr(res, 'params'))
  528. }
  531. #' Class mvIC
  532. #'
  533. #' Class stores mvIC score, method and parameter values
  534. #'
  535. #' @name mvIC-class
  536. #' @rdname mvIC-class
  537. #' @exportClass mvIC
  538. setClass("mvIC", representation(nparamsMethod = "character", params="data.frame"), contains="numeric")
  542. # Print mvIC object
  543. #
  544. # Print mvIC object
  545. #
  546. # @param x mvIC object
  547. # @export
  548. setMethod("print", "mvIC", function( x ){
  550. cat("\t\tMultivariate IC score\n\n")
  551. cat(paste(" Samples:\t", x@params$n, "\n"))
  552. cat(paste(" Responses:\t", x@params$p, "\n"))
  553. cat(paste(" Coef param:\t", round(x@params$m, digits=1), "\n"))
  554. cat(paste(" Cov param:\t", round(x@params$df_cov, digits=1), "\n"))
  555. cat(paste(" Regression:\t", x@nparamsMethod), "\n")
  556. cat(" Shrink method:", x@params$shrink.method, "\n")
  557. cat(paste(" lambda:\t", format(x@params$lambda, digits=3), "\n"))
  558. cat(" Criterion:\t", x@params$criterion, "\n")
  559. cat(" Score:\t", as.numeric(x), "\n\n")
  560. })
  564. # Show mvIC object
  565. #
  566. # Show mvIC object
  567. #
  568. # @param object mvIC object
  569. # @export
  570. setMethod("show", "mvIC", function( object ){
  571. print( object )
  572. })
  579. #' Class mvIC_result
  580. #'
  581. #' Class stores result of \code{mvForwardStepwise}
  582. #'
  583. #' @name mvIC_result-class
  584. #' @rdname mvIC_result-class
  585. #' @exportClass mvIC_result
  586. # setClass("mvIC_result", representation(formula = "formula", settings="data.frame", trace="data.frame"))
  587. setClass("mvIC_result", contains="list")
  591. # Print mvIC_result object
  592. #
  593. # Print mvIC_result object
  594. #
  595. # @param x mvIC_result object
  596. # @export
  597. setMethod("print", "mvIC_result", function( x ){
  599. cat("\t\tMultivariate IC forward stepwise regression\n\n")
  600. cat(" Samples:\t", x$settings$n, '\n')
  601. cat(" Responses:\t", x$settings$p, '\n')
  602. cat(" Shrink method:", x$settings$shrink.method, '\n')
  603. cat(" Criterion:\t", x$settings$criterion, '\n')
  604. cat(" Iterations:\t", max(x$trace$iter), "\n\n")
  605. cat(' Best model:', paste(as.character(x$formula), collapse=" "), '\n\n')
  606. })
  610. # Show mvIC_result object
  611. #
  612. # Show mvIC_result object
  613. #
  614. # @param object mvIC_result object
  615. # @export
  616. setMethod("show", "mvIC_result", function( object ){
  617. print( object )
  618. })
  620. #' @importFrom lme4 findbars
  621. .isMixedModelFormula = function(formula){
  622. !is.null(findbars(as.formula(formula)))
  623. }

evalCriterion.R at commit 15f0b63, no license · at the source

Overview

Authors: Swadha Singh1,2, Marina Iskhakova1,2, Tova Y. Lambert1,2, Aditi Valada1,2, Neda Shokrian1,2, Viviana Evans1,2, Jaroslav Bendl1,3,4, Pavan K. Auluck5, Stefano Marenco5, Minghui Wang4,6,7, Bin Zhang4,6,7,8, Gabriel E. Hoffman1,3,4,9, Kiran Girdhar1,3,4, Panos Roussos1,2,3,4,9,10, Schahram Akbarian1,2,4
ORCID iDs: Swadha Singh, Tova Y. Lambert, Jaroslav Bendl, Pavan K. Auluck, Stefano Marenco, Minghui Wang, Gabriel E. Hoffman, Kiran Girdhar, Panos Roussos, Schahram Akbarian
  1. Department of Psychiatry, Icahn School of Medicine at Mount Sinai,New York, NY USA
  2. Friedman Brain Institute, Icahn School of Medicine at Mount Sinai,New York, NY USA
  3. Center for Disease Neurogenomics, Icahn School of Medicine at Mount Sinai,New York, NY USA
  4. Department of Genetics and Genomic Sciences, Icahn School of Medicine at Mount Sinai,New York, NY USA
  5. Human Brain Collection Core, National Institute of Mental Health–Intramural Research Program,Bethesda, MD USA
  6. Mount Sinai Center for Transformative Disease Modeling, Icahn School of Medicine at Mount Sinai,New York, NY USA
  7. Icahn Genomics Institute, Icahn School of Medicine at Mount Sinai,New York, NY USA
  8. Department of Pharmacological Sciences, Icahn School of Medicine at Mount Sinai,New York, NY USA
  9. Mental Illness Research Education and Clinical Center (MIRECC), James J. Peters VA Medical Center,Bronx, NY USA
  10. Center for Precision Medicine and Translational Therapeutics, James J. Peters VA Medical Center,Bronx, NY USA
Institutions: Icahn School of Medicine at Mount Sinai (United States); National Institute of Mental Health (United States); James J. Peters VA Medical Center (United States)
Journal: Nature communications, volume 17, issue 1, article 4922
Dates: received 29 October 2024; accepted 16 March 2026; published online 6 April 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-71281-7 · PMID 41942462 · PMCID PMC13234368 · OpenAlex W7150785853
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), schizophrenia / psychosis (population), bipolar (population), cellular / molecular (subfield)
Methods: Smoothing, state filtering, decompositions, Statistics
Keywords: Computational neuroscience, Cellular neuroscience, Schizophrenia
MeSH: Dopaminergic Neurons*, Down-Regulation*, Mesencephalon*, Schizophrenia*, Bipolar Disorder, Female, Genetic Predisposition to Disease, Humans, Male, Nuclear Receptor Subfamily 4, Group A, Member 2, Transcription, Genetic (* major topic)
Topic: Schizophrenia research and treatment (Psychiatry and Mental health, Medicine), according to OpenAlex
Funding: R01DA047880 U01DA048279
Citations: not cited yet (Europe PMC); 66 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

Its files are read in the Code ↔ Paper reader above, with 3 matches between paragraphs and lines of code.

codeocean:2247653

License: none: the authors keep all their rights
State: cannot be verified, verified on 29 September 2026
Evidence: found in the paper
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 29 September 2026: cannot be verified
  • 29 September 2026: cannot be verified
At the source: codeocean.com/capsule/2247653/tree

GabrielHoffman/mvIC

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 15f0b63089645c840624e71da289aafb172bb2da, 30 August 2022
Languages: R (12), JavaScript (3)
Size: 80 files, 15 scripts
Software Heritage: not archived
Found in: the references
Holds: README, environment (DESCRIPTION), tests, documentation, 4 notebooks
Not found: license file, CITATION.cff, continuous integration
Tools: data.table (2 files), ggplot2 (2 files), lme4 (2 files), tidyverse (2 files), edgeR (1 file), limma (1 file), reshape2 (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
16 files
At the source: github.com/GabrielHoffman/mvIC

cran.r-project.org/package=remacor

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: the link answers
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
At the source: CRAN.R-project.org/package=remaCor

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41467-026-71281-7.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 3 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 15 scripts, each with its path and the digest of its content;
  • 3 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • no repository, dataset or request procedure was recognized in it

Read it in the paper: doi.org/10.1038/s41467-026-71281-7.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 15 authors, 3 keywords, 11 MeSH terms, 1 funder, 63 references.

Cite

This paper

Singh, S., Iskhakova, M., Lambert, T. Y., Valada, A., Shokrian, N., Evans, V., Bendl, J., Auluck, P. K., Marenco, S., Wang, M., Zhang, B., Hoffman, G. E., Girdhar, K., Roussos, P., & Akbarian, S. (2026). Downregulated transcription in chromosomal domains of midbrain dopamine neurons linked to schizophrenia. Nature communications, 17(1), 4922. https://doi.org/10.1038/s41467-026-71281-7

BibTeX

@article{singh2026downregulated,
author = {Singh, Swadha and Iskhakova, Marina and Lambert, Tova Y. and Valada, Aditi and Shokrian, Neda and Evans, Viviana and Bendl, Jaroslav and Auluck, Pavan K. and Marenco, Stefano and Wang, Minghui and Zhang, Bin and Hoffman, Gabriel E. and Girdhar, Kiran and Roussos, Panos and Akbarian, Schahram},
title = {{Downregulated transcription in chromosomal domains of midbrain dopamine neurons linked to schizophrenia}},
journal = {Nature communications},
year = {2026},
month = apr,
volume = {17},
number = {1},
pages = {4922},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-71281-7},
url = {https://doi.org/10.1038/s41467-026-71281-7},
pmid = {41942462},
pmcid = {PMC13234368}
}

RIS

TY - JOUR
AU - Singh, Swadha
AU - Iskhakova, Marina
AU - Lambert, Tova Y.
AU - Valada, Aditi
AU - Shokrian, Neda
AU - Evans, Viviana
AU - Bendl, Jaroslav
AU - Auluck, Pavan K.
AU - Marenco, Stefano
AU - Wang, Minghui
AU - Zhang, Bin
AU - Hoffman, Gabriel E.
AU - Girdhar, Kiran
AU - Roussos, Panos
AU - Akbarian, Schahram
TI - Downregulated transcription in chromosomal domains of midbrain dopamine neurons linked to schizophrenia
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/04/06
VL - 17
IS - 1
SP - 4922
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-71281-7
UR - https://doi.org/10.1038/s41467-026-71281-7
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-71281-7",
"type": "article-journal",
"title": "Downregulated transcription in chromosomal domains of midbrain dopamine neurons linked to schizophrenia",
"container-title": "Nature communications",
"author": [
{
"family": "Singh",
"given": "Swadha"
},
{
"family": "Iskhakova",
"given": "Marina"
},
{
"family": "Lambert",
"given": "Tova Y."
},
{
"family": "Valada",
"given": "Aditi"
},
{
"family": "Shokrian",
"given": "Neda"
},
{
"family": "Evans",
"given": "Viviana"
},
{
"family": "Bendl",
"given": "Jaroslav"
},
{
"family": "Auluck",
"given": "Pavan K."
},
{
"family": "Marenco",
"given": "Stefano"
},
{
"family": "Wang",
"given": "Minghui"
},
{
"family": "Zhang",
"given": "Bin"
},
{
"family": "Hoffman",
"given": "Gabriel E."
},
{
"family": "Girdhar",
"given": "Kiran"
},
{
"family": "Roussos",
"given": "Panos"
},
{
"family": "Akbarian",
"given": "Schahram"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "4922",
"DOI": "10.1038/s41467-026-71281-7",
"PMID": "41942462",
"PMCID": "PMC13234368",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-71281-7",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
6
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41588-026-02646-3 [code]
Co-expression-based models improve eQTL predictions for transcriptome-wide association studies and highlight new schizophrenia-associated genes.
Journal: Nature genetics
In common: limma, reshape2, data.table, 1 other tool, schizophrenia / psychosis, cellular / molecular, 8 references
[2] doi:10.1038/s41467-026-75722-1 [code]
Single-nucleus analysis of the adult human olfactory epithelium uncovers shared neurogenesis programs with the brain.
Journal: Nature communications
In common: reshape2, ggplot2, tidyverse, 2 references, 2 authors
[3] doi:10.1093/nar/gkag788 [code]
Neuronal activity-driven 3D chromatin dynamics in cortical pyramidal neurons depend on SATB2.
Journal: Nucleic acids research
In common: edgeR, limma, reshape2, 3 other tools, cellular / molecular, 4 references
[4] doi:10.1101/gr.281113.125 [code]
Single-nucleus multiomic profiling of the aging mouse substantia nigra reveals conserved gene alterations linked to Parkinson's disease.
Journal: Genome research
In common: edgeR, limma, reshape2, 3 other tools, cellular / molecular, 3 references
[5] doi:10.1038/s41467-026-71542-5 [code]
Astrocyte fatty acid metabolism as a driver of risk for major depressive disorder.
Journal: Nature communications
In common: edgeR, limma, lme4, 4 other tools, cellular / molecular, 2 references
[6] doi:10.1111/adb.70179 [code]
Transcriptional Response to Chronic Long-Access Fentanyl Self-Administration in Rat Habenula and Amygdala.
Journal: Addiction biology
In common: edgeR, limma, lme4, 4 other tools, cellular / molecular, 2 references
[7] doi:10.1038/s41467-026-74753-y [code]
A human-specific microRNA controls the timing of excitatory synaptogenesis.
Journal: Nature communications
In common: edgeR, limma, lme4, 4 other tools, cellular / molecular, 2 references
[8] doi:10.1038/s41467-026-76675-1 [code]
Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: edgeR, limma, reshape2, 3 other tools, 3 references
[9] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: edgeR, limma, lme4, 4 other tools, 2 references
[10] doi:10.1038/s41467-026-74541-8 [code]
SECmeres outperform extracellular vesicles as potential blood RNA biomarkers for Alzheimer's disease.
Journal: Nature communications
In common: 1 reference, 2 authors

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.