http://blog.gmane.org/gmane.comp.lang.r.ecology hourly 1 1901-01-01T00:00+00:00 http://gmane.org/img/gmane-25t.png http://gmane.org http://comments.gmane.org/gmane.comp.lang.r.ecology/3069 <pre>Hi all, I trying develop a code to implement a recommendation of Wilson et al (1991). In they text formula: "... LnA' = fitted mean ln abundance; sigma = fitted standar deviation of ln abundance. In ranked-abundance terms: lnAi = lnA' + sigma * fi^-1 * (S-i+0.5)/S Where S = number of species; Ai = abundance of ith out of S species; fi^-1 = inverse cumulative distribution function of a norma distribution, i.e. the ln abundance at which the area under the normal curve is the value indicated. ..." I developed this code: rank.lognormal&lt;-function(x){   S &lt;- length(x);   xlog &lt;- log(x);   p &lt;- ppoints(xlog);   mulog &lt;- mean(xlog);   sdlog &lt;- sd(xlog);   fi &lt;- function(y){ qnorm(y, mulog, sdlog)  };   sapply(1:S, FUN=function(i){ exp(1) ^ (mulog + sdlog * fi(p[i]) * (S-i+0.5)/S) }); } bci&lt;-c(25, 24, 22, 21, 18, 17, 15, 14, 14, 13, 13, 12, 12, 11, 11,        10, 9, 8, 7, 6, 6, 6, 6, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4,        3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,        2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1); sum(rank.lognormal(bci)) plot(log(rank.lognormal(bci))) But the sum of fitted rank don't match with the sum of rank abundance. The curve don't like lognormal too. I know rad.lognormal of vegan package, but I'd like follow Wilson recommendation and understand how this procedure works. Thank you. Best, Mario Wilson et al 1991 Methods for fitting dominance/diversity curves - Journal of Vegetation Science. http://onlinelibrary.wiley.com/doi/10.2307/3235896/abstract ................................................... Master student in Ecology Institute of Biology, Dept. Plant Biology, Ecology Lab. State University of Campinas - UNICAMP Campinas, São Paulo, Brazil </pre> Mario José 2012-05-25T13:11:47 http://comments.gmane.org/gmane.comp.lang.r.ecology/3066 <pre>Hello, I'm planning on using a regression model to describe seed set of plants (my response) using some sort of predictor based on temperature. I have a number of temperature variables calculated from the same set of data (hourly temperatures for the growing season, converted to variables such as average temperature, maximum temperature, minimum temperature, degree-days above zero Celsius, degree days above ten Celsius, etc...), and I want to decide which one should be included in my model. I know that I would ideally select one based on "prior knowledge" of the system (e.g. so-called "planned comparisons" or choosing a temperature threshold that is known to be important for the development of seeds), but not much is known about this system. I've been warned against testing the significance of multiple predictors using p-values, unless I use Bonferroni correction (or some equivalent). Unfortunately, using Bonferroni correction would result in something like p = 0.05/7 (for seven different temperature variables); a rather small value for detecting anything! I was wondering whether it would be appropriate to instead use likelihood-based techniques (direct comparisons of log-likelihoods or AIC scores) to compare a series of models using each of the alternative predictors in turn, and choose the most relevant temperature variable (i.e. predictor) based on that. Thoughts on the validity of this approach? Would any adjustments have to be made for multiple comparisons if I used this strategy? Jason Straka University of Victoria [[alternative HTML version deleted]] </pre> J Straka 2012-05-24T22:00:53 http://comments.gmane.org/gmane.comp.lang.r.ecology/3063 <pre>Dear all to find relation between non-normal response with independents variable i use this code: model1&lt;-gam(Clupeidae~s(depth)+s(temperature)+s(salinity),poisson) the result shown is: There were 50 or more warnings (use warnings() to see the first 50) Warning messages: 1: In dpois(y, mu, log = TRUE) : non-integer x = 2.079542 Error in cat(list(...), file, sep, fill, labels, append) : argument 2 (type 'list') cannot be handled by 'cat' what is meaning it? do i allow use it for continue? thanks </pre> Mahnaz Rabbaniha 2012-05-24T07:11:06 http://comments.gmane.org/gmane.comp.lang.r.ecology/3055 <pre>Hello R-sig-ecology group, I was wondering if anyone is aware of an example where PCA is used as a predictive model? A community analysis example might be to predict the PC values of a sample given its community composition. I had thrown this question up on a statistics forum ( http://stats.stackexchange.com/questions/28916/can-empirical-orthogonal-function-eof-analysis-be-used-as-a-predictive-model) but have gotten hardly any response. I imagined that there are some folks here that would have some insight into this problem. Many thanks, Marc [[alternative HTML version deleted]] </pre> Marc Taylor 2012-05-23T07:17:35 http://comments.gmane.org/gmane.comp.lang.r.ecology/3053 <pre>Dear all for finding the relation between this variable i use gam ( in base of non-normality) [1] "temperature" "salinity" "depth" "Clupeidae" model&lt;-gam(Clupeidae~s(temperature)+s(salinity)+s(depth)) the result: Parametric coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 0.40156 0.05467 7.345 2.12e-10 *** Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: Edf Ref.df F p-value s(temperature) 4.321 5.326 1.470 0.207 s(salinity) 2.893 3.569 1.540 0.204 s(depth) 1.000 1.000 0.015 0.903 R-sq.(adj) = 0.1 Deviance explained = 18.9% GCV score = 0.28203 Scale est. = 0.25109 n = 84 this is my question;in base of R-sq. or others parameters,is this model robust ifor analyses? also in others section this parameter is less than0.1 thanks all </pre> Mahnaz Rabbaniha 2012-05-21T16:33:00 http://comments.gmane.org/gmane.comp.lang.r.ecology/3052 <pre>Dear list, I am working on a dendrochronological investigation with the use of a linear mixed-effect model. I identified significant main and interactions effects in parameter estimation of fixed-effects Has anyone a hint on how I can generate or receive results and statements like: main/interaction effect of predictor A/AB explained X% of variation of growth as response variable. thanks CH P.S.: thanks also for all replies concerning my previous questions. this is a great community. </pre> C Hess 2012-05-19T13:26:47 http://comments.gmane.org/gmane.comp.lang.r.ecology/3051 <pre>In species distribution modeling where one uses a large sample of background points to capture background variation in presence\pseudo-absence or use\available models (0\1 response) it is frequently recommended that one weight the data so the sum of the absence weights is equal to the sum of presence weights so that the model isn?t swamped by an overwhelming and arbitrary number of background points. My concern is I don't know that glm in the R stats package handles this correctly. When weights are set in glm they are ignored in the stepAIC and I'm wondering if this should be the case. For example if I have a data set with 75 presence and 75 absence I fit glm. If I repeat the 75 absences 10 times say and then weigh the absences so total weight of absence equals total weight of presence fitting any individual glm to the two sets give almost the same answer (in coefficients, null and residual deviance but NOT in likelihood extracted using logLik). Then when I use the step command in R to select the best from a series of models I get quite different answers for two reasons. Step uses the likelihood which is not the same in the two models (and I'm not sure why given that they have the same deviance). The second reason is that the step function in R uses AIC = - 2*log L + k * edf where edf is effective degrees of freedom. It seems that if you're setting the weights then the the step function should adjust the degrees of freedom it uses accordingly. I've tried adjusting the stepwise regression by multiplying the penalty term in step by sum(data$Weights)/length(data$Weights) where I generally set the penalty to 2 for AIC or log(nrow(data)) for BIC. Indeed when I use this approach I get much more similar predictions from the two models but because the likelihood is different I don't always get exactly the same predictions from the two models. I was wondering if there is a statistical explanation for why the answers should be different and if not, why I can't get the correct answer. I tried the same thing on STATA this afternoon and it was interesting to note that STATA did give the same likelihood for the two models and the same answer for stepwise selection. I'm aware that GLM in the standard stats package in R is intended to use weights only to specify an integer number of binomial trials under the assumption that the response is proportion of success but had read that the only consequence of using non-integer weights with a Bernoulli response was that starting values are difficult to choose and convergence to the maximum likelihood isn't guaranteed. I was willing to accept this possibility but it's seeming that the issues with using glm in this way might run more deep. Another option would be to consider the presence\background as strata with different sampling probabilities and use the survey package svyglm but I'm not sure how to do the automatic model selection in that case. I've noticed that many papers and software programs use R code with non integer weights set in this or a similar way: BIOMOD uses the glm with weights and stepwise model selection without any consideration of the above issues, Phillips and Elith (Logistic Methods for Resource Selection Functions and Presence-Only Species Distribution Models), Guisan, Edwards and Hastie (Generalized linear and generalized additive models in studies of species distributions: setting the scene) though they use a custom stepwise selection based on the deviance rather than the likelihood. Anyway, if anyone has any insight on how to appropriately set weights in glm and use stepwise model selection, I'd really appreciate your recommendations. Marian Talbert [[alternative HTML version deleted]] </pre> Marian K Talbert 2012-05-18T14:26:28 http://comments.gmane.org/gmane.comp.lang.r.ecology/3050 <pre>Dear all, I'm now using the package adehabitatLT to handle with animals movements. Despite things are going ok I cannot superimpose the various animals pathways that plot.ltraj (or simple plot) function generates. My big question is: Is it possible to do such procedure? Alternatively: Does some of you can suggest some way to do this task? Thanks in advance, Cheers, </pre> Conrado Galdino 2012-05-17T19:19:07 http://comments.gmane.org/gmane.comp.lang.r.ecology/3048 <pre>Dear all i have done GAM for data ( they are non- normal and i find regression between Sillaginidae with several hydrological factors),i done : [1] "depth" "temperature" "salinity" "Sillaginidae" &gt; pairs(sc,panel=function(x,y){points(x,y);lines(lowess(x,y))}) but i take this message: Error in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots) : A term has fewer unique covariate combinations than specified maximum degrees of freedom please help me what can i do? or what is meaning this sentence? thanks Mahnaz </pre> Mahnaz Rabbaniha 2012-05-17T16:17:43 http://comments.gmane.org/gmane.comp.lang.r.ecology/3046 <pre>Greetings, I have previously tried the R mixed model mailing list but unfortunately I remain stuck with this question. I certainly hope you will be able to shed some light on this particular problem. I am interested in ecology investigations, however this is just practice data to try and accomplish the modelling parameters. My sincere thanks for an assistance. I have previously attempted to replicate the UN (unstructured), CS (compound symmetry), and AR(1) (Autoregression) covariance structures used in SAS PROC MIXED with the lme function in this nlme package. However, my efforts have fallen short on replicating the *Variance Components (VC)* structure. I have read that it is also known as a diagonal structure. Below I have copied over all the models I have tried and their output with no success. Perhaps someone here will see my error or something I have overlooked. I have attached the data for this particular model. Thanks to all, I certainly cannot thank this help list enough. If you need any further information/clarification, please ask. Cheers, Charles library(nlme) dat=read.table("C:/subset.csv" ,sep=",",header=TRUE, na.strings=".") attach(dat) dat34=dat[Group %in% c("3", "4"),] attach(dat34) liver34=within(dat34, { Group=factor(Group) Event_name=factor(Event_name) Died=factor(Died) ID=factor(ID) }) attach(liver34) ##What is should be from SAS #CV var #Type 3 Tests of Fixed Effects #Effect NumDF DenDF F Value Pr &gt; F #Group 1 22 0.25 0.6244 #Died 1 22 6.55 0.0179 #Group*Died 1 22 4.43 0.0470 fit.1=lme(var~Group*Died, random=~1|ID, data=dat34) anova(fit.1, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 101 227.58700 &lt;.0001 #Group 1 22 0.18320 0.6728 #Died 1 22 3.63388 0.0698 #Group:Died 1 22 3.04103 0.0951 fit.2=lme(var~Group*Died, data=dat34, random=~1|ID/Died) anova(fit.2, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 101 77.99004 &lt;.0001 #Group 1 22 1.46275 0.2393 #Died 1 22 5.84535 0.0243 #Group:Died 1 22 3.04103 0.0951 fit.3=lme(var~Group*Died, random=list(ID=pdSymm(~Event_name)), data=dat34) anova(fit.3, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 101 273.10918 &lt;.0001 #Group 1 22 0.69692 0.4128 #Died 1 22 1.43316 0.2440 #Group:Died 1 22 5.74399 0.0255 fit.4=lme(var~Group*Died, random=list(ID=pdSymm(~Group)), data=dat34) anova(fit.4, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 101 235.13889 &lt;.0001 #Group 1 22 0.15878 0.6941 #Died 1 22 3.83253 0.0631 #Group:Died 1 22 3.01222 0.0966 fit.5=lme(var~Group*Died, random=list(ID=pdSymm(~Group)), data=dat34, weights=varIdent(form=~1|Event_name)) anova(fit.5, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 101 277.16705 &lt;.0001 #Group 1 22 0.23901 0.6298 #Died 1 22 3.99283 0.0582 #Group:Died 1 22 3.23135 0.0860 fit.6=lme(var~Group*Died, random=list(ID=pdSymm(~Group)), data=dat34, weights=varIdent(form=~1|Event_name)) anova(fit.6, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 101 277.16705 &lt;.0001 #Group 1 22 0.23901 0.6298 #Died 1 22 3.99283 0.0582 #Group:Died 1 22 3.23135 0.0860 fit.7=lme(var~(Group*Died), random=list(ID=pdCompSymm(~Died)), data=dat34) anova(fit.7, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 101 85.83799 &lt;.0001 #Group 1 22 1.60624 0.2183 #Died 1 22 4.71795 0.0409 #Group:Died 1 22 2.65379 0.1175 fit.8=lme(var~(Group*Died), data=dat34, random=~1|ID, corr=corSymm()) anova(fit.8, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 101 119.54403 &lt;.0001 #Group 1 22 4.58972 0.0435 #Died 1 22 8.01715 0.0097 #Group:Died 1 22 5.27470 0.0315 fit.9=lme(var~(Group*Died), data=dat34, random=list(ID=pdDiag(~Group*Died)), corr=corSymm(, ~1|ID)) # Error in lme.formula(var ~ (Group * Died), data = dat34, random = list(ID = pdDiag(~Group * : # nlminb problem, convergence error code = 1 # message = iteration limit reached without convergence (9) fit.10=lme(var~(Group*Died), data=dat34, random=list(ID=pdDiag(~Group*Died)), corr=NULL) anova(fit.10, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 101 93.90211 &lt;.0001 #Group 1 22 1.75311 0.1991 #Died 1 22 6.84379 0.0158 #Group:Died 1 22 3.11458 0.0915 fit.11=lme(var~Group*Died, data=dat34, random=list(ID=pdDiag(~Group*Died))) anova(fit.11, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 101 93.90211 &lt;.0001 #Group 1 22 1.75311 0.1991 #Died 1 22 6.84379 0.0158 #Group:Died 1 22 3.11458 0.0915 fit.12=lme(var~Group*Died, data=dat34, random=list(ID=pdDiag(~Event_name))) anova(fit.12, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 101 87.33040 &lt;.0001 #Group 1 22 1.09661 0.3064 #Died 1 22 5.46329 0.0289 #Group:Died 1 22 2.94589 0.1001 fit.13=lme(var~Group*Died, data=dat34, random=list(ID=pdDiag(~Group))) anova(fit.13, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 101 77.99004 &lt;.0001 #Group 1 22 1.46275 0.2393 #Died 1 22 5.84535 0.0243 #Group:Died 1 22 3.04103 0.0951 fit.14=lme(var~Group*Died, data=dat34, random=list(ID=pdDiag(~Died))) anova(fit.14, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 101 85.83800 &lt;.0001 #Group 1 22 1.60624 0.2183 #Died 1 22 4.71795 0.0409 #Group:Died 1 22 2.65379 0.1175 fit.15=lme(var~Group*Died, data=dat34, random=~1|ID, corr=corCompSymm()) anova(fit.15, type="marginal", adjustSigma=F) #same as fit.13 fit.16=lme(var~Group*Died, data=dat34, random=~1|ID/Event_name) anova(fit.16, type="marginal", adjustSigma=F) #same as fit.13 ####################### _______________________________________________ R-sig-ecology mailing list R-sig-ecology-0bNBQ1PAWB4BXFe83j6qeQ&lt; at &gt;public.gmane.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology </pre> Charles Determan Jr 2012-05-16T16:19:01 http://comments.gmane.org/gmane.comp.lang.r.ecology/3039 <pre>Rich: one approach is to fit the data to a gamma distribution, and take the mean and variance from that. I've found Bolker's book, Ecological Models &amp; Data in R, to be helpful in understanding the gamma distribution. If you're not familiar with Newman et al. Estimating mean and variance for environmental studies with below detection limit observations Water Resources Bulletin 25, pp. 905-916, it is mandatory reading. Below is some R code provided by a statistician, Jim Baldwin, for estimating [NO3] mean and variance in an experiment involving three tree species and two levels of plant diversity. NO3 concentrations have many of the properties you describe for metal concentrations. The code has provisions for dealing with left-censored data, but the exact details escape me at the moment. --S eth # Estimating gamma distribution parameters given the data is # left-censored at 0.04. # Rather than estimate the shape and scale parameters and then # determining the mean and variance, the mean and variance are # found directly by re-parametrizing the gamma density function. # Get original data setwd("c:\\documents and settings\\jbaldwin\\desktop\\seth") g = read.csv("Task1404_Treats_3-4 NO3 min set.csv") # Define function returning the log of the likelihood logL = function(z,x) { shape = z[1]^2/z[2] # Shape parameter from mean (z[1]) and variance (z[2]) scale = z[2]/z[1] # Scale parameter from mean (z[1]) and variance (z[2]) # Non-censored data values xx = x[(!is.na(x)) &amp; (x &gt; 0.04)] # Number of censored data values n = length(x[(!is.na(x)) &amp; (x &lt;= 0.04)]) # Return log of the likelihood n*pgamma(0.04,shape=shape,scale=scale,log.p=TRUE) + sum(dgamma(xx,shape=shape,scale=scale,log=TRUE)) } # MLE function for mean and variance mle.censored = function(species,treatment) { # Extract data for specified species and treatment x = g$NO3[g$KEYSPECI==species &amp; g$TREATMNT==treatment] # Get maximum likelihood estimates mle.g = optim(c(1,1),logL,method = "L-BFGS-B",lower=0.0001,upper=Inf, hessian=TRUE,x=x,control=list(fnscale=-1)) # Calculate estimate of covariance matrix mle.g$VarCovMatrix = -solve(mle.g$hessian) # Save a copy of the data mle.g$x = x # Return the data frame with all of the results mle.g } # Get estimates for each combination of species and treatment mle13 = mle.censored(1,3) mle14 = mle.censored(1,4) mle23 = mle.censored(2,3) mle24 = mle.censored(2,4) mle33 = mle.censored(3,3) mle34 = mle.censored(3,4) # Collect all terms # Maximum likelihood estimates of population means mu.mle = c(mle13$par[1],mle14$par[1],mle23$par[1],mle24$par[1], mle33$par[1],mle34$par[1]) # Maximum likelihood estimates of population variances var.mle = c(mle13$par[2],mle14$par[2],mle23$par[2],mle24$par[2], mle33$par[2],mle34$par[2]) # Estimate of standard error of estimates of population means mu.se = c(mle13$VarCov[1,1],mle14$VarCov[1,1],mle23$VarCov[1,1], mle24$VarCov[1,1],mle33$VarCov[1,1],mle34$VarCov[1,1])^0.5 # Estimate of standard error of estimates of population variances var.se = c(mle13$VarCov[2,2],mle14$VarCov[2,2],mle23$VarCov[2,2], mle24$VarCov[2,2],mle33$VarCov[2,2],mle34$VarCov[2,2])^0.5 # Species and treatment id id = c("s1 t3","s1 t4","s2 t3","s2 t4","s3 t3","s3 t4") # Plot means and approximate 95% confidence intervals plot(c(min(mu-1.96*mu.se),max(mu+1.96*mu.se)),c(0.5,6.5),type="n", xlab="Mean NO3 concentration",ylab="",axes=FALSE) axis(1) axis(2,c(1:6),id,las=1) box() for (i in c(1:6)) { points(mu[i],i,pch=16) lines(c(mu[i]-1.96*mu.se[i],mu[i]+1.96*mu.se[i]),c(i,i)) } This electronic message contains information generated by the USDA solely for the intended recipients. Any unauthorized interception of this message or the use or disclosure of the information it contains may violate the law and subject the violator to civil or criminal penalties. If you believe you have received this message in error, please notify the sender and delete the email immediately. [[alternative HTML version deleted]] </pre> Bigelow, Seth W -FS 2012-05-16T12:59:44 http://comments.gmane.org/gmane.comp.lang.r.ecology/3032 <pre> The water chemistry data of metal concentrations are not normally distributed (based on Q-Q plots) and are not improved by transformation (log10, sqrt, cubic root). For the 30 metal species the percentage of zeros ranges from none (10 metals) to 48.6; average 5.6. Most metals are at very low concentrations with infrequent spikes which might be very high. Those with fewer zeros are not a concern, but I'd like your thoughts on 1) at what percentage do the number of zeros become a concern and 2) how to characterize and model these data. Rich </pre> Rich Shepard 2012-05-15T18:15:39 http://comments.gmane.org/gmane.comp.lang.r.ecology/3031 <pre>Dear all, I have one response variable and several explanatory variables, and I want to know the contribution of each explanatory variable to the response variable. Path analysis and partial correlation analysis seem to be suitable for my aim. But I have three questions as follows: 1) My response data is count data. I have analyzed the relationship between response data and explanatory data using Poisson regression, and have found out these significant factors. Before regression, I have stanstardized explanatory table using scale() function. Can I identify the contribution of each significant factor according to the partial regression coefficient directly? 2) Can path analysis and partial correlation analysis be based on Poisson regression (i.e. calculate the path coefficient or partial correlation coefficent without changing the link function of "glm") 3) I am not sure the difference of these two method, so I don't know which method is more suitable for my analysis. Thanks very much for any suggestion. Guojun May 14th, 2012 [[alternative HTML version deleted]] </pre> lgj200306 2012-05-14T19:36:06 http://comments.gmane.org/gmane.comp.lang.r.ecology/3028 <pre>Hello all, When I perform a cross validation with the function "crossval" of the package cocorresp (for a co-correspondence analysis), I get the following error: Error in La.svd(x, nu, nv) : error code 1 from Lapack routine 'dgesdd' I think the problem is not in the function, as other people has had the same error for other types of model (see, for example, http://comments.gmane.org/gmane.comp.lang.r.ecology/2608). From this and other posts, I suppose it is a convergence problem. I am working with presence-absence species matrices, and the problem persists even if I transform the matrices with Hellinger transformation (although not at the same site). Does anyone know how to solve, or go around, this problem? Any advice would be most welcome. Thanks, Duarte </pre> Duarte Viana 2012-05-14T13:06:31 http://comments.gmane.org/gmane.comp.lang.r.ecology/3025 <pre>Dear all, I have two species-by-site matrices and I'd like to combine them in one matrix. The rows (n=20) of each matrix present plants 1-20. So, the new matrix should have 40 rows. As these matrices have some shared species and some exclusive species, the new matrix should have the species from the two matrix, but when certain species that doesn't occur for example in Matrix 1, its value will be 0. I've tried the functions merge, cbind and rbind but I wasn't able to do it. Please, see the following example to make it clear ;-) By using for example the following matrices: matrix1 &lt;- matrix(c(1,2,3, 11,12,13), nrow = 2, ncol=3, byrow=TRUE, dimnames = list(c("plant1", "plant2"), c("sp1", "sp2", "sp3"))) matrix2 &lt;- matrix(c(1,2,3, 0,21,9), nrow = 2, ncol=3, byrow=TRUE, dimnames = list(c("plant3", "plant4"), c("sp1", "sp2", "sp4"))) I'd like to combine them like this: new.matrix &lt;- matrix(c(1, 2, 3, 0, 11, 12, 13, 0, 1, 2, 0, 3, 0,21,0,9), nrow = 4, ncol=4, byrow=TRUE, dimnames = list(c("plant1", "plant2", "plant3", "plant4"), c("sp1", "sp2", "sp3", "sp4"))) Once I have 200 columns and 80 rows, I'm pretty sure that R can help me to do it automatically ;-) All the best, Thiago. </pre> Thiago Gonçalves-Souza 2012-05-14T02:50:14 http://comments.gmane.org/gmane.comp.lang.r.ecology/3022 <pre>Hi,i would like use CA code in my data it is a metrics with 16 column ( abundance of fish larva) and 84 raw, this raw in three groups: Coralline,non-coralline and creek habitat . the data are non-normal and i have to use non-parametric analyses ,i want to find similarity or dissimilarity between these habitat. i used the nmds in past but the stress wasn't good fitness please guide me thanks Mahnaz [[alternative HTML version deleted]] </pre> Mahnaz Rabbaniha 2012-05-11T04:09:32 http://comments.gmane.org/gmane.comp.lang.r.ecology/3008 <pre> Dear all, I have tried to use the function bca() on PCA results (dudi.pca) : pca.ade4&lt;-dudi.pca(batrapair[,2:20], center=T, scale=T, scannf=T); 2; betweenclass&lt;-bca(pca.ade4,fac, scannf = TRUE);2; and I get this error: Error in if (nf &lt;= 0) nf &lt;- 2 : missing value where TRUE / FALSE is required If you have any idea to solve this please let me know. I attached my script and here is what batrapair looks like SCIENTIFIC_NAME BIO_1 BIO_2 BIO_3 BIO_4... 15513 Batrachoseps incognitus 119 159 59 3844 15585 Batrachoseps incognitus 138 130 64 2449 15649 Batrachoseps incognitus 132 149 63 3093 15700 Batrachoseps incognitus 128 144 64 2872 15701 Batrachoseps incognitus 130 141 65 2719 19094 Batrachoseps minor 135 164 58 3954 19095 Batrachoseps minor 137 162 59 3791 19096 Batrachoseps minor 139 146 61 3147... Thank you very much, Julian _______________________________________________ R-sig-ecology mailing list R-sig-ecology-0bNBQ1PAWB4BXFe83j6qeQ&lt; at &gt;public.gmane.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology </pre> Julian WITTISCHE 2012-05-09T22:09:42 http://comments.gmane.org/gmane.comp.lang.r.ecology/3006 <pre>Dear all, I have two questions about data normality. I used stepwise multiple regression to determine which variables contributed to tree growth, and want to built a model to explain tree growth. Sample size is about 50 tree species, I think it is not a large sample size, and some variables are not normal distribution. 1. Do I have to transform them to normal distributions before I perform multiple regression? 2. Two variables can not transform to normal distributions although I used some methods (e.g log, sqrt, boxcoxfit), what should I do for the two variables? Thank you and best wishes! Yong Shen [[alternative HTML version deleted]] </pre> Yong Shen 2012-05-06T08:54:42 http://comments.gmane.org/gmane.comp.lang.r.ecology/3003 <pre>As implemented in the literature, Bray-Curtis dissimilarity measures are used to look at pairwise comparisons. So, how dissimilar are 2 plots in their composition, for example. I'm wondering if anyone has run across or implemented a method to look at BC dissimilarity or otherwise that has two different twists to it 1) Negative values are possible. 2) The index is for an n-wise measure of dissimilarity. For example, instead of asking how dissimilar 2 plots are, what is the dissimilarity between 3 plots. I've been trying to play with vegdist to see if I can implement this, and have been rolling my own functions that work well with just presence-abscence data. But I'm unclear of how to do so for continuous data. Thanks! -Jarrett </pre> Jarrett Byrnes 2012-05-04T21:46:44 http://comments.gmane.org/gmane.comp.lang.r.ecology/3000 <pre>Hello R users, I may be thinking about this all wrong. If I am, please let me know. My question is about bootstrapping gam models. I want to know if there is a way to produce bootstrapped smoothers in a gam model that has more than one significant explanatory term. I know how to achieve this when the model only has one term. The code that I am using to do so is below. It produces a nice graph of y values in response to my explanatory variable (x), along with bootstrapped smoothers in a grey color and the smoother produced by the gam package in black. In my mind this is a very useful graph to help determine the "realness" of the fitted smoother that gam produces. I should note a thanks to Charles Geyer who had code on his STATS 5601 website that has helped me immensely thus far. x &lt;- SI1 ## my explanatory variable y &lt;- NO3.sp3 ## my dependent/response variable model &lt;- gam(y ~ s(x, bs = "cr")) ## I know I can use other smoother types. I am still experimenting with this aspect. plot(x, y) curve(predict(model, newdata = data.frame(x = x)), add = TRUE) n &lt;- length(x) nboot &lt;- 100 for (i in 1:nboot) { k.star &lt;- sample(n, replace = TRUE) x.star &lt;- x[k.star] y.star &lt;- y[k.star] model.star &lt;- gam(y.star ~ s(x.star, bs = "cr")) curve(predict(model.star, data.frame(x.star = x)), add = TRUE, col = "grey") } points(x, y, pch = 16) curve(predict(model, newdata = data.frame(x = x)), add = TRUE, lwd = 2) But now lets say I have a model with two explanatory variables ( x and z) and for the sake of this discussion, both terms are significant (i.e., the data is better explained by having them in the model). The model in this case would be: new.model &lt;- gam(y ~ s(x, bs = "cr")+ s(z, bs = "cr")) I know by using "plot (new.model, resid = T, pch = 16)" I can see the smoothers and 95% CI produced by the gam package with the data values overlayed. Is there a way to produce bootstrapped smoothers for each term from the same model (new.model) so that I can visualize them in the same way as the plot function allows me? I thought of having separate single-term models for each explanatory variable (Y~s(x) and y~s(z)), but doing so is not correct since removing one term from a model causes the parameter values of the other term to change. Any suggestions on the bootstrapping code that would produce such a graphical output would be appreciated, especially since I am new to the bootstrapping coding. I hope that is clear and thanks in advance for any help that anyone can offer. </pre> Basil Iannone 2012-05-04T01:36:06 http://comments.gmane.org/gmane.comp.lang.r.ecology/2998 <pre>Hello, I am attempting to run zero-inflated and hurdle models using the pscl() package, but keep encountering problems with running my fully-parameterized model and I am not sure why. My data consist of measurements of seed set (the response, which is zero-inflated and overdispersed) for two treatments and a control (B,F, and O) at three elevations (low, mid, and high). At each elevation I also measured average temperature (AvgTemp) and abundance of insects (N). My goal is to construct a model that best describes seed set using the four predictors above (and interactions, if necessary). I find that I have no problems running a partial model with only my two categorical predictors as follows (just an example using hurdle(),but results are the same for zeroinfl()): H1B &lt;- hurdle(AVInt ~ Elev + Treatment | Elev + Treatment, dist ="negbin", link = "logit", data = erytseeds) It’s when I try to add the two continuous predictors that things begin to get messy; H4B &lt;- hurdle(AVInt ~ Elev*Treatment*AvgTemp*N | Elev*Treatment*AvgTemp*N) This results in the error: Warning message: In sqrt(diag(object$vcov)) : NaNs produced Error in optim(fn = loglikfun, gr = gradfun, par = c(start$count, start$zero, : non-finite value supplied by optim Even when I run the model using one continuous predictor and one categorical predictor, or both continuous predictors, I end up with NAs in the summary table of the count, but not in the zero-inflated portion of the model. For example: Call: zeroinfl(AVInt ~ AvgTemp+N | AvgTemp+N, data = erytseeds, dist ="negbin", link = "logit") Pearson residuals: Min 1Q Median 3Q Max -0.5807 -0.4566 -0.2791 -0.2791 5.5813 Count model coefficients (negbin with log link): Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -3.770120 0.452221 -8.337 &lt;2e-16 *** AvgTemp 0.404550 0.173458 2.332 0.0197 * N 0.001211 NA NA NA Log(theta) 1.308387 NA NA NA Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -4.814e+01 8.115e+00 -5.932 2.99e-09 *** AvgTemp 4.039e+00 7.158e-01 5.642 1.68e-08 *** N 3.010e-03 4.863e-04 6.189 6.04e-10 *** --- Theta = 3.7002 Number of iterations in BFGS optimization: 31 Log-likelihood: -456.4 on 7 Df Call: zeroinfl(AVInt ~ Treatment+N | Treatment+N, data = erytseeds, dist ="negbin", link = "logit") Pearson residuals: Min 1Q Median 3Q Max -0.6608 -0.4684 -0.3565 -0.2956 6.1730 Count model coefficients (negbin with log link): Estimate Std. Error z value Pr(&gt;|z|) (Intercept) 0.9555134 NA NA NA TreatmentF 0.1407240 0.1844811 0.763 0.446 TreatmentO 0.0430166 0.1647207 0.261 0.794 N 0.0009843 NA NA NA Log(theta) 1.3049021 NA NA NA Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(&gt;|z|) (Intercept) 0.0685048 1.0056156 0.068 0.94569 TreatmentF -0.3822091 0.3631982 -1.052 0.29264 TreatmentO -1.1624582 0.3551596 -3.273 0.00106 ** N 0.0009120 0.0004773 1.911 0.05603 . --- Theta = 3.6873 Number of iterations in BFGS optimization: 31 Log-likelihood: -458.1 on 9 Df Reading other posts, I know that my values are all &gt;= 0, and there are no missing values that could be causing this, but there ARE lots of zeroes. My sample size is decent (about 35 observations per treatment, per elevation). So my questions are essentially as follows: 1) Why am I getting all these NaNs, and why can’t I run the full (saturated) model, including continuous as well as categorical predictors and interactions? How can I begin the process of model simplification if I can’t get the full model to run? 2) If this is a result of the many zeroes in my data set, is there a way to work around this? Any advice would be greatly appreciated. Jason Straka University of Victoria </pre> J Straka 2012-05-03T21:56:56 Search the mailing list at Gmane query http://search.gmane.org/?group=$group=gmane.comp.lang.r.ecology