[R-gregmisc-users] SF.net SVN: r-gregmisc: [952] trunk/PathwayModeling/thesispaper
Brought to you by:
warnes
Menu
▾
▴
[R-gregmisc-users] SF.net SVN: r-gregmisc: [952] trunk/PathwayModeling/thesispaper
|
From: <r_b...@us...> - 2006-04-03 15:34:25
|
Revision: 952 Author: r_burrows Date: 2006-04-03 08:34:00 -0700 (Mon, 03 Apr 2006) ViewCVS: http://svn.sourceforge.net/r-gregmisc/?rev=952&view=rev Log Message: ----------- moved figures, tables to separate pages Modified Paths: -------------- trunk/PathwayModeling/thesispaper/introduction.Snw trunk/PathwayModeling/thesispaper/methods.Snw trunk/PathwayModeling/thesispaper/paper.Snw trunk/PathwayModeling/thesispaper/results.Snw Modified: trunk/PathwayModeling/thesispaper/introduction.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/introduction.Snw 2006-03-29 00:44:36 UTC (rev 951) +++ trunk/PathwayModeling/thesispaper/introduction.Snw 2006-04-03 15:34:00 UTC (rev 952) @@ -23,12 +23,6 @@ that explain hyperglycemic toxicity in diabetes \cite{Oates02,Brownlee01}. One mechanism is illustrated in Figure~\ref{glycolysis}. -\begin{figure} - \centering - \includegraphics[scale=0.9]{figures/glycolysis} - \caption{Production of methylglyoxal in hyperglycemia} - \label{glycolysis} -\end{figure} Some of the damage to tissue components seen in diabetes is thought to be the result of reactions with glycolytic intermediates such as methylglyoxal. Nishikawa et al. \cite{Nishikawa00} have Modified: trunk/PathwayModeling/thesispaper/methods.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/methods.Snw 2006-03-29 00:44:36 UTC (rev 951) +++ trunk/PathwayModeling/thesispaper/methods.Snw 2006-04-03 15:34:00 UTC (rev 952) @@ -16,10 +16,8 @@ # trim the Gillespie output rawdata <- read.table('data/gillespie.out',header=T)[2528:7160,] # sampling times: -# times12 <- c(16:19,20.5,21,seq(24,39,3)) times16 <- c(16:19,20.5,21,seq(22,40,2)) -# times25 <- c(16:19,20.5,21:40) -# generating the MCMC input files, e.g., input16.dat +# generating the MCMC input file input16.dat data <- getPaperData(rawdata,20,times16,npoints=3, varVector=c(50,50,50,50,50)) rates <- getPaperRates(data,1:12,13:48) temp <- matrix(0,nrow=16,ncol=11) @@ -60,23 +58,6 @@ For the 5-reaction sequence of reactions the perturbation of $R1$ at $time = 20$ results in the time courses plotted in Figure~\ref{pulse}. -<<echo=F,eval=F>>=2 -source("R/pulse.R") -rawdata <- read.table("data/rawdata.dat",header=T) -attach(rawdata) -pdf(file="figures/tempDir/pulse.pdf",width=9,height=5) -pulse() -dev.off() -detach(rawdata) -@ - \begin{figure} - \centering - \includegraphics[scale=0.5]{figures/tempDir/pulse} - \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} For each time point there are 5 values for the reactant concentrations and 5 values for the estimated reaction rates. The data at each time point was generated by sampling 3 values for @@ -141,32 +122,6 @@ \deriv{R5}{t} &= v_5 = \frac{d_7R4}{d_8 + R4} - d_9R5 \end{align*} -\begin{table} - \begin{center} - \caption{Rate constants used in the Gillespie simulator} - \begin{tabular}{c|r} - rate constant & value \\ - \hline - $k_1$ & 1.600 \\ - $k_2$ & 0.540 \\ - $k_3$ & 19.500 \\ - $k_4$ & 2.125 \\ - $k_5$ & 0.190 \\ - $k_6$ & 8.460 \\ - $k_7$ & 2.077 \\ - $k_8$ & 0.090 \\ - $k_9$ & 4.560 \\ - $k_{10}$ & 1.400 \\ - $k_{11}$ & 0.106 \\ - $k_{12}$ & 3.670 \\ - $k_{13}$ & 1.640 \\ - $k_{14}$ & 0.400 \\ - \hline - \end{tabular} - \label{ki} - \end{center} -\end{table} - \subsection{Statistical Models} The statistical model of the parameters is @@ -237,24 +192,6 @@ The all-components algorithm (Figure~\ref{allComp}) changes all the components simultaneously by sampling from a multivariate normal distribution centered at the current point. - \begin{figure} - \centering - \subfloat[1-component]{\includegraphics[scale=0.50]{figures/1comp} - \label{1comp}} - \qquad - \subfloat[all-components]{\includegraphics[scale=0.50]{figures/allcomp} - \label{allComp}} - \qquad - \subfloat[NKC]{\includegraphics[scale=0.50]{figures/nkc}\label{nkc}} - \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.} - \label{1comp_vs_all} - \end{figure} The component-wise Metropolis algorithm has the advantage of simplicity but may move very slowly if the components are highly correlated. The all-components Metropolis avoids the problems with Modified: trunk/PathwayModeling/thesispaper/paper.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/paper.Snw 2006-03-29 00:44:36 UTC (rev 951) +++ trunk/PathwayModeling/thesispaper/paper.Snw 2006-04-03 15:34:00 UTC (rev 952) @@ -1,5 +1,5 @@ -\documentclass{elsart} -% \documentclass[doublespacing]{elsart} +% \documentclass{elsart} +\documentclass[doublespacing]{elsart} \usepackage{graphicx} \usepackage[nogin]{Sweave} \usepackage[numbers]{natbib} \bibliographystyle{plainnat} @@ -26,7 +26,8 @@ \begin{frontmatter} \title{Statistical Modeling of Biochemical Pathways} -\author{Robert B. Burrows} +% \author{Robert B. Burrows}\footnote{corresponding author:\\Robert Burrows\\1506 Chopmist Hill Road\\North Scituate, RI 02857\\401-647-5290\\rb...@ne...} +\author{Robert B. Burrows}\footnote{corresponding author:\\Robert Burrows\\[-4mm]1506 Chopmist Hill Road\\[-4mm]North Scituate, RI 02857\\[-4mm]401-647-5290\\[-4mm]rb...@ne...} \address{New England Biometrics, North Scituate, RI} \author{Gregory R. Warnes} \address{Pfizer, Inc., Groton, CT} @@ -62,4 +63,42 @@ \bibliography{./refs} +\SweaveInput{table1} + +\clearpage + +\SweaveInput{table2} + +\clearpage + +\SweaveInput{figure1} + +\clearpage + +\SweaveInput{figure2} + +\clearpage + +\SweaveInput{figure3} + +\clearpage + +\SweaveInput{figure4} + +\clearpage + +\SweaveInput{figure5} + +\clearpage + +\SweaveInput{figure6} + +\clearpage + +\SweaveInput{figure7} + +\clearpage + +\SweaveInput{figure8} + \end{document} Modified: trunk/PathwayModeling/thesispaper/results.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/results.Snw 2006-03-29 00:44:36 UTC (rev 951) +++ trunk/PathwayModeling/thesispaper/results.Snw 2006-04-03 15:34:00 UTC (rev 952) @@ -9,70 +9,8 @@ (Figure~\ref{densities}), were skewed to the right and had thicker tails than normal distributions. -<<fig4,echo=F,eval=F>>=3 -source("R/paperSSQ.R") -source("R/getModes.R") -source("R/do.optim.R") -source("R/plotDensities.R") -output16 <- read.table("data/output16.AllComp.thinned") -input16 <- read.table("data/input16.dat", header=T) -attach(input16) -temp<-output16[sample(1:20000,5000),] -temp2<-numeric(5000) -for (i in 1:5000) temp2[i]<-paperSSQ(temp[i,]) -temp<-temp[order(temp2),] -write.table(temp[1:4750,],file="figures/tempDir/Dsamp.dat",row.names=F,col.names=F) -mcmcML <- as.numeric(temp[1,-10]) -write.table(mcmcML,file="figures/tempDir/mcmcML.dat",row.names=F,col.names=F) -modes<-getModes(output16) -temp<-do.optim(start=modes[-10],scail=modes[-10]) -temp<-do.optim(start=temp$par,scail=temp$par) -optimML<-temp$par -write.table(optimML,file="figures/tempDir/optimML.dat",row.names=F,col.names=F) -SD <- numeric(9) -for (i in 1:9) SD[i] <- sd(output16[,i]) -pdf(file="figures/tempDir/densities.pdf",width=6,height=6) -par(mfrow=c(3,3)) -plotDensities(output16) -dev.off() -par(mfrow=c(1,1)) -detach(input16) -@ -\begin{figure} - \centering - \includegraphics[scale=0.8]{figures/tempDir/densities} - \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. \textit{Red vertical lines indicate the parameters values that - minimize the mean squared residuals.}} - \label{densities} -\end{figure} The effect of the number of data points on the parameter distributions can be seen in Figure~\ref{converged}. -<<fig5,echo=F,eval=F>>=4 -source("R/plotDensity.R") -source("R/plotConverged.R") -output12 <- read.table("data/output12.AllComp.thinned") -output16 <- read.table("data/output16.AllComp.thinned") -output25 <- read.table("data/output25.AllComp.thinned") -pdf(file="figures/tempDir/converged.pdf",width=6,height=6) -par(mfrow=c(3,3)) -plotConverged() -dev.off() -par(mfrow=c(1,1)) -@ -\begin{figure} - \centering - \includegraphics[scale=0.9]{figures/tempDir/converged} - \caption{\textit{Posterior distributions from different numbers} of data - points for the all-components algorithm. (---------) prior - distribution; ({\color{red} - - - -}) 12 points; ({\color{green} - $\cdots\cdots$}) 16 points; ({\color{blue} -- $\cdot$ -- $\cdot$ --}) 25 - points} - \label{converged} -\end{figure} 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 @@ -81,32 +19,6 @@ distributions. There is correlation between some pairs of parameters, e.g., $d_1$ -- $d_2$, but no evidence of multi-modality. -<<fig6, echo=F,eval=F>>=5 -source("R/paperPairs.R") -mcmcML <-scan("figures/tempDir/mcmcML.dat") -output16 <- read.table("data/output16.AllComp.thinned") -mcmc.cor <- function(x,y) { - par(usr=c(0,1,0,1)) - r <- cor(x,y) - txt <- format(r,digits=2) - text(0.5,0.5,txt,cex=1.4) - } -pdf(file="figures/tempDir/scatterPlot.pdf",width=8,height=8) -par(pch='.') -paperPairs(output16[sample(1:20000,2000),-10],labels=c('d1','d2','d3','d4','d5','d6','d7','d8','d9'),lower.panel=mcmc.cor) -dev.off() -par(pch=1) -@ -\begin{figure} - \centering - \includegraphics[scale=0.6]{figures/tempDir/scatterPlot} - \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 red lines indicate the maximum likelihood estimates of the parameters.} - \label{scatter} -\end{figure} - The value of the probability density for inference was assessed graphically. The probability density was used to find the maximum likelihood estimate for the model parameters and the 95\% confidence @@ -114,53 +26,6 @@ is shown in Figure~\ref{fits}. Quantitative measures of the fits for all the algorithms are given in Table~\ref{MSq}. -<<fig7, echo=F, eval=F>>=6 -source("R/paperSSQ.R") -source("R/getModes.R") -source("R/do.optim.R") -source("R/fitPaper.R") -source("R/get95CI.R") -source("R/plotVi.R") -mcmcML <-scan("figures/tempDir/mcmcML.dat") -optimML <-scan("figures/tempDir/optimML.dat") -output16 <- read.table("data/output16.AllComp.thinned") -input16 <- read.table("data/input16.dat",header=T) -attach(input16) -Dsamp <- as.matrix(read.table("figures/tempDir/Dsamp.dat")) -pdf(file="figures/tempDir/V%dfitted.pdf",onefile=FALSE,width=10,height=5.5) -plotVi(Dsamp,input16) -dev.off() -detach(input16) -@ -\begin{figure} - \centering - \subfloat[v2]{\includegraphics[scale=0.33]{figures/tempDir/V1fitted}} - \subfloat[v3]{\includegraphics[scale=0.33]{figures/tempDir/V2fitted}}\\ - \subfloat[v4]{\includegraphics[scale=0.33]{figures/tempDir/V3fitted}} - \subfloat[v5]{\includegraphics[scale=0.33]{figures/tempDir/V4fitted}} - \caption{Curves fit to the 16-point data with the all-components - algorithm. The green curves are drawn with the maximum likelihood estimates for the parameters found with the L-BFGS-B algorithm \cite{Byrd95} as implemented in the \textit{R optim()} function. } - \label{fits} -\end{figure} - -\begin{table}[h] - \begin{center} - \caption{Mean squared residuals} - \begin{tabular}{c||c|c|c||c|c|c} - & - \multicolumn{3}{c||}{mean residual SSQ $\times 10^{-4}$} & \multicolumn{3}{c}{$R_{adj}^2$}\\ - \cline{2-7} - algorithm & 12 pt. & 16 pt. & 25 pt. & - 12 pt. & 16 pt. & 25 pt.\\ - \hline - 1-comp & 1.34 & 0.83 & 1.06 & 0.87 & 0.92 & 0.86\\ - all-comp & 0.74 & 0.93 & 0.74 & 0.93 & 0.90 & 0.90 \\ - NKC & 0.86 & 0.73 & 0.71 & 0.92 & 0.92 & 0.91 - \end{tabular} - \label{MSq} - \end{center} -\end{table} - Rates of convergence are illustrated in Figure~\ref{SSQvsIter}. The mean sums of squared residuals (SSQ) are plotted vs. the number of likelihood evaluations for the three algorithms. The relatively @@ -173,31 +38,6 @@ This is not a problem with the all-component Metropolis and NKC algorithms since they update all the parameters at each iteration. -<<fig8, echo=F,eval=F>>=7 -source("R/getMSDvsEvals.R") -source("R/plotMSD16.R") -source('R/paperSSQ.R') -input16 <- read.table('data/input16.dat',header=T) -attach(input16) -OneComp.MSD.dat <- read.table("data/1comp.MSD.dat") -AllComp.MSD.dat <- read.table("data/AllComp.MSD.dat") -for (i in 1:10) { -file <- sub('#',as.character(i),'NKC.MSD.dat#') -assign(file, read.table(paste("data/",file,sep=""))) -} -pdf(file="figures/tempDir/MSR16.pdf",width=9,height=6) -plotMSD16() -dev.off() -detach(input16) -@ -\begin{figure} - \centering - \includegraphics[scale=0.5]{figures/tempDir/MSR16} - \caption[Mean squared residuals vs. number of likelihood - evaluations]{Mean squared residuals vs. number of likelihood - evaluations for the 5-reaction model} - \label{SSQvsIter} -\end{figure} The actual convergence times on a 3.2GHz P4 machine with 1GB of memory are 1--2 hours for the 1-component algorithm and about 1 minute for the all-components and Normal Kernel Coupler algorithms. This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
