r-gregmisc-users — Discussion list for R gregmisc packages (gplot, gmodels, gdata, gtools)

[R-gregmisc-users] SF.net SVN: r-gregmisc: [951] trunk/PathwayModeling/thesispaper/refs.bib

[<r_b...@us...>]

Revision: 951
Author:   r_burrows
Date:     2006-03-28 16:44:36 -0800 (Tue, 28 Mar 2006)
ViewCVS:  http://svn.sourceforge.net/r-gregmisc/?rev=951&view=rev

Log Message:
-----------
added Byrd95 reference

Modified Paths:
--------------
    trunk/PathwayModeling/thesispaper/refs.bib
Modified: trunk/PathwayModeling/thesispaper/refs.bib
===================================================================
--- trunk/PathwayModeling/thesispaper/refs.bib	2006-03-29 00:37:26 UTC (rev 950)
+++ trunk/PathwayModeling/thesispaper/refs.bib	2006-03-29 00:44:36 UTC (rev 951)
@@ -207,3 +207,12 @@
     "\url{http://www.sourceforge.net/projects/hydra-mcmc/} (17 Jan
     06)",
 }
+@article{Byrd95,
+    author=	"R. H. Byrd and P. Lu and J. Nocedal and C. Zhu",
+    title=	"A limited memory algorithm for bound constrained optimization",
+    journal=	"SIAM Journal for Scientific Computing",
+    year=	"1995",
+    volume=	"16",
+    number=	"",
+    pages=	"1190--1208"
+}


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[R-gregmisc-users] SF.net SVN: r-gregmisc: [950] trunk/PathwayModeling/thesispaper

[<r_b...@us...>]

Revision: 950
Author:   r_burrows
Date:     2006-03-28 16:37:26 -0800 (Tue, 28 Mar 2006)
ViewCVS:  http://svn.sourceforge.net/r-gregmisc/?rev=950&view=rev

Log Message:
-----------
incorporated Greg's edits

Modified Paths:
--------------
    trunk/PathwayModeling/thesispaper/R/plotMSD16.R
    trunk/PathwayModeling/thesispaper/discussion.Snw
    trunk/PathwayModeling/thesispaper/introduction.Snw
    trunk/PathwayModeling/thesispaper/methods.Snw
    trunk/PathwayModeling/thesispaper/paper.Snw
    trunk/PathwayModeling/thesispaper/results.Snw
Modified: trunk/PathwayModeling/thesispaper/R/plotMSD16.R
===================================================================
--- trunk/PathwayModeling/thesispaper/R/plotMSD16.R	2006-03-25 14:42:14 UTC (rev 949)
+++ trunk/PathwayModeling/thesispaper/R/plotMSD16.R	2006-03-29 00:37:26 UTC (rev 950)
@@ -9,11 +9,14 @@
     }
     temp3 <- temp3/10
     ymax <- max(temp1[,2],temp2[,2],temp3[,2])
-    plot(log10(temp1[,1]),temp1[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),xlab="log10(#evals)",ylab="mean squared residual")
+#    plot(log10(temp1[,1]),temp1[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),xlab="log10(#evals)",ylab="mean squared residual")
+    plot(temp1[,1],temp1[,2],type='l',log="x",ylim=c(5000,ymax),xlim=c(1,10^7),xlab="#evals",ylab="mean squared residual")
     par(new=T)
-    plot(log10(temp2[,1]),temp2[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),col=2,lty=2,xlab="",ylab="")
+#    plot(log10(temp2[,1]),temp2[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),col=2,lty=2,xlab="",ylab="")
+    plot(temp2[,1],temp2[,2],type='l',log="x",ylim=c(5000,ymax),xlim=c(1,10^7),col=2,lty=2,xlab="",ylab="")
     par(new=T)
-    plot(log10(temp3[,1]),temp3[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),col=3,lty=3,xlab="",ylab="")
-    legend(4.5,ymax,c("1-comp","all-comp","NKC"),lty=1:3,col=1:3)
+#    plot(log10(temp3[,1]),temp3[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),col=3,lty=3,xlab="",ylab="")
+    plot(temp3[,1],temp3[,2],type='l',log="x",ylim=c(5000,ymax),xlim=c(1,10^7),col=3,lty=3,xlab="",ylab="")
+    legend(10^4.5,ymax,c("1-comp","all-comp","NKC"),lty=1:3,col=1:3)
 }
 

Modified: trunk/PathwayModeling/thesispaper/discussion.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/discussion.Snw	2006-03-25 14:42:14 UTC (rev 949)
+++ trunk/PathwayModeling/thesispaper/discussion.Snw	2006-03-29 00:37:26 UTC (rev 950)
@@ -6,9 +6,9 @@
 
 We are interested in assessing the usefulness of Markov chain Monte Carlo 
 (MCMC) methods for the fitting of models of metabolic pathways. Fitting 
-these models can be difficult because the pathways are described by sets 
+metabolic pathway models can be difficult because pathways are described by sets 
 of reaction rate equations with overlapping sets of parameters. We have 
-found that the MCMC algorithms handled this situation well, and produced 
+found that the MCMC algorithms handle this situation well, and produce 
 reasonable joint probability densities for the model parameters. 
 This output can be used 
 for the estimation of confidence intervals for the parameters and the 
@@ -37,32 +37,22 @@
 
  All three algorithms produced very similar results but the execution times 
 varied considerably depending on the size of the model and the algorithm. 
-The all-component Metropolis and NKC algorithms converged fairly easily 
-but the 1-component algorithm was quite slow. This is undoubtably due to 
+The all-component Metropolis and NKC algorithms converged quickly 
+but the 1-component algorithm was quite slow to converge. This is undoubtably due to 
 the high correlations between some parameters which is evident in 
 Figure~\ref{scatter}. The 1-component algorithm has trouble sampling from 
 this density because it is restricted to movements parallel to the axes. 
 This is less of a problem for the all-components and NKC algorithms since 
 they can move diagonally.
 
-The similarity of the results obtained with the different algorithms is 
-apparent from Table~\ref{MSq}. The mean squared residuals do not vary 
+Each of the algorithms produced essentially identical results. The mean squared residuals do not vary 
 much, nor do the calculated values for $R_{adj}^2$. These last values are 
 indicative of a useful fit of the models to the data. These results
 also demonstrate the usefulness of MCMC simulation. The derived
 probability densities can be used to find good estimates for the
 parameters as well as reveal the correlations between the parameters.
-The result is a well-fitted model with realistic confidence intervals (Figure~\ref{fitted}).
+The result is a well-fitted model with realistic confidence intervals (Figure~\ref{fits}).
 
-It was of interest to see if the calculated parameter values could be used 
-to infer values for the underlying rate constants by fitting a simple 
-linear regression model. For this, expressions for the parameters as 
-polynomial functions of the rate constants were found by systematically 
-varying the rate constants according to a face-centered cubic experimental 
-design. It was found that the parameters are rather complex functions of 
-the rate constants, too complex to allow useful inferences concerning the 
-rate constants to be made given the parameter values using this method.
-
 In this paper we have developed models of linear sequences of irreversible 
 biochemical reactions and shown that several Markov chain methods can be 
 used to generate useful fits of the models to simulated data. We also have 

Modified: trunk/PathwayModeling/thesispaper/introduction.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/introduction.Snw	2006-03-25 14:42:14 UTC (rev 949)
+++ trunk/PathwayModeling/thesispaper/introduction.Snw	2006-03-29 00:37:26 UTC (rev 950)
@@ -33,32 +33,31 @@
 to be the result of reactions with glycolytic intermediates such as
 methylglyoxal. Nishikawa et al. \cite{Nishikawa00} have
 hypothesized that this is caused by the increased catabolism of glucose in diabetes
