[R] Help with modFit of FME package 2

[Paola Lecca]

* Apologies for multiple posting *
I attached to my previous e-mail a .r file, and it was not permitted by the
rules of the mailing lis. Again, please receive my sincere apologies for
this.

I re-send again the e-mail with .txt attachemnt in the hope someone an help
me to solve my problem.


I'm trying to fit a set an ODE to an experimental time series. In the
attachment you find the R code I wrote using modFit and modCost of FME
package and the file of the time series.

When I run summary(Fit) I obtain this error message, and the values of the
parameters are equal to the initial guesses I gave  to them.

The problem is not due to the fact that I have only one equation (I tried
also with more equations, but I still obtain this error).

I would appreciate if someone could help me in understanding the reason of
the error and in fixing it.

Thanks for your attention,
Paola Lecca.

Here the error:

> summary(Fit)

Parameters:
             Estimate Std. Error t value Pr(>|t|)
pro1_strength        1         NA      NA       NA

Residual standard error: 2.124 on 10 degrees of freedom
Error in cov2cor(x$cov.unscaled) : 'V' is not a square numeric matrix
In addition: Warning message:
In summary.modFit(Fit) : Cannot estimate covariance; system is singular




______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.




-- 
*Paola Lecca, PhD*
*The Microsoft Research - University of Trento*
*Centre for Computational and Systems Biology*
*Piazza Manci 17 38123 Povo/Trento, Italy*
*Phome: +39 0461282843*
*Fax: +39 0461282814*
-------------- next part --------------
time	pp1_mrna	
0	0
2	2.754
4	2.958
6	4.058
8	3.41
10	3.459
12	2.453
14	1.234
16	2.385
18	3.691
20	3.252
-------------- next part --------------
require(deSolve)
require(FME)


##################################################################################
#                                      PART 1                                    #
##################################################################################

# Differential equations
model_1_part_1 <- function(t, S, parameters)
                 {
                  with(as.list(parameters), {
                  #
                  cod1 = pro1_strength
                  #
                  pp1_mrna_degradation_rate <- 1
                  ###############################################################
                  #
                  v1 = cod1
                  v2 = pp1_mrna_degradation_rate * S[1]
                  #
                  #################################################################
                  #
                  dS1 = v1 - v2
                  #
                  ####################################################################
                  #
                  list(c(dS1))
                 })
                }


# Parameters
parms_part_1 <- c(pro1_strength = 1.0)

# Initial values of the species concentration
S <- c(pp1_mrna = 0)

times <- seq(0, 20, by = 2)

# Solve the system
ode_solutions_part_1 <- ode(S, times, model_1_part_1, parms = parms_part_1)
ode_solutions_part_1

summary(ode_solutions_part_1)

## Default plot method
plot(ode_solutions_part_1)


########################################################################################
########################################################################################

# Estimate of the parameters
experiment <- read.table("./wild_pp1_mrna.txt", header=TRUE)

rw <- dim(experiment)[1]

names <- array("", rw)
for (i in 1:rw)
    {
     names[i] <- "pp1_mrna"
    }
names
      
observed_data_part_1 <- data.frame(name = names,
                   time = experiment[,1], val = experiment[,2])

observed_data_part_1

ode_solutions_part_1

Cost_function <- function (pars) {
 out <- ode_solutions_part_1
 cost <- modCost(model = out, obs = observed_data_part_1, y = "val")
 cost
 }

Cost_function(parms)

# Fit the model to the observed data
Fit <- modFit(f = Cost_function, p = parms_part_1)
Fit


# Summary of the fit
summary(Fit)

# Model coefficients
coef(Fit)

# Deviance of the fit
deviance(Fit)