User manual for fraping package

User Manual

fin

Introduction

Introducción

Download packages

En la paquetería fraping se implementan las funciones necesarias para realizar el análisis de datos de FRAP, sin embargo, para que esta funcione correctamente debe instalarse también la paquetería itz. Ambos son repositorios de GitHub y su instalación se lleva a cabo mediante la función install_github de la paquetería devtools de R.

devtools::install_github("artitzco/itz")
devtools::install_github("artitzco/fraping")

De esta forma, las librerías quedarán instaladas en su equipo, por lo que cada que inicie un nuevo proyecto, para el análisis de datos de FRAP, únicamente tendrá que importar las librerías correspondientes:

library(itz)
library(fraping)

La función newDataFrap crea una base para un nuevo grupo de datos de FRAP, la función sintetiza la información y devuelve una lista que permite al usuario manipular sus datos con mayor facilidad. Los parámetros que admite son el nombre del grupo y la intensidad del fotoblanqueo en porcentaje. Para más detalles ir a newDataFrap.

Para evitar futuros problemas, el nombre del grupo

  1. No debe incluir espacios,
  2. Debe inicial con una letra,
  3. Sólo se admiten los caracteres especiales “.” y “_“.
A <- newDataFrap("Neurites", 80)
B <- newDataFrap("Secondary.Neurites", 80)
C <- newDataFrap("Dendrites", 80)

Import data

Una vez inicializados los grupos de datos es hora de añadir la información. Una forma de hacerlo es por medio de la función 'newDataFrap'$addDir, la cual requiere que el usuario organice las tablas de FRAP en tres ficheros de su equipo: el primero con los datos de recuperación, el segundo con las tablas de control y el tercero con la información del fondo. Los parámetros que admite la función 'newDataFrap'$addDir son la duración del experimento en segundos, el índice del fotograma en el que se aplica el fotoblanqueo, la ubicación de los ficheros de datos, creados previamente, (recuperación, control y fondo, en el mismo orden). Los datos del fondo pueden ser omitidos, sin embargo, debe tener en cuenta que esto repercute directamente en la estandarización de los datos.

A$addDir(480, 4,
         "data/GFP/GFP Secondary Neurite",
         "data/GFP/GFP Secondary Neurite Control",
         "data/GFP/GFP Secondary Neurite Back")

B$addDir(480, 4,
         "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,
         "data/LIFEACT GFP/LIFEACT GFP Dendrite",
         "data/LIFEACT GFP/LIFEACT GFP Dendrite Control",
         "data/LIFEACT GFP/LIFEACT GFP Dendrite Back")

Al momento de crear un nuevo grupo de datos, por medio de la función newDataFrap, el programa asigna al grupo un color aleatorio que servirá para diferenciarlo del resto a la hora de realizar las gráficas. En caso de que el usuario quiera reasignar dicho color puede hacerlo por medio de la función 'newDataFrap'$setColor, la cual recibe como parámetro el valor del color en formato hex, también puede elegir uno del vector colors() o asignarlo por medio de la función color de la paquetería itz.

A$setColor("#05753D")
B$setColor("#0000A6")
C$setColor("#E014F7")

Antes de realizar un ajuste paramétrico es posible que el usuario prefiera realizar uno o más gráficos de los datos, para hacer una inspección visual de los mismos, esto se logra con la función plotRecover, recibe como parámetros las bases creadas con newDataFrap, tantas como se requiera, además de diversos parámetros que permitirán personalizar los gráficos. Para más detalles ir a plotRecover.

plotRecover(A, B, C, plot.lines=T)

plotRecover(A, B, C, plot.shadow=T, plot.mean=T)

plotRecover(A, B, C, plot.points=T)

Choose a probability model

El modelo paramétrico propuesto para describir el comportamiento de la curva de curva de recuperación de la fluorescencia a partir del momento del fotoblanqueo, es el siguiente:

\[F(t)=f_{min}+\alpha\,(1-f_{min})\,p(t-t_0|\Theta),\] donde \(p\) es una función de distribución acumulada parametrizada por \(\Theta\), \(t_0\) es el tiempo del fotoblanqueo y \(f_{min}\) es el valor mínimo de la fluorescencia después del fotoblanqueo, también se define el valor teórico \(f_{max}\) como el valor máximo que puede alcanzar la fluorescencia después del fotoblanqueo, y está dado por: \[f_{max}= f_{min}+\alpha(1-f_{min}).\]