-which leads to an increase of superoxide. Superoxide is known to inhibit the enzyme
+which leads to an increase of superoxide anion ($O_2^-$). Superoxide is known to inhibit the enzyme
 glyceraldehyde-3-phosphate dehydrogenase (G3PDH) by about 50\% under
 these conditions. The inhibition of this
 enzyme then increases the concentration of glyceraldehyde-3-phosphate
 (G3P) and dihydroxyacetone phosphate
 (DHAP) which breaks down to form methylglyoxal, an active aldehyde which readily reacts with and inactivates proteins. To test this
 hypothesis we need to be able to answer
-questions such as: Is the observed inhibition of this enzyme by 50\% sufficient
+questions such as: Is the observed inhibition of glyceraldehyde-3-phosphate dehydrogenase by 50\% sufficient
 to explain the observed increase in methylglyoxal concentration? If
 not, then what other mechanisms might be operative? What is the most
 effective way to modulate these reactions?
 We address these questions by developing analytical methods that will allow us to characterize these types of biochemical systems at the molecular level.
 
-We will use a Bayesian modeling approach with fitting via Markov chain
-Monte Carlo (MCMC) simulation
+We have used a Bayesian modeling approach with fitting via Markov chain
+Monte Carlo (MCMC) methods
 \cite{Metropolis53,Gilks95} to estimate the kinetic parameters of the enzymes in
-a metabolic pathway. With this method we can combine known information
-about the pathway with newly-generated experimental data to produce
+a metabolic pathway. With this approach we can combine known information
+about the pathway with newly-generated experimental data to produce the
 joint
-probability densities for all of the parameters. This joint
+probability density of all of the parameters. This joint
 probability density can then be used to answer various questions
 about the behavior of the system.
 
-We will explore the utility of the the Bayesian approach using a
-simulation technique: First we will generate data for a model with
-known parameters, and then we will fit this data using Bayesian
+We have explored the utility of the the Bayesian approach using a
+simulation technique: First we generated data for hypothetical sequences of enzymatic reactions, and then we fit models to this data using Bayesian
 methods.
 The biochemical data is generated using Gillespie's
 stochastic simulation algorithm \cite{Gillespie77}. This algorithm is
@@ -75,13 +74,12 @@
 systems.
 
 We simulate perturbation experiments. In these
-experiments, the concentration of a reactant in a pathway is
+experiments, the concentration of the first reactant in a pathway is
 transiently increased and the concentrations of all the reactants is
 monitored as they return to a steady state. This type of
 experiment is quite feasible with cultured animal cells or suspensions
 of microbial cells such as yeast, and it yields information on the
-behavior of the pathway as a single system. The MCMC methods described
-here are very useful for analyzing this type of data.
+behavior of the pathway as a single system.
 
 MCMC simulation was found to be a very useful method for fitting
 systems of equations to data. It is capable of fitting several

Modified: trunk/PathwayModeling/thesispaper/methods.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/methods.Snw	2006-03-25 14:42:14 UTC (rev 949)
+++ trunk/PathwayModeling/thesispaper/methods.Snw	2006-03-29 00:37:26 UTC (rev 950)
@@ -4,11 +4,9 @@
 $Id$
 \end{verbatim}
 
-  \subsection{The Data}
+  \subsection{The Pathway}
   
-    The data is generated with Gillespie's
-    algorithm for stochastic simulation of chemical reactions
-    \cite{Gillespie77}. In this paper we present the results for a 
+    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. 
 <<data, echo=F, eval=F>>=1
@@ -22,7 +20,7 @@
 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
-data <- getPaperData(rawdata,20,times16,npoints=3, varVector=c(50,70,60,40,40))
+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)
 for (i in 1:16) temp[i,] <- (rates[3*i,]+rates[3*i-1,]+rates[3*i-2,])/3
@@ -43,18 +41,22 @@
       R5 &\yields^{k_{14}}& sink
     \end{chemarray*}
 
+  \subsection{The Experiment}
+
     To estimate the kinetic parameters of the enzymes in a sequence of
     reactions, we need to find expressions of the form
     \[
     \text{reaction velocity} = f(\mathit{[substrate],[enzyme]})
     \]
-    These expressions can be found by perturbing a reaction network
+    These expressions are found by perturbing a reaction network
     and observing the metabolite concentrations as the system relaxes
-    back to a steady state. The metabolite concentrations are
-    measured, and the velocities are calculated as the
-    slopes of the curves of metabolite concentration as functions of
+    back to a steady state. Specifically, the Gillespie simulation was run with the initial values for the reactant concentrations until a steady state was achieved. The steady state appeared to have been reached by \textit{time} = 15. At \textit{time} = 20 the concentration of $R1$ was increased to 10000. At the indicated time points the metabolite concentrations are
+    measured, and each reaction velocity is calculated as the
+    slope of the curve of metabolite concentration as a function of
     time.
 
+    \subsection{The data}
+
     For the 5-reaction sequence of reactions 
     the perturbation of $R1$ at $time = 20$ results in the time
     courses plotted in Figure~\ref{pulse}.
@@ -79,12 +81,12 @@
     concentrations and 5 values for the estimated reaction rates.
     The data at each time point was generated by sampling 3 values for
     each reactant from a normal distribution centered on the Gillespie
-    output and with a standard deviation of 50. Three sets of time
+    output and with a standard deviation of 50, i.e., the measurement error is 50. Three sets of time
     points containing 12,16, and 25 points were used.
     The reaction velocities were estimated with the
     \textit{smooth.spline()} function in \textit{R}~\cite{R}. 
 
-  \subsection{Biochemical Models}
+  \subsection{The Biochemical Model}
 
     Individual reactions were fit using the Michaelis-Menten
     equation \cite{Michaelis13} for individual enzymes. This is a
@@ -183,10 +185,9 @@
     \]
     where the $\mu_i$ are the estimates of the parameter values before
     any data are collected. With this distribution, the $d_i$s will
-    have a mode of $\mu_i$ and large values will be unlikely. In this situation, $\mu_i \equiv 50$ for
+    have a mode of $\mu_i$ and large values will be unlikely. For this situation, we selected $\mu_i \equiv 50$ for
     the rate constants and $\sigma \equiv 2$ for the data standard
-    deviation seem appropriate, and are used in the
-    priors.
+    deviation.
   
     The likelihood is the probability of the
     estimated reaction velocities $v_j$ given the model.
@@ -217,7 +218,7 @@
 
     \subsubsection{MCMC sampling algorithms}
 
-    The performance of a MCMC simulation
+    The efficiency of a MCMC simulation
     \cite{Metropolis53,Hastings70} depends heavily on the
     method used to find candidate points to add to the Markov chain
     [$q(\cdot |\cdot)$]. 
@@ -230,7 +231,7 @@
     The component-wise Metropolis and the all-components Metropolis
     algorithms both operate on single Markov chains. The
     component-wise algorithm (Figure~\ref{1comp}) generates candidate points by
-    changing the value of only one component (dimension) of the current point at
+    changing the value of only one component (dimension) of the current state at
     each iteration by sampling from a univariate normal distribution
     centered at the current point.
     The all-components algorithm (Figure~\ref{allComp}) changes all the components
@@ -257,7 +258,7 @@
     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
-    correlation by using a covariance matrix and updating two or more
+    correlation, if an appropriate covariance matrix is supplied, by updating all
     components simultaneously but may perform poorly in high
     dimensions. A problem with both of these algorithms
     is that if the component distributions have two or more modes
@@ -268,7 +269,7 @@
     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 kernel density estimate is generated using the entire set of
