(* 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[     15143,        869]
NotebookOptionsPosition[     14896,        856]
NotebookOutlinePosition[     15240,        871]
CellTagsIndexPosition[     15197,        868]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{
Cell["\<\
(* BACKTRN2-09 - This program is used to calculate back transformations

  (transform back to proportion organisms responding to stimulus) after

  probit, logit, or complementary log-log transformation of bioassay

  data. This program is for analyzing dose-mortality or time-mortality 

  data that have independent observations at each dose or time. Use

  BACKTRAN for analyzing serial time-mortality data where observations

  are made on the same group of organisms over time. Dose may be

  transformed as log10(dose). The user may also choose to calculate

  residuals of both untransformed and transformed data. The program 

  reads input from the file BT2input. 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[\"BT2input\"];



filename = data[[1]];

title1 = data[[2]];

title2 = data[[3]];

intstart = data[[4]];

intend = data[[5]];

ints = data[[6]];

tranx = data[[7]];

trany = data[[8]];

a = data[[9]];

b = data[[10]];

resids = data[[11]];

nodoses = data[[12]];

contdead = data[[13]];

expdead = data[[14]];

doses = data[[15]];

ncont = data[[16]];

inexp = data[[17]];



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



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



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



WriteString[out,title1]

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

WriteString[out,title2]

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

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



(* Descriptions of input data. *)



(*   intstart   Lowest dose at which observations were made.  Note: a

                control dose should not be included.  Control data are

                not used to calculate the regression line, but are used

                only to adjust for control mortality. *)



(*   intend     Highest dose at which observations were made. *)



(*   ints       Size of intervals desired for output data set.  For

                example, if dose goes from 0 to 1 mg you may want output

                 in 0.1 mg intervals so set ints equal to 0.1. *)



(*   tranx      Set to yes if log10 of dose was used in analysis. 

                Otherwise, set to no.  Enclose in quotation marks. *)



(*   trany      Indicate whether probit, logit, or cll transformation of

                y was used in analysis.  Enclose in quotation marks. *)



(*   a          Intercept from analysis.  *)



(*   b          Slope from analysis. *)



(*   resids     Indicate whether residuals should be calculated and

                printed to the output file (yes or no in quotation 

                marks). *)



(*   nodoses    Number of doses for which data will be entered. *)



(*   contdead   The observed number of insects dead in the control. *)



(*   expdead    The observed number of insects dead in the treatment at

                each dose, enclosed in brackets. *)



(*   doses      The doses at which observations were made, enclosed in

                brackets. *)



(*   ncont      The initial number of insects in the control.  If the 

                number of insects dead in the treatment at each dose 

                has already been corrected for control mortality, set

                ncont equal to any integer greater than 0. *)



(*  inexp      The number of insects tested at each dose. *)





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



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

WriteString[out,\"Starting dose = \"]

Write[out,intstart]

WriteString[out,\"Ending dose = \"]

Write[out,intend]

WriteString[out,\"Size of intervals = \"]

Write[out,ints]

WriteString[out,\"Logarithmic transformation of dose = \"]

Write[out,tranx]

WriteString[out,\"Transformation = \"]

Write[out,trany]

WriteString[out,\"Intercept = \"]

Write[out,a]

WriteString[out,\"Slope = \"]

Write[out,b]

WriteString[out,\"Residuals desired = \"]

Write[out,resids]

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

ndist=NormalDistribution[0,1];



Do[If[resids==\"no\",Goto[no]];



(* This section of the program calculates residuals. *)



WriteString[out,\"Input data for calculating residuals: \\n\\n\"];

WriteString[out,\"Number of doses = \"];

Write[out,nodoses];

WriteString[out,\"Number dead in control = \"];

Write[out,contdead];

WriteString[out,\"Number dead at each dose = \"];

Write[out,expdead];

WriteString[out,\"Doses = \"];

Write[out,doses];

WriteString[out,\"Number of controls = \"];

Write[out,ncont];

WriteString[out,\"Number treated at each dose = \"];

Write[out,inexp];

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



(* Apply log10 transformation to dose. *)



If[tranx==\"yes\",doses=Log[10,doses] //N];



(* Set up vectors containing the doses (t), the number dead in the

control (dcont), the number dead at each dose in the treatment (dexp),

the number tested at each dose (nexp), and the number of doses (m).  

This is a renaming of variables that is due to an earlier version of the

program being changed. *)



t=doses;

dcont=contdead;

dexp=expdead;

nexp=inexp;

m=nodoses;



(* Calculate the proportion of insects that die at each dose (pexpinit),

and the proportion that die in the control (pcont). *)



pexpinit=dexp/nexp //N;

pcont=dcont/ncont;



(* Apply Abbott's formula to correct for control mortality, and 

   correct for 0 or all dead at a dose. *)

                       

pexp=Table[(pexpinit[[i]]-pcont)/(1-pcont),{i,m}];

up3=Table[If[pexp[[i]]<=0,(0.5)/(nexp[[i]]+1),pexp[[i]]],{i,m}];

up=Table[If[up3[[i]]>=1.,(dexp[[i]]+0.5)/(nexp[[i]]+1.),

  up3[[i]]],{i,m}];



(* Calculate the standard error for each residual for calculating

   standardized residuals. *)



seresid=Sqrt[up*(1-up)/(nexp)];



(* Use the regression equation to calculate predicted probit for

   each time interval or dose *)



tz=a+b*t;



(* Calculate the predicted proportion kill at each x depending on

   the transformation used *)



If[trany == \"probit\",tp=Map[CDF[ndist,#]&,tz]];

If[trany == \"logit\",tp=Exp[tz]/(1+Exp[tz])];

If[trany == \"cll\",tp=1-Exp[-Exp[tz]]];



(* Calculate the probit corresponding to the observed proportion

   kill. *)



Do[If[trany != \"probit\",Goto[other]];

y=Map[Quantile[ndist,#]&,up];

uz=Table[y[[i]],{i,m}];

Goto[endtrans];

Label[other];

If[trany != \"logit\",Goto[cll],uz=Table[Log[up[[i]]/

  (1-up[[i]])],{i,m}]];

Goto[endtrans];

Label[cll];

uz=Table[Log[-Log[1-up[[i]]]],{i,m}];

Label[endtrans];,{i2,1}];



(* Calculate residuals for proportion kill and probits.  Residuals

   are calculated as observed minus predicted. *)



pres=pexp-tp;

zres=uz-tz;



(* Calculate standardized residual for proportion kill. *)



stdresid=pres/seresid;



(* Convert log(dose) to dose, if necessary. *)



If[tranx==\"yes\",t=10^t];



(* If residuals are very small, set it equal to 0.0000 so that

   scientific notation will not be used and the data can be easily

   imported into a graphics package.  The point will essentially look

   like 0 in a graph, no matter how much smaller than 0.0001 it is. *)



presout2=Table[If[pres[[i]]>0. && pres[[i]]<0.00005,0.0000,pres[[i]]],

  {i,m}];

zresout2=Table[If[zres[[i]]>0. && zres[[i]]<0.00005,0.0000,zres[[i]]],

  {i,m}];

sresout2=Table[If[stdresid[[i]]>0. && stdresid[[i]]<0.00005,0.0000,

  stdresid[[i]]],{i,m}];

presout=Table[If[presout2[[i]]<0. && presout2[[i]]>-0.00005,0.0000,

  presout2[[i]]],{i,m}];

zresout=Table[If[zresout2[[i]]<0. && zresout2[[i]]>-0.00005,0.0000,

  zresout2[[i]]],{i,m}];

sresout=Table[If[sresout2[[i]]<0. && sresout2[[i]]>-0.00005,0.0000,

  sresout2[[i]]],{i,m}];



(* Write output for proportion kill to file. *)



Write[out,

\"This section of the output contains the observed proportion killed, the 

proportion killed backtransformed from predictions using the regression 

equation, the residual (observed minus predicted) proportion killed, and

the standardized residual (residuals divided by their standard errors).\"];

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

Write[out,\"        Dose   Untransformed Proportion Killed           Residual \
     Std.

                                                                     Residual\
\\n\"];

Write[out,\"                     Observed    Predicted\\n\"];

Do[Write[out,PaddedForm[SequenceForm[\"  \",t[[i]],\"      \",

  pexp[[i]],\"   \",tp[[i]],\"              \",presout[[i]],\"

  \",sresout[[i]]], {8,4}]],{i,m}];



(* Write output for probits to file. *)



If[trany != \"probit\",Goto[other1]];

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

Write[out,

\"This section of the output contains the probit corresponding to the 

observed, corrected proportion killed; the probit predicted using the 

regression equation; and the residual (observed minus predicted) \
probits.\\n\\n\"];

Goto[otherend];

Label[other1];

If[trany != \"logit\",Goto[other2]];

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

Write[out,

\"This section of the output contains the logit corresponding to the 

observed, corrected proportion killed; the logit predicted using the 

regression equation; and the residual (observed minus predicted) \
logits.\\n\\n\"];

Goto[otherend];

Label[other2];

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

Write[out,

\"This section of the output contains the CLL transform corresponding to the 

observed, corrected proportion killed; the CLL transform predicted using the 

regression equation; and the residual (observed minus predicted) CLL 

transforms.\\n\\n\"];

Label[otherend];

Write[out,\"            Dose     Transformed Proportion Killed         Residual

\"];

Write[out,\"\\n                          Observed      Predicted\\n\"];

Do[Write[out,PaddedForm[SequenceForm[\"      \",t[[i]],\"      \",

  uz[[i]],\"     \",tz[[i]],\"          \",zresout[[i]]],{8,4}]],{i,m}];



Label[no];

,{i5,1}]



(* This section does the backtransformation at specified intervals. *)



(* Calculate number of doses *)



inc=((intend-intstart)/ints)+1;



(* Set up a vector containing the doses (x) *)



x=Table[i,{i,intstart,intend,ints}];



(*  Apply log10 transformation to dose, if applicable. *)



If[tranx==\"yes\",logx=Log[10,x] //N];



(* Use the regression equation to calculate predicted y for each dose

   *)



If[tranx==\"yes\",y=a+b*logx];

If[tranx==\"no\",y=a+b*x];



(* Calculate the proportion kill at each x depending on the

   transformation used *)



If[trany == \"probit\",p=Map[CDF[ndist,#]&,y]];

If[trany == \"logit\",p=Exp[y]/(1+Exp[y])];

If[trany == \"cll\",p=1-Exp[-Exp[y]]];



(* If p is very small set it equal to 0.0001 so that scientific

   notation will not be used and the data can be easily imported

   into a graphics package.  The point will essentially look like

   0 in a graph, no matter how much smaller than 0.0001 it is. *)



pout=Table[If[p[[i]]<0.0001,0.0001,p[[i]]],{i,inc}];



If[trany == \"probit\",title=\"  Probit\"];

If[trany == \"logit\",title=\"   Logit\"];

If[trany == \"cll\",title=\"CLL transform\"];



(* Write data to output file. *)



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

WriteString[out,\"Backtransformed output from the regression equation\"];

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

Write[out,\"            Dose\",\"          \",title,\"        Prop. Kill\\n\"];

Do[Write[out,PaddedForm[SequenceForm[x[[i]],\"    \",

  y[[i]],\"    \",pout[[i]]],{14,4}]],{i,inc}];



(* If resids=no, can't create graphs. *)



Do[If[resids==\"no\",Goto[no2]];



(* Graph output data and write to file. *)



(* Plot obs. vs. pred. % kill *)



g1=ListPlot[Table[{t[[i]],pexp[[i]]},{i,m}],DisplayFunction ->Identity];

g2=ListPlot[Table[{x[[i]],pout[[i]]},{i,inc}],Joined->True,

DisplayFunction ->Identity];

Show[g1,g2,DisplayFunction ->(Print[#]; &),Frame->True,PlotLabel->

\"Observed and Calculated Proportion kill vs Dose\", DisplayFunction->

$DisplayFunction,AxesOrigin->{0,0}];



(* Plot residual % kill vs. dose *)



g3=ListPlot[Table[{t[[i]],presout[[i]]},{i,m}],DisplayFunction

  ->Identity];

Show[g3,DisplayFunction ->(Print[#]; &),Frame->True,PlotLabel->

\"Residual of proportion kill vs dose\",DisplayFunction->

$DisplayFunction,AxesOrigin->{0,0}];



(* Plot standardized residual % kill vs. dose *)



g4=ListPlot[Table[{t[[i]],sresout[[i]]},{i,m}],DisplayFunction

  ->Identity];

Show[g4,DisplayFunction ->(Print[#]; &),Frame->True,PlotLabel->

\"Standardized residual of proportion kill vs dose\",DisplayFunction->

$DisplayFunction,AxesOrigin->{0,0}];



(* Plot observed and predicted probits vs. dose *)



g5=ListPlot[Table[{t[[i]],uz[[i]]},{i,m}],DisplayFunction

  ->Identity];

g6=ListPlot[Table[{x[[i]],y[[i]]},{i,inc}],Joined->True,

DisplayFunction ->Identity];

If[trany != \"probit\",Goto[other3]];

Show[g5,g6,DisplayFunction ->(Print[#]; &),Frame->True,PlotLabel->

\"Observed and Calculated Probits vs Dose\",DisplayFunction->

$DisplayFunction,AxesOrigin->{0,0}];

Goto[endoth3];

Label[other3];

If[trany != \"logit\",Goto[other4]];

Show[g5,g6,DisplayFunction ->(Print[#]; &),Frame->True,PlotLabel->

\"Observed and Calculated Logits vs Dose\",DisplayFunction->

$DisplayFunction,AxesOrigin->{0,0}];

Goto[endoth3];

Label[other4];

Show[g5,g6,DisplayFunction ->(Print[#]; &),Frame->True,PlotLabel->

\"Observed and Calculated CLL transforms vs Dose\",DisplayFunction->

$DisplayFunction,AxesOrigin->

{0,0}];

Label[endoth3];



(* Plot residual probits vs. dose *)



g7=ListPlot[Table[{t[[i]],zresout[[i]]},{i,m}],DisplayFunction

  ->Identity];

If[trany != \"probit\",Goto[other7]];

Show[g7,DisplayFunction ->(Print[#]; &),Frame->True,PlotLabel->

\"Residual probits vs dose\",DisplayFunction->$DisplayFunction,

AxesOrigin->{0,0}];

Goto[endoth7];

Label[other7];

If[trany != \"logit\",Goto[other8]];

Show[g7,DisplayFunction ->(Print[#]; &),Frame->True,PlotLabel->

\"Residual logits vs Dose\",DisplayFunction->$DisplayFunction,

AxesOrigin->{0,0}];

Goto[endoth7];

Label[other8];

Show[g7,DisplayFunction ->(Print[#]; &),Frame->True,PlotLabel->

\"Residual CLL transforms vs dose\",DisplayFunction->$DisplayFunction,

AxesOrigin->{0,0}];

Label[endoth7];



Label[no2];

,{i6,1}]



Close[out]\
\>", "Input"]
},
WindowSize->{707, 787},
WindowMargins->{{72, Automatic}, {38, Automatic}},
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, 14335, 834, 14107, "Input"]
}
]
*)

(* End of internal cache information *)

