am carrying out research to look at differences in fish density and fish species richness when using two different underwater visual census methods. My data was originally count data but then converted to fish density and so decided to use an offset using the following formula:   model1 &lt;- glm(g_den ~ method + site + depth + offset(log(area)), poisson)  My 3 predictor variables are method, site and depth which I ordered as factors when I input them.  My response variables are grouper species richness, grouper density and the same for other fish groups.   Despite using the offset i keep getting the same error, could anyone advise me why this is?  In dpois(y, mu, log = TRUE) : non-integer x = 0.037500
78313
I am carrying out research to look at differences in fish density and fish species richness when using two different underwater visual census methods. My data was originally count data but then typically this is changed to fish density but I have still decided to use a Poisson GLM, which I hope is right.   model1 &lt;- glm(g_den ~ method + site + depth, poisson)   My 3 predictor variables are method, site and depth which I ordered as factors when I input them.   My response variables are grouper species richness, grouper density and the same for other fish groups. I am aware that density is not an integer and it is numerical data e.g 1.34849. I am now however getting this error:      In dpois(y, mu, log = TRUE) : non-integer x = 0.037500   I've been reading up and many people suggest using an offset, is this the most advisable thing to do?
12043
I have built a binomial regression model using the equation: logit_1 &lt;- glm(past_adl_binary~., family = binomial(link = &quot;cloglog&quot;),                 data = newdata)  The outcome of the model is as follows:  Can anyone help me interpret the outcome of this model? I understand that all of the coefficients are significant apart from sunangle. Can I use this equation (100*exp(parameter estimate)-1) to calculate the increase/ decrease in percent probability of my y-variable (1 or 0) changing based on my predictor variables?  For example, for habitat_categoryONDEEP: 100*(exp(-0.173)-1) = -15.886.  Would this suggest that a habitat category of ONDEEP decreases the probability of my y-variable being 1 by 15.9%, compared to the reference category?
122822
I am using a binomial regression model for presence/absence, with 20 independent variables to test. The data has x and y coordinates and I would like to understand how can I take into account the spatial autocorrelation.  I already studied the correlation between the variables and run the same model for 1000 different samples (I have a big dataset that allows me to do this) to understand the distribution of each parameter and check for variables that might be introducing problems in my model.   glm_model &lt;- glm(PA ~ Var1 + Var2 + Var3 + Var4 + Var5,family=binomial(link=logit))  However I believe I also need to account for spatial autocorrelation. I saw that there is a package that might help me (spdep), however I am not sure I completely understand if I can use my model or not. My question is what are my options ?
88765
I have a question about the interpretation of the coefficients of an interaction between continuous and categorical variable. here is my model:  model_glm3=glm(cog~lg_hag+race+pdg+sex+as.factor(educa)+(lg_hag:as.factor(educa)),                 data=base_708)  Coefficients:                          Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)               21.4836     2.0698  10.380  &lt; 2e-16 *** lg_hag                     8.5691     3.7688   2.274  0.02334 *   raceblack                 -8.4715     1.7482  -4.846 1.61e-06 *** racemexican               -3.0483     1.7073  -1.785  0.07469 .   racemulti/other           -4.6002     2.3098  -1.992  0.04687 *   pdg                        2.8038     0.4268   6.570 1.10e-10 *** sexfemale                  4.5691     1.1203   4.078 5.15e-05 *** as.factor(educa)2         13.8266     2.6362   5.245 2.17e-07 *** as.factor(educa)3         21.7913     2.4424   8.922  &lt; 2e-16 *** as.factor(educa)4         19.0179     2.5219   7.541 1.74e-13 *** as.factor(educa)5         23.7470     2.7406   8.665  &lt; 2e-16 *** lg_hag:as.factor(educa)2 -21.2224     6.5904  -3.220  0.00135 **  lg_hag:as.factor(educa)3 -19.8083     6.1255  -3.234  0.00129 **  lg_hag:as.factor(educa)4  -8.5502     6.6018  -1.295  0.19577     lg_hag:as.factor(educa)5 -17.2230     6.3711  -2.703  0.00706 ***   lets say the equation of the model is:  E[cog] = a + b1(lg_hag) + b2(educa2*lg_hag) + b3(educa3*lg_hag) + b4(educa4*lg_hag)  + b5(pdg, centered) + other covars, where  b1 = difference in cog  with higher lg_hag among lowest education (coded as 1) b1 + b2 = difference in cog with higher lg_hag among middle education (coded as 2) b1 + b3 = difference in cog with higher lg_hag among high education (coded as 3) b1 + b3 = difference in cog with higher lg_hag among very high education (coded as 4) b5 = difference in cog with each unit increase in pdg   My question is: if my interpretation is right, how to construct confidence intervals for each effect estimate of interactions (e.g: b1+b2) from the confidence intervals of b1 and b2.
19133
I recently discovered how to model exposures over time using the log of (e.g.) time as an offset in a Poisson regression.  I understood that the offset corresponds to having time as covariate with coefficient 1.  I'd like to better understand the difference between using time as an offset or as a normal covariate (therefore estimating the coefficient). In which situation should I want to use one method or the other?  UPGRADE: I don't know if it is interesting, but I ran a validation of the two methods using randomly split data repeated 500 times and I noticed that using the offset method leads to a larger test error.
49240
I'm working on a VAR model and am doing a log transformation of the raw data.   x = log(x)   After differencing, running various tests and running a VAR(2) model, I add the predicted differences cumulatively to the original time series (logged) like this:  prd.diffs &lt;- as.ts(prd$mean$d.x1) U &lt;- tail(undiff,2) cumu &lt;- cumsum(prd.diffs) + U[1]   However, while the output now looks as it should, it is still in logged values.  My question is as follows:  How does one return the logged output to original values? I'm under the impression that using exp() will not give the correct results, and having tried some variations it hasn't worked for me.   Does the fact that it is a VAR($p$) model (and possibly what lag it has) influence the method, and how should that be implemented to get the correct output?  P.S. I realise that the above code does not constitute a replicable example as I didn't deem that necessary for the question, but should that be desired I can add some more code.
102173
I built a ZI model and it is producing predicted values that are from a very small range when compared to the observed values. Plus it does not produce any zeros. See the fitted vs. observed graph below.     My model structure is given below. My response is a count variable and the predictor is continuous water level measurements (scaled).   m.u.meanwl.gu &lt;- zeroinfl(GU ~ MeanWLScaled | MeanWLScaled, data = uluabat_wl_fg, dist = "negbin", link = "logit")   I need help to figure out what's wrong and fix it. Thanks in advance.
37989
I have following simple X and Y vectors:   &gt; X [1] 1.000 0.063 0.031 0.012 0.005 0.000 &gt; Y [1] 1.000 1.000 1.000 0.961 0.884 0.000 &gt;  &gt; plot(X,Y)     I want to do regression using log of X. To avoid getting log(0), I try to put +1 or +0.1 or +0.00001 or +0.000000000000001 :   &gt; summary(lm(Y~log(X))) Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :    NA/NaN/Inf in 'x' &gt; summary(lm(Y~log(1+X)))  Call: lm(formula = Y ~ log(1 + X))  Residuals:        1        2        3        4        5        6  -0.03429  0.22189  0.23428  0.20282  0.12864 -0.75334   Coefficients:             Estimate Std. Error t value Pr(&gt;|t|)   (Intercept)   0.7533     0.1976   3.812   0.0189 * log(1 + X)    0.4053     0.6949   0.583   0.5910   --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 0.4273 on 4 degrees of freedom Multiple R-squared:  0.07838,   Adjusted R-squared:  -0.152  F-statistic: 0.3402 on 1 and 4 DF,  p-value: 0.591  &gt; summary(lm(Y~log(0.1+X)))  Call: lm(formula = Y ~ log(0.1 + X))  Residuals:        1        2        3        4        5        6  -0.08099  0.20207  0.23447  0.21870  0.15126 -0.72550   Coefficients:              Estimate Std. Error t value Pr(&gt;|t|)   (Intercept)    1.0669     0.3941   2.707   0.0537 . log(0.1 + X)   0.1482     0.2030   0.730   0.5058   --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 0.4182 on 4 degrees of freedom Multiple R-squared:  0.1176,    Adjusted R-squared:  -0.103  F-statistic: 0.5331 on 1 and 4 DF,  p-value: 0.5058  &gt; summary(lm(Y~log(0.00001+X)))  Call: lm(formula = Y ~ log(1e-05 + X))  Residuals:        1        2        3        4        5        6  -0.24072  0.02087  0.08796  0.13872  0.14445 -0.15128   Coefficients:                Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)     1.24072    0.12046  10.300 0.000501 *** log(1e-05 + X)  0.09463    0.02087   4.534 0.010547 *   --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 0.1797 on 4 degrees of freedom Multiple R-squared:  0.8371,    Adjusted R-squared:  0.7964  F-statistic: 20.56 on 1 and 4 DF,  p-value: 0.01055  &gt;  &gt; summary(lm(Y~log(0.000000000000001+X)))  Call: lm(formula = Y ~ log(1e-15 + X))  Residuals:         1         2         3         4         5         6  -0.065506  0.019244  0.040983  0.031077 -0.019085 -0.006714   Coefficients:                Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)     1.06551    0.02202   48.38 1.09e-06 *** log(1e-15 + X)  0.03066    0.00152   20.17 3.57e-05 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 0.04392 on 4 degrees of freedom Multiple R-squared:  0.9903,    Adjusted R-squared:  0.9878  F-statistic: 406.9 on 1 and 4 DF,  p-value: 3.565e-05   The output is different in all cases. What is the correct value to put to avoid log(0) in regression? What is the correct method for such situations.  Edit: my main aim is to improve prediction of the regression model by adding log term, i.e.: lm(Y ~ X + log(X))
35890
I am trying to generate two continuous predictors and one continuous outcome from mvrnorm. Now, I am trying to set up both standardized slopes at 0.3. This is my R code. n.tr=10 mu = c(3, 3, 3) Sigma = cbind(c(0.91,0,0.3),               c(0 0.91,0.3),               c(0.3,0.3,0.91) dat =  as.data.frame(mvrnorm(n=n.tr, mu=mu, Sigma=Sigma,                       empirical=TRUE)) colnames(dat) = c('z', 'w','y')  I want both of the standardized slopes to be 0.3, and 0 correlation between z and w. I also want the variance of y to be 1, so the residual variance of y is 0.91. If i want the diagonal of the matrix to be residual variance, how can I set up the matrix?
53093
I'm working on a multiple logistic regression in R using glm. The predictor variables are continuous and categorical. An extract of the summary of the model shows the following:  Coefficients:                Estimate Std. Error z value Pr(&gt;|z|) (Intercept)   2.451e+00  2.439e+00   1.005   0.3150 Age           5.747e-02  3.466e-02   1.658   0.0973 . BMI          -7.750e-02  7.090e-02  -1.093   0.2743 ... --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   Confidence intervals:                    2.5 %       97.5 % (Intercept)  0.10969506 1.863217e+03 Age          0.99565783 1.142627e+00 BMI          0.80089276 1.064256e+00 ...   Odd ratios:                   Estimate Std. Error   z value Pr(&gt;|z|) (Intercept)  1.159642e+01  11.464683 2.7310435 1.370327 Age          1.059155e+00   1.035269 5.2491658 1.102195 B            9.254228e-01   1.073477 0.3351730 1.315670 ...   The first output shows that $Age$ is significant. However, the confidence interval for $Age$ includes the value 1 and the odds ratio for $Age$ is very close to 1. What does the significant p-value from the first output mean? Is $Age$ a predictor of the outcome or not?
99542
I am running a GLM, using the function glm.nb (pscl package) trying to figure how what could influence a particular trait in several locations and years. The output as follow (with slight modification and removing the things beyond this question)            Estimate Std. Error z value Pr(&gt;|z|)               state2        -0.163190   0.807307  -0.202  0.00081               state3        -0.530588   0.783023  -0.678  0.00302            state4         0.942953   0.737328   1.279  0.00094               year2         -1.759102   0.733214  -2.399  0.01643            year3         -0.633870   0.662903  -0.956  0.33897       I understand the output as such: compared to year1, year 2 is statistically different, but year3 is not different from year1. but what about year2 and year3?  For the state factor, state2, state3 and state4 are all statistically different than state1. so I need a post-hoc test between state2 state3 and state4, am I right?  thanks
62050
I have data from a survey experiment in which respondents were randomly assigned to one of four groups:  &gt; summary(df$Group)        Control     Treatment1     Treatment2     Treatment3              59             63             62             66    While the three treatment groups do vary slightly in the stimulus applied, the main distinction that I care about is between the control and treatment groups. So I defined a dummy variable Control:  &gt; summary(df$Control)      TRUE FALSE         59   191    In the survey, respondents were asked (among other things) to choose which of two things they preferred:   &gt; summary(df$Prefer)       A   B  NA's      152  93   5    Then, after receiving some stimulus as determined by their treatment group (and none if they were in the control group), respondents were asked to choose between the same two things:  &gt; summary(df$Choice)   A    B  149  101    I want to know if the being in one of the three treatment groups had an effect on the choice that respondents made in this last question. My hypothesis is that respondents who received a treatment are more likely to choose A than B.       Given that I am working with categorical data, I have decided to use a logit regression (feel free to chime in if you think that's incorrect). Since respondents were randomly assigned, I am under the impression that I shouldn't necessarily need to control for other variables (e.g. demographics), so I have left those out for this question. My first model was simply the following:  &gt; x0 &lt;- glm(Product ~ Control + Prefer, data=df, family=binomial(link="logit")) &gt; summary(x0)  Call: glm(formula = Choice ~ Control + Prefer, family = binomial(link = "logit"),      data = df)  Deviance Residuals:      Min       1Q   Median       3Q      Max   -1.8366  -0.5850  -0.5850   0.7663   1.9235    Coefficients:                     Estimate Std. Error z value Pr(&gt;|z|)     (Intercept)           1.4819     0.3829   3.871 0.000109 *** ControlFALSE         -0.4068     0.3760  -1.082 0.279224     PreferA              -2.7538     0.3269  -8.424  &lt; 2e-16 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  (Dispersion parameter for binomial family taken to be 1)      Null deviance: 328.95  on 244  degrees of freedom Residual deviance: 239.69  on 242  degrees of freedom   (5 observations deleted due to missingness) AIC: 245.69  Number of Fisher Scoring iterations: 4   I am under the impression that the intercept being statistically significant is not something that holds interpretable meaning. I thought perhaps that I should include an interaction term as follows:  &gt; x1 &lt;- glm(Choice ~ Control + Prefer + Control:Prefer, data=df, family=binomial(link="logit")) &gt; summary(x1)  Call: glm(formula = Product ~ Control + Prefer + Control:Prefer, family = binomial(link = "logit"),      data = df)  Deviance Residuals:      Min       1Q   Median       3Q      Max   -2.5211  -0.6424  -0.5003   0.8519   2.0688    Coefficients:                                  Estimate Std. Error z value Pr(&gt;|z|)     (Intercept)                         3.135      1.021   3.070  0.00214 **  ControlFALSE                       -2.309      1.054  -2.190  0.02853 *   PreferA                            -5.150      1.152  -4.472 7.75e-06 *** ControlFALSE:PreferA                2.850      1.204   2.367  0.01795 *   --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  (Dispersion parameter for binomial family taken to be 1)      Null deviance: 328.95  on 244  degrees of freedom Residual deviance: 231.27  on 241  degrees of freedom   (5 observations deleted due to missingness) AIC: 239.27  Number of Fisher Scoring iterations: 5   Now the respondents status as in a treatment group has the expected effect. Was this a valid set of steps? How can I interpret the interaction term ControlFALSE:PreferA? Are the other coefficients still the log odds?
70806
Can anyone tell me how to interpret the 'residuals vs fitted', 'normal q-q', 'scale-location', and  'residuals vs leverage' plots?  I am fitting a binomial GLM, saving it and then plotting it.
33007
I have count data (demand/offer analysis with counting number of customers, depending on - possibly - many factors). I tried a linear regression with normal errors, but my QQ-plot is not really good. I tried a log transformation of the answer: once again, bad QQ-plot.  So now, I'm trying a regression with Poisson Errors. With a model with all significant variables, I get:  Null deviance: 12593.2  on 53  degrees of freedom Residual deviance:  1161.3  on 37  degrees of freedom AIC: 1573.7  Number of Fisher Scoring iterations: 5   Residual deviance is larger than residual degrees of freedom: I have overdispersion.  How can I know if I need to use quasipoisson? What's the goal of quasipoisson in this case? I read this advise in "The R Book" by Crawley, but I don't see the point nor a large improvement in my case.
19801
I'm trying to predict a binary outcome using 50 continuous explanatory variables (the range of most of the variables is $-\infty$ to $\infty$). My data set has almost 24,000 rows. When I run glm in R, I get:  Warning messages:   1: glm.fit: algorithm did not converge   2: glm.fit: fitted probabilities numerically 0 or 1 occurred    I've read the other responses that suggest perfect separation might be occurring, but I'm confident that isn't the case in my data (though quasi-complete separation could exist; how can I test to see if that's the case?). If I remove some variables, the "did not converge" error might go away. But that's not always what happens.  I tried using the same variables in a bayesglm function and got the same errors.  What steps would you take to figure out exactly what's going on here? How do you figure out which variables are causing the problems?
41895
I try to reproduce with optim the results from a simple linear regression fitted with glm or even nls R functions. The parameters estimates are the same but the residual variance estimate and the standard errors of the other parameters are not the same particularly when the sample size is low. I suppose that this is due differences in the way the residual standard error is calculated between Maximum Likelihood and Least square approaches (dividing by n or by n-k+1 see bellow in the example). I understand from my readings on the web that optimization is not a simple task but I was wondering if it would be possible to reproduce in a simple way the standard error estimates from glm while using optim.   Simulate a small dataset  set.seed(1) n = 4 # very small sample size ! b0 &lt;- 5 b1 &lt;- 2 sigma &lt;- 5 x &lt;- runif(n, 1, 100) y =  b0 + b1*x + rnorm(n, 0, sigma)    Estimate with optim  negLL &lt;- function(beta, y, x) {     b0 &lt;- beta[1]     b1 &lt;- beta[2]     sigma &lt;- beta[3]     yhat &lt;- b0 + b1*x     likelihood &lt;- dnorm(y, yhat, sigma)     return(-sum(log(likelihood))) }  res &lt;- optim(starting.values, negLL, y = y, x = x, hessian=TRUE) estimates &lt;- res$par     # Parameters estimates se &lt;- sqrt(diag(solve(res$hessian))) # Standard errors of the estimates cbind(estimates,se)       &gt; cbind(estimates,se)       estimates         se b0     9.016513 5.70999880 b1     1.931119 0.09731153 sigma  4.717216 1.66753138   Comparison with glm and nls  &gt; m &lt;- glm(y ~ x) &gt; summary(m)$coefficients             Estimate Std. Error   t value    Pr(&gt;|t|) (Intercept) 9.016113  8.0759837  1.116411 0.380380963 x           1.931130  0.1376334 14.030973 0.005041162 &gt; sqrt(summary(m)$dispersion) # residuals standard error [1] 6.671833 &gt;  &gt; summary(nls( y ~ b0 + b1*x, start=list(b0 = 5, b1= 2)))  Formula: y ~ b0 + b1 * x  Parameters:    Estimate Std. Error t value Pr(&gt;|t|)    b0   9.0161     8.0760   1.116  0.38038    b1   1.9311     0.1376  14.031  0.00504 ** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 6.672 on 2 degrees of freedom   I can reproduce the different residual standard error estimates like this :   &gt; # optim / Maximum Likelihood estimate &gt; sqrt(sum(resid(m)^2)/n) [1] 4.717698 &gt;  &gt; # Least squares estimate (glm and nls estimates) &gt; k &lt;- 3 # number of parameters &gt; sqrt(sum(resid(m)^2)/(n-k+1)) [1] 6.671833
62061
This is a fairly complex question so I will attempt to ask it in a fairly basic manner.   I have data on the abundance of 99 different species of estuarine macroinvertebrate species and the sediment mud content (0 - 100 %) in which each observation was obtained. I have a total of 1402 observations for each species (i.e. a massive dataset).   Here is a subset of the raw data for one species to give you an idea of the data I'm working with (if I had 10 reputation points I'd upload a plot of real raw data):  Abundance: 10,14,10,3,3,3,3,4,5,5,0,0,0,0,0,0,0,0,0,0,0,0,0,6,6,6,0,0,0,0,12,0,0,0,34,0,0 Mud %:     0.9,4,2,10,13,14,6,5,5,7,22,27,34,37,47,58,54,70,54,80,90,65,56,7,8,34,67,54,32,1,57,45,49,4,78,65,45,35   The primary aim of my research is to determine an "optimum mud % range" (e.g. 15 - 45 %) and "distribution mud % range" (e.g. 0 - 80 %) for each of the 99 invertebrate species.  As you can see the abundance data for the above species contains a significant number of zero values. Although this significantly skews any sort of model that I run on the data (i.e. GLM, GAM), even if I model the non-zero data only, the model for certain species does not fit the data at all well.  So, my question is: what would be the best, most robust way to determine an "optimum" and "distribution" mud range for each species, given that responses vary significantly between species?     Just to clarify - the above data is a hypothetical example to give you an idea of how messy the abundance (that is count) data can be for a given species.  Regarding the poisson regression approach: I'm considering conducting a two-step GLM or GAM approach for each species; Step (1) uses logistic regression to model the "probability of presence"  for a given species over the mud gradient - using presence/absence data. This obviously takes into account the zero counts; and Step (2) models the "maximum abundance" over the mud gradient - using presence only count data. This step gives me an idea of the species typical response to mud where they DO occur. What are your thoughts on this approach?  I have R code for both steps for one particular species. Heres the code:       ## BINARY  aa1&lt;-glm(Bin~Mud,dist=binomial,data=Antho) xmin &lt;- ceiling(min(Antho$Mud))     xmax &lt;- floor(max(Antho$Mud)) Mudnew &lt;- seq(from=xmin, to=xmax, by=0.1) pred.dat &lt;- data.frame(Mudnew) names(pred.dat) &lt;- "Mud" pred.aa1 &lt;- data.frame(predict.glm(aa1, pred.dat, se.fit=TRUE, type="response")) pred.aa1.comb &lt;- data.frame(pred.dat, pred.aa1) names(pred.aa1.comb) plot(fit ~ Mud, data=pred.aa1.comb, type="l", lwd=2, col=1, ylab="Probability of presence", xlab="Mud content (%)", ylim=c(0,1))   ## Maximum abundance   aa2&lt;-glm(Maxabund~Mud,family=Gamma,data=antho)  xmin &lt;- ceiling(min(antho$Mud))      xmax &lt;- floor(max(antho$Mud))  Mudnew &lt;- seq(from=xmin, to=xmax, by=0.1)  pred.dat &lt;- data.frame(Mudnew)  names(pred.dat) &lt;- "Mud"  pred.aa2 &lt;- data.frame(predict.glm(aa2, pred.dat, se.fit=TRUE, type="response"))  pred.aa2.comb &lt;- data.frame(pred.dat, pred.aa2)  names(pred.aa2.comb)  plot(fit ~ Mud, data=pred.aa2.comb, type="l", lwd=2, col=1, ylab="Maximum abundance per 0.0132 m2", xlab="Mud content (%)")  AIC(aa2)   My question is: for step (2); will the model code need to be altered (i.e. family=) depending on the shape of each species abundance data, if so, would I just need to inspect a scatter plot of the raw presence only abundance data to confirm the use of a certain function? and how would the code be written for a certain species exhibiting a certain response/functional form?