+    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
@@ -277,10 +278,10 @@
   \section{Convergence}
     
     The MCMC simulations are run until the Markov chains have reached
-    stable distributions as assessed by \textit{mcgibbsit()}
+    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 $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, which are used
-    in this paper, are $q=0.025$, $r=0.0125$, and $s=0.95$.
+    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.

Modified: trunk/PathwayModeling/thesispaper/paper.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/paper.Snw	2006-03-25 14:42:14 UTC (rev 949)
+++ trunk/PathwayModeling/thesispaper/paper.Snw	2006-03-29 00:37:26 UTC (rev 950)
@@ -34,16 +34,16 @@
 \address{University of Rhode Island, Kingston, RI}
   
 \begin{abstract}
-\textit{We have examined} the usefulness of \textit{Bayesian statistical methods} for the modeling
+We have examined the usefulness of Bayesian statistical methods for the modeling
 of biochemical reactions. With simulated data, it is
-shown that mechanistic models can \textit{effectively} be fit to sequences of enzymatic
-reactions. These methods have the advantages of being relatively easy
+shown that these methods can effectively fit mechanistic models of sequences of enzymatic
+reactions to experimental data. These methods have the advantages of being relatively easy
 to use and producing probability distributions for the model
-parameters rather than point estimates, \textit{allowing more informative statistical inferences to be drawn}.
+parameters rather than point estimates, allowing more informative inferences to be drawn.
 
 Three Markov chain Monte Carlo algorithms are used to fit models to
 data from a
-sequence of 4 enzymatic reactions. These algorithms 
+sequence of 4 enzymatic reactions. The algorithms 
 are evaluated with respect to
 the goodness-of-fit of the fitted models and the time to completion.
 It is shown that the algorithms produce essentially the same

Modified: trunk/PathwayModeling/thesispaper/results.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/results.Snw	2006-03-25 14:42:14 UTC (rev 949)
+++ trunk/PathwayModeling/thesispaper/results.Snw	2006-03-29 00:37:26 UTC (rev 950)
@@ -4,21 +4,10 @@
 $Id$
 \end{verbatim}
 
-The joint probability density for the parameters of a model of a
-metabolic pathway has been estimated using MCMC simulation. Three
-algorithms were used for the simulations and they converged to similar
-distributions. 
-
+All three algorithms converged to similar distributions and produced essentially the same metabolic pathway model.
 The marginal distributions for the parameters were mound-shaped 
-(Figure~\ref{densities}). They
-generally were skewed to the right and had thicker
-tails than normal distributions. The modes of the distributions 
-turned out to be poor estimates of the maximum likelihood values for 
-the parameters. The maximum likelihood values were all larger than 
-the modes, and the difference increased as the position of the reaction 
-in the sequence increased. Correlated pairs of parameters 
-(e.g., $d_1$ --$d_2$ and $d_3$ --$d_4$, Figure~\ref{scatter}) were similar relative distances from their 
-modes. The significance of all this is not clear at this point.
+(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")
@@ -84,7 +73,7 @@
   points}
   \label{converged}
 \end{figure}
-There is some narrowing of the distributions as the number of points
+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$.
@@ -114,11 +103,11 @@
   \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).}
+  triangle). The red lines indicate the maximum likelihood estimates of the parameters.}
   \label{scatter}
 \end{figure}
 
-The usefulness of the probability density for inference was assessed
+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
 intervals. The fits of the resulting models to the 16-point data set
@@ -203,12 +192,7 @@
 @
 \begin{figure}
   \centering
