[R] test regression against given slope for reduced major axis regression (RMA)

[Patrick Drechsler]

Hi,

for testing if the slope of experimental data differs from a
given slope I'm using the function
"test_regression_against_slope" (see below).

I am now confronted with the problem that I have data which
requires a modelII regression (also called reduced major axes
regression (RMA) or geometric mean regression). For this I use
the function "modelII" (see below).

What would be a good way of adapting
"test_regression_against_slope" for use with RMA regression?

The question I am trying to answer is: "Does the slope acquired
from experimental data differ significantly from theoretical
predictions?"

Any feedback on this is highly appreciated. And if you need more
info do not hesitate to ask.

Kind Regards

Patrick



*test_regression_against_slope*

--8<------------------------schnipp------------------------->8---
test_regression_against_slope <- function(x,y,slope_2)
{
### TEST_REGRESSION_AGAINST_SLOPE tests a measured regression against a
### given regression.
###
### INPUT: 
###    
### x and y: raw data
### slope_2: the slope you would like to test against (ie 1/3)
###    
### OUTPUT:
###    
### pvalue: the P-value...
### upperlimit95 and lowerlimit95 give the 95 percent confidence
### intervall (two-tailed).
###
### see Sokal and Rohlf, p. 465/471
  
  n <- length(x)
  mydf <- n-2

  ## least square fit:
  x2 <- (x-mean(x))^2
  y2 <- (y-mean(y))^2

  ## regression (pedestrian solution):
  xy <- (x-mean(x))*(y-mean(y))
  slope1 <- sum(xy)/sum(x2)
  intercept_a <- mean(y) - slope1 * mean(x)

  ## model data y_hat:
  y_hat <- intercept_a + slope1 * x
  ## least squares of model data:
  y_hat2 <- (y - y_hat)^2

  s2yx <- sum(y_hat2) / (n-2)
  sb <- sqrt(s2yx/sum(x2))
  ts <- (slope1 - slope_2) / sb
  pvalue <-  2*(pt(abs(ts), df, lower.tail=FALSE))

  ## 0.95 for one-tailed 0.975 for two-tailed t-distribution with
  ## alpha<-5%:
  tval <- qt(.975, df=mydf)
  ts2 <- tval*sb

  lowerlimit95 <- slope1 - ts2
  upperlimit95 <- slope1 + ts2

  list(pvalue = pvalue,
       lowerlimit95 = lowerlimit95,
       upperlimit95 = upperlimit95)
}
--8<------------------------schnapp------------------------->8---




*modelII*

--8<------------------------schnipp------------------------->8---
modelII <- function(XjArray,YjArray){
###  ============================================================
###
### Purpose:
###
### Calculates MODEL II Regression paramaters. Also called "reduced
### major axis regression" (Prentice 1987) or "geometric mean
### regression" (Webb et al. 1981).
###
### Input:
###
### Two one dimensional arrays XjArray and YjArray containing the X
### and Y vectors.
###
### XjArray = [0 0.9 1.8 2.6 3.3 4.4 5.2 6.1 6.5 7.4]
### YjArray = [5.9 5.4 4.4 4.6 3.5 3.7 2.8 2.8 2.4 1.5]
###
### Output:
###
### A list with the following:
###
### sumXjYj As Double
### sumXj As Double, sumYj As Double
### sumXjSquared As Double, sumYjSquared As Double
### n As Long
### varXj, varYj
### output(7)
###
###  ============================================================

  sumXjYj <- 0
  sumXj <- 0
  sumYj <- 0
  n <- 0
  n <- length(XjArray)
  sumXjSquared <- 0
  sumYjSquared <- 0
  covariancexy <- 0
  
  for(i in 1:n){
    sumXjYj <- sumXjYj + XjArray[i] * YjArray[i]
    sumXj <- sumXj + XjArray[i]
    sumYj <- sumYj + YjArray[i]
    sumXjSquared <- sumXjSquared + XjArray[i]^2
    sumYjSquared <- sumYjSquared + YjArray[i]^2
  }
  
  ## Mean of X and Y vectors
  meanyj <- sumYj / n
  meanxj <- sumXj / n
  
  ## Create covariance
  for(i in 1:n){
    covariancexy <- covariancexy + ((XjArray[i] - meanxj) * (YjArray[i] - meanyj))
  }

  covariancexy <- covariancexy / n
  
  
  ## get variance of X and Y (SD)
  varXj <- (n * sumXjSquared-sumXj^2)/(n*(n - 1))
  varYj <- (n * sumYjSquared-sumYj^2)/(n*(n - 1))
  sdxij <- (sumXjSquared)-(sumXj^2/n)
  sdxik <- (sumYjSquared)-(sumYj^2/n)

  ## make beta    'sgn function to return sign with magnitude of 1
  betacoeff <- sign(covariancexy) * ((varYj^0.5) / (varXj^0.5))
  ##     'make intercept
  Intercept <- meanyj - meanxj * betacoeff
  
  ## Make R the pearson produce moment correlation coefficient
  if (varYj==0 | varXj==0){
    corrCoeff <- 0
  }else{
    corrCoeff <- (sumXjYj - ((sumXj * sumYj) / n)) / ((sdxij * sdxik)^0.5)
  }
  
  ## Make sample variances of betacoefficient and intercept
  variancebeta <- (varYj / varXj) * ((1 - (corrCoeff ^ 2)) / n)
  varianceintercept <- (varYj / n) * (1 - corrCoeff) * (2 + ((meanxj ^ 2) * ((1 + corrCoeff) / varXj)))
  sdbeta <- variancebeta^0.5
  sdintercept <- varianceintercept^0.5

  list(betacoeff=betacoeff, # <- Steigung
       Intercept=Intercept,
       sdbeta=sdbeta, # standard deviation
       sdintercept=sdintercept,
       meanxj=meanxj,
       meanyj=meanyj,
       corrCoeff=corrCoeff) # <- pearson correlation koeffizient ;
  
}
--8<------------------------schnapp------------------------->8---

-- 
Snoopy (on being house-trained with a rolled-up newspaper): 
It does tend however to give one a rather distorted view of the press!