User Manual for fraping
Gabriela Itzetl Medina-Ruíz 1,3 , Arturo Itzcóatl Medina-Ruiz 2 and Julio Morán 1 * 1. División de Neurociencias, Instituto de Fisiología Celular, Universidad Nacional Autónoma de México, Coyoacán, 04510, Ciudad de México 2. Facultad de Ciencias, Universidad Nacional Autónoma de México, Coyoacán, 04510, Ciudad de México, México 3. Posgrado en Ciencias Biológicas, Unidad de Posgrado, Edificio D, 1° Piso, Circuito de Posgrados, Ciudad Universitaria, Coyoacán, 04510, Ciudad de México, México
June 12, 2023
Introduction
Fluorescence Recovery After Photobleaching (FRAP), a widely used technique in biology and biomedicine, provides a unique opportunity to investigate the mobility and dynamics of molecules in living cells. It offers valuable insights into the recovery of fluorophore-labeled molecules after photobleaching.
This manual presents an advanced methodology for analyzing fluorescence recovery in the FRAP technique. Its main objective is to provide researchers with a comprehensive tool to understand and quantify molecule dynamics in living cells using FRAP.
The content of this manual covers the following key points:
Fundamentals of FRAP technique: A detailed description of the fundamental principles of FRAP technique is provided, including basic concepts of photobleaching, fluorescence recovery, and the use of fluorophores in FRAP experiments.
Models of fluorescence recovery: Different mathematical and statistical models are introduced to describe fluorescence recovery after photobleaching. The theoretical foundations of each model are explained, and examples of their application in FRAP data analysis are provided.
Curve fitting and parameter estimation: A methodology for fitting fluorescence recovery curves to experimental data obtained through FRAP is presented. Optimization techniques are discussed, and practical tools for estimating recovery model parameters are provided.
Results analysis and quality assessment: The assessment of curve fitting quality for fluorescence recovery curves is addressed using statistical methods and result visualization. Criteria for determining goodness of fit are discussed, and graphical tools for analyzing and comparing results are presented.
Data simulation and model validation: The process of stochastic simulation of fluorescence recovery curves based on estimated parameters is explained. The importance of simulation in FRAP analysis is discussed, and tools for validating the recovery model through comparison with real data are provided.
In summary, this manual offers a comprehensive and practical guide for analyzing fluorescence recovery in the FRAP technique. Researchers will find a valuable reference in this manual for understanding theoretical foundations, applying statistical models, curve fitting, result evaluation, and conducting simulations that enable rigorous and quantitative analysis of data obtained through FRAP.
Get start
In the fraping
package, the necessary functions for FRAP data analysis are implemented.
However, to use it properly, it is also necessary to install the itz package. Both packages
are repositories on GitHub, and their
installation can be done using the install_github function
from the devtools
package in R.
devtools::install_github("artitzco/itz")
devtools::install_github("artitzco/fraping")This way, the libraries will be installed on your system, so every time you start a new project for FRAP data analysis, you only need to import the corresponding libraries:
library(itz)
library(fraping)Getting Started (optional)
If you want to start using the fraping package right away, you can
download a sample dataset using the downloadData function.
This function will download a file with FRAP data into the “data” folder
within your working directory, which you can use as a starting point to
practice the analysis methods described in this manual. Please note that
this download is optional, and you can use your own FRAP data if you
have it available.
To use the sample data, run the following code in the R command
console: fraping::downloadData(). After executing this
function, in addition to the “data” folder, you will find a test script
called “test.fraping.R” in the root directory. You can use this script
as a starting point to explore and analyze FRAP data.
The newDataFrap function creates a base for a new FRAP
data group. The function synthesizes the information and returns a list
that allows the user to manipulate the data more easily. The supported
parameters are the group name and the photobleaching intensity in
percentage. For more details, refer to newDataFrap.
To avoid future issues, the group name should follow the following recommendations:
- It should not contain spaces.
- It should start with a letter.
- Only the special characters “.” and “_” are allowed.
A <- newDataFrap("Neurites", 80)
B <- newDataFrap("Secondary.Neurites", 80)
C <- newDataFrap("Dendrites", 80)Remember to follow these recommendations when creating new FRAP data groups. This will ensure consistency and proper handling of the data in future analyses.
Import data FRAP
Once the data groups are initialized, it’s time to add the
information. One way to do this is through the addDir
function of the newDataFrap object. This function requires
the user to organize the FRAP tables into three files on their computer:
the first one with the recovery data, the second one with the control
tables, and the third one with the background information.
The addDir function in the fraping package accepts the
following parameters to add the information from the FRAP data
files:
seconds: Duration of the experiment in seconds.photo: Frame index where photobleaching is applied.areaColumn: Name of the column containing the area data.inputColumn: Name of the column containing the intensity data.dirR: Path to the recovery data file.dirC: Path to the control data file.dirB: Path to the background data file (optional).
It’s important to note that the background data is optional. If you decide to omit the background data file, you should consider that this will have a direct impact on the data standardization during the analysis.
Here’s an example of how to use the addDir function:
A$addDir(480, 4, "Area", "IntDen",
"data/GFP/GFP Secondary Neurite",
"data/GFP/GFP Secondary Neurite Control",
"data/GFP/GFP Secondary Neurite Back")
B$addDir(480, 4, "Area", "IntDen",
"data/LIFEACT GFP/LIFEACT GFP Secondary Neurite",
"data/LIFEACT GFP/LIFEACT GFP Secondary Neurite Control",
"data/LIFEACT GFP/LIFEACT GFP Secondary Neurite Back")
C$addDir(480, 4, "Area", "IntDen",
"data/LIFEACT GFP/LIFEACT GFP Dendrite",
"data/LIFEACT GFP/LIFEACT GFP Dendrite Control",
"data/LIFEACT GFP/LIFEACT GFP Dendrite Back")When creating a new data group using the newDataFrap
function, a random color is automatically assigned to the group. This
color is used to differentiate it from others in the plots. If you want
to assign a specific color to the group, you can use the
setColor function of the newDataFrap object.
This function accepts the color value in hex format as a
parameter, or you can choose a predefined color from the
grDevices::colors() vector or use the color
function from the itz
package.
A$setColor("#05753D")
B$setColor("#0000A6")
C$setColor("#E014F7")Before performing a parametric fitting, you might prefer to visualize
the data through plots for visual inspection. To achieve this, you can
use the plotRecover function, which takes the bases created
with newDataFrap as parameters. You can customize the plots
using various parameters. For more details about this function, refer to
plotRecover.
plotRecover(A, B, C, plot.lines=T)
plotRecover(A, B, C, plot.shadow=T, plot.mean=T)```{r
out.width=“75%”, fig.align=‘center’, echo=FALSE} graph(“img/002.jpg”)#>>><<<
```r
plotRecover(A, B, C, plot.points=T)
Choose a probability model
The proposed parametric model to describe the behavior of the fluorescence recovery curve from the photobleaching moment is as follows:
\[F(t)=f_{min}+\alpha\,(1-f_{min})\,p(t-t_0|\Theta),\]
where \(p\) is a cumulative distribution function parameterized by \(\Theta\), \(t_0\) is the photobleaching time, and \(f_{min}\) is the minimum fluorescence value after photobleaching. The theoretical maximum value of fluorescence after photobleaching, \(f_{max}\), is defined as:
\[f_{max}= f_{min}+\alpha(1-f_{min}).\]
It can be observed that the function \(F(t)\) is defined in the time interval \(t\in[t_0, \infty)\). Similarly, the function \(F^{^{AB}}(t)\) is defined to represent the fluorescence recovery curve after photobleaching, and it is defined in the time interval \(t\in[0, \infty)\):
\[F^{^{AB}}(t)=f_{min}+\alpha\,(1-f_{min})\,p(t|\Theta).\]
More information about the construction of the model can be found in
the appendix of this document. In order for the user
to perform the parametric fitting for subsequent data analysis, the user
only needs to choose the most appropriate probability model, \(p\). The newFit function
allows the declaration of a new probability model. The function takes
parameters such as the name of the model, the cumulative distribution
function, the names of the probability function parameters, and a list
with the value ranges for those parameters. For more details, refer to
newFit.
Exp<-newFit("Exponential", pexp,"rate", list(c(0,1)))
Gam<-newFit("Gamma", pgamma, c("shape","rate"), list(c(0,5), c(0,5)))
Wei<-newFit("Weibull", pweibull, c("shape","scale"), list(c(0,5),c(0,2000)))One of the most important criteria for choosing between different probability models is the mean squared error. Models with lower error should be preferred. The mean squared error is given by the following expression:
\[E_i =\frac{1}{m_i}\sum_{j=1}^{m_i} [f_i(t_j)-F(t_j)]^2,\]
where \(\{f_1,..,f_n\}\) is a group
of fluorescence data and \(\{t_1,t_2,\cdots,t_{m_i}\}\) are all the
values greater than the photobleaching time quantified for \(f_i\). The compareFit function
compares the mean squared error of different probability models fitted
to a group of data by creating a Q-Q plot and the histogram of the error
for each model. The function takes parameters such as the data group to
be compared, a list of the models to be compared, and several
customizable parameters for the plot. For more details, refer to compareFit.
compareFit(B, fit=l(Exp, Gam, Wei), col.lines=c("red","yellow"), lty.lines=c(2,1), lwd.lines=2, lwd.mean=2)
compareFit(C, fit=l(Exp, Gam, Wei), col
.lines=c("red","yellow"), lty.lines=c(2,1), lwd.lines=2, lwd.mean=2)
For more information about choosing a good probability model, refer to the appendix.
Fit a parametric model
For this particular case, the Weibull model yields better results.
Therefore, the Wei function is used to perform the
corresponding fitting of the data. The Wei function,
applied to a specific data group, returns a list with the information
obtained from the parametric fitting, including the estimation of the
fluorescence recovery function \(F(t)\)
and other relevant data. For more details, refer to newFit.
To plot the parametric fitting of the fluorescence recovery, the
plotFit function is used. It takes as parameters the data
groups to be fitted and the probability model chosen for each data
group. There are also optional parameters that allow customization and
formatting of each plot. For more details, refer to plotFit. If the chosen probability model is the same
for all data groups, it is sufficient to specify it once as shown
below:
plotFit(B, C, fit=Wei, plot.lines=T)
plotFit(B, C, fit=Wei, plot.shadow=T, plot.mean=T, alp.shadow=0.2, ylim=c(0.2, 1), xdigits=0)
Using the plotFit function, future estimates of the
fluorescence recovery curve can be plotted by modifying the range of the
graph. The plots obtained from the plotRecover and
plotFit functions can be stacked by changing the optional
parameter new.plot to FALSE.
plotFit(B, fit=Wei, plot.shadow=T, alp.shadow=0.3, xlim=c(0,500), ylim=c(0.2, 1))
plotRecover(B, AB=T, plot.lines=T, new.plot=F)
Another use for the plotFit function is to visually
compare two or more probability models on the same data group.
plotRecover(B, plot.shadow =T, ylim=c(0.2, 1.2), xdigits=0,lwd.border=1, col="#808080")
plotFit(B, B, fit=l(Wei, Exp), plot.lines=T, new.plot=F, displacement=T, col=c("blue2","red2"))
Compare data sets
Conventional Analysis
You can compare two groups of data based on the variables associated
with the parametric model using the compareParam function.
This function takes parameters such as the two groups of data, the
probability model that fits each group of data, and the variable to be
compared. You can also customize the graph by specifying additional
values. For more details, refer to compareParam. You can obtain the names of the
variables as follows:
colnames(Wei(B)$table)## [1] "shape" "scale" "alpha" "fmax" "UF" "MF" "error"The compareParam function generates a Q-Q plot of the
chosen variable. By setting the parameter new.plot to
FALSE, you can stack the graphs. You can obtain the details
of the t-student test for the difference of means with a 95% confidence
level by setting return=TRUE (the graph will be omitted).
The confidence level and alternative hypothesis can be modified using
the alternative and conf.level parameters,
respectively. For more details, refer to compareParam.
compareParam(B, C, fit=Wei, param="shape")
This allows for a conventional analysis of fluorescence recovery, where the immobile fraction (UF) and mobile fraction (MF) are compared. Recall that:
\[UF=\frac{1-f_{max}}{1-f_{min}}\] \[MF=\frac{f_{max}-f_{min}}{1-f_{min}}\] \[UF+MF=1\]
compareParam(B, C, fit=Wei, param="UF", xlim=c(0,1), ylim=c(0,1),lty.lines=c(2,1), col.lines=c("white","red"), lwd.lines=2)
compareParam(B, C, fit=Wei, param="MF", new.plot=F, lwd.lines=2)
You can also use the compareParam function to compare
two probability models fitted to the same group of data based on the
error produced by each fit. Note that this graph is already included in
the mosaic produced by the compareFit function seen previously.
compareParam(B, B, fit=l(Exp,Wei), param="error", col.lines=c("red", "yellow"), lwd.lines=2, lty.lines=c(2,1), ydigits=2, xdigits=2)
compareParam(B, B, fit=l(Exp,Wei), param="error", return=T)#>>><<<##
## Welch Two Sample t-test
##
## data: error_Exponential_Secondary.Neurites and error_Weibull_Secondary.Neurites
## t = 1.2307, df = 21.424, p-value = 0.2318
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.09746458 0.38089484
## sample estimates:
## mean of x mean of y
## 0.7222331 0.5805180Mean Behavior Analysis
The analysis of mean behavior involves evaluating the average recovery of fluorescence over time. The mean curve is calculated as follows:
\[\overline F(t)=\frac{1}{n}\sum_{i=1}^{n}F^{^{AB}}_i(t),\]
where \(n\) is the size of the data set, and \(F^{AB}_i\) represents the recovery
curve of fluorescence after photobleaching associated with the quantification of fluorescence \(f_i\).
plotFit(B, fit=Wei, plot.lines=T, plot.mean=T, col.mean="purple",lwd.mean=2, xdigits=0, ydigits=2, alp.lines=0.5)
The compareMean function allows you to compare the mean
curves of two groups of data using a continuous t-student test for the
difference of means over a period of time. The function takes parameters
such as the groups of data to be compared, the probability model fitted
to the data, and additional values to customize the graph.
The resulting graph represents the p-value of the t-student test,
with the reference line indicating the significance level. By setting
return=TRUE, you can obtain the details of the test stored
in a list (without generating the graph). The confidence level and
alternative hypothesis can be modified using the
alternative and conf.level parameters,
respectively. For more details, refer to compareMean.
compareMean(B, C, fit=Wei, lty.lines=c(2,1), col.lines=c("red","purple"), ydigits=2)
By setting p.value=FALSE, you can calculate the
difference of means between the two groups of data, with the reference
line located at zero. In this case, by setting return=TRUE,
you obtain a vector with the difference of means between the two
curves.
compareMean(B, C, fit=Wei, p.value=F, lty.lines=c(2,1), col.lines=c("red","magenta"), ydigits=3)
Simulation
The implementation of simulation methods for fluorescence recovery analysis is one of the main contributions of this work. Simulation allows you to increase the sample size without altering its central tendency and dispersion measures. This is especially useful in cases where the sample size is not sufficiently large. Remember that hypothesis tests that rely on limit processes are particularly sensitive to sample size. Additionally, the simulation process is informed by an appropriate choice of the probability model, as explained later.
The simulation of the fluorescence recovery model \(F(t)\) is achieved by simulating the
parameters \(\Theta\) of the
probability model \(p(t|\Theta)\).
Refer to the appendix for more information. The
simFit function allows you to simulate data fitted to a
probability model. It takes a group of data fitted to a probability
model as a mandatory parameter, and optional parameters such as the size
of the simulated sample and a simulation seed. For more details about
the simulation seed, refer to the appendix.
For the case of the Weibull fit, the data simulation is performed as
follows: simFit(Wei(B)) and simFit(Wei(C)).
For more specific information about the data generated by this function,
refer to simFit.
Furthermore, the plotFit, compareParam, and
compareMean functions have the simulation method
implemented. To access this feature, you need to change the
simulated parameter to TRUE. We can replicate
the previously made graphs now with simulated data:
plotFit(B, C, fit=Wei, simulated=T, plot.lines=T)
plotFit(B, C, fit=Wei, simulated=T, plot.shadow=T, plot.mean=T, alp.shadow=0.2, ylim=c(0.2, 1), xdigits=0)
plotFit(B, fit=Wei, simulated=T, plot.shadow=T, alp.shadow=0.3, xlim=c(0,500), ylim=c(0.2, 1))
plotRecover(B, AB=T, plot.lines=T, new.plot=F)
plotRecover(B, AB=T, plot.shadow=T,ylim=c(0.2, 1.2), xdigits=0,lwd.border=1, col="#808080")
plotFit(B, B, fit=l(Wei, Exp), simulated=T, plot.lines=T, new.plot=F, displacement=T, col=c("blue2","red2"))
compareParam(B, C, fit=Wei, simulated=T, param="shape")
compareParam(B, C, fit=Wei, param="UF", simulated=T, xlim=c(0,1), ylim=c(0,1),
lty.lines=c(2,1), col.lines=c("white","red"), lwd.lines=2)
compareParam(B, C, fit=Wei, param="MF", simulated=T, new.plot=F, lwd.lines=2)
compareParam(B, C, fit=Wei, simulated=T, param="UF", return=T)##
## Welch Two Sample t-test
##
## data: UF_WeibullSim_Secondary.Neurites and UF_WeibullSim_Dendrites
## t = -1.4645, df = 97.842, p-value = 0.1463
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.0891179 0.0134351
## sample estimates:
## mean of x mean of y
## 0.07055279 0.10839419compareParam(B, C, fit=Wei, simulated=T, param="MF", return=T)##
## Welch Two Sample t-test
##
## data: MF_WeibullSim_Secondary.Neurites and MF_WeibullSim_Dendrites
## t = 1.4645, df = 97.842, p-value = 0.1463
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.0134351 0.0891179
## sample estimates:
## mean of x mean of y
## 0.9294472 0.8916058plotFit(B, fit=Wei, simulated=T, plot.lines=T, plot.mean=T, col.mean="purple",lwd.mean=2, xdigits=0, ydigits=2, alp.lines=0.5)
compareMean(B, C, fit=Wei, simulated=T, lty.lines=c(2,1), col.lines=c("red","purple"), ydigits=2)
compareMean(B, C, fit=Wei, simulated=T, p.value=F, lty.lines=c(2,1), col.lines=c("red","magenta"), ydigits=3)
Appendix
The Fluorescence Recovery Model
The quantification of fluorescence in the FRAP technique follows a specific recovery pattern. Initially, the fluorescence remains constant, but at the time of photobleaching (also known as \(t_0\)), the fluorescence intensity decreases to a minimum value called \(f_{min}\). From that point on, a gradual process of fluorescence recovery begins, characterized by two fundamental conditions: firstly, the recovery does not decrease as time progresses, and secondly, the fluorescence reaches a maximum limit (\(f_{max}\)) that cannot exceed the initial value before photobleaching. This evolution is graphically represented for a better understanding of the process.
In this sense, the proposed model for fluorescence recovery after photobleaching must comply with the aforementioned characteristics. The general model for fluorescence recovery after photobleaching is as follows:
\[F(t)=f_{min}+\alpha\,(1-f_{min})\,p(t-t_0|\Theta)\]
where \(p\) is a continuous cumulative distribution function (CDF) parametrized by \(\Theta\).
Observations:-
The function \(F(t)\) is defined in the interval \([t_0,\infty)\).
-
The value \(f_{max}\) is given by \[f_{max} =\lim_{t\to\infty}F(t)\]
-
The value \(\alpha\) is bounded by the values \[\alpha\in[0,1]\]
-
The cumulative distribution function \(p\) must be associated with a random variable with support in \((0,\infty)\). Remember that the CDF is defined as follows:
\[p(t|\Theta)=\int_0^td(x|\Theta)\,\mbox{d}x\] where \(d\) is a parametric density function.
Similarly, the function \(F^{^{AB}}(t)\) is defined to represent the fluorescence recovery curve after photobleaching, which is defined in the time interval \(t\in[0, \infty)\):
\[F^{^{AB}}(t)=f_{min}+\alpha\,(1-f_{min})\,p(t|\Theta).\]
Let’s prove that \(F\) satisfies the characteristics of fluorescence recovery after photobleaching:
- \(F(t_0)=f_{min}\)
- \(F(t)\) is non-decreasing
- \(f_{max}=f_{min}+\alpha\,(1-f_{min})\) and \(f_{max}\in[f_{min}, 1]\)
Evaluating:
\[\begin{aligned} F(t_0) &=f_{min}+\alpha\,(1-f_{min})\,p(t_0-t_0|\Theta)\\ &=f_{min}+\alpha\,(1-f_{min})\,p(0|\Theta)\\ &=f_{min}+\alpha\,(1-f_{min})\int_0^0d(x|\Theta)\,\mbox{d}x\\ &=f_{min}+\alpha\,(1-f_{min})\,(0)\\ &=f_{min} \end{aligned}\] \[\therefore F(t_0)=f_{min}\]
Consider \(t_1,t_2\in\mathbb R^+\) such that \(t_1\leq t_2\), then
$$ \[\begin{aligned} t_1-t_0\leq &\;t_2-t_0\\ \Rightarrow p(t_1-t_0|\Theta )\leq &\;p(t_2-t_0|\Theta)\\ f_{min}+\alpha\,(1-f_{min})\,p(t_1-t_0|\Theta)\leq &\;f_{min}+\alpha\,(1-f_{min})\,p(t_2-t_0|\Theta) \end{aligned}\]\[ \]F(t_1)F(t_2)$$
The above proves that the parametric model \(F\) is non-decreasing over time.
We calculate the limit of \(F(t)\) as \(t\) tends to infinity:
\[\begin{aligned} f_{max}=\lim_{t\to\infty}F(t) &=\lim_{t\to\infty}f_{min}+\alpha\,(1-f_{min})\,p(t-t_0|\Theta)\\ &=f_{min}+\alpha\,(1-f_{min})\lim_{t\to\infty}\,p(t-t_0|\Theta)\\ &=f_{min}+\alpha\,(1-f_{min})\lim_{t\to\infty}\int_0^{t-t_0}d(x|\Theta)\,\mbox{d}x\\ &=f_{min}+\alpha\,(1-f_{min})\int_0^{\infty}d(x|\Theta)\,\mbox{d}x\\ &= f_{min}+\alpha\,(1-f_{min})\,(1) \end{aligned}\]
\[\therefore f_{max}=f_{min}+\alpha\,(1-f_{min})\] On the other hand, we have \(\alpha\in[0,1]\), so
\[\begin{aligned} &0\leq \alpha\leq 1\\ \Rightarrow &0\leq \alpha\,(1-f_{min})\leq 1-f_{min}\\ \Rightarrow &f_{min}\leq f_{min}+\alpha\,(1-f_{min})\leq f_{min}+1-f_{min}\\ \Rightarrow &f_{min}\leq f_{min}+\alpha\,(1-f_{min})\leq 1\\ \Rightarrow &f_{min}\leq f_{max}\leq 1\\ \end{aligned}\] \[\therefore f_{max}\in[f_{min},1]\]Given a set of fluorescence data \(f\) and a cumulative distribution function \(p(\bullet|\Theta)\), the parameters \((\alpha,\Theta)\) are estimated using the least squares method, which involves minimizing the mean square error between the curves \(f\) and \(F\). The mean square error function is as follows:
\[\begin{aligned}E(\alpha,\Theta)&=\sum_{j=1}^m\left[f(t_j)-F(t_j)\right]^2\\&=\sum_{j=1}^m\left[f(t_j)-f_{min}-\alpha\,(1-f_{min})\,p(t_j-t_0|\Theta)\right]^2\end{aligned}\] where \(\{t_1,t_2,\cdots,t_m\}\) are all the times greater than \(t_0\) quantified for \(f\).
Choosing a Probability Model
The cumulative distribution function \(p(t|\Theta)\) forms the basis for
constructing the parametric model of fluorescence recovery \(F(t|\Theta)\). The choice of the function
\(p\) is at the discretion of the user,
but it is essential that \(p\) has
support in \((0,\infty)\). Another
aspect that the user must take care of is assigning the intervals where
the parameters \(\Theta\) take values.
This can be done using the newFit function.
Example:
Suppose the user has decided to use a Pareto probability function, then \[p(t|\,\mbox{shape},\mbox{scale})=1-\left(\frac{\mbox{scale}}{t+\mbox{scale}}\right)^{\mbox{shape}}\] \[t>0,\;\;\mbox{scale}>0\;\;\&\;\;\mbox{shape}>0\] We implement the Pareto cumulative distribution function and declare its associated parametric fit:
ppareto <- function(t, shape, scale){
ppar <- NULL
for (t in t)
ppar <- c(ppar, 1-(scale/(t+scale))^shape)
ppar
}
Par <- newFit("Pareto", ppareto, c("shape","scale"), list(c(0,a), c(0,b)))The above code sets \(\mbox{shape}\in(0,a)\) and \(\mbox{scale}\in(0,b)\). After fitting the
model Par to a data set, the user needs to verify that the
shape and scale values are fully contained
within the intervals \((0,a)\) and
\((0,b)\), respectively. Otherwise, the
user will have to modify the values of \(a\) and/or \(b\) and repeat the process. This can be
achieved by visually inspecting the variables returned by the fit: the
codes Par(B)$table$shape and
Par(B)$table$scale should be executed and the results
inspected, where B is a FRAP data set obtained using the
newDataFrap function.
The functions created by newFit aim to find the best
values of \(\Theta\) that minimize the
mean square error between the model \(F(t|\Theta)\) and the actual fluorescence
recovery data. Therefore, the main criterion for choosing the most
suitable probability model \(p\) should
be to prefer the model that produces the lowest mean square error. This
can be achieved, as seen earlier, using the compareFit
function.
compareFit(B, fit=l(Exp, Gam, Wei, Par), col.lines=c("red","yellow"), lty.lines=c(2,1), lwd.lines=2, lwd.mean=2)
On the other hand, a second criterion for evaluating the viability of a probability model \(p\) is to examine how the simulations of \(F\) behave. If there are discrepancies between the actual fluorescence data and the simulations of \(F\), then it is reasonable to discard the \(p\) model. This would imply that there is no causal relationship between the parameters \(\Theta\) and the fluorescence recovery data. As will be seen in the next appendix, simulations of \(f_{max}\) are particularly sensitive to variations in \(\Theta\).
plotFit(B, B, fit=l(NULL,Exp), col=c("#
916F91","#FF4D66"), plot.lines=c(F,T), plot.mean=c(T,F), lwd.mean=2, plot.shadow=c(T,F), lwd.border=1, alp.border=0.4, alp.shadow=0.3, ylim=c(0.2, 1), xdigits=0, simulated=T, seed=1992, main="Exponential simulation", cex.main=1.5)
plotFit(B, B, fit=l(NULL,Gam), col=c("#916F91","#4DB34D"), plot.lines=c(F,T), plot.mean=c(T,F), lwd.mean=2, plot.shadow=c(T,F), lwd.border=1, alp.border=0.4, alp.shadow=0.3, ylim=c(0.2, 1), xdigits=0, simulated=T, seed=1992, main="Gamma simulation", cex.main=1.5)
plotFit(B, B, fit=l(NULL,Wei), col=c("#916F91","#664DFF"), plot.lines=c(F,T), plot.mean=c(T,F), lwd.mean=2, plot.shadow=c(T,F), lwd.border=1, alp.border=0.4, alp.shadow=0.3, ylim=c(0.2, 1), xdigits=0, simulated=T, seed=1992, main="Weibull simulation", cex.main=1.5)
plotFit(B, B, fit=l(NULL,Par), col=c("#916F91","#B34DB3"), plot.lines=c(F,T), plot.mean=c(T,F), lwd.mean=2, plot.shadow=c(T,F), lwd.border=1, alp.border=0.4, alp.shadow=0.3, ylim=c(0.2, 1), xdigits=0, simulated=T, seed=1992, main="Pareto simulation", cex.main=1.5)
Theory of Stochastic Simulation
Stochastic simulation is a process by which the behavior of a random variable is reproduced. The basis of simulation follows from the following theorem:
Let \(X\) be a random variable with cumulative distribution function \(\mathbb F_X\), and let \(U\) be a random variable with distribution \(\mbox{Uniform}(0,1)\). Then the random variable \(\mathbb F^{-1}_X(U)\) has the same distribution as \(X\), i.e.,
\[\mathbb F_X = \mathbb F_{\mathbb F^{-1}_X(U)}\]
Proof:
-
$$
\[\begin{aligned}\mathbb F_{\mathbb
F^{-1}_X(U)}
&=\mathbb P\left[\mathbb F^{-1}_X(U)\leq t\right]
=\mathbb P\left[\mathbb F_X\left(\mathbb F^{-1}_X(U)\right)\leq \mathbb
F_X\left(t
\right)\right]\\
&=\mathbb P\left[U\leq \mathbb F_X\left(t
\right)\right]
=\int_0^{\mathbb F_X\left(t
\right)}f_U(u)\,\mbox{d}u
=\int_0^{\mathbb F_X\left(t
\right)}1\,\mbox{d}u\\
&=\left.u\right|_0^{\mathbb F_X
\left(t
\right)}={\mathbb F_X\left(t
\right)}-0={\mathbb F_X\left(t
\right)}
\end{aligned}\]
\[ \]=F_{F^{-1}_X(U)}$$
Note that since \(U\) is a uniform variable, \(f_U(u)=1,\;\; 0<u<1\). (\(f_U\) is the probability density function of \(U\))
The above theorem allows us to simulate any random variable \(X\) as long as its cumulative distribution function \(F_X\) is known. Let \(\{u_1,u_2,\cdots,u_n\}\) be a random sample of observed values from the uniform variable \(\mbox{Uniform}(0,1)\). Then \(\{\mathbb F_X^{-1}(u_1),\mathbb F_X^{-1}(u_2),\cdots,\mathbb F_X^{-1}(u_n)\}\) is a sample of simulated values from variable \(X\). In computational terms, to perform the simulation process, generating a random uniform sample is accomplished using pseudorandom number generation algorithms.
There will be situations where it is necessary to simulate the behavior of a data sample \(x_1,x_2,\cdots,x_n\) from an unknown distribution, for which the cumulative distribution function is not available. In these cases, the empirical cumulative distribution function (ecdf) is used for simulation. The ecdf is defined as follows:
\[\widehat{\mathbb F}_{x_1,x_2,\cdots,x_n}(t)=\frac{1}{n}\sum_{i=1}^n\mathbb {I}_{x_i\leq t}\]
where \(\mathbb{I}_{x_i\leq t}=\left\lbrace\begin{aligned} 1,\; x_i\leq t\\ 0,\; x_i> t\end{aligned}\right.\)