-%  \subfloat[12 points]{
-%    \includegraphics[scale=0.5]{figures/MSR12}}
-%  \subfloat[16 points]{
     \includegraphics[scale=0.5]{figures/tempDir/MSR16}
-%  \subfloat[25 points]{
-%    \includegraphics[scale=0.5]{figures/MSR25}}
   \caption[Mean squared residuals vs. number of likelihood
   evaluations]{Mean squared residuals vs. number of likelihood
     evaluations for the 5-reaction model}


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[R-gregmisc-users] SF.net SVN: r-gregmisc: [949] trunk/PathwayModeling/thesispaper/R/plotVi.R

[<r_b...@us...>]

Revision: 949
Author:   r_burrows
Date:     2006-03-25 06:42:14 -0800 (Sat, 25 Mar 2006)
ViewCVS:  http://svn.sourceforge.net/r-gregmisc/?rev=949&view=rev

Log Message:
-----------
added ML curves

Modified Paths:
--------------
    trunk/PathwayModeling/thesispaper/R/plotVi.R
Modified: trunk/PathwayModeling/thesispaper/R/plotVi.R
===================================================================
--- trunk/PathwayModeling/thesispaper/R/plotVi.R	2006-03-25 14:41:09 UTC (rev 948)
+++ trunk/PathwayModeling/thesispaper/R/plotVi.R	2006-03-25 14:42:14 UTC (rev 949)
@@ -17,7 +17,7 @@
     par(new=T)
     plot(times,temp2[[1]],ylim=c(low,high),type="l",xlab="",ylab="",col=2)
     abline(h=0)
-    legend(32,2000,c("data","fitted curve","95% CI"),lty=c(0,1,2),pch=c(1,NA,NA),col=c(1,2,2))
+    legend(32,2000,c("data","fitted curve","95% CI","ML"),lty=c(0,1,2,1),pch=c(1,NA,NA,NA),col=c(1,2,2,3))
     for (i in 3:5) {
         CI <- get95CI(i,input,1,Dsamp)
         low <-min(CI,-400)


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[R-gregmisc-users] SF.net SVN: r-gregmisc: [948] trunk/PathwayModeling/thesispaper

[<r_b...@us...>]

Revision: 948
Author:   r_burrows
Date:     2006-03-25 06:41:09 -0800 (Sat, 25 Mar 2006)
ViewCVS:  http://svn.sourceforge.net/r-gregmisc/?rev=948&view=rev

Log Message:
-----------
minor edits

Modified Paths:
--------------
    trunk/PathwayModeling/thesispaper/discussion.Snw
    trunk/PathwayModeling/thesispaper/introduction.Snw
    trunk/PathwayModeling/thesispaper/methods.Snw
    trunk/PathwayModeling/thesispaper/results.Snw
Modified: trunk/PathwayModeling/thesispaper/discussion.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/discussion.Snw	2006-03-20 17:31:46 UTC (rev 947)
+++ trunk/PathwayModeling/thesispaper/discussion.Snw	2006-03-25 14:41:09 UTC (rev 948)
@@ -9,11 +9,10 @@
 these models can be difficult because the pathways are described by sets 
 of reaction rate equations with overlapping sets of parameters. We have 
 found that the MCMC algorithms handled this situation well, and produced 
-reasonable joint probability densities for the model parameters. In 
-addition to point estimates for the parameters, this output can be used 
+reasonable joint probability densities for the model parameters. 
+This output can be used 
 for the estimation of confidence intervals for the parameters and the 
-detection of correlations and multimodality. This additional information 
-cannot be obtained with any other method in routine use.
+detection of correlations and multimodality. Thus MCMC compares favorably with maximum likelyhood methods that produce point estimates of the parameters and nonlinear regression methods that find approximations of the parameters and their variances.
 
 The models used here have the form of the Hill function, 
 $\frac{x^n}{\theta^n+x^n}$, with an exponent of 1. This form was chosen 
@@ -36,7 +35,7 @@
 Normal Kernel Coupler, is more computationally intensive but it can sample 
 from complex densities more efficiently than the other two algorithms.
 
- All four algorithms produced very similar results but the execution times 
+ All three algorithms produced very similar results but the execution times 
 varied considerably depending on the size of the model and the algorithm. 
 The all-component Metropolis and NKC algorithms converged fairly easily 
 but the 1-component algorithm was quite slow. This is undoubtably due to 
@@ -53,7 +52,7 @@
 also demonstrate the usefulness of MCMC simulation. The derived
 probability densities can be used to find good estimates for the
 parameters as well as reveal the correlations between the parameters.
-The result is a well-fitted model with realistic confidence intervals.
+The result is a well-fitted model with realistic confidence intervals (Figure~\ref{fitted}).
 
 It was of interest to see if the calculated parameter values could be used 
 to infer values for the underlying rate constants by fitting a simple 

Modified: trunk/PathwayModeling/thesispaper/introduction.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/introduction.Snw	2006-03-20 17:31:46 UTC (rev 947)
+++ trunk/PathwayModeling/thesispaper/introduction.Snw	2006-03-25 14:41:09 UTC (rev 948)
@@ -5,20 +5,21 @@
 \end{verbatim}
 
 We know a great deal about the individual components of living
-organisms but much less about how the properties of these components
-affect the systems of which they are a part. Our ignorance in this
+organisms but much less about the properties of the systems
+ of which they are a part. Our ignorance in this
 regard prevents us from developing a deeper theoretical understanding
 of living systems and from developing more effective methods for
 altering the system behavior in beneficial ways.
 
 In order to gain some insight into the properties of biological
-systems, we would like to use quantitative models of sequences of
-reactions.
+systems, we would like to use quantitative models that describe the processes at the molecular level.
  Specifically, we would like to develop methods that can be used to
-deduce the kinetic parameters of the enzymes in a pathway. With this
+deduce the kinetic parameters of enzyme-catalyzed reactions. With this
 information, it should be possible to identify the most effective
 control points in a pathway and modify the pathway behavior in
-predictable ways. For example, there are thought to be four mechanisms
+predictable ways.
+
+An example of how a model of a biochemical pathway could be used is in the elucidation of a mechanism for the pathology seen with diabetes. There are thought to be four mechanisms
 that explain hyperglycemic toxicity in diabetes
 \cite{Oates02,Brownlee01}. One mechanism is illustrated in
 Figure~\ref{glycolysis}.
@@ -37,14 +38,13 @@
 these conditions. The inhibition of this
 enzyme then increases the concentration of glyceraldehyde-3-phosphate
 (G3P) and dihydroxyacetone phosphate
-(DHAP) which breaks down to form methylglyoxal. To test this
+(DHAP) which breaks down to form methylglyoxal, an active aldehyde which readily reacts with and inactivates proteins. To test this
 hypothesis we need to be able to answer
 questions such as: Is the observed inhibition of this enzyme by 50\% sufficient
 to explain the observed increase in methylglyoxal concentration? If
 not, then what other mechanisms might be operative? What is the most
 effective way to modulate these reactions?
-We address these questions by developing analytical methods that will allow us to characterize
-the intact systems at the molecular level.
+We address these questions by developing analytical methods that will allow us to characterize these types of biochemical systems at the molecular level.
 
 We will use a Bayesian modeling approach with fitting via Markov chain
 Monte Carlo (MCMC) simulation
@@ -83,13 +83,6 @@
 behavior of the pathway as a single system. The MCMC methods described
 here are very useful for analyzing this type of data.
 
-The MCMC results are compared to
-parameter estimates found by least squares fitting of the models to
-the data. The modes of the MCMC parameter distributions are
-essentially identical to the least squares parameter values, and this
-is taken as evidence that the MCMC simulation is producing valid
-estimates.
-
 MCMC simulation was found to be a very useful method for fitting
 systems of equations to data. It is capable of fitting several
 equations simultaneously, and does not require an excessive number of

Modified: trunk/PathwayModeling/thesispaper/methods.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/methods.Snw	2006-03-20 17:31:46 UTC (rev 947)
+++ trunk/PathwayModeling/thesispaper/methods.Snw	2006-03-25 14:41:09 UTC (rev 948)
@@ -11,7 +11,7 @@
     \cite{Gillespie77}. 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. 
-<<data, echo=F, eval=T>>=1
+<<data, echo=F, eval=F>>=1
 # gillespie.out - output from the Gillespie simulation
 source('R/getPaperData.R')
 source('R/getPaperRates.R')
@@ -58,8 +58,7 @@
     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=T>>=2
+<<echo=F,eval=F>>=2
 source("R/pulse.R")
 rawdata <- read.table("data/rawdata.dat",header=T)
 attach(rawdata)
@@ -78,7 +77,6 @@
     \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
     each reactant from a normal distribution centered on the Gillespie
     output and with a standard deviation of 50. Three sets of time
@@ -222,16 +220,9 @@
     The performance of a MCMC simulation
     \cite{Metropolis53,Hastings70} depends heavily on the
     method used to find candidate points to add to the Markov chain
-    [$q(\cdot |\cdot)$]. A
-    standard proposal distribution for the continuous problem is a
-    multivariate normal distribution. However,
-    the size of the space to be sampled grows exponentially with the
-    number of dimensions and if the joint density is multimodal,
-    finding all the modes can be a problem. For this reason proposal
-    distributions have been developed that result in more efficient
-    sampling than a simple one-component-at-a-time sampler.
-    
-    Three algorithms are implemented in the Hydra
+    [$q(\cdot |\cdot)$]. 
+    For this reason different methods have been developed to increase the efficiency of sampling from the sample space. We have used 
+    three algorithms from the Hydra
     library \cite{Hydra}:
      component-wise Metropolis, all-components Metropolis, 
     and Normal Kernel Coupler \cite{Warnes00}.

Modified: trunk/PathwayModeling/thesispaper/results.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/results.Snw	2006-03-20 17:31:46 UTC (rev 947)
+++ trunk/PathwayModeling/thesispaper/results.Snw	2006-03-25 14:41:09 UTC (rev 948)
@@ -12,15 +12,15 @@
 The marginal distributions for the parameters were mound-shaped 
 (Figure~\ref{densities}). They
 generally were skewed to the right and had thicker
-tails than normal distributions. /textit{The modes of the distributions 
+tails than normal distributions. The modes of the distributions 
 turned out to be poor estimates of the maximum likelihood values for 
 the parameters. The maximum likelihood values were all larger than 
 the modes, and the difference increased as the position of the reaction 
 in the sequence increased. Correlated pairs of parameters 
 (e.g., $d_1$ --$d_2$ and $d_3$ --$d_4$, Figure~\ref{scatter}) were similar relative distances from their 
-modes. The significance of all this is not clear at this point.}
+modes. The significance of all this is not clear at this point.
 
-<<fig4,echo=F,eval=T>>=3
+<<fig4,echo=F,eval=F>>=3
 source("R/paperSSQ.R")
 source("R/getModes.R")
 source("R/do.optim.R")
@@ -62,8 +62,7 @@
 \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=T>>=4
+<<fig5,echo=F,eval=F>>=4
 source("R/plotDensity.R")
 source("R/plotConverged.R")
 output12 <- read.table("data/output12.AllComp.thinned")
@@ -86,15 +85,14 @@
   \label{converged}
 \end{figure}
 There is some narrowing of the distributions as the number of points
-increases from 16 to 25 but it is not pronounced. \textit{Overall, the 
+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$}.
-
+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.
 
-<<fig6, echo=F,eval=T>>=5
+<<fig6, echo=F,eval=F>>=5
 source("R/paperPairs.R")
 mcmcML <-scan("figures/tempDir/mcmcML.dat")
 output16 <- read.table("data/output16.AllComp.thinned")
@@ -127,7 +125,7 @@
 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=T>>=6
+<<fig7, echo=F, eval=F>>=6
 source("R/paperSSQ.R")
 source("R/getModes.R")
 source("R/do.optim.R")
@@ -147,12 +145,12 @@
 @
 \begin{figure}
   \centering
-  \subfloat[v2]{\includegraphics[scale=0.35]{figures/tempDir/V1fitted}}
-  \subfloat[v3]{\includegraphics[scale=0.35]{figures/tempDir/V2fitted}}\\
-  \subfloat[v4]{\includegraphics[scale=0.35]{figures/tempDir/V3fitted}}
-  \subfloat[v5]{\includegraphics[scale=0.35]{figures/tempDir/V4fitted}}
+  \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}
+  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}
 