Puede observarse que la función \(F(t)\) está definida en el intervalo de tiempo \(t\in[t_0, \infty)\). De forma análoga, se define la función \(F^{^{AB}}(t)\) para representar a la curva de recuperación de la fluorescencia después del fotoblanqueo la cual está definida en el intervalo de tiempo \(t\in[0, \infty)\):

\[F^{^{AB}}(t)=f_{min}+\alpha\,(1-f_{min})\,p(t|\Theta).\]

En el apéndice de este documento se puede encontrar más información acerca de la construcción del modelo. Para que el usuario pueda realizar el ajuste paramétrico, para el posterior análisis de datos sus datos, tendrá que ocuparse únicamente de elegir el modelo de probabilidad, \(p\), más adecuado. Mediante la función newFit es posible declarar un nuevo modelo de probabilidad, la función recibe como parámetros el nombre que se le quiera dar al modelo, la función de distribución acumulada, el nombre de los parámetros de la función de probabilidad y una lista con los rangos de valores para dichos parámetros. Para más detalles ir a 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)))

Uno de los criterios más importantes para elegir entre distintos modelos de probabilidad es mediante el error medio cuadrático, aquellos modelos que tengan menor error deben ser mayormente preferidos. El error medio cuadrático está dado por la siguiente expresión: \[E_i =\frac{1}{m_i}\sum_{j=1}^{m_i} [f_i(t_j)-F(t_j)]^2,\] donde \(\{f_1,..,f_n\}\) es un grupo de datos fluorescencia y \(\{t_1,t_2,\cdots,t_{m_i}\}\) son todos aquellos valores mayores al tiempo de fotoblanqueo cuantificados para \(f_i\). La función compareFit comparar el error medio cuadrático de diferentes modelos de probabilidad ajustados a un grupo de datos, por medio de una gráfica Q-Q y el histograma del error de cada uno de los modelos. La función recibe como parámetros al grupo de datos en cuestión, una lista de los modelos que se requieran comparar, así como una serie de parámetros que permiten personalizar el gráfico. Para más detalles ir a 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)

Para más información acerca de la elección de un buen modelo de probabilidad ir al apéndice.

Fit a parametric model

Para este caso, en particular, se tiene que el modelo Weibull arroja mejores resultados, por lo que se recurre a la función Wei, para realizar el ajuste correspondiente de los datos, la función Wei, aplicada en un grupo de datos específico, devuelve una lista con la información obtenida por el ajuste paramétrico, así como la estimación de la función de recuperación de la fluorescencia \(F(t)\), entre otros datos relevantes, para más detalles ir a newFit.

Para realizar el gráfico del ajuste paramétrico de la recuperación de la fluorescencia se utiliza la función plotFit la cual recibe como parámetros, en primer lugar, a los grupos de datos por ajustar y posteriormente el modelo de probabilidad que se haya elegido para cada grupo de datos, entre parámetros opcionales que permiten personalizar y dar formato a cada gráfica, ir a plotFit para saber más. Si el modelo de probabilidad elegido es el mismo para todos los grupos de datos, basta con especificarlo una única vez como se ve a continuación:

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)

Con la función plotFit se pueden graficar estimaciones futuras, de la curva de recuperación de la fluorescencia, con solo modificar el rango del gráfico. Las gráficas obtenidas por las funciones plotRecover y plotFit pueden ser apiladas cambiando el parámetro opcional new.plot a 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)

Otro uso para la función plotFit es hacer una comparación visual entre dos o más modelos de probabilidad sobre un mismo grupo de datos.

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

Por medio de la función compareParam es posible comparar dos grupos de datos por medio de las variables asociadas al ajuste del modelo paramétrico. Los parámetros que recibe esta función son los dos grupos de datos, el modelo de probabilidad que ajusta a cada grupo de datos, la variable que se desea comparar, entre otros valores que permiten personalizar el gráfico. para saber más ir a compareParam. Los nombres de las variables pueden obtenerse de la siguiente forma:

colnames(Wei(B)$table)
## [1] "shape" "scale" "alpha" "fmax"  "UF"    "MF"    "error"

