[R] Comparing slopes in two linear models

Arthur Weiss arthurw at lncc.br
Thu Feb 12 13:00:36 CET 2009


Hi everyone,

I have a data frame (d), wich has the results of mosquitoes trapping in
three different places.
I suspect that one of these places (Local=='Palm') is biased by low
numbers and will yield slower slopes in the variance-mean regression over
the areas. I wonder if these slopes are diferents.

I've looked trought the support list for methods for comparing slopes and
found the three explained above:

#d - whole datafrmae
d$Local <- as.factor(d$Local)
dHT <- subset(d, Local!='Palm')
dP <- subset(d, Local!='Palm')


1 - Single dataframe procedure:

fit1 <- lm(log(varCDI) ~ Local + I(log(CDI)) + Local:I(log(CDI)), data=d)
fit2 <- lm(log(varCDI) ~ Local + I(log(CDI)) + Local:I(log(CDI)),
data=dHT)
a1 <- anova(fit1)
a2 <- anova(fit2)

a1
Local:I(log(CDI))   2   1.183   0.592    7.8741 0.0005574 ***
Residuals         152  11.421   0.075

a2
Local:I(log(CDI))   1  0.061   0.061   0.8156 0.36859
Residuals         102  7.643   0.075


>>>> Concluding that the places are different


2- proposed by Spencer Graves:

Spencer <- function(df1, df2){
 fit1 <- lm( log(df1$varCDI)~log(df1$CDI) )
 fit2 <- lm( log(df2$varCDI)~log(df2$CDI) )
 s1 <- summary(fit1)$coefficients
 s2 <- summary(fit2)$coefficients
 db <- (s2[2,1]-s1[2,1])

#### TEST DEGREES OF FREEDOM
 df <- (fit1$df.residual+fit2$df.residual)

### Pooled Variance ####

 sd <- sqrt(s2[2,2]^2+s1[2,2]^2)

 td <- db/sd
 2*pt(-abs(td), df)
}
Spencer(dP, dHT)
[1] 0.0001717730




3- I wrote one based on Zar (1999) and Kleinbaum (1987) procedure:

Zar <- function(df1, df2){
 fit1 <- lm( log(df1$varCDI)~log(df1$CDI) )
 fit2 <- lm( log(df2$varCDI)~log(df2$CDI) )
 s1 <- summary(fit1)$coefficients
 s2 <- summary(fit2)$coefficients
 db <- (s1[2,1]-s2[2,1])

#### TEST DEGREES OF FREEDOM
 n1 <- fit1$df.residual
 n2 <- fit2$df.residual
 n <- n1 + n2

### Pooled Variance ####

 ssx1 <- var(df1$CDI) * (n1+1)
 ssx2 <- var(df2$CDI) * (n2+1)
 RSS1 <- anova(fit1)[2,2]
 RSS2 <- anova(fit2)[2,2]
 sspvm <- ( RSS1 + RSS2 ) / n
 sd1_2 <- sqrt( sspvm/ssx1 + sspvm/ssx2 )

 td <- db/sd1_2
 2*pt(-abs(td), n)
}
Zar(dP, dHT)
[1] 0.4040566


Is there anything wrong with the function I wrote?
Assuming that there is nothing wrong with the function I wrote, I'd like
to know why R dos not follow this procedure.
Why does the variation over the X-axis is not considered for the pooled
variance?
Perhaps a simpler question: how is calculated the coeficient standart
error wich is printed by 'summary'?

Thanks,
Arthur




More information about the R-help mailing list