@@ -176,7 +174,7 @@
 
 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. \textit{The relatively 
+likelihood evaluations for the three algorithms. The relatively 
 long convergence time for the 1-component Metropolis algorithm is a result 
 of the high correlations between some of the parameters (Figure~\ref{scatter}).
 With the 1-component Metropolis algorithm, one parameter is
@@ -184,9 +182,9 @@
 searches the parameter space by moving parallel to the axes and thus 
 it can sample a density with a diagonal axis of symmetry very slowly. 
 This is not a problem with the all-component Metropolis and NKC 
-algorithms since they update all the parameters at each iteration}.
+algorithms since they update all the parameters at each iteration.
 
-<<fig8, echo=F,eval=T>>=7
+<<fig8, echo=F,eval=F>>=7
 source("R/getMSDvsEvals.R")
 source("R/plotMSD16.R")
 source('R/paperSSQ.R')


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[R-gregmisc-users] SF.net SVN: r-gregmisc: [947] trunk/PathwayModeling/thesispaper/R/plotMSD16.R

[<r_b...@us...>]

Revision: 947
Author:   r_burrows
Date:     2006-03-20 09:31:46 -0800 (Mon, 20 Mar 2006)
ViewCVS:  http://svn.sourceforge.net/r-gregmisc/?rev=947&view=rev

Log Message:
-----------
changed positioning of legend

Modified Paths:
--------------
    trunk/PathwayModeling/thesispaper/R/plotMSD16.R
Modified: trunk/PathwayModeling/thesispaper/R/plotMSD16.R
===================================================================
--- trunk/PathwayModeling/thesispaper/R/plotMSD16.R	2006-03-20 17:29:20 UTC (rev 946)
+++ trunk/PathwayModeling/thesispaper/R/plotMSD16.R	2006-03-20 17:31:46 UTC (rev 947)
@@ -8,11 +8,12 @@
 	temp3 <- temp3 + getMSDvsEvals(eval(as.name(file)),2,20)
     }
     temp3 <- temp3/10
-    plot(log10(temp1[,1]),temp1[,2],type='l',ylim=c(5000,max(temp1[,2],temp2[,2],temp3[,2])),xlim=c(0,7),xlab="log10(#evals)",ylab="mean squared residual")
+    ymax <- max(temp1[,2],temp2[,2],temp3[,2])
+    plot(log10(temp1[,1]),temp1[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),xlab="log10(#evals)",ylab="mean squared residual")
     par(new=T)
-    plot(log10(temp2[,1]),temp2[,2],type='l',ylim=c(5000,max(temp1[,2],temp2[,2],temp3[,2])),xlim=c(0,7),col=2,lty=2,xlab="",ylab="")
+    plot(log10(temp2[,1]),temp2[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),col=2,lty=2,xlab="",ylab="")
     par(new=T)
-    plot(log10(temp3[,1]),temp3[,2],type='l',ylim=c(5000,max(temp1[,2],temp2[,2],temp3[,2])),xlim=c(0,7),col=3,lty=3,xlab="",ylab="")
-    legend(4.5,14000,c("1-comp","all-comp","NKC"),lty=1:3,col=1:3)
+    plot(log10(temp3[,1]),temp3[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),col=3,lty=3,xlab="",ylab="")
+    legend(4.5,ymax,c("1-comp","all-comp","NKC"),lty=1:3,col=1:3)
 }
 


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[R-gregmisc-users] SF.net SVN: r-gregmisc: [946] trunk/PathwayModeling/thesispaper

[<r_b...@us...>]

Revision: 946
Author:   r_burrows
Date:     2006-03-20 09:29:20 -0800 (Mon, 20 Mar 2006)
ViewCVS:  http://svn.sourceforge.net/r-gregmisc/?rev=946&view=rev

Log Message:
-----------
set eval=T for code chunks

Modified Paths:
--------------
    trunk/PathwayModeling/thesispaper/methods.Snw
    trunk/PathwayModeling/thesispaper/results.Snw
Modified: trunk/PathwayModeling/thesispaper/methods.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/methods.Snw	2006-03-20 17:10:18 UTC (rev 945)
+++ trunk/PathwayModeling/thesispaper/methods.Snw	2006-03-20 17:29:20 UTC (rev 946)
@@ -11,7 +11,7 @@
     \cite{Gillespie77}. 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. 
-<<data, echo=F, eval=F>>=1
+<<data, echo=F, eval=T>>=1
 # gillespie.out - output from the Gillespie simulation
 source('R/getPaperData.R')
 source('R/getPaperRates.R')
@@ -59,7 +59,7 @@
     the perturbation of $R1$ at $time = 20$ results in the time
     courses plotted in Figure~\ref{pulse}.
 
-<<echo=F,eval=F>>=2
+<<echo=F,eval=T>>=2
 source("R/pulse.R")
 rawdata <- read.table("data/rawdata.dat",header=T)
 attach(rawdata)

Modified: trunk/PathwayModeling/thesispaper/results.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/results.Snw	2006-03-20 17:10:18 UTC (rev 945)
+++ trunk/PathwayModeling/thesispaper/results.Snw	2006-03-20 17:29:20 UTC (rev 946)
@@ -20,7 +20,7 @@
 (e.g., $d_1$ --$d_2$ and $d_3$ --$d_4$, Figure~\ref{scatter}) were similar relative distances from their 
 modes. The significance of all this is not clear at this point.}
 
-<<fig4,echo=F,eval=F>>=3
+<<fig4,echo=F,eval=T>>=3
 source("R/paperSSQ.R")
 source("R/getModes.R")
 source("R/do.optim.R")
@@ -63,7 +63,7 @@
 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
+<<fig5,echo=F,eval=T>>=4
 source("R/plotDensity.R")
 source("R/plotConverged.R")
 output12 <- read.table("data/output12.AllComp.thinned")
@@ -94,7 +94,7 @@
 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
+<<fig6, echo=F,eval=T>>=5
 source("R/paperPairs.R")
 mcmcML <-scan("figures/tempDir/mcmcML.dat")
 output16 <- read.table("data/output16.AllComp.thinned")
@@ -127,7 +127,7 @@
 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
+<<fig7, echo=F, eval=T>>=6
 source("R/paperSSQ.R")
 source("R/getModes.R")
 source("R/do.optim.R")
@@ -186,7 +186,7 @@
 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
+<<fig8, echo=F,eval=T>>=7
 source("R/getMSDvsEvals.R")
 source("R/plotMSD16.R")
 source('R/paperSSQ.R')


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[R-gregmisc-users] SF.net SVN: r-gregmisc: [945] trunk/PathwayModeling/thesispaper/elsart.cls

[<r_b...@us...>]

Revision: 945
Author:   r_burrows
Date:     2006-03-20 09:10:18 -0800 (Mon, 20 Mar 2006)
ViewCVS:  http://svn.sourceforge.net/r-gregmisc/?rev=945&view=rev

Log Message:
-----------
added ID keyword

Modified Paths:
--------------
    trunk/PathwayModeling/thesispaper/elsart.cls

Property Changed:
----------------
    trunk/PathwayModeling/thesispaper/elsart.cls
Modified: trunk/PathwayModeling/thesispaper/elsart.cls
===================================================================
--- trunk/PathwayModeling/thesispaper/elsart.cls	2006-03-20 11:51:45 UTC (rev 944)
+++ trunk/PathwayModeling/thesispaper/elsart.cls	2006-03-20 17:10:18 UTC (rev 945)
@@ -22,7 +22,8 @@
  \def\RCSdate{#4}%
 }
 \readRCS
-$Header$
+%$Header$
+$Id$
 \def\@shortjnl{\relax}
  \def\esp@fileversion{2.18}
  \def\esp@filedate{2001/01/05}


Property changes on: trunk/PathwayModeling/thesispaper/elsart.cls
___________________________________________________________________
Name: svn:keywords
   - Author Date Id Revision
   + Id


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[R-gregmisc-users] SF.net SVN: r-gregmisc: [944] trunk/PathwayModeling/thesispaper/R/paperPairs.R

[<r_b...@us...>]

Revision: 944
Author:   r_burrows
Date:     2006-03-20 03:51:45 -0800 (Mon, 20 Mar 2006)
ViewCVS:  http://svn.sourceforge.net/r-gregmisc/?rev=944&view=rev

Log Message:
-----------
initial commit

Added Paths:
-----------
    trunk/PathwayModeling/thesispaper/R/paperPairs.R
Added: trunk/PathwayModeling/thesispaper/R/paperPairs.R
===================================================================
--- trunk/PathwayModeling/thesispaper/R/paperPairs.R	                        (rev 0)
+++ trunk/PathwayModeling/thesispaper/R/paperPairs.R	2006-03-20 11:51:45 UTC (rev 944)
@@ -0,0 +1,125 @@
+paperPairs <- function(x, ...) UseMethod("pairs")
+
+pairs.formula <-
+function(formula, data = NULL, ..., subset, na.action = stats::na.pass)
+{
+    m <- match.call(expand.dots = FALSE)
+    if(is.matrix(eval(m$data, parent.frame())))
+        m$data <- as.data.frame(data)
+    m$... <- NULL
+    m$na.action <- na.action # force in even if  default
+    m[[1]] <- as.name("model.frame")
+    mf <- eval(m, parent.frame())
+    pairs(mf, ...)
+}
+
+#################################################
+## some of the changes are from code
+## Copyright 1999 Dr. Jens Oehlschlaegel-Akiyoshi
+## Others are by BDR and MM
+#################################################
+
+pairs.default <-
+function (x, labels, panel = points, ...,
+          lower.panel = panel, upper.panel = panel,
+          diag.panel = NULL, text.panel = textPanel,
+          label.pos = 0.5 + has.diag/3,
+          cex.labels = NULL, font.labels = 1,
+          row1attop = TRUE, gap = 1)
+{
+    textPanel <-
+        function(x = 0.5, y = 0.5, txt, cex, font)
+            text(x, y, txt, cex = cex, font = font)
+
+    localAxis <- function(side, xpd, bg, col=NULL, main, oma, ...) {
+      ## Explicitly ignore any color argument passed in as
+      ## it was most likely meant for the data points and
+      ## not for the axis.
+      axis(side, xpd = NA, ...)
+    }
+
+    localPlot <- function(..., main, oma, font.main, cex.main) plot(...)
+    localLowerPanel <- function(..., main, oma, font.main, cex.main)
+        lower.panel(...)
+    localUpperPanel <- function(..., main, oma, font.main, cex.main)
+        upper.panel(...)
+
+    dots <- list(...); nmdots <- names(dots)
+    if (!is.matrix(x)) x <- data.matrix(x)
+    if (!is.numeric(x)) stop("non-numeric argument to 'pairs'")
+    panel <- match.fun(panel)
+    if((has.lower <- !is.null(lower.panel)) && !missing(lower.panel))
+        lower.panel <- match.fun(lower.panel)
+    if((has.upper <- !is.null(upper.panel)) && !missing(upper.panel))
+        upper.panel <- match.fun(upper.panel)
+    if((has.diag  <- !is.null( diag.panel)) && !missing( diag.panel))
+        diag.panel <- match.fun( diag.panel)
+
+    if(row1attop) {
+        tmp <- lower.panel; lower.panel <- upper.panel; upper.panel <- tmp
+        tmp <- has.lower; has.lower <- has.upper; has.upper <- tmp
+    }
+
+    nc <- ncol(x)
+    if (nc < 2) stop("only one column in the argument to 'pairs'")
+    has.labs <- TRUE
+    if (missing(labels)) {
+        labels <- colnames(x)
+        if (is.null(labels)) labels <- paste("var", 1:nc)
+    }
+    else if(is.null(labels)) has.labs <- FALSE
+    oma <- if("oma" %in% nmdots) dots$oma else NULL
+    main <- if("main" %in% nmdots) dots$main else NULL
+    if (is.null(oma)) {
+        oma <- c(4, 4, 4, 4)
+        if (!is.null(main)) oma[3] <- 6
+    }
+    opar <- par(mfrow = c(nc, nc), mar = rep.int(gap/2, 4), oma = oma)
+    on.exit(par(opar))
+
+    for (i in if(row1attop) 1:nc else nc:1)
+        for (j in 1:nc) {
+            localPlot(x[, j], x[, i], xlab = "", ylab = "",
+                      axes = FALSE, type = "n", ...)
+	    if (i<j) {
+	        abline(v=mcmcML[j],col=2)
+	        abline(h=mcmcML[i],col=2)
+	    }
+            if(i == j || (i < j && has.lower) || (i > j && has.upper) ) {
+                box()
+                if(i == 1  && (!(j %% 2) || !has.upper || !has.lower ))
+                    localAxis(1 + 2*row1attop, ...)
+                if(i == nc && (  j %% 2  || !has.upper || !has.lower ))
+                    localAxis(3 - 2*row1attop, ...)
+                if(j == 1  && (!(i %% 2) || !has.upper || !has.lower ))
+                    localAxis(2, ...)
+                if(j == nc && (  i %% 2  || !has.upper || !has.lower ))
+                    localAxis(4, ...)
+                mfg <- par("mfg")
+                if(i == j) {
+                    if (has.diag) diag.panel(as.vector(x[, i]))
+                    if (has.labs) {
+                        par(usr = c(0, 1, 0, 1))
+                        if(is.null(cex.labels)) {
+                            l.wid <- strwidth(labels, "user")
+                            cex.labels <- max(0.8, min(2, .9 / max(l.wid)))
+                        }
+                        text.panel(0.5, label.pos, labels[i],
+                                   cex = cex.labels, font = font.labels)
+                    }
+                } else if(i < j)
+                    localLowerPanel(as.vector(x[, j]), as.vector(x[, i]), ...)
+                else
+                    localUpperPanel(as.vector(x[, j]), as.vector(x[, i]), ...)
+                if (any(par("mfg") != mfg))
+                    stop("the 'panel' function made a new plot")
+            } else par(new = FALSE)
+
+        }
+    if (!is.null(main)) {
+        font.main <- if("font.main" %in% nmdots) dots$font.main else par("font.main")
+        cex.main <- if("cex.main" %in% nmdots) dots$cex.main else par("cex.main")
+        mtext(main, 3, 3, TRUE, 0.5, cex = cex.main, font = font.main)
+    }
+    invisible(NULL)
+}


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[R-gregmisc-users] SF.net SVN: r-gregmisc: [943] trunk/PathwayModeling/thesispaper/R

[<r_b...@us...>]

Revision: 943
Author:   r_burrows
Date:     2006-03-20 03:49:41 -0800 (Mon, 20 Mar 2006)
ViewCVS:  http://svn.sourceforge.net/r-gregmisc/?rev=943&view=rev