La función compareParam realiza un gráfico Q-Q de la variable elegida, cambiando el parámetro new.plot=FALSE es posible apilar los gráficos. Los detalles de la prueba \(t\mbox{-}student\) para la diferencia de medias al \(95\%\) de confianza se obtienen cambiando el parámetro return=TRUE (en consecuencia, se omitirá el gráfico), el nivel de confianza y la hipótesis alternativa pueden se modificados mediante los parámetros alternative y conf.level respectivamente, para más detalles ir a compareParam.

compareParam(B, C, fit=Wei, param="shape")

De esta forma es posible realizar un análisis convencional de la recuperación de la fluorescencia, es decir, comparar la Fracción Inmóvil (UF) y la Fracción Móvil (MF). Recuerde que \[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)

compareParam(B, C, fit=Wei, param="UF", return=T)#>>><<<
## 
##  Welch Two Sample t-test
## 
## data:  UF_Weibull_Secondary.Neurites and UF_Weibull_Dendrites
## t = -0.51751, df = 4.2562, p-value = 0.6305
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2846008  0.1933821
## sample estimates:
##  mean of x  mean of y 
## 0.08361768 0.12922699
compareParam(B, C, fit=Wei, param="MF", return=T)#>>><<<
## 
##  Welch Two Sample t-test
## 
## data:  MF_Weibull_Secondary.Neurites and MF_Weibull_Dendrites
## t = 0.51751, df = 4.2562, p-value = 0.6305
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1933821  0.2846008
## sample estimates:
## mean of x mean of y 
## 0.9163823 0.8707730

También, con la función compareParam es posible comparar dos modelos de probabilidad, ajustados a un mismo grupo de datos, con respecto al error que produce cada ajuste, puede notar que esta gráfica ya se encuentra incluida dentro del mosaico que arroja la función compareFit vista anteriormente.

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.5805180

Mean curve analysis

El análisis de la curva media consiste en evaluar el comportamiento medio de la recuperación de la fluorescencia a lo largo del tiempo, la curva media está dada por:

\[\overline F(t)=\frac{1}{n}\sum_{i=1}^{n}F^{^{AB}}_i(t),\]

donde \(n\) es el tamaño del conjunto de datos en cuestión y \(F^{AB}_i\) es la curva de recuperación de la fluorescencia después del fotoblanqueo asociada a la cuantificación de la fluorescencia \(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)

La función compareMean permite contrastar las curvas medias de dos grupos de datos por medio de una prueba \(t\mbox{-}student\) para la diferencia de medias aplicada de manera continua en un periodo de tiempo. La función recibe como parámetros a los grupos de datos por comparar, el modelo probabilidad que ajusta a los datos, entre otros valores que permiten personalizar la gráfica.

El gráfico que devuelve representa al \(p\mbox{-}value\) de la prueba \(t\mbox{-}student\), la línea de referencia indica el nivel de significancia. Cambiando el parámetro return=TRUE se obtienen los detalles de la prueba almacenados en una lista (omitiendo la realización de la gráfica), el nivel de confianza y la hipótesis alternativa pueden ser modificados mediante los parámetros alternative y conf.level respectivamente, para más detalles ir a compareMean.

compareMean(B, C, fit=Wei, lty.lines=c(2,1), col.lines=c("red","purple"), ydigits=2)

Cambiando el parámetro p.value=FALSE se calculará la diferencia de medias entre los dos grupos de datos, la línea de referencia estará situada en cero. En este caso al cambiar el parámetro return=TRUE se obtiene un vector con la diferencia de medias entre ambas curvas.

compareMean(B, C, fit=Wei, p.value=F, lty.lines=c(2,1), col.lines=c("red","magenta"), ydigits=3)

Simulation

La implementación de métodos de simulación para el análisis de la recuperación de la fluorescencia es una de las principales aportaciones de este trabajo. La simulación permite incrementar el tamaño de una muestra sin alterar sus medidas de tendencia centra ni dispersión, esto es especialmente útil para los casos en los que no se cuenta con un tamaño de muestra suficientemente grande, recuerde que las pruebas de hipótesis que dependen de procesos límite son particularmente sensibles al tamaño de la muestra. Por otro lado, el proceso de simulación se retroalimenta con una elección adecuada del modelo de probabilidad, esto se explica más adelante.

