(* Content-type: application/vnd.wolfram.mathematica *)

(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)

(* CreatedBy='Mathematica 9.0' *)

(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[       157,          7]
NotebookDataLength[      5443,        339]
NotebookOptionsPosition[      5198,        326]
NotebookOutlinePosition[      5542,        341]
CellTagsIndexPosition[      5499,        338]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{
Cell["\<\
(* RELPOT-09 - This program is used to calculate relative potency of two

   stimuli. The bioassay data may have been transformed using probit,

   logit, or complementary log-log transformations. The independent

   variable (time or dose) may have been transformed using log10. 

   Note that the same transformations must have been used for the two

   regression lines. That is, if you used the probit transformation of 

   the dependent variable and the log10 transformation of the 

   independent variable for one data set, you must have used those same

   transformations for the other data set. The program reads input from 

   the file RPinput. The user should make no changes to the program 

   below. *)



(* This program was developed by:



   Jim Throne

   USDA, ARS, SJVASC

   9611 South Riverbend Avenue

   Parlier CA 93648

   (559) 596-2702

   james.throne@ars.usda.gov



   Please notify me of problems or suggestions.  3 June 2013 *)



Off[]



(* Read the input data set. *)



data = ReadList[\"RPinput\"];



filename = data[[1]];

title1 = data[[2]];

title2 = data[[3]];

a1 = data[[4]];

a2 = data[[5]];

b1 = data[[6]];

b2 = data[[7]];

sea1 = data[[8]];

sea2 = data[[9]];

seb1 = data[[10]];

seb2 = data[[11]];

cov1 = data[[12]];

cov2 = data[[13]];

theta1 = data[[14]];

theta2 = data[[15]];

logx = data[[16]];



(* The next line opens the output file for writing. *)



out=OpenWrite[filename,FormatType->OutputForm,PageWidth->76];



(* The next lines put descriptions on the output file. *)



WriteString[out,title1]

WriteString[out,\" \\n\"]

WriteString[out,title2]

WriteString[out,\" \\n\"]

WriteString[out,\" \"]

WriteString[out,\" \\n\\n\"]



(* Descriptions of input data. *)



(*   a1, a2       Intercepts for the two regression lines.  *)



(*   b1, b2       Slopes for the two regression lines. *)



(*   sea1, sea2   Standard errors for the intercepts for the two

                  regression lines.  *)



(*   seb1, seb2   Standard errors for the slopes for the two

                  regression lines. *)



(*   cov1, cov2   Covariances of slope and intercept for the two

                  regression lines. *)



(* theta1,theta2  Lethal time or dose values to be compared.  For

                  example, one pesticide may take 5 days to kill

                  50% of the insects while another pesticide takes

                  7 days.  We would compare theta1 equal 5 with theta2

                  equal 7.  If log transformation of time or dose was

                  used in the analysis, then input theta as logs. *)



(* The user should not make any changes to the program below. *)



WriteString[out,\"Input data: \\n\\n\"]

WriteString[out,\"Intercept1 = \"]

Write[out,a1]

WriteString[out,\"Intercept2 = \"]

Write[out,a2]

WriteString[out,\"Slope1 = \"]

Write[out,b1]

WriteString[out,\"Slope2 = \"]

Write[out,b2]

WriteString[out,\"SE Intercept1 = \"]

Write[out,sea1]

WriteString[out,\"SE Intercept2 = \"]

Write[out,sea2]

WriteString[out,\"SE Slope1 = \"]

Write[out,seb1]

WriteString[out,\"SE Slope2 = \"]

Write[out,seb2]

WriteString[out,\"Covariance of Intercept1 and Slope1 = \"]

Write[out,cov1]

WriteString[out,\"Covariance of Intercept2 and Slope2 = \"]

Write[out,cov2]

WriteString[out,\"Lethal time or dose 1 = \"]

Write[out,theta1]

WriteString[out,\"Lethal time or dose 2 = \"]

Write[out,theta2]

WriteString[out,\"Logarithmic transformation of x used? = \"]

Write[out,logx]

WriteString[out,\" \\n\\n\"]



(* Calculate the variance of theta1 and theta2. *)



vart1=(1/(b1^2))*((sea1^2)+(2*theta1*cov1)+((theta1^2)*(seb1^2)));

vart2=(1/(b2^2))*((sea2^2)+(2*theta2*cov2)+((theta2^2)*(seb2^2)));



Do[

If[logx==\"yes\",Goto[yes]];



(* Calculate the variance of theta1 divided by theta2. *)



vart1t2=(((1/theta2)^2)*vart1)+(((-theta1/(theta2^2))^2)*vart2);



(* Calculate the confidence limits for theta1 divided by theta2. *)



ratio=theta1/theta2;

lcl95=(theta1/theta2)-(1.96*Sqrt[vart1t2]);

ucl95=(theta1/theta2)+(1.96*Sqrt[vart1t2]);



Goto[end];

Label[yes];



(* Calculate the variance of theta1 divided by theta2 for logx. *)



ratio=10^(theta1-theta2);

sigma=Sqrt[vart1+vart2];

lcl95=10^(theta1-theta2-(1.96*sigma));

ucl95=10^(theta1-theta2+(1.96*sigma));



Label[end],{i,1}]



(* Write confidence limits to file. *)



Write[out,\"theta1 divided by theta2 equals:\"]

Write[out,ratio]

WriteString[out,\" \\n\\n\"]

Write[out,\"The 95% confidence limits for theta1 divided by theta2 are:\"]

Write[out,lcl95,\" - \",ucl95]

WriteString[out,\" \\n\\n\"]



Close[out]\
\>", "Input"]
},
WindowSize->{707, 787},
WindowMargins->{{24, Automatic}, {Automatic, 30}},
FrontEndVersion->"9.0 for Microsoft Windows (64-bit) (January 25, 2013)",
StyleDefinitions->"Default.nb"
]
(* End of Notebook Content *)

(* Internal cache information *)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[557, 20, 4637, 304, 5165, "Input"]
}
]
*)

(* End of internal cache information *)