Log Message:
-----------
minor update

Modified Paths:
--------------
    trunk/PathwayModeling/thesispaper/R/plotDensities.R
    trunk/PathwayModeling/thesispaper/R/plotVi.R
Modified: trunk/PathwayModeling/thesispaper/R/plotDensities.R
===================================================================
--- trunk/PathwayModeling/thesispaper/R/plotDensities.R	2006-03-20 11:41:09 UTC (rev 942)
+++ trunk/PathwayModeling/thesispaper/R/plotDensities.R	2006-03-20 11:49:41 UTC (rev 943)
@@ -6,5 +6,6 @@
         hist(output[,i],breaks=30,xlim=c(low,high),xlab=paste("d",i,sep=""),main="")
         par(new=T)
         plot(seq(low,high,step),dnorm(seq(low,high,step),modes[i],SD[i]),type='l',            xlim=c(low,high),xlab="",yaxt='n',ylab="")
+	abline(v=mcmcML[i],col=2)
     }
 }

Modified: trunk/PathwayModeling/thesispaper/R/plotVi.R
===================================================================
--- trunk/PathwayModeling/thesispaper/R/plotVi.R	2006-03-20 11:41:09 UTC (rev 942)
+++ trunk/PathwayModeling/thesispaper/R/plotVi.R	2006-03-20 11:49:41 UTC (rev 943)
@@ -3,7 +3,8 @@
 #    attach(as.data.frame(input))
     times <- input[,1]
     CI <- get95CI(2,input,1,Dsamp)
-    temp <- fitPaper(paramsML$par,input)[[1]]
+    temp <- fitPaper(optimML,input)
+    temp2 <- fitPaper(mcmcML,input)
     low <-min(CI,input[,3])
     high <- max(CI,input[,3])
     plot(times,input[,3],ylim=c(low,high),xlab="time",ylab="v2")
@@ -12,12 +13,13 @@
     par(new=T)
     plot(times,CI[,2],ylim=c(low,high),type="l",lty=2,xlab="",ylab="",col=2)
     par(new=T)
-    plot(times,temp,ylim=c(low,high),type="l",xlab="",ylab="",col=2)
+    plot(times,temp[[1]],ylim=c(low,high),type="l",xlab="",ylab="",col=3)
+    par(new=T)
+    plot(times,temp2[[1]],ylim=c(low,high),type="l",xlab="",ylab="",col=2)
     abline(h=0)
     legend(32,2000,c("data","fitted curve","95% CI"),lty=c(0,1,2),pch=c(1,NA,NA),col=c(1,2,2))
     for (i in 3:5) {
         CI <- get95CI(i,input,1,Dsamp)
-        temp <- fitPaper(paramsML$par,input)[[i-1]]
         low <-min(CI,-400)
         high <- max(CI,800)
         plot(times,input[,i+1],ylim=c(low,high),
@@ -27,7 +29,9 @@
         par(new=T)
         plot(times,CI[,2],ylim=c(low,high),type="l",lty=2,xlab="",ylab="",col=2)
         par(new=T)
-        plot(times,temp,ylim=c(low,high),type="l",xlab="",ylab="",col=2)
+        plot(times,temp[[i-1]],ylim=c(low,high),type="l",xlab="",ylab="",col=3)
+        par(new=T)
+        plot(times,temp2[[i-1]],ylim=c(low,high),type="l",xlab="",ylab="",col=2)
         abline(h=0)
     }
 #    detach(as.data.frame(input))


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[R-gregmisc-users] SF.net SVN: r-gregmisc: [942] trunk/PathwayModeling/thesispaper

[<r_b...@us...>]

Revision: 942
Author:   r_burrows
Date:     2006-03-20 03:41:09 -0800 (Mon, 20 Mar 2006)
ViewCVS:  http://svn.sourceforge.net/r-gregmisc/?rev=942&view=rev

Log Message:
-----------
minor revisions

Modified Paths:
--------------
    trunk/PathwayModeling/thesispaper/introduction.Snw
    trunk/PathwayModeling/thesispaper/paper.Snw
    trunk/PathwayModeling/thesispaper/results.Snw
Modified: trunk/PathwayModeling/thesispaper/introduction.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/introduction.Snw	2006-03-20 01:13:55 UTC (rev 941)
+++ trunk/PathwayModeling/thesispaper/introduction.Snw	2006-03-20 11:41:09 UTC (rev 942)
@@ -12,7 +12,7 @@
 altering the system behavior in beneficial ways.
 
 In order to gain some insight into the properties of biological
-systems, we would like to prepare quantitative models of sequences of
+systems, we would like to use quantitative models of sequences of
 reactions.
  Specifically, we would like to develop methods that can be used to
 deduce the kinetic parameters of the enzymes in a pathway. With this
@@ -58,7 +58,7 @@
 
 We will explore the utility of the the Bayesian approach using a
 simulation technique: First we will generate data for a model with
-known parameters, and then we will fit this data using our Bayesian
+known parameters, and then we will fit this data using Bayesian
 methods.
 The biochemical data is generated using Gillespie's
 stochastic simulation algorithm \cite{Gillespie77}. This algorithm is

Modified: trunk/PathwayModeling/thesispaper/paper.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/paper.Snw	2006-03-20 01:13:55 UTC (rev 941)
+++ trunk/PathwayModeling/thesispaper/paper.Snw	2006-03-20 11:41:09 UTC (rev 942)
@@ -34,16 +34,16 @@
 \address{University of Rhode Island, Kingston, RI}
   
 \begin{abstract}
-The usefulness of Markov chain Monte Carlo methods for the modeling
-of biochemical reactions is examined. With simulated data, it is
-shown that mechanistic models can be fit to sequences of enzymatic
+\textit{We have examined} the usefulness of \textit{Bayesian statistical methods} for the modeling
+of biochemical reactions. With simulated data, it is
+shown that mechanistic models can \textit{effectively} be fit to sequences of enzymatic
 reactions. These methods have the advantages of being relatively easy
 to use and producing probability distributions for the model
-parameters rather than point estimates.
+parameters rather than point estimates, \textit{allowing more informative statistical inferences to be drawn}.
 
 Three Markov chain Monte Carlo algorithms are used to fit models to
 data from a
-sequence of 4 enzymatic reactions. The algorithms 
+sequence of 4 enzymatic reactions. These algorithms 
 are evaluated with respect to
 the goodness-of-fit of the fitted models and the time to completion.
 It is shown that the algorithms produce essentially the same

Modified: trunk/PathwayModeling/thesispaper/results.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/results.Snw	2006-03-20 01:13:55 UTC (rev 941)
+++ trunk/PathwayModeling/thesispaper/results.Snw	2006-03-20 11:41:09 UTC (rev 942)
@@ -9,16 +9,37 @@
 algorithms were used for the simulations and they converged to similar
 distributions. 
 
-The marginal distributions for the parameters were mound-shaped. They
+The marginal distributions for the parameters were mound-shaped 
+(Figure~\ref{densities}). They
 generally were skewed to the right and had thicker
-tails than normal distributions as illustrated in
-Figure~\ref{densities}.
+tails than normal distributions. /textit{The modes of the distributions 
+turned out to be poor estimates of the maximum likelihood values for 
+the parameters. The maximum likelihood values were all larger than 
+the modes, and the difference increased as the position of the reaction 
+in the sequence increased. Correlated pairs of parameters 
+(e.g., $d_1$ --$d_2$ and $d_3$ --$d_4$, Figure~\ref{scatter}) were similar relative distances from their 
+modes. The significance of all this is not clear at this point.}
 
-<<fig4,echo=F,eval=T>>=3
+<<fig4,echo=F,eval=F>>=3
+source("R/paperSSQ.R")
+source("R/getModes.R")
+source("R/do.optim.R")
 source("R/plotDensities.R")
-source("R/getModes.R")
 output16 <- read.table("data/output16.AllComp.thinned")
-modes <- getModes(output16)
+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)
@@ -26,6 +47,7 @@
 plotDensities(output16)
 dev.off()
 par(mfrow=c(1,1))
