[R-gregmisc-users] SF.net SVN: r-gregmisc: [1090] trunk/PathwayModeling/thesispaper
Brought to you by:
warnes
From: <wa...@us...> - 2007-04-17 15:00:17
|
Revision: 1090 http://svn.sourceforge.net/r-gregmisc/?rev=1090&view=rev Author: warnes Date: 2007-04-17 08:00:06 -0700 (Tue, 17 Apr 2007) Log Message: ----------- Finalize revisions for IET Systems Biology Modified Paths: -------------- trunk/PathwayModeling/thesispaper/R/plotConverged_nocolor.R trunk/PathwayModeling/thesispaper/R/plotDensity_nocolor.R trunk/PathwayModeling/thesispaper/bibliography.Snw trunk/PathwayModeling/thesispaper/discussion.Snw trunk/PathwayModeling/thesispaper/figure1.Snw trunk/PathwayModeling/thesispaper/figure2.Snw trunk/PathwayModeling/thesispaper/figure3.Snw trunk/PathwayModeling/thesispaper/figure4.Snw trunk/PathwayModeling/thesispaper/figure5.Snw trunk/PathwayModeling/thesispaper/figure6.Snw trunk/PathwayModeling/thesispaper/figure7.Snw trunk/PathwayModeling/thesispaper/figure8.Snw trunk/PathwayModeling/thesispaper/methods.Snw trunk/PathwayModeling/thesispaper/paper.Snw trunk/PathwayModeling/thesispaper/paper.pdf trunk/PathwayModeling/thesispaper/priorPlot.Snw trunk/PathwayModeling/thesispaper/results.Snw Modified: trunk/PathwayModeling/thesispaper/R/plotConverged_nocolor.R =================================================================== --- trunk/PathwayModeling/thesispaper/R/plotConverged_nocolor.R 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/R/plotConverged_nocolor.R 2007-04-17 15:00:06 UTC (rev 1090) @@ -7,19 +7,17 @@ maxY <- 0.0012 # for (i in 1:8) { - for (i in c(1,3,8)) { + for (i in c(1,4)) { p <- paste("d",i,sep="") plotDensity_nocolor(output12[,i],p,5000,Ylim=maxY,lineType=2) - par(new=T) - plotDensity_nocolor(output16[,i],p,5000,Ylim=maxY,lineType=3) - par(new=T) - plotDensity_nocolor(output25[,i],p,5000,Ylim=maxY,lineType=4) + plotDensity_nocolor(output16[,i],p,5000,Ylim=maxY,lineType=3, add=TRUE) + plotDensity_nocolor(output25[,i],p,5000,Ylim=maxY,lineType=4, add=TRUE) } -# plotDensity_nocolor(output12[,9],'d9',5,Ylim=5.7,lineType=2) -# par(new=T) -# plotDensity_nocolor(output16[,9],'d9',5,Ylim=5.7,lineType=3) -# par(new=T) -# plotDensity_nocolor(output25[,9],'d9',5,Ylim=5.7,lineType=4) + # prior mode is 5, but fits are all much less than 5 so scale down so range + # is [0,5] + plotDensity_nocolor(output12[,9],'d9',5,Ylim=5.7,lineType=2, Xlim=5) + plotDensity_nocolor(output16[,9],'d9',5,Ylim=5.7,lineType=3, Xlim=5, add=TRUE) + plotDensity_nocolor(output25[,9],'d9',5,Ylim=5.7,lineType=4, Xlim=5, add=TRUE) par(oldpar) Modified: trunk/PathwayModeling/thesispaper/R/plotDensity_nocolor.R =================================================================== --- trunk/PathwayModeling/thesispaper/R/plotDensity_nocolor.R 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/R/plotDensity_nocolor.R 2007-04-17 15:00:06 UTC (rev 1090) @@ -2,7 +2,9 @@ parameter="di", prior, Ylim=0, - lineType=1 + Xlim=5*prior, + lineType=1, + add=FALSE ) { @@ -10,13 +12,18 @@ deltaX<-(max(dens$x)-min(dens$x))/length(dens$x) sumY<-sum(dens$y) maxY<-max(dens$y, Ylim) - X<-seq(1,5*prior, 0.05*prior) - - plot(X, 3*dchisq(3*X/prior, 5)/prior, type='l', - xlim=c(0, 5*prior), ylim=c(0, maxY), - ylab="density", xlab=parameter, - lty=1, lwd=2) + X <- seq(0,Xlim, 0.05*prior) + if(!add) + plot(X, + 3*dchisq(3*X/prior, 5)/prior, + type='l', + xlim=c(0,Xlim), + ylim=c(0, maxY), + ylab="density", + xlab=parameter, + lty=1, lwd=2) + lines(dens$x, dens$y, lty=lineType, lwd=2) } Modified: trunk/PathwayModeling/thesispaper/bibliography.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/bibliography.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/bibliography.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -58,7 +58,4 @@ Byrd, R.H., Lu, P., Nocedal, J., and Zhu, C., 1995, A limited memory algorithm for bound constrained optimization, \emph{SIAM J. Sci. Comput.}, \textbf{16}, 1190--1208. - \bibitem{Geyer}%16 - Geyer, CJ, (no date), \emph{Burn-In is Unnecessary}, web site, - http://www.stat.umn.edu/~charlie/mcmc/burn.html \end{thebibliography} Modified: trunk/PathwayModeling/thesispaper/discussion.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/discussion.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/discussion.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -25,7 +25,7 @@ of the system and use the experimental data to refine the estimates (Figure~\ref{converged}). Thus the model fitting procedure described here lends itself to iterative experimentation where the experimental -results, even if they consit of a single datum, can be used to update +results, even if they consist of a single datum, can be used to update the prior for the next experiment. The models used here have the form of the Hill function, @@ -53,7 +53,7 @@ times varied considerably depending on the size of the model. The all-component Metropolis and NKC algorithms converged quickly while the 1-component algorithm was quite slow to converge. This is -undoubtably due to the high correlations between some parameters (see +undoubtedly due to the high correlations between some parameters (see Figure~\ref{scatter}). The 1-component algorithm has trouble sampling from this density because it is restricted to moves parallel to the coordinate axes. This is less of a problem for the all-components and Modified: trunk/PathwayModeling/thesispaper/figure1.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/figure1.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/figure1.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -1,6 +1,6 @@ -\begin{figure} +\begin{figure*} \centering \includegraphics[scale=0.9]{figures/glycolysis} \caption{Production of methylglyoxal in hyperglycemia} \label{glycolysis} -\end{figure} +\end{figure*} Modified: trunk/PathwayModeling/thesispaper/figure2.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/figure2.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/figure2.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -21,12 +21,18 @@ N <- dev.off() detach(rawdata.trimmed) @ - \begin{figure} + +\begin{figure*} + \centering % \includegraphics[scale=0.5]{figures/tempDir/pulse} \includegraphics[scale=0.5]{figures/tempDir/pulse_nocolor} - \caption[Reactant concentrations following a pulse of R1]{Reactant concentrations following a pulse of R1 at - $time=20$ for the sequence of reactions $R1\rightarrow - R2\rightarrow R3\rightarrow R4\rightarrow R5\rightarrow sink$.} + \caption[Reactant concentrations following a pulse of R1]{ + Reactant concentrations following a pulse of R1 at + $time=20$ for the sequence of reactions $R1\rightarrow + R2\rightarrow R3\rightarrow R4\rightarrow R5\rightarrow sink$. + } \label{pulse} - \end{figure} + +\end{figure*} + Modified: trunk/PathwayModeling/thesispaper/figure3.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/figure3.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/figure3.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -1,4 +1,5 @@ - \begin{figure} +\begin{figure*} + \centering \subfloat[1-component]{\includegraphics[scale=0.45]{figures/1comp} \label{1comp}} @@ -8,12 +9,14 @@ \qquad \subfloat[NKC]{\includegraphics[scale=0.45]{figures/nkc}\label{nkc}} \qquad - \caption[Movement of Markov chains with the - Metropolis algorithms]{Movement of Markov chains with the component-wise and - all-components Metropolis algorithms. Movement in a - 2-dimensional space is illustrated, so each point has 2 - components. The dotted lines are contours of equal probability - density for the proposal distributions and the dashed lines are - probabilty contours of the target distribution.} + \caption[Movement of Markov chains with the Metropolis algorithms]{ + Movement of Markov chains with the component-wise and + all-components Metropolis algorithms. Movement in a + 2-dimensional space is illustrated, so each point has 2 + components. The dotted lines are contours of equal probability + density for the proposal distributions and the dashed lines are + probability contours of the target distribution. + } \label{1comp_vs_all} - \end{figure} + +\end{figure*} Modified: trunk/PathwayModeling/thesispaper/figure4.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/figure4.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/figure4.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -36,17 +36,24 @@ par(mfrow=c(1,1)) detach(input16) @ -\begin{figure} + +\begin{figure*} + \centering % \includegraphics[scale=0.8]{figures/tempDir/densities} \includegraphics[scale=0.8]{figures/tempDir/densities_nocolor} - \caption[Histograms of the marginal probability distributions]{Histograms of the marginal probability distributions for - the 5-reaction model generated with the all-components Metropolis - algorithm and the 16-point dataset. The curves are normal densities with means equal to the - medians of the distributions and variances equal to the variances of - the distributions. Broken + \caption[Histograms of the marginal probability distributions]{ + + Histograms of the marginal probability distributions for + the 5-reaction model generated with the all-components Metropolis + algorithm and the 16-point dataset. The curves are normal + densities with means equal to the medians of the distributions and + variances equal to the variances of the distributions. Broken % Red -vertical lines indicate the parameters values that - minimize the mean squared residuals.} + vertical lines indicate the parameters values that minimize the + mean squared residuals. + } \label{densities} -\end{figure} + +\end{figure*} + Modified: trunk/PathwayModeling/thesispaper/figure5.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/figure5.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/figure5.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -1,4 +1,4 @@ -<<fig5,echo=F,eval=F>>=4 +<<fig5,echo=F,eval=T>>=4 # source("R/plotDensity.R") # source("R/plotConverged.R") source("R/plotDensity_nocolor.R") @@ -9,7 +9,7 @@ # pdf(file="figures/tempDir/converged.pdf",width=6,height=6) # postscript(file="figures/tempDir/converged.eps",width=6,height=6) -pdf(file="figures/tempDir/converged_nocolor.pdf",width=2*6,height=2*2) +pdf(file="figures/tempDir/converged_nocolor.pdf",width=6,height=2) # plotConverged() plotConverged_nocolor() N <- dev.off() @@ -22,11 +22,11 @@ par(mfrow=c(1,1)) @ -\begin{figure} +\begin{figure*} \centering \includegraphics[width=1.0\textwidth]{figures/tempDir/converged_nocolor} \caption{% - Posterior distributions for selected parameters fitte1d using the % + Posterior distributions for selected parameters fitted using the % all-components algorithm. For each parameter, the model was fitted % with 4 different time point granularities across the fixed time interval: % --------- prior distribution (0 points), % @@ -35,4 +35,4 @@ -- $\cdot$ -- $\cdot$ -- 25 points % } \label{converged} -\end{figure} +\end{figure*} Modified: trunk/PathwayModeling/thesispaper/figure6.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/figure6.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/figure6.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -20,14 +20,14 @@ N <- dev.off() par(pch=1) @ -\begin{figure} +\begin{figure*} \centering % \includegraphics[width=1.0\textwidth]{figures/tempDir/scatterPlot} \includegraphics[width=1.0\textwidth]{figures/tempDir/scatterPlot_nocolor} \caption[Bivariate scatter plots of the parameter distributions]{Bivariate scatter plots of the parameter distributions for the 5-reaction model found with the all-components Metropolis algorithm (upper triangle); correlation coefficients (lower - triangle). The veritcal and horizontal lines indicate the maximum + triangle). The vertical and horizontal lines indicate the maximum likelihood estimates of the parameters.} \label{scatter} -\end{figure} +\end{figure*} Modified: trunk/PathwayModeling/thesispaper/figure7.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/figure7.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/figure7.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -23,7 +23,7 @@ N <- dev.off() detach(input16) @ -\begin{figure} +\begin{figure*} \centering % \subfloat[v2]{\includegraphics[scale=0.30]{figures/tempDir/V1fitted}} % \subfloat[v3]{\includegraphics[scale=0.30]{figures/tempDir/V2fitted}}\\ @@ -42,4 +42,4 @@ function. } \label{fits} -\end{figure} +\end{figure*} Modified: trunk/PathwayModeling/thesispaper/figure8.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/figure8.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/figure8.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -30,7 +30,7 @@ N <- dev.off() detach(input16) @ -\begin{figure} +\begin{figure*} \centering % \includegraphics[scale=0.5]{figures/tempDir/MSR16} \includegraphics[scale=0.5]{figures/tempDir/MSR16_nocolor} @@ -38,4 +38,4 @@ evaluations]{Mean squared residuals vs. number of likelihood evaluations for the 5-reaction model} \label{SSQvsIter} -\end{figure} +\end{figure*} Modified: trunk/PathwayModeling/thesispaper/methods.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/methods.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/methods.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -6,9 +6,11 @@ \subsection{The Pathway} - While we examined several pathway structures, in this paper we - present the results for a 5-reaction sequence that is sampled at 12, - 16, or 25 time points with 3 replicates at each time point. The + We examined pathways containing 3, 5, 7, and 9 reactions with both + linear and branched structures using two sets of parameter values. + As essentially identical modeling performance was obtained for each, + we present only the results for a 5-reaction sequence, sampled at + 12, 16, or 25 time points with 3 replicates at each time point. The units for time and concentration are arbitrary; for the purpose of illustration we will use minutes as the unit of time and micromoles/liter ($\mu$M) as the unit of concentration. @@ -112,7 +114,7 @@ half-maximal velocity} \end{eqnarray*} This form of the equation is very useful because $v$ and $S$ are - usually measureable and $V_{max}$ and $K_m$ can be obtained by + usually measurable and $V_{max}$ and $K_m$ can be obtained by fitting equation 3 to the data. In contrast, modeling $v$ as a function of the individual rate constants is less useful because that requires the measurement of the concentration of the @@ -236,38 +238,39 @@ The Normal Kernel Coupler (NKC) algorithm operates on multiple chains, so that there are several ``current'' states. The NKC uses - a normal kernel density estimate as the proposal distribution for - choosing candidate states. The estimate is generated using the - entire set of current points. Since the algorithm uses the entire - set of current points it can move over areas of low probability in - the parameter space, especially if the user takes care to seed - each mode with a few points in the starting set. + a normal kernel density estimate based on the entire set of + current states as the proposal distribution for choosing candidate + states. Since the algorithm uses the entire set of current points + it can move over areas of low probability in the parameter space, + especially if the user takes care to seed the starting set with a + few points in each mode. \subsection{Computation} \subsubsection{Tuning} All three MCMC algorithms require an appropriate covariance matrix - to move efficiently through the paramter space. The optimal - covaraince matrix for each method is a scaled version of the true - covarance matrix. To obtain this covariance matrix, we performed + to move efficiently through the parameter space. The optimal + covariance matrix for each method is a scaled version of the true + covariance matrix. To obtain this covariance matrix, we performed several runs of the all-components algorithm. After each run, two calculations were performed. First, an estimate of the posterior covariance matrix was computed from the MCMC output to be used in the next run. Second, \textit{mcgibbsit} (see \ref{convergence} below) was used to estimate the necessary length of the next - run. Only three runs were required to achieve both an acceptible + run. Only three runs were required to achieve both an acceptable estimate according to \textit{mcgibbsit} and a covariance matrix with good performance for all three algorithms. \subsubsection{Execution} To enable a ``fair'' comparison between the methods, each MCMC - algorithm was started at the same starting point (jittered for - the multiple starting states of the NKC) and the proposal covariance - was set to an appropriately scaled version of the posterior - covariance matrix described above. With this covariance matrix, each - of the algorithms demonstrated an appropriate acceptance rate. + algorithm was started at the same starting point (jittered to + provide the multiple starting states of the NKC) and the proposal + covariance was set to an appropriately scaled version of the + posterior covariance matrix described above. With this covariance + matrix, each of the algorithms demonstrated an appropriate + acceptance rate. \subsubsection{Convergence} \label{convergence} @@ -275,10 +278,10 @@ The MCMC simulations were run until the Markov chains have reached stable distributions as assessed by the \textit{mcgibbsit} algorithm \cite{Warnes00}. \textit{mcgibbsit} calculates the - number of iterations necessary to estimate a user-specfied quantile + number of iterations necessary to estimate a user-specified quantile $q$ to within $\pm r$ with probability $s$, i.e., \textit{mcgibbsit} indicates when the MCMC sampler has run long - enough to provide good confidence interval estimates. The defaults, + enough to provide good credible interval estimates. The defaults, which are used in this paper, are $q=0.025$, $r=0.0125$, and $s=0.95$. These values generate estimates of the 2.5\% and 97.5\% quantiles to $\pm$ 0.0125 quantiles with probability 0.95. @@ -290,7 +293,7 @@ be heavily influenced by the initial starting point, and can result in bias in the computed estimates. We adopted an alternative approach to avoiding bias in our estimates; each MCMC sampler was - run until \textit{mcgibbsit} was satistified that the quantile + run until \textit{mcgibbsit} was satisfied that the quantile estimates were properly estimated. One necessary condition for \textit{mcgibbsit} to declare that these quantiles are properly estimated is stability of these estimates. This will only occur @@ -299,7 +302,7 @@ posterior distribution. Likewise, the high serial correlation of most MCMC samples has - induced many MCMC papers to recommend \emph{thinning} the ouput of + induced many MCMC papers to recommend \emph{thinning} the output of the MCMC sampler by discarding all but 1 out of,say, every 10 or 100 generated values. While this does not improve in the quality of estimates obtained from the MCMC run itself, it can dramatically Modified: trunk/PathwayModeling/thesispaper/paper.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/paper.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/paper.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -15,6 +15,7 @@ \usepackage{dchem} %\usepackage{color} \usepackage{setspace} +\usepackage{endfloat} \bibliographystyle{vancouver} \newcommand{\deriv}[2]{\ensuremath{\frac{\mathrm{d} #1}{\mathrm{d} #2}}} Modified: trunk/PathwayModeling/thesispaper/paper.pdf =================================================================== (Binary files differ) Modified: trunk/PathwayModeling/thesispaper/priorPlot.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/priorPlot.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/priorPlot.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -1,4 +1,6 @@ -\begin{figure} + +\begin{figure*} + \centering <<priorPlot,echo=F,eval=T,fig=T,width=6,height=3>>= @@ -17,7 +19,11 @@ mtext(text="3", side=1, at=3, line=0.5) @ -\caption{Prior density (scaled $\chi^2 <- 5$) for $\mu_i$ and $\sigma$ - model parameters.} +\caption{ + Prior density (scaled $\chi^2 <- 5$) for $\mu_i$ and $\sigma$ + model parameters. +} \label{priorPlot} -\end{figure} + +\end{figure*} + Modified: trunk/PathwayModeling/thesispaper/results.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/results.Snw 2007-04-17 03:13:14 UTC (rev 1089) +++ trunk/PathwayModeling/thesispaper/results.Snw 2007-04-17 15:00:06 UTC (rev 1090) @@ -15,14 +15,16 @@ The effect of the number of data points within a fixed time window on the parameter distributions can be seen in Figure~\ref{converged}. +% \SweaveInput{figure5} -We see some improved precision as the number of points -increases from 16 to 25 but it is not pronounced. Overall, the -reduction in width varied from 18-fold for $d_9$ and 9-fold for $d_1$ to -1.5-fold for $d_4$. -Figure~\ref{scatter} is an example of a bivariate scatterplot of the -distributions. There is correlation between some pairs of parameters, e.g., -$d_1$ -- $d_2$, but no evidence of multi-modality. +% +We see some improved precision as the number of +points increases from 16 to 25 but it is not pronounced. Overall, the +reduction in width varied from 18-fold for $d_9$ and 9-fold for $d_1$ +to 1.5-fold for $d_4$. Figure~\ref{scatter} is an example of a +bivariate scatterplot of the distributions. There is strong +correlation between some pairs of parameters, e.g., $d_1$ -- $d_2$, +but no evidence of multi-modality. \SweaveInput{figure6} The value of the posterior probability density for inference was @@ -30,7 +32,7 @@ maximum posterior density estimate for the model parameters and the 95\% credible intervals. The fits of the resulting models to the 16-point data set is shown in Figure~\ref{fits}. Quantitative measures -of the fits for +of the fits for % \SweaveInput{figure7} % This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |