[R] qpcR propagate problem

[Robert Keams Valenzuela]

I am using the propagate function of the qpcR package to estimate the
standard error for an expression.  This expression is simple enough
that I am able to calculate the first-order propagation of error,
which is what the documentation on propagate states that it does.
However, the results are not consistent.  Can someone help me figure
out where the error is?  I should note that I have applied the same
script to a simpler expression(i.e., y/x), but with different partial
derivatives, and the results are consistent (i.e., my calculation
equals the output of propagate).

The script I am using is as follows:
#######################################################################
rm(list=ls())
j=1

r_mat=matrix(0,9,j)
act=matrix(0,1,j)
est=matrix(0,1,j)

for (i in 1:j){

    r=matrix(rnorm(9),9,1)
    r_mat[,i]=r

    a=r[1]
    b=r[2]
    c=r[3]
    var_a=abs(r[4])
    var_b=abs(r[5])
    var_c=abs(r[6])
    var_ab=abs(r[7])
    var_ac=abs(r[8])
    var_bc=abs(r[9])

    #f=(a*b+c)/(a*c+b)
    #partial derivatives
    pfpa=((a*c+b)*b-(a*b+c)*c)/(a*c+b)^2
    pfpb=((a*c+b)*a-(a*b-c))/(a*c+b)^2
    pfpc=((a*c+b)-(a*b+c)*a)/(a*c+b)^2
    #first-order propagation of error
    var_f=(pfpa^2)*var_a + (pfpb^2)*var_b + (pfpc^2)*var_c +
2*(pfpa*pfpb*var_ab + pfpa*pfpc*var_ac + pfpb*pfpc*var_bc)
    act_stderr=sqrt(var_f)

    print(act_stderr)
    act[1,i]=act_stderr

    #*****Estimate standard error
    EXPR <- expression((aa1*bb1+cc1)/(aa1*cc1+bb1))
    aa1=c(a,sqrt(var_a))
    bb1=c(b,sqrt(var_b))
    cc1=c(c,sqrt(var_c))
    DF <- cbind(aa1,bb1,cc1)
    RES1 <- propagate(EXPR, DF, type = "stat", use.cov = TRUE,
                      do.sim = FALSE, verbose = FALSE, plot = TRUE)
    est_stderr=RES1$summary$Prop[2]

    print(est_stderr)
    est[1,i]=est_stderr

}
plot(act,est, xlim=c(0,200), ylim=c(0,200))
lm_mod <- lm(est[1,] ~ act[1,])
abline(lm_mod, col="red") # regression line (y~x)
lm_coef <- round(coef(lm_mod), 3) # extract coefficients
mtext(bquote(y == .(lm_coef[2])*x + .(lm_coef[1])), adj=1, padj=0) #
display equation

data=rbind(r_mat,act,est)
act
est
#######################################################################
Thank you in advance for any input.
rkv1 at email.arizona.edu