La simulación del modelo de recuperación de la fluorescencia \(F(t)\) se logra gracias a la simulación de los parámetros \(\Theta\) del modelo de probabilidad \(p(t|\Theta)\), ir al apéndice para saber más. La función simFit permite simular los datos ajustados a un modelo de probabilidad, recibe como parámetro obligatorio a un grupo de datos ajustados a un modelo de probabilidad y como parámetros opcionales el tamaño de la muestra simulada y una semilla de simulación, para más detalles sobre la semilla de simulación ir al apéndice.

Para el caso del ajuste Weibull la simulación de los datos se realiza de la siguiente forma simFit(Wei(B)) y simFit(Wei(C)), para más especificaciones sobre los datos que arroja esta función ir simFit. Por otro lado, las funciones plotFit, compareParam y compareMean tienen implementado el método de simulación, para acceder a esta función es necesario cambiar el parámetro simulated=FALSE. Podemos replicar las gráficas hechas anteriormente ahora con datos simulados:

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.10839419
compareParam(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.8916058
plotFit(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

La cuantificación de la fluorescencia, denotada por \(f\), durante el proceso de la técnica de FRAP se describe a continuación: en un inicio es constante pues no hay variación de la fluorescencia, al momento del fotoblanqueo, denominado \(t_0\), la fluorescencia cae hasta un valor mínimo \(f_{min}\) y a partir de este momento comienza un proceso gradual de recuperación de la fluorescencia que cumple necesariamente con dos condiciones, la primera es que la recuperación es no decreciente a lo largo del tiempo, y la segunda es que la fluorescencia alcanza el valor límite \(f_{max}\) que no puede ser mayor que el valor inicial antes del botoblanqueo. Lo anterior se representa en la siguiente gráfica.

En ese sentido, el modelo propuesto para la recuperación de la fluorescencia debe cumplir con las características anteriores, el modelo general para la recuperación de la fluorescencia después del fotoblanqueo es el siguiente:

\[F(t)=f_{min}+\alpha\,(1-f_{min})\,p(t-t_0|\Theta)\]

En donde \(p\) es una función de distribución acumulada continua (CDF), parametrizada por \(\Theta\).

Observaciones

  • La función \(F(t)\) está definida en el intervalo \([t_0,\infty)\)

  • El valor \(f_{max}\) está dado por \[f_{max} =\lim_{t\to\infty}F(t)\]

  • El valor \(\alpha\) está acotado por los valores \[\alpha\in[0,1]\]

  • La función de distribución acumulada \(p\) debe estar asociada a una variable aleatoria que con soporte en \((0,\infty)\), recuerde que la CDF se define de la siguiente manera:

    \[p(t|\Theta)=\int_0^td(x|\Theta)\,\mbox{d}x\] donde \(d\) es una función de densidad paramétrica.

De forma análoga, se define la función \(F^{^{AB}}(t)\) para representar a la curva de recuperación de la fluorescencia después del fotoblanqueo la cual está definida en el intervalo de tiempo \(t\in[0, \infty)\):

\[F^{^{AB}}(t)=f_{min}+\alpha\,(1-f_{min})\,p(t|\Theta).\]

Probemos que \(F\) cumple con las características de la recuperación de la fluorescencia después del fotoblanqueo:
  1. \(F(t_0)=f_{min}\)

    2. Se evalua

    \[\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}\]

  2. \(F(t)\) es no decreciente

    3. Considere a \(t_1,t_2\in\mathbb R^+\) tales que \(t_1\leq t_2\), entonces

    \[\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}\] \[\therefore F(t_1)\leq F(t_2)\]

    Lo anterior prueba que el modelo paramétrico \(F\) es no decreciente a lo largo del tiempo.

  3. \(f_{max}=f_{min}+\alpha\,(1-f_{min})\) and \(f_{max}\in[f_{min}, 1]\)

    4. Se calcula el límite de \(F(t)\) cuando \(t\) tiende a infinito:

    \[\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})\] Por otro lado se tiene que \(\alpha\in[0,1]\), entonces

    \[\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]\]

