Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to compare slopes in R

I am performing an ANCOVA so as to test what is the relationship between body size (covariate, logLCC) and different head measures (response variable, logLP) in each sex (cathegorical variable, sexo). I got the slopes for each sex in the lm and I would like to compare them to 1. More specifically, I would like to know if the slopes are significantly higher or less than 1, or if they are equal to 1, as this would have different biological meanings in their allometric relationships. Here is my code:

#Modelling my lm#
> lm.logLP.sexo.adu<-lm(logLP~logLCC*sexo, data=ADU)
> anova(lm.logLP.sexo.adu)
Analysis of Variance Table

Response: logLP
         Df Sum Sq Mean Sq  F value    Pr(>F)    
logLCC        1 3.8727  3.8727 3407.208 < 2.2e-16 ***
sexo          1 0.6926  0.6926  609.386 < 2.2e-16 ***
logLCC:sexo   1 0.0396  0.0396   34.829 7.563e-09 ***
Residuals   409 0.4649  0.0011                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

#Obtaining slopes#
> lm.logLP.sexo.adu$coefficients
 (Intercept)       logLCC        sexoM logLCC:sexoM 
  -0.1008891    0.6725818   -1.0058962    0.2633595 
> lm.logLP.sexo.adu1<-lstrends(lm.logLP.sexo.adu,"sexo",var="logLCC")
> lm.logLP.sexo.adu1
 sexo logLCC.trend         SE  df  lower.CL  upper.CL
 H       0.6725818 0.03020017 409 0.6132149 0.7319487
 M       0.9359413 0.03285353 409 0.8713585 1.0005241

Confidence level used: 0.95

#Comparing slopes#
> pairs(lm.logLP.sexo.adu1)
 contrast   estimate         SE  df t.ratio p.value
 H - M    -0.2633595 0.04462515 409  -5.902  <.0001

#Checking whether the slopes are different than 1#
#Computes Summary with statistics      
> s1<-summary(lm.logLP.sexo.adu)
> s1

Call:
lm(formula = logLP ~ logLCC * sexo, data = ADU)

Residuals:
 Min       1Q   Median       3Q      Max 
-0.13728 -0.02202 -0.00109  0.01880  0.12468 

Coefficients:
         Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.10089    0.12497  -0.807     0.42    
logLCC        0.67258    0.03020  22.271  < 2e-16 ***
sexoM        -1.00590    0.18700  -5.379 1.26e-07 ***
logLCC:sexoM  0.26336    0.04463   5.902 7.56e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Residual standard error: 0.03371 on 409 degrees of freedom
Multiple R-squared:  0.9083,    Adjusted R-squared:  0.9076 
F-statistic:  1350 on 3 and 409 DF,  p-value: < 2.2e-16

#Computes t-student H0: intercept=1. The estimation of coefficients and   their s.d. are in s1$coefficients
> t1<-(1-s1$coefficients[2,1])/s1$coefficients[2,2]
#Calculates two tailed probability
> pval<- 2 * pt(abs(t1), df = df.residual(lm.logLP.sexo.adu), lower.tail = FALSE)
> print(pval)
[1] 3.037231e-24

I saw this whole process in several threads here. But all that I can understand is that my slopes are just different from 1. How could I check that they are greater or smaller than 1?

EDITED

Solved!

#performs one-side test H0=slope bigger than 1
pval<-pt(t1, df = df.residual(lm.logLP.sexo.adu), lower.tail = FALSE)
#performs one-side test H0=slope smaller than 1
pval<-pt(t1, df = df.residual(lm.logLP.sexo.adu), lower.tail = TRUE)

Also, tests should be performed in single-sex models.

like image 207
Alicia Avatar asked Nov 08 '22 13:11

Alicia


1 Answers

How could I check that they are greater or smaller than 1?

As in this post, this post, and as your in question, you can make Wald test which you compute by

t1<-(1-s1$coefficients[2,1])/s1$coefficients[2,2]

Alternatively, use the vcov and coef function to make the code more readable

fit <- lm.logLP.sexo.adu
t1<-(1-coef(fit)[1])/vcov(fit)[1, 1]

The Wald test gives you t-statistics which can be used to make both a two-sided or one-sided test. Thus, you can drop the abs and set the lower.tail argument according to which tail you want to test in.

like image 172
Benjamin Christoffersen Avatar answered Nov 14 '22 21:11

Benjamin Christoffersen