<html><head></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; ">Dear Lennart,<div><br></div><div>You did terrific job converting ProbABEL in quite professional-looking product!</div><div><div><br></div><div>I have just committed few changes related to documentation and adding '--mmscore' option for binary traits; also changed version to 0.2-0 as you suggested. Seems all worked all right, though I had to merge few files. </div><div><br></div><div>One questions: what is the right way to clean up local source before committing? It seems that a number of newly generated files (Makefile's,srd/.deps ...) hang around after running 'make clean' (naturally). I can just force them to be ignored by commit, but wonder if there may be a more proper way to do this?</div><div><br></div><div>best wishes,</div><div>Yurii<br><div>
<div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><span class="Apple-style-span" style="border-collapse: separate; color: rgb(0, 0, 0); font-family: Helvetica; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: 2; text-align: -webkit-auto; text-indent: 0px; text-transform: none; white-space: normal; widows: 2; word-spacing: 0px; -webkit-border-horizontal-spacing: 0px; -webkit-border-vertical-spacing: 0px; -webkit-text-decorations-in-effect: none; -webkit-text-size-adjust: auto; -webkit-text-stroke-width: 0px; font-size: medium; "><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div><div>-------------------------------------------------------</div><div>Yurii Aulchenko, PhD, Dr. Habil.</div><div>Independent researcher and consultant</div><div>yurii [dot] aulchenko [at] gmail [dot] com</div></div></div></span></div>
</div>
<br><div><div>On Aug 23, 2011, at 9:48 PM, <a href="mailto:noreply@r-forge.wu-wien.ac.at">noreply@r-forge.wu-wien.ac.at</a> wrote:</div><br class="Apple-interchange-newline"><blockquote type="cite"><div>Author: yurii<br>Date: 2011-08-23 21:48:36 +0200 (Tue, 23 Aug 2011)<br>New Revision: 762<br><br>Modified:<br> pkg/ProbABEL/doc/CHANGES.LOG<br> pkg/ProbABEL/doc/ProbABEL_manual.tex<br> pkg/ProbABEL/src/main.cpp<br> pkg/ProbABEL/src/mematri1.h<br> pkg/ProbABEL/src/reg1.h<br> pkg/ProbABEL/src/usage.h<br>Log:<br>update manual; changes.log; add '--mmscore' option for binary traits<br><br>Modified: pkg/ProbABEL/doc/CHANGES.LOG<br>===================================================================<br>--- pkg/ProbABEL/doc/CHANGES.LOG<span class="Apple-tab-span" style="white-space:pre"> </span>2011-08-22 08:32:43 UTC (rev 761)<br>+++ pkg/ProbABEL/doc/CHANGES.LOG<span class="Apple-tab-span" style="white-space:pre"> </span>2011-08-23 19:48:36 UTC (rev 762)<br>@@ -1,3 +1,54 @@<br>+***** v.0.2-0 (2011.08.23)<br>+<br>+* ProbABEL can now (experimentally!) analyze binary traits accounting for relationship <br>+structure (thanks to Yurii and Nicola Pirastu). This adds '--mmscore' <br>+option for logistic regression. Important: on the contrast to the <br>+'palinear' the matrix which should be supplied to palogist with <br>+--mmscore should contain an inverse of the correlation matrix (not <br>+the inverse of var-cov matrix, as is the case for 'palinear' with <br>+--mmscore). This matrix can be obtained through (GenABEL notation):<br>+<br>+h2.object$InvSigma * h.object2$h2an$estimate[length(h2$h2an$estimate)]<br>+<br>+where h2.object is the object returned by 'polygenic()'. Documentation <br>+explains this procedure in more details; a wrapper function to prepare <br>+and export correct objects for ProbABEL is planned. (thanks to Yurii) <br>+<br>+- in probabel.pl the location for the probabel config file is now set<br>+ to /etc/, the default location it will end up. Before it still<br>+ pointed to the location the probabel.pl file was in. However, a<br>+ better fix would be to somehow put in the actual value of the<br>+ --sysconfdir option to ./configure.<br>+- The PDF of the LaTeX documentation is now only generated if the<br>+ pdflatex binary can be found. So now we also build the documentation<br>+ using autotools.<br>+<br>+<br>+* ProbABEL uses autoconf and automake now (thanks to Lennart). Now the <br>+compilation process is like this: first, from the ProbABEL directory, <br>+generate documentation by<br>+<br>+cd doc<br>+make<br>+make clean_doc<br>+cd ..<br>+<br>+Then, run<br>+<br>+./configure<br>+make<br>+make install <br>+<br>+To see options, run<br>+<br>+./configure --help<br>+<br>+In the ProbABEL autotools integration branch:<br>+- Removed the old Makefile<br>+- Added configure.ac that servers as input for auto(re)conf<br>+- Added the Makefile.am files needed by automake to generate the final Makefiles.<br>+<br>+<br> ***** v.0.1-9e (2011.05.15)<br> Previous version wouldn't compile because of missing frerror. file.<br> Fixes bug #1339 (the bug fix was done a month ago, but no package had<br><br>Modified: pkg/ProbABEL/doc/ProbABEL_manual.tex<br>===================================================================<br>--- pkg/ProbABEL/doc/ProbABEL_manual.tex<span class="Apple-tab-span" style="white-space:pre"> </span>2011-08-22 08:32:43 UTC (rev 761)<br>+++ pkg/ProbABEL/doc/ProbABEL_manual.tex<span class="Apple-tab-span" style="white-space:pre"> </span>2011-08-23 19:48:36 UTC (rev 762)<br>@@ -1,6 +1,9 @@<br> \title{ProbABEL manual}<br>-\author{Yurii Aulchenko, Maksim Struchalin, Lennart Karssen \\<br>-<span class="Apple-tab-span" style="white-space:pre"> </span>Erasmus MC Rotterdam<br>+\author{<br>+Maksim Struchalin$^{1}$, Lennart Karssen$^{1}$, Yurii Aulchenko$^{1,2}$ \\ <br>+\\<br>+$^{1}$Erasmus MC Rotterdam \\ <br>+$^{2}$Institute of Cytology and Genetics SD RAS <br> }<br> \date{\today}<br><br>@@ -28,6 +31,10 @@<br><br> \makeindex<br><br>+\newcommand{\PA}{\texttt{ProbABEL-package}}<br>+\newcommand{\GA}{\texttt{GenABEL-package}}<br>+\newcommand{\DA}{\texttt{DatABEL-package}}<br>+<br> \begin{document}<br> \maketitle<br> \tableofcontents<br>@@ -76,16 +83,16 @@<br> performing association analysis by means of regression of the<br> outcome of interest onto estimated genotypic probabilities.<br><br>-The \texttt{ProbABEL} package was designed to perform such regression<br>+The \PA{} package was designed to perform such regression<br> in a fast, memory-efficient and consequently genome-wide feasible manner.<br>-Currently, \texttt{ProbABEL} implements linear, logistic regression,<br>+Currently, \PA{} implements linear, logistic regression,<br> and Cox proportional hazards models. The corresponding analysis<br> programs are called \texttt{palinear}, \texttt{palogist},<br> and \texttt{pacoxph}.<br><br><br> \section{Input files}<br>-\texttt{ProbABEL} takes three files as input: a file containing SNP<br>+\PA{} takes three files as input: a file containing SNP<br> information (e.g.~the MLINFO file of MACH), a file with genome- or<br> chromosome-wide predictor information (e.g.~the MLDOSE or MLPROB file of MACH),<br> and a file containing the phenotype of interest and covariates.<br>@@ -93,6 +100,13 @@<br> Optionally, the map information can be supplied (e.g.~the "legend"<br> files of HapMap).<br><br>+The dose/probability file may be supplied in filevector format<br>+in which case \PA{} will operate much faster, and<br>+in low-RAM mode (approx. $\approx$ 128 MB). See the R libraries \GA{} and<br>+\DA{} on how to convert MACH and IMPUTE files to<br>+filevector format (functions: \texttt{mach2databel()} and<br>+\texttt{impute2databel()}, respectively).<br>+<br> \subsection{SNP information file}<br> \label{ssec:infoin}<br> In the simplest scenario, the SNP information file is an MLINFO<br>@@ -103,7 +117,7 @@<br> ``Rsq'', the proportion of variance decrease after imputations).<br><br> Actually,<br>-for \texttt{ProbABEL}, it does not matter what is written in this file --<br>+for \PA{}, it does not matter what is written in this file --<br> this information is just brought forward to the output. However,<br> \textbf{it is critical} that the number of columns is seven and the number<br> of lines in the file is equal to the number of SNPs in the<br>@@ -239,7 +253,7 @@<br> the program).<br><br> There are in total 11 command line options you can specify to the<br>-\texttt{ProbABEL} analysis functions \texttt{palinear} or<br>+\PA{} analysis functions \texttt{palinear} or<br> \texttt{palogist}. If you run either program without any argument, you<br> will get a short explanation to command line options:<br> \begin{verbatim}<br>@@ -350,25 +364,25 @@<br> first column contains id names, the second the trait. The others are<br> covariates.<br><br>-An example of how a polygenic object estimated by GenABEL can be used<br>+An example of how a polygenic object estimated by \GA{} can be used<br> with ProbABEL is provided in \texttt{../example/mmscore.R}<br><br> Though technically \texttt{--mmscore} allows for inclusion of multiple<br> covariates, these should be kept to minimum as this is a score test. We<br> suggest that any covariates explaining an essential proportion of<br>-variance should be fit as part of \texttt{GenABEL}'s<br>+variance should be fit as part of \GA{}'s<br> \texttt{polygenic} procedure.<br><br> %Option \texttt{--interaction\_only} is like \texttt{--interaction} but does not<br> %include in the model the main effect of the covariate, which is acting in<br> %interaction with SNP. This option is useful when running \texttt{--mmscore},<br> %in whch case the main effect should normally estimated in the polygenic<br>-%model and only the interaction term in the \texttt{ProbABEL} analysis.<br>+%model and only the interaction term in the \PA{} analysis.<br><br> \subsection{Running multiple analyses at once: \texttt{probabel.pl}}<br><br> The Perl script \texttt{bin/probabel.pl\_example} represents a handy<br>-wraper for \texttt{ProbABEL} functions. To start using it the<br>+wraper for \PA{} functions. To start using it the<br> configuration file \texttt{bin/probabel\_config.cfg\_example} needs to<br> be edited. The configuration file consists of five columns. Each<br> column except the first is a pattern for files produced by<br>@@ -383,7 +397,7 @@<br> for the ``mlinfo'' files. Probably you also have to change the variable<br> \texttt{\$config} in the script to point to the full path of the<br> configuration file and the variable \texttt{@anprog} to point full<br>-path to the \texttt{ProbABEL} scripts.<br>+path to the \PA{} scripts.<br><br><br> \section{Output file format}<br>@@ -412,7 +426,7 @@<br> At this point, unless you have complete measurements on all<br> subjects, you should feel alarmed if the number here is exactly the<br> number of people in the file -- this probably indicates you did not code<br>-missing values according to \texttt{ProbABEL} format ('NA', 'NaN', or 'N').<br>+missing values according to \PA{} format ('NA', 'NaN', or 'N').<br><br> The next column, nine (``Mean\_predictor\_allele''), gives the estimated<br> frequency of the predictor allele in subjects with complete phenotypic data.<br>@@ -433,14 +447,24 @@<br><br> \section{Memory use and performance}<br> Maximum likelihood regression is implemented in<br>-\texttt{ProbABEL}. With 6,000 people and 2.5 millions SNPs, a<br>+\PA{}. With 6,000 people and 2.5 millions SNPs, a<br> genome-wide scan is completed in less that an hour for a linear model<br> with 1-2 covariates and overnight for logistic regression or the Cox<br>-proportional hazards model.<br>+proportional hazards model (figures for a PC bought back in 2007).<br><br>-Memory is an issue with \texttt{ProbABEL} -- large chromosomes,<br>+Memory may be an issue with \PA{} if you use <br>+MACH text dose/probability files, e.g. for large chromosomes,<br> such as chromosome one consumed up to 5 GB of RAM with 6,000 people.<br><br>+We suggest that dose/probability file is to be supplied in filevector format<br>+in which case \PA{} will operate about 2-3 times faster, and<br>+in low-RAM mode (approx.~128 MB). See the R libraries \GA{} and<br>+\DA{} on how to convert MACH and IMPUTE files to<br>+filevector format (functions: \texttt{mach2databel()} and<br>+\texttt{impute2databel()}, respectively).<br>+<br>+When '--mmscore' option is used, the analysis may take quite some time. <br>+<br> \section{Methodology}<br> \label{sec:methodology}<br> \subsection{Analysis of population-based data}<br>@@ -598,7 +622,7 @@<br> \index{Cox proportional hazards model}<br> \index{proportional hazards model}<br> \index{regression!Cox proportional hazards}<br>-in \texttt{ProbABEL} is entirely based on the code of \texttt{R}<br>+in \PA{} is entirely based on the code of \texttt{R}<br> library \texttt{survival} developed by Thomas Lumley<br> (function \texttt{coxfit2}), and is therefore not described here.<br><br>@@ -742,14 +766,14 @@<br><br> \section{How to cite}<br><br>-If you used \texttt{ProbABEL} for<br>-your analysis please give a link to the \texttt{ABEL} home<br>+If you used \PA{} for<br>+your analysis please give a link to the \texttt{GenABEL project} home<br> page<br><br> \begin{quote}<br> \url{<a href="http://www.genabel.org/">http://www.genabel.org/</a>}<br> \end{quote}<br>-and cite the \texttt{ProbABEL} paper to give us some credit:<br>+and cite the \PA{} paper to give us some credit:<br> \begin{quote}<br> Aulchenko YS, Struchalin MV, van Duijn CM.\\<br> \emph{ProbABEL package for genome-wide association analysis of imputed data.}\\<br>@@ -757,15 +781,15 @@<br> \end{quote}<br> A proper reference may look like<br> \begin{quote}<br>-For the analysis of imputed data, we used the \texttt{ProbABEL} package<br>-from the GenABEL suite of programs (Aulchenko \emph{et al.}, 2010).<br>+For the analysis of imputed data, we used the \PA{} <br>+from the \texttt{GenABEL} suite of programs (Aulchenko \emph{et al.}, 2010).<br> \end{quote}<br><br> If you have used the Cox proportional hazard model, please mention the<br> \texttt{R} package \texttt{survival} by Thomas Lumley. Additionally<br> to the above citation, please tell that<br> \begin{quote}<br>-The Cox proportional hazards model implemented in \texttt{ProbABEL}<br>+The Cox proportional hazards model implemented in \PA{}<br> makes use of the source code of the \texttt{R} package ''\texttt{survival}''<br> as implemented by T. Lumley.<br> \end{quote}<br><br>Modified: pkg/ProbABEL/src/main.cpp<br>===================================================================<br>--- pkg/ProbABEL/src/main.cpp<span class="Apple-tab-span" style="white-space:pre"> </span>2011-08-22 08:32:43 UTC (rev 761)<br>+++ pkg/ProbABEL/src/main.cpp<span class="Apple-tab-span" style="white-space:pre"> </span>2011-08-23 19:48:36 UTC (rev 762)<br>@@ -303,9 +303,12 @@<br><br> <span class="Apple-tab-span" style="white-space:pre"> </span>mematrix<double> invvarmatrix;<br><br>+/*<br>+ * now should be possible... delete this part later when everything works<br> #if LOGISTIC<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><br> <span class="Apple-tab-span" style="white-space:pre"> </span>if(inverse_filename != NULL) {std::cerr<<"ERROR: mmscore is forbidden for logistic regression\n";exit(1);}<br> #endif<br>+*/<br><br> #if COXPH<br> <span class="Apple-tab-span" style="white-space:pre"> </span>if(inverse_filename != NULL) {std::cerr<<"ERROR: mmscore is forbidden for cox regression\n";exit(1);}<br><br>Modified: pkg/ProbABEL/src/mematri1.h<br>===================================================================<br>--- pkg/ProbABEL/src/mematri1.h<span class="Apple-tab-span" style="white-space:pre"> </span>2011-08-22 08:32:43 UTC (rev 761)<br>+++ pkg/ProbABEL/src/mematri1.h<span class="Apple-tab-span" style="white-space:pre"> </span>2011-08-23 19:48:36 UTC (rev 762)<br>@@ -368,7 +368,7 @@<br> {<br> <span class="Apple-tab-span" style="white-space:pre"> </span>if (M.ncol != M.nrow) <br> <span class="Apple-tab-span" style="white-space:pre"> </span>{<br>-<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stderr,"invert: only suare matrices possible\n");<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stderr,"invert: only square matrices possible\n");<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>exit(1);<br> <span class="Apple-tab-span" style="white-space:pre"> </span>}<br> <span class="Apple-tab-span" style="white-space:pre"> </span>if (M.ncol == 1) <br>@@ -376,12 +376,19 @@<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>mematrix<DT> temp(1,1);<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>temp[0] = 1./M[0];<br> <span class="Apple-tab-span" style="white-space:pre"> </span>}<br>+<span class="Apple-tab-span" style="white-space:pre"> </span>/*<br> <span class="Apple-tab-span" style="white-space:pre"> </span>for (int i=0;i<M.ncol;i++) <br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if (M.data[i*M.ncol+i]==0) <br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>{<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stderr,"invert: zero elements in diagonal\n");<br>-<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>exit(1);<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>mematrix<DT> temp = M;<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>for (int i = 0; i < M.ncol; i++)<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>for (int j = 0; j < M.ncol; j++)<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>temp.put(NAN,i,j);<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>return temp;<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//exit(1);<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>}<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>*/<br> <span class="Apple-tab-span" style="white-space:pre"> </span>int actualsize = M.ncol;<br> <span class="Apple-tab-span" style="white-space:pre"> </span>int maxsize = M.ncol;<br> <span class="Apple-tab-span" style="white-space:pre"> </span>mematrix<DT> temp = M;<br><br>Modified: pkg/ProbABEL/src/reg1.h<br>===================================================================<br>--- pkg/ProbABEL/src/reg1.h<span class="Apple-tab-span" style="white-space:pre"> </span>2011-08-22 08:32:43 UTC (rev 761)<br>+++ pkg/ProbABEL/src/reg1.h<span class="Apple-tab-span" style="white-space:pre"> </span>2011-08-23 19:48:36 UTC (rev 762)<br>@@ -626,7 +626,7 @@<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//Oct 26, 2009<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>residuals.reinit((rdata.X).nrow,1);<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>sigma2=-1.;<br>-<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>loglik=-9.999e+32;<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>loglik=-9.999e+32; // should actually be MAX of the corresponding type<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>niter=-1;<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>chi2_score=-1.;<br> <span class="Apple-tab-span" style="white-space:pre"> </span>}<br>@@ -635,10 +635,36 @@<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>delete beta;<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>delete sebeta;<br> <span class="Apple-tab-span" style="white-space:pre"> </span>}<br>-<span class="Apple-tab-span" style="white-space:pre"> </span>void estimate(regdata &rdatain, int verbose, int maxiter, double eps, double tol_chol, int model, int interaction, int ngpreds, mematrix<double> & invvarmatrix, int robust, int nullmodel=0)<br>+<br>+<span class="Apple-tab-span" style="white-space:pre"> </span>void estimate(regdata &rdatain, int verbose, int maxiter, double eps, double tol_chol,<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>int model, int interaction, int ngpreds, mematrix<double> & invvarmatrixin, int robust,<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>int nullmodel=0)<br> <span class="Apple-tab-span" style="white-space:pre"> </span>{<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>// on the contrast to the 'linear' 'invvarmatrix' contains<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>// inverse of correlation matrix (not the inverse of var-cov matrix)<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>// h2.object$InvSigma * h.object2$h2an$estimate[length(h2$h2an$estimate)]<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>// the inverse of var-cov matrix scaled by total variance<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>regdata rdata = rdatain.get_unmasked_data();<br><br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>// a lot of code duplicated between linear and logistic...<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>// e.g. a piece below...<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>mematrix<double> invvarmatrix;<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if(invvarmatrixin.nrow!=0 && invvarmatrixin.ncol!=0)<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>{<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>invvarmatrix.reinit(rdata.nids,rdata.nids);<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>int i1=0, j1=0;<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>for (int i=0;i<rdata.nids;i++)<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if (rdatain.masked_data[i]==0) {<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>j1 = 0;<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>for (int j=0;j<rdata.nids;j++) if (rdatain.masked_data[j]==0) {<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>invvarmatrix.put(invvarmatrixin.get(i,j),i1,j1);<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>j1++;<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>}<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>i1++;<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>}<br>+<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>}<br>+<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>mematrix<double> X = apply_model(rdata.X,model, interaction, ngpreds, false, nullmodel);<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>int length_beta = X.ncol;<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>beta.reinit(length_beta,1);<br>@@ -668,10 +694,20 @@<br><br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>mematrix<double> tX = transpose(X);<br><br>-<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if(invvarmatrix.nrow!=0 && invvarmatrix.ncol!=0) tX = X*invvarmatrix;<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if(invvarmatrix.nrow!=0 && invvarmatrix.ncol!=0)<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>{<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>tX = tX*invvarmatrix;<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>}<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>/*<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stdout,"\n");<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stdout,"X %f %f %f\n",X.get(0,0),X.get(0,1),X.get(0,2));<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if (X.ncol==4) fprintf(stdout,"X[4] %f\n",X.get(0,3));<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stdout,"Inv %f %f %f\n",invvarmatrix.get(0,0),invvarmatrix.get(0,1),invvarmatrix.get(0,2));<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if (X.ncol==4) fprintf(stdout,"X[4] %f\n",invvarmatrix.get(0,3));<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stdout,"tXInv %f %f %f\n",tX.get(0,0),tX.get(1,0),tX.get(2,0));<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if (X.ncol==4) fprintf(stdout,"X[4] %f\n",tX.get(3,0));<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>*/<br><br>-<br>-<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>double N;<br><br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>niter=0;<br>@@ -700,21 +736,29 @@<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>N = tXWX.get(0,0);<br><br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if (verbose) {printf("tXWX:\n");tXWX.print();}<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//printf("tXWX:\n");tXWX.print();<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>// use cholesky to invert<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>tXWX_i = tXWX;<br>-<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>cholesky2_mm(tXWX_i,tol_chol);<br>-<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if (verbose) {printf("chole tXWX:\n");tXWX_i.print();}<br>-<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>chinv2_mm(tXWX_i);<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//cholesky2_mm(tXWX_i,tol_chol);<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//if (verbose) {printf("chole tXWX:\n");tXWX_i.print();}<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//printf("chole tXWX:\n");tXWX_i.print();<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//chinv2_mm(tXWX_i);<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>// was before<br>-<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>// tXWX_i=invert(tXWX);<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>tXWX_i=invert(tXWX);<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if (verbose) {printf("tXWX-1:\n");tXWX_i.print();}<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//fprintf(stdout,"*** tXWX_i\n");tXWX_i.print();<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>mematrix<double> tmp1 = productMatrDiag(tX,W);<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>mematrix<double> tXWz = tmp1*z;<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if (verbose) {printf("tXWz:\n");tXWz.print();}<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>beta = tXWX_i*tXWz;<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//fprintf(stdout,"*** res: %f %f %f\n",residuals[0],residuals[1],residuals[2]);<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//mematrix<double> txres = tx * residuals;<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//fprintf(stdout,"*** txres\n");txres.print();<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//beta = txwx_i* txres;<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if (verbose) {printf("beta:\n");beta.print();}<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//printf("beta:\n");beta.print();<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>// compute likelihood<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>prevlik = loglik;<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>loglik=0.;<br>@@ -780,6 +824,9 @@<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>}<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>}<br> <span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>if (verbose) {printf("sebeta (%d):\n",sebeta.nrow);sebeta.print();}<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//printf("sebeta (%d):\n",beta.nrow);beta.print();<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//printf("sebeta (%d):\n",sebeta.nrow);sebeta.print();<br>+<span class="Apple-tab-span" style="white-space:pre"> </span><span class="Apple-tab-span" style="white-space:pre"> </span>//exit(1);<br> <span class="Apple-tab-span" style="white-space:pre"> </span>}<br> <span class="Apple-tab-span" style="white-space:pre"> </span>// just a stupid copy from linear_reg<br> <span class="Apple-tab-span" style="white-space:pre"> </span>void score(mematrix<double> &resid,regdata &rdata,int verbose, double tol_chol, int model, int interaction, int ngpreds, mematrix<double> & invvarmatrix, int nullmodel=0)<br><br>Modified: pkg/ProbABEL/src/usage.h<br>===================================================================<br>--- pkg/ProbABEL/src/usage.h<span class="Apple-tab-span" style="white-space:pre"> </span>2011-08-22 08:32:43 UTC (rev 761)<br>+++ pkg/ProbABEL/src/usage.h<span class="Apple-tab-span" style="white-space:pre"> </span>2011-08-23 19:48:36 UTC (rev 762)<br>@@ -22,7 +22,7 @@<br> <span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stdout,"\t --allcov : report estimates for all covariates (large outputs!)\n");<br> <span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stdout,"\t --interaction: Which covariate to use for interaction with SNP analysys (default is no ineraction, 0)\n");<br> <span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stdout,"\t --interaction_only: like previos but without covariate acting in interaction with SNP (default is no ineraction, 0)\n");<br>-<span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stdout,"\t --mmscore : score test in samples of related individuals. File with inverse of variance-covariance matrix as input parameter\n");<br>+<span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stdout,"\t --mmscore : score test in samples of related individuals. File with inverse of variance-covariance matrix (for palinear) or inverse correlation (for palogist) as input parameter\n");<br> <span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stdout,"\t --robust : report robust (aka sandwich, aka Hubert-White) standard errors\n");<br> <span class="Apple-tab-span" style="white-space:pre"> </span>fprintf(stdout,"\t --help : print help\n");<br> <span class="Apple-tab-span" style="white-space:pre"> </span>exit(exit_code);<br><br>_______________________________________________<br>Genabel-commits mailing list<br><a href="mailto:Genabel-commits@lists.r-forge.r-project.org">Genabel-commits@lists.r-forge.r-project.org</a><br>https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/genabel-commits<br></div></blockquote></div><br></div></div></body></html>