+detach(input16)
 @
 \begin{figure}
   \centering
@@ -34,13 +56,14 @@
   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.}
+  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=T>>=4
+<<fig5,echo=F,eval=F>>=4
 source("R/plotDensity.R")
 source("R/plotConverged.R")
 output12 <- read.table("data/output12.AllComp.thinned")
@@ -55,7 +78,7 @@
 \begin{figure}
   \centering
   \includegraphics[scale=0.9]{figures/tempDir/converged}
-  \caption{Parameter distributions as functions of the number of data
+  \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
@@ -63,12 +86,17 @@
   \label{converged}
 \end{figure}
 There is some narrowing of the distributions as the number of points
-increases from 16 to 25 but it is not pronounced. 
+increases from 16 to 25 but it is not pronounced. \textit{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.
 
-<<fig6, echo=F,eval=T>>=5
+<<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))
@@ -78,7 +106,7 @@
     }
 pdf(file="figures/tempDir/scatterPlot.pdf",width=8,height=8)
 par(pch='.')
-pairs(output16[sample(1:20000,2000),],labels=c('d1','d2','d3','d4','d5','d6','d7','d8','d9'),lower.panel=mcmc.cor)
+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)
 @
@@ -99,25 +127,19 @@
 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=T>>=6
+<<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)
-temp<-output16[sample(1:20000,5000),]
-temp2<-numeric(5000)
-for (i in 1:5000) temp2[i]<-paperSSQ(temp[i,])
-temp<-temp[order(temp2),]
-Dsamp<-as.matrix(temp[1:4750,])
-modes<-getModes(output16)
-temp<-do.optim(start=modes[-10],scail=modes[-10])
-temp<-do.optim(start=temp$par,scail=temp$par)
-paramsML<-temp
+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()
@@ -154,13 +176,17 @@
 
 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.
-For the 1-component Metropolis algorithm, one parameter is
-updated for each evaluation of the likelihood method whereas the
-all-component Metropolis and NKC algorithms update all the parameters
-for each likelihood evaluation.
+likelihood evaluations for the three algorithms. \textit{The relatively 
+long convergence time for the 1-component Metropolis algorithm is a result 
+of the high correlations between some of the parameters (Figure~\ref{scatter}).
+With the 1-component Metropolis algorithm, one parameter is
+updated for each evaluation of the likelihood method, i.e., the algorithm 
+searches the parameter space by moving parallel to the axes and thus 
+it can sample a density with a diagonal axis of symmetry very slowly. 
+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=T>>=7
+<<fig8, echo=F,eval=F>>=7
 source("R/getMSDvsEvals.R")
 source("R/plotMSD16.R")
 source('R/paperSSQ.R')
@@ -172,7 +198,7 @@
 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)
+pdf(file="figures/tempDir/MSR16.pdf",width=9,height=6)
 plotMSD16()
 dev.off()
 detach(input16)


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[R-gregmisc-users] SF.net SVN: r-gregmisc: [941] trunk/PathwayModeling/thesispaper/methods.Snw

[<r_b...@us...>]

Revision: 941
Author:   r_burrows
Date:     2006-03-19 17:13:55 -0800 (Sun, 19 Mar 2006)
ViewCVS:  http://svn.sourceforge.net/r-gregmisc/?rev=941&view=rev

Log Message:
-----------
minor edits

Modified Paths:
--------------
    trunk/PathwayModeling/thesispaper/methods.Snw
Modified: trunk/PathwayModeling/thesispaper/methods.Snw
===================================================================
--- trunk/PathwayModeling/thesispaper/methods.Snw	2006-03-14 18:00:18 UTC (rev 940)
+++ trunk/PathwayModeling/thesispaper/methods.Snw	2006-03-20 01:13:55 UTC (rev 941)
@@ -1,4 +1,4 @@
-\section{METHODS}
+\section{Methods}
 
 \begin{verbatim}
 $Id$
@@ -59,7 +59,7 @@
     the perturbation of $R1$ at $time = 20$ results in the time
     courses plotted in Figure~\ref{pulse}.
 
-<<echo=F,eval=T>>=2
+<<echo=F,eval=F>>=2
 source("R/pulse.R")
 rawdata <- read.table("data/rawdata.dat",header=T)
 attach(rawdata)
@@ -215,7 +215,7 @@
       v_5 &\sim N\left(\frac{d_7R4}{d_8+R4} - d_9R5, \;\; \sigma^2\right)
     \end{align*}
     In this case there are 10 parameters to be estimated: the 9
-    coefficients $d_1$ -- $d_5$ and $\sigma^2$.
+    coefficients $d_1$ -- $d_9$ and $\sigma^2$.
 
     \subsubsection{MCMC sampling algorithms}
 


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[R-gregmisc-users] SF.net SVN: r-gregmisc: [940] trunk/gdata

[<nj...@us...>]

Revision: 940
Author:   nj7w
Date:     2006-03-14 10:00:18 -0800 (Tue, 14 Mar 2006)
ViewCVS:  http://svn.sourceforge.net/r-gregmisc/?rev=940&view=rev

Log Message:
-----------
Fixed R CMD check errors and added trim.default to NAMESPACE

Modified Paths:
--------------
    trunk/gdata/NAMESPACE
    trunk/gdata/R/trim.R
    trunk/gdata/man/resample.Rd
Modified: trunk/gdata/NAMESPACE
===================================================================
--- trunk/gdata/NAMESPACE	2006-03-13 19:50:48 UTC (rev 939)
+++ trunk/gdata/NAMESPACE	2006-03-14 18:00:18 UTC (rev 940)
@@ -33,5 +33,8 @@
 S3method(nobs,data.frame)
 S3method(nobs,default)
 S3method(nobs,lm)
+S3method(trim, character)
+S3method(trim, default)
+S3method(trim, factor)
 S3method(reorder,factor)
 

Modified: trunk/gdata/R/trim.R
===================================================================
--- trunk/gdata/R/trim.R	2006-03-13 19:50:48 UTC (rev 939)
+++ trunk/gdata/R/trim.R	2006-03-14 18:00:18 UTC (rev 940)
@@ -1,6 +1,9 @@
 # $Id$
 
 trim <- function(s)
+  UseMethod("trim",s)
+
+trim.default <- function(s)
   {
     s <- sub("^ +","",s)
     s <- sub(" +$","",s)
@@ -9,12 +12,12 @@
 
 trim.character <- function(s)
 {
-  return(trim(s))
+  return(trim.default(s))
 }
 
 trim.factor <- function(s)
 {
-  levels(s) <- trim(levels(s))
+  levels(s) <- trim.default(levels(s))
   return(s)
 }
 

Modified: trunk/gdata/man/resample.Rd
===================================================================
--- trunk/gdata/man/resample.Rd	2006-03-13 19:50:48 UTC (rev 939)
+++ trunk/gdata/man/resample.Rd	2006-03-14 18:00:18 UTC (rev 940)
@@ -57,5 +57,4 @@
       sample(x, n, ...)
   }
 }
-\keyword{ ~kwd1 }% at least one, from doc/KEYWORDS
-\keyword{ ~kwd2 }% __ONLY ONE__ keyword per line
+\keyword{misc}% at least one, from doc/KEYWORDS


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.