Dados un conjunto de datos de fluorescencia \(f\) y una función de distribución acumulada \(p(\bullet|\Theta)\), entonces los parámetros \((\alpha,\Theta)\) son estimados por medio del metodo de mínimos cuadrados, que consiste en minimizar el error medio cuadrático producido entre las curvas \(f\) y \(F\), la función del error medio cuadrático es la siguiente:

\[\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}\] donde \(\{t_1,t_2,\cdots,t_m\}\) son todos aquellos tiempos mayores a \(t_0\) cuantificados para \(f\).

Choosing a probability model

La función de distribución acumulada \(p(t|\Theta)\) es la base para construir el modelo paramétrico de recuperación de la fluorescencia \(F(t|\Theta)\). La elección de la función \(p\) queda a libre criterio del usuario, sin embargo, es indispensable que \(p\) tenga soporte en \((0,\infty)\), otro aspecto del que debe encargarse el usuario es la asignación de los intervalos donde toman valores los parámetros \(\Theta\), esto se logra por medio de la función newFit.

Ejemplo:

Suponga que el usuario ha decidido utilizar una función de probabilidad Pareto, entonces \[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\] Se implementa la función de distribución acumulada Pareto y se declara su ajuste paramétrico asociado:

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)))

El código anterior establece que \(\mbox{shape}\in(0,a)\) y \(\mbox{scale}\in(0,b)\). Después de ajustar el modelo Par a un grupo de datos el usuario debe verificar que los valores shape y scale se encuentren completamente contenidos en el interior de los intervalos \((0,a)\) y \((0,b)\) respectivamente, de lo contrario, tendrá que modificar los valores de \(a\) y/o \(b\) y repetir el procedimiento. Lo anterior se consigue mediante una inspección visual de las variables devueltas por el ajuste: se deben ejecutar los códigos Par(B)$table$shape y Par(B)$table$scale e inspeccionar los resultados, donde B es un conjunto de datos de FRAP obtenido por medio de la función newDataFrap.

Las funciones creadas por newFit tienen el propósito de encontrar los mejores valores \(\Theta\) tales que minimicen el error cuadrático entre el modelo \(F(t|\Theta)\) y los datos reales de recuperación de la fluorescencia. Por lo que el principal criterio para elegir el modelo de probabilidad \(p\) más adecuado debe ser: preferir el modelo que produzca el menor error cuadrático medio, esto se logra, como se vio anterior mente, por medio de la función compareFit.

compareFit(B, fit=l(Exp, Gam, Wei, Par), col.lines=c("red","yellow"), lty.lines=c(2,1), lwd.lines=2, lwd.mean=2)

Por otro lado, un segundo criterio para evaluar la viabilidad de un modelo de probabilidad \(p\) es viendo cómo se comportan las simulaciones de \(F\). Si existen discrepancias entre los datos reales de fluorescencia y las simulaciones de \(F\) entonces lo más razonable es descartar el modelo \(p\). Ya que esto implicaría que no hay una relación causal entre los parámetros \(\Theta\) y los datos de recuperación de la fluorescencia y, como se verá en el siguiente apéndice, las simulaciones de \(f_{max}\) son particularmente sensibles a las variaciones de \(\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)

Stochastic simulation theory

La simulación estocástica es un proceso mediante el cual se pretende reproducir el comportamiento de una variable aleatoria, la base de la simulación se sigue del siguiente teorema:

Sea \(X\) es una variable aleatoria y \(\mathbb F_X\) su función de distribución acumulada, sea \(U\) una variable aleatoria con distribución \(\mbox{Uniform}(0,1)\), entonces la variable aleatoria \(\mathbb F^{-1}_X(U)\) tiene la misma distribución que \(X\), es decir:

\[\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}\] \[\therefore {\mathbb F_X\left(t \right)}=\mathbb F_{\mathbb F^{-1}_X(U)}\]

    Nótese que como \(U\) es una variable uniforme, se tiene que \(f_U(u)=1,\;\; 0<u<1\). (\(f_U\) es la función de densidad de \(U\))

El teorema anterior permite simular cualquier variable aleatoria \(X\) siempre que se conozca su función de distribución acumulada \(F_X\). Sea \(\{u_1,u_2,\cdots,u_n\}\) una muestra aleatoria de valores observados de la variable aleatoria \(\mbox{Uniform}(0,1)\), entonces \(\{\mathbb F_X^{-1}(u_1),\mathbb F_X^{-1}(u_2),\cdots,\mathbb F_X^{-1}(u_n)\}\) es una muestra de valores simulados de la variable \(X\). En términos computacionales, para realizar el proceso de simulación de simulación, la generación de una muestra aleatoria uniforme se lleva a cabo por medio de algoritmos de generación de números pseudoaleatorios.

