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. r...@email.arizona.edu ______________________________________________ R-help@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.