0
$\begingroup$

I am trying to find the best fit between an species dataset and prevailing climatic conditions, in order to be able to predict the environmental conditions from the species dataset (paleoclimate research).

I have 15 species(sp1-sp15), expressed as relative amounts (some are 0). I have done some data exploration in excel, and have seen that I get good fits using the ratio sp2/(sp1+sp2).

I have done multiple linear regression on this dataset, but none give a correlation as good as the ratio I stumbled on. The ultimate goal is of course to test whether this ratio I stumbled upon is yields the best model to explain changes in the climate variable.

I would thus like to create the ratios: sum(spj)/sum(spj), where spj refers to one of the species. Here, both the numerator and denominator can be any possible combination of all species.

Can I generate a formula model to include a ratio? Can I then select the best model using regbsubsets?

EDIT Based on answers below, this will not be possible for all 15 variables. How would it be possibel for 5 variables (I have seen this published before, using R).

EDIT On stackoverflow (where I posted the same question), I got the comment to use the function (I), to write expressions inside formulae. However, if I use

leapsMAT<-regsubsets(x1 ~ I(1+ (y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8 + y9 + y10 + y11 + y12+ y13 + y14 +y15 )/(y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8 + y9 + y10 + y11 + y12+ y13 + y14 + y15)), force.in= FTFALSE,nvmax=15,data=my.data, nbest=1)

the regsubsets doesn't recognize the variables as variables, of course. How do I implement the I(function), or is there another way to deal with this?

At the moment I have obtained the same linear regression using two methods:

library(leaps)
attach(my.data)
FTFALSE<-c(FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE)
leapsMAT<-regsubsets(x1 ~ 1+ y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8 + y9 + y10 + y11 + y12+ y13 + y14 + y15 , force.in= FTFALSE,nvmax=15,data=my.data, nbest=1)

And, using a longer method specified here (implemented from somewhere else on the internet):

allModelsList <- apply(regMat, 1, function(x) as.formula (paste(c("x1 ~ 1", namevar2[x]),collapse=" + " )))
allModelsList
warnings()

#Calculating the model
my.data
allModelsResults <- lapply(allModelsList, function(x) lm(x, data=my.data))
allModelsResults

dfCoefNum   <- ldply(allModelsResults, function(x) as.data.frame(t(coef(x))))

dfStdErrors <- ldply(allModelsResults, function(x) as.data.frame(t(coef(summary(x))[, "Std. Error"])))


dftValues   <- ldply(allModelsResults, function(x) as.data.frame(t(coef(summary(x))[, "t value"])))

dfpValues   <- ldply(allModelsResults, function(x) as.data.frame(t(coef(summary(x))[, "Pr(>|t|)"]))) 
dfpValues
(warnings)
# rename DFs so we know what the column contains
names(dfStdErrors) <- paste("se", names(dfStdErrors), sep=".")
names(dftValues) <- paste("t", names(dftValues), sep=".")
names(dfpValues) <- paste("p", names(dfpValues), sep=".")

# p-value for overall model fit
calcPval <- function(x){
    fstat <- summary(x)$fstatistic
    pVal <- pf(fstat[1], fstat[2], fstat[3], lower.tail = FALSE)
    return(pVal)
}

# Before creating ONE data frame with all important entries,
# we need to compute some more indices 
NoOfCoef <- unlist(apply(regMat, 1, sum))
R2       <- unlist(lapply(allModelsResults, function(x)
                          summary(x)$r.squared))
    adjR2    <- unlist(lapply(allModelsResults, function(x)
                              summary(x)$adj.r.squared))
RMSE     <- unlist(lapply(allModelsResults, function(x)
                          summary(x)$sigma))
fstats   <- unlist(lapply(allModelsResults, calcPval))



# now we can combine all the data into one data frame
results <- data.frame( model = as.character(allModelsList),
                       NoOfCoef = NoOfCoef,
                       dfCoefNum,
                       dfStdErrors,
                       dftValues,
                       dfpValues,
                       R2 = R2,
                       adjR2 = adjR2,
                       RMSE = RMSE,
                       pF = fstats  )
results[1:20,]
# round the results
results[,-c(1,2)] <- round(results[,-c(1,2)], 3)
results

model.maxRadj<-which(results$adjR2 == max(results$adjR2), arr.ind = TRUE)
maxRadj<-results[model.maxRadj,]

Many thanks in advance! Please let me know if more information is required.

$\endgroup$
3
$\begingroup$

You wrote

I would thus like to create the ratios: sum(spj)/sum(spj), where spj refers to one of the species. Here, both the numerator and denominator can be any possible combination of all species.

No. You don't want to do this. 15 species may be combined in $2^{15}$ ways. That is in your numerator and denominator. You will then have $2^{30}$ ratios to sort through. That can't be a good strategy.

In addition, your dependent variable is "prevailing environmental conditions", it is unlikely that this is a continuous variable. A single condition might be continuous (e.g. temperature, humidity) but over all conditions are not likely to be. Even just combining two conditions yields noncontinuous data: Hot dry, cool wet etc. Therefore, I don't think you want linear regression; you probably want multinomial logistic.

Third, you write

The ultimate goal is of course to test whether this ratio I stumbled upon is yields the best model to explain changes in the climate variable.

OK, I am not a climate scientist, but the relatively abundance of species cannot explain climate although the reverse is possible. It isn't hot because there are more reptiles; it may be that there are more reptiles because it is hot.

$\endgroup$
  • $\begingroup$ Thank you for your reply. You indicate that finding the best ratio is not a good strategy, but is it impossible? I will explore the multinomial logistic option, although our goal is to predict continuous variables. $\endgroup$ – Deborah Jun 12 '13 at 11:40
  • 3
    $\begingroup$ It would mean adding $2^{30}$ explanatory variables, which is approximately one bilion ($10^9$) variables. This is to all intends and purposes impossible $\endgroup$ – Maarten Buis Jun 12 '13 at 11:47
  • $\begingroup$ Ok, thank you for this clear answer. Would it be possible for a dataset with less species? I have been searching on this problem for a number of days now, it would be good to learn how to do it, for future reference. $\endgroup$ – Deborah Jun 12 '13 at 11:56
  • $\begingroup$ This technique, examining the correlations of all possible ratios with an environmental parameter has been published in the palaeoclimate field before, I know of an example with 5 species. Can somebody indicate how I could write a function that could do this? $\endgroup$ – Deborah Jun 12 '13 at 12:21
  • $\begingroup$ What continuous variable do you want to predict? Also, just because it's published doesn't mean it's a good thing to do. If you generated completely random data, this method would still find highly significant results: It's a way of generating type I errors $\endgroup$ – Peter Flom - Reinstate Monica Jun 12 '13 at 12:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.