Habrá situaciones en las que se requiera simular el comportamiento de una muestra de datos \(x_1,x_2,\cdots,x_n\) de alguna distribución desconocida, por lo que no se dispone de función de distribución acumulada, en estos caso para realizar la simulación se utiliza a la función de distribución acumulada empírica (ecdf), la cual se define de la siguiente manera:

\[\widehat{\mathbb F}_{x_1,x_2,\cdots,x_n}(t)=\frac{1}{n}\sum_{i=1}^n\mathbb {I}_{x_i\leq t}\]

donde \(\mathbb{I}_{x_i\leq t}=\left\lbrace\begin{aligned} 1,\; x_i\leq t\\ 0,\; x_i> t\end{aligned}\right.\)

Simulation of a FRAP curves

Cosidere el conjunto de datos de cuantificación de fluorescencia denotados por \(\{f_1,f_2,...,f_n\}\) y el conjunto de curvas de FRAP asociadas \(\{F^{^{AB}}_1,F^{^{AB}}_2,...,F^{^{AB}}_n\}\), el proceso de simulación permite generar un conjunto de \(N_{sim}\) nuevas curvas de FRAP \(\{\mathcal F^{^{AB}}_1,\mathcal F^{^{AB}}_2,...,\mathcal F^{^{AB}}_{N_{sim}}\}\). Cada curva \(\mathcal F^{^{AB}}_j\) no se encuentra directamente asociada a ningún dato de florescencia \(f_i\), en cambio, \(\mathcal F^{^{AB}}_j\) es producto de un proceso se simulación estocastica sobre el conjunto de parámetros \(\{(\alpha_1, \Theta_1),(\alpha_2, \Theta_2),\cdots,(\alpha_n, \Theta_n)\}\), donde \((\alpha_i, \Theta_i)\) son aquellos valores que se obtuvieron por medio de la optimización por mínimos cuadrados entre \(f_i\) y \(F_i\).

Incluir imagen

De manera más específica, el proceso de simulación de las curvas de FRAP es el siguiente:

Por medio de la función de distribución acumulada empírica (vista en la sección anterior de este apéndice) de \(\{\Theta_1,\Theta_2,\cdots,\Theta_n\}\) se generan los valores \(\{ \Theta^{sim}_1,\Theta^{sim}_2,\cdots,\Theta^{sim}_m\}\), es decir:

\[\Theta_j^{sim}=\widehat{\mathbb F}_{\Theta_1,\Theta_2,\cdots,\Theta_n}^{-1}(U_j)\]

\[\mathcal F^{^{AB}}_j(t)=f_{min}+\alpha_j^{sim}\,(1-f_{min})\,p(t|\Theta_j^{sim}).\]

fraping-doc

newDataFrap

plotRecover

newFit

compareFit

plotFit

compareParam

compareMean

simFit

inicio

Comentarios

Publicar un comentario

Entradas más populares de este blog

Probabilidad de que tiemble un 19 de septiembre

Imagen
Probabilidad de que tiemble un 19 de septiembre Probabilidad de que tiemble un 19 de septiembre @artitzco CDMX, septiembre 2022 Resumen La probabilidad de que en México ocurran más de \(2\) sismos de magnitud \(7.0\) con fecha del 19 de septiembre es de \(0.0018\) \((0.1832\%)\) . Una cifra aparentemente pequeña si se le considera a éste como un evento aislado, sin embargo, esta misma probabilidad se replica de forma idéntica para 365 días del año. En realidad, lo ocurrido el 19 de septiembre es tan sólo una observación de la variable aleatoria que mide el número máximo de sismos que se registran en una misma fecha. La probabilidad de que el número máximo de sismos que se registran en una misma fecha sea igual a \(3\) es de \(0.4506\) \((45.06\%)\) y esto puede ocurrir para \(k\in[0,3]\) fechas distintas, con un \(95\%\) de confianza. La probabilidad de que ocurra lo anterior para al ...
Leer más