[R-sig-eco] PhD position in ecology and statistics at NTNU, Norway
(the deadline for this is the 31st of October, sorry about the short notice! If you apply, please mention that you saw this on r-sig-ecology) We have a vacancy for a position as a PhD research fellow at the Department of Mathematical Sciences, for someone interested in statistical modelling and applying it to problems of biodiversity and sustainability. I can promise the chance to work with a cross-disciplinary team: there will also be a student working on the ecological modelling, and we also have a very active group of students working with citizen science data. The successful applicant will work as part of an inter-disciplinary team to develop methods to aid in the understanding of the effects of water resource management on biodiversity. They will join a team working on statistical methods to apply ecological models to messy data. The project will focus on developing the statistical models to integrate different data from different sources (e.g. different surveys, and reports from citizen scientists) to model the biodiversity of Norwegian rivers and use these models to predict the impacts of resource management decisions. Duties of the position • make a workflow to bring data in from a variety of sources • develop models, based on recent work at NTNU, to model the distributions of riverine species • work with other members of the project to look at the effects of water resource management, such as the effects of dams. Required selection criteria The PhD-position's main objective is to qualify for work in research positions. The qualification requirement is that you have completed a master’s degree or second degree (equivalent to 120 credits) with a strong academic background in statistics or equivalent education with a grade of B or better in terms of NTNU’s grading scale. If you do not have letter grades from previous studies, you must have an equally good academic foundation. If you are unable to meet these criteria you may be considered only if you can document that you are particularly suitable for education leading to a PhD degree. The appointment is to be made in accordance with the regulations in force concerning State Employees and Civil Servants and national guidelines for appointment as PhD, post doctor and research assistant. Required selection criteria: • ability to work in a cross-disciplinary team • interest in applying statistical methods to ecological problems • good knowledge of statistical modeling methods Preferred selection criteria • experience with Bayesian methods • skills in hierarchical modelling • good written and oral English language skills More details here (including how to apply, and the formalities): https://www.jobbnorge.no/en/available-jobs/job/194344/phd-research-fellow-at-the-department-of-mathematical-sciences If any potential applicant have questions, they should feel free to ask me. Bob -- Bob O'Hara Institutt for matematiske fag NTNU 7491 Trondheim Norway Mobile: +47 915 54 416 Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] CI for stratified Cox PH model
Is there any reason you can#'t use confint()? e.g. test1 <- list(time=c(4,3,1,1,2,2,3), status=c(1,1,1,0,1,1,0), x=c(0,2,1,1,1,0,0), sex=c(0,0,0,0,1,1,1)) confint(coxph(Surv(time, status) ~ x + strata(sex), test1)) Another trick is to use predict(), with new data at all of hte levels you are interested in: pred1 <- list(x=c(0:2, 0:2), sex=c(0,0,0,1,1,1)) predict(coxph(Surv(time, status) ~ x + strata(sex), test1), newdata=pred1, type="lp" , se.fit=TRUE) (it doesn't provide CIs, but mean +/-1.96s.e. should work OK for the linear predictor scale) Bob On 05/06/18 18:47, Bertolo, Andrea wrote: Hi everyone, I have a doubt about the way to calculate 95% CI for coefficients in the stratified Cox proportional hazard models and your help is welcomed . Say that I have a variable of interest Imi and a stratifying variable UV in a interaction model (since the interaction between Imi and UV is of interest for me and the interaction model has a better fit than the no-interaction model: library(survival) model1 <- coxph(Surv(start,stop, Status.time) ~ Imi + Imi:UV + strata(UV) + cluster(ID), weights = NB_Event, data=Data.unfold) Whereas it is pretty straightforward to calculate the coefficients (and associated HR) for each combination of Imi and UV, I am not sure about how to calculate the associated CI (note that, of course, I got the CI for the estimate of "Imi:UV" from the output of model1). Is it correct to calculate separately a model for each UV level and use the CI for the Imi variable to get the CI for the two levels of UV (see below) ? # 2 models (one per UV level) data.Low <- subset(Data.unfold,UV=="low") model2.1 <- coxph(Surv(start,stop, Status.time) ~ Imi + cluster(ID), weights = NB_Event, data=data.Low) data.High <- subset(Data.unfold,UV=="high") model2.2 <- coxph(Surv(start,stop, Status.time) ~ Imi + cluster(ID), weights = NB_Event, data=data.High) Alternatively, is there a way to get the CI directly from the output of the stratified model ? Many thanks Andrea Bertolo Université du Québec à Trois-Rivières 3351, bd des Forges C.P.500, Trois-Rivières (Québec) Canada G9A 5H7 ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Institutt for matematiske fag NTNU 7491 Trondheim Norway Mobile: +49 1515 888 5440 Homepage: http://www.ntnu.edu/employees/bob.ohara Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] correlating three time-series in R
You can pass the columns to ccf() directly: df <- data.frame(x=rnorm(6), y=rnorm(6)) ccf(df$x, df$y) print(ccf(df$x, df$y)) You should probably also check the time series task view: <https://cran.r-project.org/web/views/TimeSeries.html>, in particular the zoo package, to see what can be done with irregular time series. But with 6 data points I'd be surprised if you have the power to detect anything that doesn't jump out when you simply plot the data. Bob On 26/07/17 11:07, Tania Bird wrote: I have three data sets of abundances through time for plants, insects and reptiles. There are 6 samples over a ten year period (all taxa sampled at the same time). I recognise this is a small data set for time series. I would like to correlate the time series to see if a) increases in abundance of one taxon are correlated to another, and b) to see if the correlation between plants:insects is greater than plants:reptiles. I thought to use the cross-correlation function in R e.g. ccf(insects, reptiles) Currently the data is in one dataframe with time as one column and abundance of each taxa is the next three columns. How do I convert the data to a time.series format as given in the R example? How can I compare the two ccf outputs? Thanks Tania Tania Bird MSc *"There is a sufficiency in the world for man's need but not for man's greed" ~ Mahatma Gandhi* [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara NOTE: this email will die at some point, so please update you records to bob.oh...@ntnu.no Institutt for matematiske fag NTNU 7491 Trondheim Norway Mobile: +49 1515 888 5440 Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] mixed effects in ordination (?)
This doesn't sound like an ordination problem: it's more similar to a standard mixed model regression. There are several ways to approach this, for example Ben Bolker has some (old) notes here: <https://rpubs.com/bbolker/3336> I'd suggest you start by running the analyses on each trait individually, as that way you can make sure you have the model sorted correctly. The step up to the multivariate model just needs some work understanding what the relevant R package wants. HTH Bob On 12/01/17 14:12, Martin Weiser wrote: Dear friends, Could you please help me with analysis? I am afraid that it involves crossed random effects in the mixed-effect constrained ordination setting, so to say. Goal: Show an effect of the species trait (single one) and treatment (four levels, quantitative scale) on parameters. Trait x treatment interaction is possible. If possible, show changes through time. Data: Individuals of 7 species, subjected to 4 treatment levels (fully factorial) - from 6 to 12 individuals in each combination. Each individual scored 4 times (same times: 1st wk, 2nd wk, 3rd wk, 4th wk). Several (10) parameters scored every time on each individual. What I did: In order to avoid multiple testing (the parameters are likely to be correlated to each other), I decided to use multivariate analysis (RDA). I am by far much more accustomed to vegan than ade4, so excuse me if I use some "veganisms". Predictors: time, trait, treatment (possibly with interactions), conditioned on individual identity to avoid treating records from the same plant as independent. Variance partitioning. Here comes the problem: how to set permutations in order to be able to report p_vals (some people just are not happy without them)? Since individuals of the same species share the same trait value, maybe the right way is to: shift observations within individual (if time is among predictors for the particular model) and permute trait value among species. Is it so? Is this treatment of the pseudoreplication at the species level (i.e. only in the significance testing, not in the ordination per se) ok? I also tried to use different approach: I averaged all params individual-wise (getting rid of temporal pseudoreplication, but also time effects), and subsequently I averaged the result within treatment x species levels. I assume that I can go for simple free permutations this way? Pity is that this way, I cannot see development in time. And another way: I averaged params for species x treatment x time groups, ignoring interdependence of records from the same individual, hoping that the effect of an individual "dissolves" in the average. Is that meaningful? If yes, what is the appropriate permutation structure in this case? But maybe I miss something and there are better ways how to deal with this problem? Any suggestions (ok: not any, just those made in an attempt to help :-) ) are appreciated. With the best regards, Martin Weiser ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara NEW EMAIL: bob.oh...@ntnu.no NOTE NEW ADDRESS!!! Institutt for matematiske fag NTNU 7491 Trondheim Norway Mobile: +49 1515 888 5440 Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] LMER: applying a random term to only one level of a factor
On 06/04/16 20:12, Aislinn Pearson wrote: Hi, I've tried googling this but haven't been very successful. Essentially, I'd like to know what is the most statistically valid way of dealing with a random term which doesn't apply to every level of fixed-effect factor. I have a mixed effect model that looks like this Disease level <- weight + Flown +(1|DateFlown) Either I flew my insects on a flight mill (which can be thought of as a 'treadmill' for flying insects) or I didn't, hence flown is a two level factor (Yes or No) and I want to understand how this affects the amount of disease in my insect. To get as many replicates as I could on a single day, I had two different banks of flight mills (A & B), each bank containing 16 individual insect treadmills. The insects were randomly assigned to one of the two sets of 16 flight mills. Previous studies tell me there are differences between these two sets of flight mills, so I would like to account for them as a random term in my model. As a practical matter, it's not worth setting a level with two levels as random: you don't gain anything in the analysis and the variance component is really poorly estimated. In practice, this might actually make things cleaner, as you will have to look a bit more at the flight mill effects, so you should get a better feel for what's happening. However, I can't run this in LMER. When I tried I got the error; Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : contrasts can be applied only to factors with 2 or more levels Which I imagine means that one of my factors (i.e. Flown) doesn't include any levels for the random term mill set (i.e. for all unflown insects the value in the mill set column is NA) Is it possible to include this form of experimental design in LMER (the package I know best) or, alternatively, nlme (which I am a lot less accustomed to using)? I can think of two ways of doing this: either set up a factor with three levels (Flight mill A, Flight mill B, Not Flown) or set the Not Flown to one of the flight mill levels. The first way feels less confusing, but you might have to set up some contrasts to estimate the differences. But hopefully your insects will cooperate nicely and make the difference between the flight ills will be much smaller than between flight mills and not flown. HTH Bob I would be really grateful if anyone has any insight. Many thanks Rothamsted Research is a company limited by guarantee, registered in England at Harpenden, Hertfordshire, AL5 2JQ under the registration number 2393175 and a not for profit charity number 802038. [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] AIC in R: back-transforming standardized model parameters (slopes)
-0.38008 0.13722 -2.77 [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology ___________ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] NMDS axes scores
houldn't be used individually for further analysis. I therefore would like to include both of my NMDS site scores as a response into a GLM model simultaneously. Unfortunately, I couldn't find any advice on how to actually do this. I found a couple of papers using NMDS scores in GLMs, but they all seem to use them individually, fitting separate models to each of the ordination axes. I'm a bit at a loss here and any advice is very much appreciated, Conny [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- -- Pokud je tento e-mail součástí obchodního jednání, Přírodovědecká fakulta Univerzity Karlovy v Praze: a) si vyhrazuje právo jednání kdykoliv ukončit a to i bez uvedení důvodu, b) stanovuje, že smlouva musí mít písemnou formu, c) vylučuje přijetí nabídky s dodatkem či odchylkou, d) stanovuje, že smlouva je uzavřena teprve výslovným dosažením shody na všech náležitostech smlouvy. ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] best choice of GLMM for seed set data
” in Tinn-R :-). -- *Mariano Devoto* [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Mehdi Abedi Department of Range Management Faculty of Natural Resources Marine Sciences Tarbiat Modares University (TMU) 46417-76489, Noor Mazandaran, IRAN mehdi.ab...@modares.ac.ir Homepage Tel: +98-122-6253101 Fax: +98-122-6253499 -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Zero infalted model- T50
On 08/07/15 09:56, Mehdi Abedi wrote: Hi, I want to calculate time to 50 percent reduction of seed viability. Mainly researcher calculate P50 using probit model for calculation T50. In some treatments seed viability decline quit fast therefore i have only value 100 and 0 values . Can i calculate T50 for this data using probit?! Could we apply Zero inflated model for this conditions with 0? All the best, Mehdi Time(weed) Seed viability 0 100 8 0 16 0 24 0 32 0 40 0 You can use a probit, but it'll be unhappy. The problem is that the curve is steep compared to the times, so you don't see it. The only real solution is to change the assays, so you observe more times. But I appreciate this might not be practical. I don't think zero inflation will help. There isn't a unique maximum likelihood estimate for this, but I think a sensible estimate would be mid way between the times where you observe 100 and the first 0. If you want to be Bayesian, this would be the posterior mean. The best confidence interval is probably the central 95% of the interval that the T50 is in (i.e. assuming the likelihood is flat). It's probably a little bit wider than it should be, but I doubt that will be a big problem. I used to get data where the data would flip between 0 and 100 over a range of doses. In the end I gave up with those problems and became an ecologist Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] mvabund SIMPER alternative
On 02/06/15 15:38, carvin90 wrote: Hi community, I am new to stats to this level in ecology. I am trying to compare DNA and RNA libraries with thousands of OTUs. I summarized taxa to get the most abundant species, but I can obtain only relative abundances. I was thinking to use SIMPER as I read in several paper to test which species differ the most per station between DNA and RNA based libraries. However I found that SIMPER is a more or less robust test. I was wondering if the manyglm was also an alternative for my question or if you suggest another function. Thank you for your help! Yes, you can use mvabund for OTUs: I've helped a colleague do this (e.g. dx.doi.org/10.1371/journal.pone.0053987). Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Association Routine?
On 02/28/2015 06:48 PM, Alexandre F. Souza wrote:
Dear friends,
I need to write a code to find data using one variable as reference. The
code I wrote, however, is not working and I can't figure it out why. Could
anyone help me?
Imagine a data set with two variables, B and C. Now I have variable A,
which is the same variable as variable B but the data are not in the same
order nor have necessarily the same extension as B (it may be a sample of
B, for example).
I want to find the values of variable C that match each line in variable A
using B as the association criterion. So the code should perform a loop in
which it would take the first line in A, search B until it finds it there,
then copy the corresponding value of C and store it in a new variable D. Do
it until all lines in A have been associated to a C value.
starting with...
df-data.frame(B=sample(letters[1:10],replace=FALSE), C=rnorm(10),
stringsAsFactors=FALSE)
A=letters[1:10]
two thoughts spring to mind:
(a) would merge() do what you want? e.g. df2 -
merge(df,data.frame(A=A), by.x=B, by.y=A), and then extract the
values of C with df2$C[df2$B==f], for example.
(b) sapply(A, function(lt, DF) DF$C[DF$B==lt], DF=df)
R's looping is generally more efficient when it's done internally, so it
will be easier for you if you understand the R mentality, in particular
vectorisation. usually if you have a for() loop, you're not writing R
code efficiently.
Bob
Here is the code I wrote:
# Considering that matrices data.ref and data.assoc have been already read,
containing the
# User-defined number of columns to be associated with A (I imagined that
more than one variable could be associated at once)
col.assoc = 20
# To assure that data will not be in a non-usable data category
ref = as.matrix(data.ref)
assoc = as.matrix(data.assoc)
# Table where results will be stored
# Number of columns = n associated variables plus one column
# Reserved to receive the initial data (example column A)
result = matrix(nrow = nrow(ref), ncol = col.assoc + 1)
# Fulfill the first column of the result table with the original reference
variable
result[,1] = ref[,1]
for (i in 1:nrow(ref)){
for (j in 1:nrow(assoc))
if (ref[i, 1] == assoc[j, 1]){
resultado[i, 2] == assoc[j, 2]
}
}
col = ncol(dados)
Any thoughts?
Thanks in advance,
Alexandre
--
Bob O'Hara
Biodiversity and Climate Research Centre
Senckenberganlage 25
D-60325 Frankfurt am Main,
Germany
Tel: +49 69 7542 1863
Mobile: +49 1515 888 5440
WWW: http://www.bik-f.de/root/index.php?page_id=219
Blog: http://blogs.nature.com/boboh
Journal of Negative Results - EEB: www.jnr-eeb.org
___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Multivariate quasi-bionomial analysis of proportion data?
On 08/02/15 12:27, Amanda Greer wrote: Hi All, I am trying to best analyse a set of foraging ecology data with 10 behaviour categories (DVs) and 3 levels of IV (season, sex, age). The time which an animal spent engaged in a behaviour was recorded and then divided by the total time spent in sight of the observer, so my data are proportional. As is typical, not all animals engaged in all behaviours and there are a large number of zeros in my dataset which is severely over-dispersed. I had initially analysed all the data using the glm function (family = quasibinomial, followed by anova. The intention was then to use the false discovery rate alpha to account for the large number of analyses. However, it was pointed out to me that a multivariate approach might be better so I have been trying to figure out (a) if it's possible to run a quasi-binomial multivariate analysis of proportion data (b) how to go about it. In the R Documentation quasi-binomial family function page ( http://artax.karlin.mff.cuni.cz/r-help/library/VGAM/html/quasibinomialff.html ) it is stated that if multivariate response = TRUE the response matrix should be binary. This seems a pretty straightforward indictment of my idea to run this analysis on my proportion data, but I am wondering why - is this just not possible, or is there a particular package that could help? If anyone could provide me with an answer or some much needed guidance on this topic I would be very grateful. Ignoring the zeroes problem for the moment, I think (quasi-)binomial distributions are a distraction: binomials are based on counts of things (see Petr Keil's post: http://www.petrkeil.com/?p=603). If you're looking at proportions of times, then it might be better to think in terms of gamma distributions, which lead to a beta distribution for the proportion of times spent doing one thing, and a Dirichlet distribution if you have several items (as you do here). Once you have to worry about the zeroes, you need to do something more, for example see this paper: http://onlinelibrary.wiley.com/doi/10./2041-210X.12122/abstract Bob Thanks, Amanda [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] SPEI index
On 12/23/2014 09:57 AM, ONKELINX, Thierry wrote: If the file is tab delimited then you need to use read.table(sep = \t) instead of read.table(sep = ). read.delim() is another option. Or save the file as a .csv and use read.csv (this makes the reading in easier if you have things like spaces in the file). Also, it's worth checking that you don't have a header row and then some empty rows. It's often worth opening up the .txt or .csv file as a text file (e.g. in Notepad) to see if there are any problems. Oh, and the gdata package has a read.xls() function that is another alternative. Bob ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest team Biometrie Kwaliteitszorg / team Biometrics Quality Assurance Kliniekstraat 25 1070 Anderlecht Belgium + 32 2 525 02 51 + 32 54 43 61 85 thierry.onkel...@inbo.be www.inbo.be To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger Brinner The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey Van: R-sig-ecology [mailto:r-sig-ecology-boun...@r-project.org] Namens Manuel Esteban Lucas Borja Verzonden: dinsdag 23 december 2014 9:13 Aan: r-sig-ecology@r-project.org Onderwerp: [R-sig-eco] SPEI index Dear All, I am starting to use R software in order to calculate an aridity index (SPEI). I have installed the spei package and the library. Then, I have created a .xls file with year, month, Precip and Temp variables. The xls file has been exported to txt format (separated by tabulations). Then to try, the only output I have obtained can be see below. Could anyone of you please help me? If you want, I can send you the txt file. All the best, Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 1 NA NA NA NA NA NA NA NA NA NA NA NA 2 NA NA NA NA NA NA NA NA NA NA NA NA Script: library(SPEI) SPEI - read.table(/Users/PALAN.txt, header=TRUE, sep=) str(SPEI) PET - thornthwaite(SPEI$Tm,40.0) spei(SPEI$PREC-PET,1) --- Manuel Esteban Lucas Borja Universidad de Castilla La Mancha Escuela T�cnica Superior de Ingenieros Agr�nomos y de Montes Departamento de Ciencia y Tecnolog�a Agroforestal y Gen�tica Campus Universitario s/n, C.P. 02071, Albacete (Spain) T�lf.; 967599200 ext. 2818 Mail: manuelesteban.lu...@uclm.es Web:http://www.uclm.es/profesorado/manuelestebanlucas/ [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology Disclaimer Bezoek onze website / Visit our websitehttps://drupal.inbo.be/nl/disclaimer-mailberichten-van-het-inbo ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] species richness, GLM and negative values
On 10/29/2014 04:27 PM, Ludovico Frate wrote: Dear all,I'am trying to fit a very simple linear model. I am analyzing the differences in the number of species (DS) found in several permanent plots in two year of observations. Firstly, I have calculated the differences per plot (i.e. number of species in Plot 1 in Time A - number of species in Plot 1 in Time B and so on for all the plots).Secondly, those differences were tested for deviation from zero by means of a linear model M2-lm(DC~1, data = gransasso)summary(M2)E2-residuals(M2)qqnorm(E2, pch = 19, col = blue); qqline(E2, col = red) The qqnorm has shown that residuals were not normally distributed, thus I need to use a GLM. However GLM (poisson family) does not work with negative values (DS has negative values).I've tried to add a constant value to these differences (i.e. +100) but the result is misleading since I am testing for deviation from zero. Do you have any suggestions? Regards,Ludovico Use the Poisson to model the number of species in each sample, so use data like this: Plot Time DS 1 A 4 1 B 7 2 A 0 2 B 1 3 A 23 3 B 7 ... Then you fit the model Mgood - glm(DS ~ Plot + Time, family=poisson()) where Plot and Time are factors (use Plot - factor(Plot), for example, if you need to). You're interested in the Time effect, which is the average difference between the numbers of species in the plots in the different times. The Plot effect controls for different plots having different numbers of species overall. If you look at summary(Mgood), the Time effect is the log of the ratio of the species richnesses in times A and B. It will be written as TimeB, which means it's log(E.B/E.A) (where E.A and E.B are the expected species richnesses at times A and B). So, for example, an estimate of 0.4 would mean that at Time B there are exp(0.4)=1.49 times more species at time B than time A. Bob Ludovico Frate PhD student (University of Molise - Italy) Environmetrics Lab http://www.distat.unimol.it/STAT/environmetrica/organico/collaboratori/ludovico-frate-1 Department of Biosciences and Territory - DiBT Universit� del Molise. Contrada Fonte Lappone, 86090 - Pesche (IS) ITALIA. Cel: ++39 767557 Fax: ++39 (0874) 404123 E-mail ludovico.fr...@unimol.it ludovicofr...@hotmail.it [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] EstimateR and Chao 1 standard deviation
On 10/17/2014 05:36 PM, Gavin Simpson wrote: On 17 October 2014 06:44, Bob O'Hara boh...@senckenberg.de mailto:boh...@senckenberg.de wrote: On 17/10/14 13:54, José M. Blanco Moreno wrote: Dear users I have been checking the values returned by the function estimateR against the formulae in the appendix to nonparametric estimators of species richness in EstimateS (http://viceroy.eeb.uconn.edu/EstimateS/EstimateSPages/EstSUsersGuide/EstimateSUsersGuide.htm#AppendixB). http://viceroy.eeb.uconn.edu/EstimateS/EstimateSPages/EstSUsersGuide/EstimateSUsersGuide.htm#AppendixB%29... and they do not match each other. Where does this line of code in vegan:::estimateR.default come from? sd.Chao1 - sqrt(a[2] * ((G^4)/4 + G^3 + (G^2)/2)) It is from any other reference that I should have under control? Unless Jari added something, the functions come from my paper from about 10 years ago: O’Hara RB (2005) Species richness estimators: how many species can dance on the head of a pin? J Anim Ecol 74: 375–386. http://doi.wiley.com/10./j.1365-2656.2005.00940.x, which refers to Chao (1987), which is cited in the estimateR documentation. EstimateS now uses a small sample correction: http://viceroy.eeb.uconn.edu/EstimateS/EstimateSPages/EstSUsersGuide/EstimateSUsersGuide.htm#Chao1AndChao2 Does this account for the discrepancy? Bob Yep, vegan does not include the small sample correction. The Chao1 computation is here: https://github.com/vegandevs/vegan/blob/master/R/estimateR.default.R#L44 but we can certainly add it if this is useful? (And I suppose it is if EstimateS has included it by default now...) Might as well. it looks like it shouldn't make a huge difference unless sample sizes are small. Bob G -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 tel:%2B49%2069%207542%201863 Mobile: +49 1515 888 5440 tel:%2B49%201515%20888%205440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org http://www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org mailto:R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Gavin Simpson, PhD -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Binomial GLMM or GAMM with random intercept and temporal correlation
On 08/26/2014 07:34 PM, SamiC wrote: Dear ecology mailing list, I am trying to model some binomial data (0/1) as a function of sex (0/1) and DistanceToFeature (continuous km’s) with an interaction between the two. My data is nested and I therefore want to include a random intercept for InidividualID and within that I want to include an AR1 correlation structure as the data is serially/temporally auto-correlated. I understand any correlation structure should be nested within the random effect. So far I have tried running the model using glmmPQL as glmmPQL(Y ~ DistanceToFeature * Sex + (1|InidividualID), correlation=corAR1(form=~1|IndividualID/ContinuousBout), family=’binomial’, data=’birds’) (note – ContinuousBout is an ID for where there are time gaps in the data). However, although this runs, am I right in understanding that I should not use PQL estimation with binomial data as it gives biased results? Does anyone know of a way I can model this? Bernoulli responses are tricky things, so I'd be a bit worried that you might be trying to fit too complex a model. But anyway If you don't mind going Bayesian, you could try INLA (R-INLA.org). Alternatively, could you simply include the previous time step's value as a covariate? I understand that this is also the case if I wish to use GAMM (as later I will be modelling a non-linear explanatory as well)? Yes - a GAM is just a GLMM with a complicated covariate structure. Additionally I will also be running a similar set up but where the data are not equally spaced in time (and therefore an AR1 structure would not apply). Can anyone give a recommendation of a modelling framework for this also. That should come from the subject matter - I can think of a couple of approaches, but they make different assumptions about how your system is behaving. Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Question on unexpected regression results
linear regression is not the method to use? Thanks for any additional suggestions. Bruce [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] deviance as a goodness of fit in GLM
On 21/07/14 19:52, Samantha PameLa wrote: Good day everybody, I'm a marine biologist student, working on my bachelor thesis and I'm stucked with a statistical doubt in the process, I hope someone here could help me. My thesis aims to understand which biological and environmental factors influences the male aggressive rate of male California sea lions. For that purpose I'm using GLM’s where the response variable is the male aggressive rate. Right now I am testing the goodness of fit of the global models and for that I'm using the deviances as a goodness of fit test. I calculated pseudoR2 (Zuur, 2009) in order to know the percentage of explanation of each candidate model. However I’m not sure how to choose the “good models”; since I am not sure over which percent of explanation indicates a “model with good fit”. For my data I am working with three different scenarios, and it seems that 20%, is a good value to could indicate the best models, but I’m not sure how to choose the value. I thank you in advance for your time and the help you can give me. Best regards, Just to follow up, AIC isn't a measure of model fit either: it's a measure of model adequacy, and can be used to compare different models. For model fit, you can (usually) compare the residual deviance to the degrees of freedom (they should be approximately equal, and you can use a chi-squared test, if you feel the need to generate a p-value). This doesn't work if your response is binomial with a small N (and psuedo R2 doesn't work terribly well for this case either). A better approach to checking if the model fits is to check to see if the residuals have any gross patterns in them, by plotting them against the fitted values and also against covariates. If you'll excuse my self-promotion, PNAS recently gave me an excuse to show some residual plots to the general public: http://www.theguardian.com/science/grrlscientist/2014/jun/04/hurricane-gender-name-bias-sexism-statistics (my explanation for why this is not just gratuitous self-promotion is that I linked to the R code to generate the plots too, so you have something to work from). Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Question about GLM post hoc and chi square (Angel)
On 16/06/14 08:49, Åström, Jens wrote: Hi, Adding to what's already been said: Just in case you're not aware of this, I think you have a typo in your model. Are you looking for glm(Count~Sex+Time+Behaviour+Sex*Time+Sex*Behaviour+Time*Behaviour,family=poisson)? By the way, I think this could also be written glm(Count~Sex+Time+Behaviour+Sex:Time+Sex:Behaviour+Time:Behaviour,family=poisson) or simply glm(Count~+Sex*Time+Sex*Behaviour+Time*Behaviour,family=poisson). The * means both main effects and all possible interactions. FWIW, glm(Count~(Sex +Time+Behaviour)^2,family=poisson) will also work: it expands to the main effects and first order interactions. (oh, and the advice I've seen is don't use Wald tests: they are conservative because they assume the other parameters are fixed at their MLEs, so ignore any uncertainty in them. Use likelihood ratio tests instead) Bob Also, you should probably look into the issue of overdispersion. Overdispersion is very common in ecological count data and basically means that you have more variation in your data than the Poisson distribution assumes. This typically leads to anti-consevative p-values, i.e. too small p-values, and needs to be accounted for. Read more about it and potential solutions here: http://glmm.wikidot.com/faq Good luck, Jens -- Message: 1 Date: Sat, 14 Jun 2014 18:55:44 -0700 (PDT) From: Angel alexander.angelo...@hotmail.com To: r-sig-ecology@r-project.org Subject: Re: [R-sig-eco] Question about GLM post hoc and chi square Message-ID: blu174-w408ebbb0f91199f57e4d90f2...@phx.gbl Content-Type: text/plain You are able to obtain Chi squared values by using a wald chi square post hoc test. To do this you can use the aod package, function wald.test. This function is perfect for generalised linear models using poisson distribution. Also, as long as you have got interaction terms in your results from GLM, you could get the wald chi square (and an associated p-value) for these terms, hence giving you the table which you are after. If this response is not clear enough, I can post some example (I am not sure of the etiquette) AA. Date: Fri, 13 Jun 2014 18:07:27 -0700 From: ml-node+s471788n757894...@n2.nabble.com To: alexander.angelo...@hotmail.com Subject: Question about GLM post hoc and chi square Dear all, I am making an analysis using a GLM using three explanatory variables and a response variable. I need to obtain a table similar to this one, http://postimg.org/image/5sau79wlt/r nevertheless, I have not been able to do it. I am having a hard time specially getting the chi square values. I would like to know how to obatin them. I also would like to know what function could help me to make ad hoc comparisons for single variables and interactions. If any of you knows how to do both estimations, I would really appreciate it. All the best!!! This is my script a=read.table(ricis3.txt,header=T) attach(a) model7=glm(Count~Sex+Time+Behaviour+Sex*Time+Sex*Behaviour+Time+Behaviour*Sex,family=poisson) summary(model7) ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Question about GLM post hoc and chi square
On 06/14/2014 03:05 AM, Luis Fernando García wrote: Dear all, I am making an analysis using a GLM using three explanatory variables and a response variable. I need to obtain a table similar to this one, http://postimg.org/image/5sau79wlt/r nevertheless, I have not been able to do it. I am having a hard time specially getting the chi square values. I would like to know how to obatin them. Use anova(). The deviance follows a chi-squared distribution (usually - if you have overdispersion it gets a bit more complicated). I also would like to know what function could help me to make ad hoc comparisons for single variables and interactions. These comparisons are called contrasts. There is a contrasts() function in R, and also a contrast package (which, I'm guessing will be of more use). Googling R contrast might help too - there seems to be plenty of material, so hopefully one or two results will be exactly what you want. Contrasts can get esoteric, so if you can find some a page with code that gives you the comparisons you want, that should help a lot. Good luck! Bob If any of you knows how to do both estimations, I would really appreciate it. All the best!!! This is my script a=read.table(ricis3.txt,header=T) attach(a) model7=glm(Count~Sex+Time+Behaviour+Sex*Time+Sex*Behaviour+Time+Behaviour*Sex,family=poisson) summary(model7) ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] mantel within mantel? multidimensional variance
Do you need to run Mantel tests at all? I don't see that you have measured any distances, as such, so I think you can do at leat as well by using standard [genrealized] linear [mixed] models. I'm not sure I understand everything, but it sounds like you run an experiment with 2 treatments, and N species, and you measure T traits. So, for each trait you could have a model y_t ~ (Treat1 + Treat2)*Species (you could run this separately for each species, if that makes things easier) and size of the Treat1:Species effect tells you how the species responds to the treatment (if you want a measure over all traits, a PCA on the coefficients should suffice). It sounds like you are then asking how these reactions, i.e. the Treat1:Species effects are related to traits which are constant within a species, i.e. a Treat1:Trait effect. So, for each y_t you could run a model y_t ~ (Treat1 + Treat2)*(Trait+Species) (or even make Species a random effect) So, now I think you just want the Treat1:Trait effects. if the Trait's are discrete classes, then you can (again) run a PCA on the effects, to visualise the effects. If you want to analyse all of the y_t's together, it gets a bit more complicated, but in principle it's the same, except you have a MANOVA-like structure. I think you could use a Seemingly Unrelated Regression approach. Bob On 06/08/2014 05:56 PM, Mgr. Martin Weiser wrote: Dear friends, I need to quantify variance within variance, and what is worse, say if it is non-random. This is the setup: in the experiment, there were 2 (partly correlated) treatments, each of six levels (but theoretically of infinite levels, so I treated them as continuous). Different species responded to them, and we measured various traits (endotraits). We used R2 from the per-species redundancy analyses (RDA) as a reaction norm: higher R2, more the species responds to the treatment (so traits were used as species in the community ecology jargon). We have twice as much RDAs as there were species (because of 2 treatments). Next, we correlated these R2 (reaction norms) with some other traits per species (exotraits). As there were 2 correlated treatments (lets say irrigation and fertilisation), and we sed the same reaction norms for correlation with different exotraits, the correlations were obviously non independent. To overcome this, we run Mantel test (matrix1= species x reaction norms to treatment, matrix2=species x exotraits). Next, we are interested in endotraits x exotraits correlation. For this, we used endotrait scores from the per-species RDAs on the treatment axis (constraining variable). Correlations are non-independent again, so we wanted an overall test. And here comes the problem: results of Mantel tests (matrix1=species x endotrait scores, matrix2=species x exotraits) are suprisingly weak. Some single correlations of endotrait scores x exotraits seem to be pretty afar from random, but the overall mantel... I think this is because single endotrait can not show better correlation than the overall reaction norm, which is based on them, so something like (pval of this mantel)*(1-pval of the correlation of reaction norm with exotraits) may be desirable for the overall test, but I simply do not know. Any advice? Best, Martin Weiser -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Help with a function
On 06/06/14 15:39, Rodrigues wrote:
Dear R users,
Iâm trying to build a function to select random samples idfâs from a
database.
So, my data frame had 2 columns and 575 rows. Follow bellow an example of my
database
Idf1 casod
121
141
151
161
173
183
193
213
251
241
261
281
293
323
333
351
363
371
483
So my function is
blinding=function(sample){
sort=sample(idf1,10,replace=F)
return(sort2)
}
Does this work? It looks like you're passing a function to the function!
Also, you're not returning something that has been calculated. In
practice, you might as well just use sample(idf1,10,replace=F).
It is pretty simple and I would like to add one more step in my choice. I
would like to link my choice to casod stats. Thus if casod==3 sample would
be random idfs could not be an idf with casod=1. Does someone can help me?
If I understand you, you only want to sample within casod==3. So this
will work:
sample(idf1[casod==3],10,replace=F)
Or, if you need a wrapper function, this should work:
blinding=function(which.casod, df, N=10)
sample(df$idf1[df$casod==which.casod], N, replace=F)
e.g. SampledData=sapply(unique(Data$casod), blinding, df=Data)
will sample for each value of casod
Bob
--
Bob O'Hara
Biodiversity and Climate Research Centre
Senckenberganlage 25
D-60325 Frankfurt am Main,
Germany
Tel: +49 69 7542 1863
Mobile: +49 1515 888 5440
WWW: http://www.bik-f.de/root/index.php?page_id=219
Blog: http://blogs.nature.com/boboh
Journal of Negative Results - EEB: www.jnr-eeb.org
[[alternative HTML version deleted]]
___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] inference from GLMM fit using glmer
+ (1 | dist) stat df p.value LRT513.43 1 2.2e-16 *** PBtest 513.43 0.009901 ** --- mc.oc.yr.noyr Parametric bootstrap test; time: 39.83 sec; samples: 100 extremes: 0; large : oc ~ bfi + lyr + (1 | dist) small : oc ~ bfi + (1 | dist) stat df p.value LRT98.323 1 2.2e-16 *** PBtest 98.323 0.009901 ** --- We calculated bootstrap confidence intervals using confint. They seem to indicate that the estimated effects were unambiguous, but we were unsure whether our sample size was adequate to apply this method. confint (oc.sr.bfi.yr, method = c(boot)) Computing bootstrap confidence intervals ... 2.5 % 97.5 % sd_(Intercept)|dist 0.13871500 0.83969172 (Intercept) 6.54440740 7.56234905 bfi -0.53389761 -0.44979598 lyr 0.02835891 0.04325712 Simple diagnostic plots (residuals vs fitted values and a normal QQ plot of residuals) of the GLMM with effects of both bfi and lyr did not reveal violations of assumptions. We would be happy to report significant effects of food availability and a significant trend across years, but are concerned that we may not be interpreting these results correctly. Any help or suggestions would be greatly appreciated. Sincerely, Eric [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Bayesian analysis with MCMC
On 01/18/2014 04:39 AM, 明 wrote: Dear all I have a model with four paramters, I want to estimate the parameter uncertainty, so Bayesian analysis with MCMC method is applied. But every sigle mcmc chain seems give quite different parameter marginal distributions. In order to get the true parameter marginal distributions, I do like this: (1) I take 100 MCMC chain, and each chain has 1 iterations. (2) caculate the different parameter marginal distributions according to the frequence of paramter in step 1 sampling. The result seems reasonable. but is it right? Looking forward for your reply, Thanks in advance Han Ming [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology Have you checked for convergence? If the chains are giving different marginal distributions, this suggests they haven't converged, so they are not sampling from the correct posterior distribution (or, at least, some chains aren't). Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] This chains contains uninitialized variables and the trap report says undefined real result
On 28/11/13 15:09, kassahun Alemu wrote:
Dear all,
When I do analysis, I have troubled with frequent pop up message:
This chains contains uninitialized variables and the trap report says
undefined real result
I am looking for your courage and help to continue my further analysis.
My guess is that some of your mu's are getting too large, so exp(mu[])
becomes infinite. In your inits you're giving the precisions values of
dgamma(1,5,10)=0.1891664 (I think you meant to use rgamma()), so the
variance is about 5. With alpha and all of the betas being positive, I
can easily see exp(mu) exploding.
I'd recommend you use BUGS to do the debugging, because it's easier to
look at the problems (e.g. here I'd use save state to see what parameter
estimates you have: I'm guessing some mu's will be large).
In models like this I usually set the inits for my beta's to be close to
zero, so they don't explode the mu's unexpectedly. You can always relax
the inits later, once it's working.
Bob
model{
for (j in 1:2160) {
y[j] ~ dpois(mu[j])
log(mu[j]) - log(e[districts[j]]) + alpha + beta1*x1[j] + beta2*x2[j] +
beta3*x3[j] + beta4* x4[j] +
beta5*x5[j] + beta6*x6[j] + beta7*x7[j] + beta8*x8[j] +
phi[districts[j]]+ delta[time[j]] + omega[time[j]]
rho[j]-exp(alpha + beta1*x1[j] + beta2*x2[j] + beta3*x3[j] + beta4* x4[j] +
beta5*x5[j] + beta6*x6[j] + beta7*x7[j] + beta8*x8[j] +
phi[districts[j]] + delta[time[j]] + omega[time[j]])
# RR
}
for (i in 1:18) {
e[i] ~ dnorm(0,sigma.e)
phi[i] ~ dnorm(0,tau.phi)
}
delta[1]-0
omega[1] ~ dnorm(0,tau.omega1)
# Adjusted RR in month 1 over all districts/years
for (t in 2:120){
delta[t] ~ dnorm(0,tau.delta)
omega[t] ~ dnorm(omega[t-1],tau.omega)
}
#PRIORS
tau.e ~ dgamma(0.001,0.001)
alpha ~ dflat()
beta1 ~ dgamma(1,1)
beta2 ~ dgamma(1,1)
beta3 ~ dgamma(1,1)
beta4 ~ dgamma(1,1)
beta5 ~ dgamma(1,1)
beta6 ~ dgamma(1,1)
beta7 ~ dgamma(1,1)
beta8 ~ dgamma(1,1)
tau.phi ~ dgamma(1,1)
tau.delta ~ dgamma(1,1)
tau.omega1 ~ dgamma(1,1)
tau.omega ~ dgamma(0.001, 0.001)
sigma.e - 1/sqrt(tau.e)
sigma.phi - sqrt(1/tau.phi)
sigma.omega1 - sqrt(1/tau.omega1)
sigma.delta - sqrt(1/tau.delta)
sigma.omega - sqrt(1/tau.omega)
}
#R script to call R2WinBUGS and plot density
# Define the dataset
library(R2WinBUGS)
library(help=BRugs)
library(help=R2WinBUGS)
data1-read.csv(file.choose(Data1.csv), header=TRUE, sep=(,))
names(data1)
districslist-as.data.frame(unique(data1$districts))
data2-as.data.frame(cbind((1:18),districslist))
names(data2)-c(D_ID,districts)
data3-merge(data1, data2,by=districts)
data3-as.data.frame(data3)
names(data3)
data6-NULL
for (i in 1:18){
data4-data3[data3$D_ID==i,]
data5-cbind(data4,seq(1,length(data4[,1]),by=1))
data6-rbind(data6,data5)
}
districts-data6[,13]
time- data6[,14]
y- data1$y
x1- data1$x1
x2- data1$x2
x3- data1$x3
x4- data1$x4
x5- data1$x5
x6- data1$x6
x7- data1$x7
x8- data1$x8
data - list(districts=districts,time=time,
y=y,x1=x1,x2=x2,x3=x3,x4=x4,x5=x5,x6=x6,x7=x7,x8=x8)
# Define inits
inits - function(){
list(tau.e=dgamma(1,5, 10), alpha = dnorm(1,1, 1), beta1 = dgamma(1,5,10),
beta2 = dgamma(1,5,10),
beta3 = dgamma(1,5,10), beta4 = dgamma(1,5,10), beta5 = dgamma(1,5,10),
beta6 = dgamma(1,5,10), beta7 = dgamma(1,5,10), beta8 = dgamma(1,5,10),
tau.phi=dgamma(1,5,10), tau.delta=dgamma(1,5,10),
tau.omega=dgamma(1,5,10), tau.omega1=dgamma(1,5,10))
}
# run the MCMC analysis
model_G- bugs(data, inits, model.file = C:/model_G.txt,
parameters = c(tau.e, alpha, beta1, beta2, beta3, beta4,
beta5, beta6, beta7, beta8, tau.phi, tau.delta, tau.omega1,
tau.omega),
n.chains = 3, n.iter = 2000, n.burnin=500,bugs.directory = C:/Program
Files/WinBUGS14/,
working.directory=NULL, WINE=NULL, debug=TRUE)
BUT, this does not work, the log file says this chain contains
uninitialized variables and trap says undefined real resul.
I cannot know where the problems are.
Thanks for your kindest help,
Sincerely,
Kassahun
[[alternative HTML version deleted]]
___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
--
Bob O'Hara
Biodiversity and Climate Research Centre
Senckenberganlage 25
D-60325 Frankfurt am Main,
Germany
Tel: +49 69 7542 1863
Mobile: +49 1515 888 5440
WWW: http://www.bik-f.de/root/index.php?page_id=219
Blog: http://blogs.nature.com/boboh
Journal of Negative Results - EEB: www.jnr-eeb.org
[[alternative HTML version deleted]]
___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Two-sided significant test for null models
On 07/23/2013 12:06 PM, Cayetano Gutiérrez Cánovas wrote: Dear all, I'd like to test the two-side significance of an observed value being lower or higher in comparison to a null distribution, but I don't know how to build this in R code. Can anybody help me with some code or any R function? Let's say 'null.dist' the null distribution, 'obs' the observed value and 'alpha' the probability of Type I error. set.seed(10) rnorm(999)-null.dist obs-1.5 alpha-0.05 pval=mean(obsnull.dist) min(pval, 1-pval)0.05 Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] interpretation of interaction between explanatory variables
It's difficult to give definitive advice about what it means, but the simplest approach is to look at what it means: plot the model and see what it looks like. Basically, the interaction says that as the variables both increase, the relationship gets less negative. The details depend on the ranges of the variables, hence plotting the results is helpful. I would write some code, but that would probably mainly show that I'm out of date with R. But you could use expand.grid() to create new data to predict, use predict() to predict it, and then plot it as several lines on one plot (e.g. plot against C.abundances, with one line for each value of C.diversity). Bob On 05/08/2013 02:56 PM, Iris Kröger wrote: Dear list members, I want to analyse the impact of a competitor community (i.e. community abundances on the one hand and community species diversity on the other hand) on mosquito larval populations of species A and B. Each variable on its own has a negative impact on mosquitoes - but when both variables are interacting, there is a positive impact... How can I interpret that? For mosquito A only C.diversity has a significant impact - but the interaction between C.abundances and C.diversity is significant? What does that mean? I used the model: lm (mosquito ~ C.abundances * C.diversity) output Mosquito A: Estimate Std. Error t value Pr(|t|) (Intercept) -0.2120 0.1159 -1.829 0.074855 . C.abundances -0.1277 0.1616 -0.790 0.434067 C.diversity -0.5787 0.1385 -4.178 0.000155 *** C.abundances:C.diversity 0.4096 0.1712 2.393 0.021476 * Output Mosquito B: Estimate Std. Error t value Pr(|t|) (Intercept) -0.2900 0.1220 -2.377 0.02233 * C.abundances -0.3856 0.1701 -2.266 0.02891 * C.diversity -0.3470 0.1458 -2.381 0.02213 * C.abundances:C.diversity 0.5367 0.1801 2.980 0.00489 ** Thanks a lot for your help! Iris [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 798 40226 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] GLM
On 03/07/2013 04:24 PM, Mahnaz Rabbaniha wrote: Dear all I want to find regression between fish larva abundance and some abiotic factor ,i used this code: glm(formula = mychto ~ po4 + No3 + Si + Tn) result: Deviance Residuals: Min 1Q Median 3Q Max -26.586 -18.262 -12.296 -2.949 226.229 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 67.421173.9781 0.9110.371 po4 -0.2887 1.6037 -0.1800.859 No3 0.9151 4.5261 0.2020.841 Si -0.1145 0.4850 -0.2360.815 Tn -1.1568 4.4818 -0.2580.798 (Dispersion parameter for gaussian family taken to be 2444.917) Null deviance: 63156 on 29 degrees of freedom Residual deviance: 61123 on 25 degrees of freedom AIC: 325.72 my question is about the acceptable this AIC, or this result with goodness of fit? thanks AIC tells you nothing about goodness of fit, it is used to compare different models. A better way to assess goodness of fit is to look at the data and the model. It's clear that you have very skewed data (look at the distances between the quartiles of the residuals), so a Gaussian model, which assumes the residuals are symmetric, is inappropriate. It might make more sense to use a Poisson distribution, although you will probably need to account for over-dispersion (there are a few ways of doing this). Also, once you've fitted a glm, you can use plot() to get some useful model-checking plots (e.g. residuals). There is more advice in the links here: http://r.789695.n4.nabble.com/Checking-the-assumptions-for-a-proper-GLM-model-tp1559502p1560503.html. Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 / +49 69 798 40226 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] patterns in weather data that could relate to pathogen prevalence
On 11/30/2012 01:09 PM, Anto Raja wrote: Hi all I am searching for a tool that would help me to identify weather patterns that influence the prevalence of a pathogen, 'Pn'. Say, I have annual prevalence data (collected in April) and I know that the prevalence of 'Pn' is affected by the weather conditions since November. I also have daily data for different weeather variables. You have a lot of weather data, but I assume you don't have so much prevalence data. So it's going to be difficult, whatever you do... The objective is to identify the relationship between the weather from Nov-Mar and the prevalence of Pn. We know that weather has an influence on Pn. The question is to find out weather from which period is relevant or what kind of weather is relevant. It could be that the first two winter months (Nov-Dec) is the decisive factor or that a certain weather situation (like 20 consecutive days of below zero conditions) occuring at any time is important or a combination of both. I have tried correlations between prevalence and monthly means for Mar, Feb-Mar, Jan-Mar and so on and nothing definite turned up. I could also do it on a weekly basis manually. But I wonder if there is a tool that uses a moving window of different sizes (say, from a min size of 1 week to a max of 4 months) and checks correlations for each of these periods. I am thinking of ARMA, but my present intention is not to forecast but only to study. Can it still be used? Or ARMA in combination with multivariate analysis to study the relative importance of each weather variable. I don't see why an ARMA model would help you, as that assumes a covariance between times (i.e. autocorrelation) in the response (i.e. prevalence). There are methods for assuming that the response has an autocorrelation, but I don't think that's your big problem. My reaction (without seeing the data, of course) is that you might be asking too much of your data to get anything meaningful out of it. Any suugestions are welcome. I have used R for basic stats analysis but never worked with time-series data or the advanced tools of data mining. So, it could also be possible I am not thinking along the right lines. Feel free to correct if I am looking in the wrong place. It sounds like you're trying to mine your data for any pattern. To be honest, if you do that, I wouldn't trust the results unless you can validate them independently: you'll find some relationship if you try enough models, but will it make biological sense? This is particularly problematic when you have correlated variables, which you will do (especially when you start sliding windows around) I'd suggest you start by using what's known of the pathogen or its host, or of similar host-pathogen systems, to develop a smaller number of hypotheses about what sort of effects are likely. Plant ecologists often use GDD5 (Growing Degree Days above 5°C), which might be a useful way of reducing the temperature data to something smaller. Of course, another temperature than 5°C might work better for you. Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 7542 1863 / +49 69 798 40226 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Spatial AND temporal correlation in linear mixed-effects models
On 07/26/2012 10:45 AM, Tscheulin Thomas wrote: Hi, I have a question regarding the nesting structure in linear mixed models of data, which is spatially and at the same time temporally correlated. So far I have tried to get around the problem by averaging out the temporal component but I would really like to keep everything in the model. Here is the experimental set-up: We have measured insect abundance in two different islands, each having 5 sites, each site having 4 plots, in which each we have 5 collection points (traps). The sampling was repeated 5 times. Two times in 2011 and 3 times in 2012. After averaging the temporal correlation out I composed the following model for the abundance: model-lme(abundance~continous_explanatory_variable,random=~1|island/site/plot/trap,method=REML) This works fine but I have problems making a model that allows the temporal component to stay in. How can I do this using the function lme? I know with lmer I could do this: model2-lmer(abundance~continous_variable+(1|island/site/distance/triplet)+(1|year/round)) but I really want to use the function lme. How can I insert multiple levels of grouping in lme? Any help is very much appreciated! Best, A bit of googling brought up this, from 10 years ago: http://tolstoy.newcastle.edu.au/R/help/02b/2068.html I don't think lme has changed that much. Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 798 40216 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Using residuals as dependent variables
On 06/21/2012 07:06 PM, Chris Mcowen wrote: Dear List, I am wondering if the methodological approach I am taking is correct and would be very grateful if people could comment and make suggestions, much appreciated. The simple answer is no it isn't: http://gking.harvard.edu/files/mist.pdf I have developed the best model ( AIC model selection) using oceanographic data ( i.e. SST, chlorophyll, NPP...x6) to predict reported fisheries catch for 52 countries. I then extract the residuals from the model and anything positive has a higher catch than would be predicted given the level of productivity in the country, with the opposite being true. What I want to do is: 1. Regress a suite of ecological and socioeconomic variables against the residuals from the oceanographic model to determine which factors cause some countries to be above and some below. I.E as trophic level increase the residuals become increasingly negative. You could simply put everything into one big model. 2. Ideally ( and I am unsure how or if it is possible) work out for each country which variables/s cause the poor fit of that country to the oceanographic model. Try using partial residuals (http://en.wikipedia.org/wiki/Partial_residual_plot). Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 798 40226 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Multiple comparisons among predictors generated from same data
On 05/25/2012 10:18 AM, Gavin Simpson wrote: On Thu, 2012-05-24 at 15:00 -0700, J Straka wrote: 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. What is the model for? Understanding so you want to interpret the coefficients directly as something meaningful or for prediction? If the latter I would say it doesn't really matter; choose the model that gives the best out-of-sample predictions (lowest error etc), or average predictions over a set of best/good models. Simply choosing the best model via some sort of selection procedure may result in a model with high variance (change the data a bit and different variables would be selected). If so, consider a regression method that applies shrinkage to the coefficients such as the lasso or the elastic net; this will lead to a small bit of bias in the estimates of the coefficients but should reduce the variance of the final model because you are considering the selection of variables as part of the model itself. If you want to interpret the model coefficients as something real then you have to be very careful doing any form of selection; the stepwise procedures and best subsets all can potentially lead to strong bias in the model coefficients. Be removing a variable from the model in effect you are saying that the sample estimate of the effect of that variable on the response is 0, not some small (statistically insignificant) value. This is a very tricky thing to get right and I'm not sure I know the right answer (or even if there is one!?). An additional complication here is that the variables are going to be correlated, so a model with all or most in it could be unstable. If a single temperature variable is enough, then I'd suggest either trying your best to pick one, or use what everyone else uses (GDD5?), so the study can be comparable. Once you have a model, it might be worth checking to see if the other variables tell a different story. If it's the same story but with different p-values, you might as well stick to the original analysis. Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 798 40226 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] PCA as a predictive model
On 05/23/2012 10:55 AM, Marc Taylor wrote: Hi Jari - one more question if you don't mind. Since the weights of the PCs are related to the the amount of variance that they explain in the original data - is it problematic to predict the PC scores with a second data set that has a different amount of variance (e.g. due to differing number of samples)? In both the 1st and 2nd data sets I have been using scaled values for the variables (mean=0 and sd=1 for each sample). Cheers, Marc I'll pretend to be Jari for a moment. :-) PCA just scales and rotates the data in cunning ways, so with the new data you need to scale and rotate it in the same way. If you scale the values first then you've already changed the scaling. What you need to do is either do PCA on the raw data or scale the new data using the mean and varianes of the old data. library(MASS) NVar=5; NObs=50 Sigma=matrix(c( 10,0.2, 0, 0,0.4, 0.2, 5,0.1, 0,0.6, 0,0.1,1.0, 0.2,0, 0, 0,0.2, 5, 0, 0.4,0.6, 0, 0,1), nrow=5) # simulate data Data=mvrnorm(NObs, rnorm(NVar), Sigma=Sigma) # Do PCA on scaled data Data.Sc=scale(Data) PC=princomp(Data.Sc) # Simulate new data NewData=mvrnorm(10, rnorm(NVar), Sigma=Sigma) # Do PCA on new data. First do it wrong... PC.wrong=predict(PC, newdata=scale(NewData)) # Now scale correctly NewData.Sc=scale(NewData, center=attr(Data.Sc, scaled:center), scale=attr(Data.Sc, scaled:scale) PC.right=predict(PC, newdata=NewData.Sc) HTH Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 798 40226 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] gam
On 05/17/2012 06:17 PM, Mahnaz Rabbaniha wrote:
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] depthtemperature salinity Sillaginidae
pairs(sc,panel=function(x,y){points(x,y);lines(lowess(x,y))})
model-gam(Sillaginidae~s(temperature)+s(salinity)+s(depth))
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?
It means that you don't have a lot of different values for one of your
covariates. You can check which with
length(unique(depth))
length(unique(temperature))
length(unique(salinity))
Unless you have a lot of data, that the sensible thing to is probably to
fit a straight line rather than a smooth curve. The other thing to do is
to define the maximum flexibility of the curve, e.g.
model-gam(Sillaginidae~s(temperature, k=6)+s(salinity)+s(depth))
I hope this helps.
Bob
--
Bob O'Hara
Biodiversity and Climate Research Centre
Senckenberganlage 25
D-60325 Frankfurt am Main,
Germany
Tel: +49 69 798 40226
Mobile: +49 1515 888 5440
WWW: http://www.bik-f.de/root/index.php?page_id=219
Blog: http://blogs.nature.com/boboh
Journal of Negative Results - EEB: www.jnr-eeb.org
___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] data normality
On 05/06/2012 10:54 AM, Yong Shen wrote: 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? No. The only area where a Normal assumption comes in is that the residuals are normally distributed. So you can happily fit the model without worrying about normality until after you've got the model. 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? Leave them as they are. Advice that makes life simpler - always the best sort. Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 798 40226 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] ICC confidence intervals and power analysis for random effects in lmer?
On 04/13/2012 04:17 AM, Bradley Carlson wrote: Thanks for the tips everyone - I'll look into MCMC sampling for the CI. As far as the power analysis goes, I'm somewhat familiar with the criticisms regarding power analysis. I think this reviewer was curious about it because there was a small sample size (small number of levels of the random effect) and the ICC point estimate was not so low as to be biologically insignificant. It would be nice to state in the paper how much larger of a sample size would have enabled us to detect an effect given the observed variation, as though this were a pilot study for planning a bigger experiment. I'll certainly bring up the suggested citations, but I would still be interested to know if there is a method available for performing a power analysis for an LRT of a random effect. Thanks again, I'd use simulation: you can write a function to simulate data and fit the model to it, so you can run that for different sample sizes. It's not too fiddly to write, and you can leave the computer to run it over the weekend, so you can run reasonable sample sizes. Bob Brad On Thu, Apr 12, 2012 at 12:18 PM, Brian Inouyebdino...@bio.fsu.edu wrote: In addition to Bob O'Hara's suggestion, here is another citation you can give the reviewer/editor, as to why retrospective power analyses are a waste of time. Hoenig, J. M. and D. M. Heisey (2001). The abuse of power: the pervasive fallacy of power calculations for data analysis. American Statistician 55(1): 19 - 24. Last year the ESA updated its author guidelines for reporting statistics, and removed a suggestion to report power analyses that had been inserted in the 1980s. -Brian Inouye Florida State University Chair, statistical ecology section of the ESA On 4/12/2012 6:00 AM, r-sig-ecology-request@r-**project.orgr-sig-ecology-requ...@r-project.orgwrote: 2) A reviewer requested a power analysis of the ability to detect a significant random effect. Any tips on how to approach that? Report the random effect and confidence intervals. Retrospective power analyses are pretty pointless (e.g. see http://beheco.oxfordjournals.** org/content/14/3/446.fullhttp://beheco.oxfordjournals.org/content/14/3/446.full), unless you're planning to repeat the experiment. Bob __**_ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/**listinfo/r-sig-ecologyhttps://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 798 40226 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Which function would fit these data (Version 2)
On 04/13/2012 06:16 AM, Beutel, Terry S wrote: Hi list members Sending this again as symbols in the previous version did not translate. Apologies if this is not too specifically R related, but I am looking fit a model to some simulated data. X is distance to sample point, y is binary outcome (present/absent). I was hoping someone can suggest a (presumably) non linear function that might meet the following criteria A. 0= y=1 (actual responses are binary) B. xmax = infinity C. xmin= 0 D. Where x = 0 then y= 0 E. Where x equals Infinity then y equals 1 F. As x approaches Infinity from xmin, the slope of the function declines monotonically toward 0. my initial reaction was to think logistic regression, but my experience with questions like this is that we'll save a lot of time if you explain the context. BTW, I'm not sure why you want the function to decrease to zero when x approaches infinity when at the same time y should equal 1 at infinity. It's not a big problem, though. Just use 1-y instead of y. Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 798 40226 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] ICC confidence intervals and power analysis for random effects in lmer?
On 04/12/2012 03:50 AM, Bradley Carlson wrote: I'm performing an analysis of behavioral variation among individual tadpoles, using individual ID as a random effect and time as a continuous fixed covariate in the lmer function in lmer4 package. I'm really interested in making inferences about the random effect (i.e. the extent of variation among individuals). I'd like to do two things that I can't seem to find straightforward answers about and I'm hoping someone can help or point me to a good resource. 1) The intraclass correlation coefficient is of particular interest to me, as it describes the proportion of variation that occurs among individuals. Ideally I'd like to report a confidence interval of the ICC but I can't find any way to calculate one, other than a function in the psychometric package that appears to only work when there are no covariates in the model (random effect only). MCMC has already been mentioned and lme4 still has its mcmcsamp() function. Failing that, you could try a parametric bootstrap, which requires a little bit of coding but simulate() makes it much easier. 2) A reviewer requested a power analysis of the ability to detect a significant random effect. Any tips on how to approach that? Report the random effect and confidence intervals. Retrospective power analyses are pretty pointless (e.g. see http://beheco.oxfordjournals.org/content/14/3/446.full), unless you're planning to repeat the experiment. Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 798 40226 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Output for interactions in models that do not include all main effects
On 04/03/2012 11:31 PM, Kristen Gorman wrote: Dear all, I have R code to run AIC including multi-model inference. I am running into a problem in calling the output from models where both parameters in an interaction are not included as main effects. Why would you want to do that? Why would you (for example) expect the average of the Rlipid slope to be zero if the slope varies with the value of RFGinit? Does this make sense? (this is the sort of thing that makes statisticians splutter into their tea when they see someone do it: it rarely makes sense. Well, unless you have nested effects - which you don't have here- where the interaction is the nested effect) if you respect marginality, there won't be a problem because the main effect is always included. If you really want to include interactions without main effects, you can either write the formula by hand, using paste(): something=Rlipid form = paste(Slipid ~ , something, + RFGinit:, something, sep=) lm(form, data = DataSet) and then work out how to get the order. Or you could try using update(): mod1 = lm(formula = Slipid ~ RFGinit*Rlipid, data = DataSet) mod2=update(mod1, . ~ . -RFGinit) HTH Bob In R, the interaction will be called depending on the parameter that was used as the only main effect in the model. So, I end up generating 2 different interactions (e.g., Rlipid:RFGinit vs RFGinit:Rlipid) that are actually the same. This becomes a problem in the remaining R code that requires weighted and summed values for the parameter and SE estimates. Thus, I would like to call the interaction consistently across models. See the following code: -- lm(formula = Slipid ~ Rlipid + RFGinit:Rlipid, data = DataSet) Residuals: Min 1Q Median 3Q Max -74.075 -19.047 7.233 20.445 45.391 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept)120.338475.30405 22.6882e-16 *** Rlipid 0.304930.23615 1.2910.202 Rlipid:RFGinit -0.020990.01773 -1.1840.241 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 30.88 on 60 degrees of freedom Multiple R-squared: 0.02721,Adjusted R-squared: -0.005221 F-statistic: 0.839 on 2 and 60 DF, p-value: 0.4372 lm(formula = Slipid ~ RFGinit + Rlipid:RFGinit, data = DataSet) Residuals: Min 1Q Median 3QMax -76.35 -21.63 7.09 22.46 45.71 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept)131.028546 8.717104 15.0312e-16 *** RFGinit -0.933483 0.742083 -1.2580.213 RFGinit:Rlipid 0.003926 0.009283 0.4230.674 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 30.9 on 60 degrees of freedom Multiple R-squared: 0.02586,Adjusted R-squared: -0.00661 F-statistic: 0.7964 on 2 and 60 DF, p-value: 0.4556 -- Is there a way to tell R to call the interaction based on alphabetical order of the 2 interaction terms and not based on the term that was used as a main effect? Thanks very much for any insight. Kristen Gorman ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 798 40226 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Problem with R2WinBUGS packages
On 01/22/2012 05:50 PM, sergio vignali wrote: Hi, I have a problem with the R2WinBUGS package, when I run the bugs() function with the argument debug=TRUE, WinBUGS executes the calculations but R is blocked, instead with the option debug=FALSE it is all ok! Is there anybody who can help me? Sergio Vignali Uwe Ligges (lig...@statistik.tu-dortmund.de) is the maintainer. But isn't this expected behaviour? The help says debug if FALSE (default), WinBUGS is closed automatically when the script has fin- ished running, otherwise WinBUGS remains open for further investigation IIRC, you just need to close WinBUGS, and it'll return to R. You only need to use debug=T if BUGS is trapping or throwing errors (and have fun understanding the traps). Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 798 40216 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] Mixed model with zero truncated Poisson distributed data
On 01/11/2012 04:58 PM, Helen Ward wrote: Hello, I am a PhD student currently trying to derive a mixed effects model. I would like to describe the relationship between age and male reproductive success in a population of greater horseshoe bats. My data consists of three columns: MaleID, Age, NumberofPups (at that age). Many of the males appear multiple times in the data set, so I believe I need to derive a mixed model with MaleID as a random variable. The data is Poisson distributed, but zero-truncated. So far I have only succeeded in making a mixed model with a poisson distribution (using glmmPQL in the MASS package), and a zero truncated poisson model (using vglm in the VGAM package), but not a mixed model capable of handling zero truncated Poisson data. It has been suggested that I could just minus 1 from each value in the NumberofPups column to make a more usual Poisson distribution, so I can ignore the zero truncated bit. I have tried this and it changes the results of the model, but is this an acceptable transformation? If not, can anyone advise me on a mixed model that can handle zero truncated Poisson data please? I also intend to post this on the R-sig-mixed-models https://stat.ethz.ch/mailman/options/r-sig-mixed-models/h.l.ward%40qmul.ac.uk list, so apologies if you've seen it twice! If an individual would have had pups, would you have seen it? I'm wondering if you could simply add the zeroes back into the data. Bob -- Bob O'Hara Biodiversity and Climate Research Centre Senckenberganlage 25 D-60325 Frankfurt am Main, Germany Tel: +49 69 798 40216 Mobile: +49 1515 888 5440 WWW: http://www.bik-f.de/root/index.php?page_id=219 Blog: http://blogs.nature.com/boboh Journal of Negative Results - EEB: www.jnr-eeb.org ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
