====== Isometry test and size correction ======

<code>
isometry.test=function(response, reference, specimen.id, expected.slope, alpha=0.05){ #1 Creates function "isometry.test", with arguments "response", "reference", "specimen.id", "expected.slope" and "alpha".
  
  ## VERIFYING THE ARGUMENTS
  
  if ((missing(response)) || (is.na(suppressWarnings(as.numeric((response)[[1]])))) || is.factor(data.frame(response)[[1]])) #2 Verifies if "response" was inserted and if its coercible to numeric.
  { 
    stop("response should be specified as a numeric vector containing the response measurement(s) for isometry test") #3 If not, interrupts the function and gives a message to the user.
  }
  
  response=as.data.frame(response) #4 Transforms "response" in a data frame, so its length can be read as number of variables.
  
  if ((missing(reference)) || (length(reference) != max(c(length(response),length(as.data.frame(t(response)))))) || (is.na(suppressWarnings(as.numeric(reference[[1]])))) || is.factor(data.frame(reference)[[1]])) #5 Verifies if "reference" was inserted, if it has the same length as response, and is coercible to numeric.
  {
    stop("reference should be specified as a numeric vector containing the reference measurements for isometry test.") #6 If not, interrupts the function and gives a message to the user.
  }
  
  if ((missing(specimen.id)) || (class(specimen.id) != "character") || (length(specimen.id) != length(reference))) #7 Verifies if "specimen.id" was inserted, if it is a character and if it has the same length as "reference".
  {
    stop("specimen.id should be specified as a character vector containing the names of specimens in data table.") #8 If not, interrupts the function and gives a message to the user.
  }
  
  if ((missing(alpha)) || (class(alpha) != "numeric") || (alpha<=0) || (alpha>=1)) #9 Verifies if "alpha" was inserted, if it is numeric and >0 and <1.
  {
    message("\talpha should be numeric and >0 or <1; herein default set to 0.05.\n") #10 If not, gives a message to the user and uses the default value set above with the argument.
  }
  
  if ((colnames(response)==as.character(rep(1:length(colnames(response)),1))) || (match(colnames(response), specimen.id, nomatch = 0) != rep(0,length(response)))) #11 Verifies if the input data frame "response" has the specimens in rows.
  {
    response=as.data.frame(t(response)) #12 If not, transposes so the specimens are in rows. Will have further implications in calculating number of response variables.
  }
  
  if ((missing(expected.slope)) || (class(expected.slope) != "numeric") || (length(expected.slope) != length(colnames(response)))) #13 Verifies if "expected.slope" was inserted, if it is numeric and has the same length as "response".
  {
    message("\texpected.slope should be a concatenation of the ratios: power of response variable/power of the reference variable; herein default set to 1 for all.\n") #14 #10 If not, gives a message to the user.
    expected.slope=rep(1,length(colnames(response))) #15 Uses the default value (set here).
  }
  
  ## FUNCTION FOR ISOMETRY TEST AND SIZE CORRECTION
  
  # Setting and configuring specific parameters
  
  ref.name=readline("\tPlease tell me the name of your reference variable: ") #16 Asks for name of reference variable, that will be used in the output.
  result=matrix(data=NA, nrow=1, ncol=length(colnames(response))) #17 Creates the empty "result" matrix, with 1 row and "response" columns.
  result=as.data.frame(result) #18 Transforms "result" into a data frame, so it will later accept characters.
  result2=matrix(data=NA, nrow=length(rownames(response)), ncol=length(colnames(response))+1) #19 Creates the empty "result2" matrix with as muchs rows as specimens, and as much columns as reference+response variables.
  reference=as.numeric(reference) #20 Coerces "reference" to numeric, in case it was read before as characters.
  result2[,1]=reference #21 Sets reference variable to the first column of "result2".
  
  if (length(response)==1) #22 Verifies if there is only one response variable.
  {
    res.name=readline("\tPlease tell me the name of your response variable: ") #23 If yes, asks the user for the name of the response variable.
    colnames(response)=res.name #24 Attributes the name of the response variable to its specific column, which will be used to construct the output matrices.
  }
  
  response.all=as.matrix(response) #25 Creates the matrix "response.all" with "response", as "response" will enter the for loop that follows.
  
  # Loop to compare each response variable to the reference variable, perform the isometry test and the specific size-correction.
  
  for (k in 1:length(colnames(response))) #26 Loop using counter k from 1 to the length of "response"
  {
    response=response.all[,k] #27 Creates object "response" with k column of "response.all"
    response=as.numeric(response) #28 Coerces "response" to numeric, in case it was read before as character.
    response=log(response) #29 Log transformation of response variable attributed to "response"
    reference2=log(reference) #30 Log transformation of reference variable attributed to "reference2", so that it overwrites "reference", similar to what was done with "response.all" and "response".
    response.m=lm(response~reference2) #31 Creates the object "response.m" with the linear model of "response" as a function of "reference2".
    f=(summary(response.m))$fstatistic #32 Creates the object "f" with the F-statistic values of the linear model.
    p=pf(f[1],f[2],f[3],lower.tail=F) #33 Creates the object "p" with the p-value from the F-statistic values above.
    
    if (p > alpha) #34 Verifies if "p" is bigger than "alpha".
    {
      result[k]="lm not significant" #35 If yes, attribute "lm not significant" to element k in "result".
    }
    else #36 If not, linear model is significant, so proceeds with calculations.
    {
      slope=(summary(response.m))$coefficients[2] #37 Creates the object "slope" and stores the slope coefficient from "response.m".
      slope.se=(summary(response.m))$coefficients[4] #38 Creates the object "slope.se" and stores the slope's standard error coefficient from "response.m".
      tvalue=abs((slope-expected.slope[k])/slope.se) #39 Creates the object "tvalue", with the t test value from the comparison between slope and expected slope.
      pvalue=(1-pt(tvalue,length(response)-2))*2 #40 Creates the object "pvalue", with the p-value corresponding to the t test value above (degrees of freedom = n-2).
      
      if (pvalue > alpha) #41 Verifies if "pvalue" is bigger than "alpha".
      {
        result[k]="Isometry" #42 If yes, the slope from response.m is not significantly different from the expected slope, so it attributes "Isometry" to element k in "result".
        result2[,k+1]=as.numeric(response.all[,k])*((mean(result2[,1])/result2[,1])^expected.slope[k]) #43 Makes isometric transformation based on expected slope, following the formula in Lleonart et al. (2000).
      }
      else #44 Verifies if "pvalue" is not bigger than alpha.
      { 
        result[k]="Allometry" #45 If it is not bigger, the slope from response.m is significantly different from the expected.slope, so it attributes "Allometry" to element k in "result".
        result2[,k+1]=as.numeric(response.all[,k])*((mean(result2[,1])/result2[,1])^slope) #46 Makes allometric transformation based on observed slope, following the formula in Lleonart et al. (2000).
      }
      
    }
    
    colnames(result)[k]=colnames(response.all)[k] #47 Establishes column names of "result" as the column names in response.all.
  }
  
  ## SETTING AND RETURNING THE OUTPUT
  
  result2[,1]=mean(result2[,1]) #48 Makes the size correction of the reference variable, by changing it in the first column of "result2" for its mean.
  colnames(result2)=c(ref.name,colnames(response.all)) #49 Establishes the column names of "result2" as the reference name given by the user, and the column names of "response.all". 
  rownames(result2)=specimen.id #50 Establishes the row names of "result2" as the specimen names given by the user.
  all=array(list(result, result2)) #51 Creates the output object "all", as array with a list that includes "result" and "result2".
  return(all) #52 Returns to the user the output array "all" with the results of the isometry test and size correction.
}

</code>