tag:blogger.com,1999:blog-224463442020-02-29T02:29:14.403-08:00Statistics Forum Functional Ecology Research SchoolA blog for discussions on statistical issues, mainly related to glm and mixed models. Most members participated in courses of Dutch PhD research schools.Tom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.comBlogger32125tag:blogger.com,1999:blog-22446344.post-64756978383198173692010-02-13T04:59:00.000-08:002010-02-13T04:59:21.838-08:00draft of new mixed models bookI just received this email through the r-sig-mixed-models group, from Douglas Bates:<br /><br /><pre wrap="">I have finally gotten the chapter drafts of the book, "lme4:<br />Mixed-effects Modeling with R", that I am writing to the point where I<br />think they can be helpful to others. I believe I have properly<br />committed them to the R-forge site but they will not appear for an<br />hour or so. They should be available under<br /><a class="moz-txt-link-freetext" href="http://lme4.r-forge.r-project.org/book/">http://lme4.R-forge.R-project.org/book/</a><br /><br />Chapters 1, 2, 4 and 5 are currently available as are the .R files for<br />those chapters. Chapters 4 and 5 will have material added to them. I<br />don't expect to change chapters 1 and 2 until I get reviewers reports.<br /><br />When pointing out errors, inconsistencies, bad writing, obtuse<br />explanations, etc., please be gentle.</pre><pre wrap=""> </pre><pre wrap=""> </pre>Tom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-31850294617082148152009-11-24T14:49:00.000-08:002009-11-24T14:49:06.127-08:00lme - lmer mixed modelsDear Kim,<br /><br />first of all, please use lmer() from library "lme4", instead of lme. I think that is currently better maintained and developed than the "nlme" library.<br />The model you suggest is a random regression per individual (obs_nr)?<br />The syntax would become<br /><br />library(lme4)<br />cocon_mixed<-lmer(MiWi~host_spec*ins_spec+(date|obs_nr),data=cocon_nash)<br /><br />It should work provided there are several observations of MiWi for each level of obs_nr, which seems to be the case given your explanation of the data.<br /><br />Cheers, TomTom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-17904574220677431672009-11-24T05:30:00.000-08:002009-11-24T05:37:14.620-08:00Mixed effects modelsHi,<br /><br />I am currently trying the exercises of Crawley on mixed effects model (ex. 7). I remade the example on page 33/34 and am currently trying tho do this analysis with my own data. I have a response variable (size=MiWi), measured at different dates (date) on the same individual (obs_nr = random variable). Further more I have insects species (ins_spec) and host species (host_spec) as fixed variables. When I run it I get the error message:<br /><br />cocon_mixed<-lme(fixed=MiWi~host_spec*ins_spec,data=cocon_nash,random=~date|obs_nr) Error in lme.formula(fixed = MiWi ~ host_spec * ins_spec, data = cocon_nash, : <br /> nlminb problem, convergence error code = 1<br /> message = iteration limit reached without convergence (9)<br /><br />What is going wrong?<br /><br />KimUnknownnoreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-68012684339131255112009-11-12T15:00:00.000-08:002009-11-12T15:03:27.149-08:00Course GLM GroningenHi,<br />I just returned from the glm course in Groningen. I will invite all participants to this blog in the coming days. Meanwhile, I already add a link to a page with several simple scripts:<br /><br /><a href="http://tomvandooren.eu/tuto.html">http://tomvandooren.eu/tuto.html</a><br /><br />As you will see, these files haven't been updated at all.<br />If things in them have become obsolete, please let me know.<br /><br />Cheers, TomTom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-55750102959199532382007-04-07T12:01:00.000-07:002007-04-07T12:09:02.381-07:00random effectsDear Tom,<br />my apologies for the long delay between you request and this answer.<br />First of all, why not include precipitation ("Precip") as a main effect too?<br /><br />Using libray "nlme" and lme():<br />Random intercepts per station:<br /> random= ~1|station<br />Don't forget the vertical bar!<br />Including random slope:<br /> random= ~Precip|station<br />or<br /> random= ~1+Precip|station<br /><br />Make sure that all stations have different names, otherwise you run the risk that stations with identical names at different distances will be considered as a single station.<br /><br />You could also consider to use the updated library now called "lme4".<br />Then the random effects are simply included as extra terms in the model statement.<br /><br />Cheers, Tom<br /><span style="text-decoration: underline;"><br /></span>Tom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-3718708441764403332007-02-21T04:32:00.000-08:002007-02-21T05:20:44.301-08:00random components of lme-modelHi everyone,<br /><br />I'm currently analysing a large dataset with measurements of nutrients in the Dutch marine systems in R. I have a variable "OffshoreDist" representing the distance from the sampling location to the coast. I also calculated the total amount of precipitation over The Netherlands and called this variable "Precip".What I want to do is link mean annual dissolved inorganic nitrogen (DIN) to the amount of precipitation. I expect that this relationship is dependent on the offshore distance. Distance classes can contain 1 or more sampling stations.<br /><br />That's why I propose these two models with different random components:<br />model1 <- lme(log(DIN) ~ Precip:factor(OffshoreDist) + factor(OffshoreDist),<br /> random= ~1station)<br />model2 <- lme(log(DIN) ~ Precip:factor(OffshoreDist) + factor(OffshoreDist),<br /> random= ~Precipstation)<br /><br />Are these the right models for random intercepts per station (model1) and random slope+intercept per station (model2)? Or should this be something like:<br /><br />random = ~ factor(OffshoreDist) -1 station (model 1)<br /> and<br />random = ~ Precip:factor(OffshoreDist)station (model 2)<br /><br />I've read Pinheiro & Bates, 2000 which is a good reference, but explanation of the formulae only by means of examples does not seem to work for me.<br /><br />Thanks in advance.<br /><br />Cheers,<br /><br />Tom Van EngelandTomVEhttp://www.blogger.com/profile/09278631261618627114noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-17688027452364653602007-02-07T01:24:00.000-08:002007-02-07T01:27:35.933-08:00moved to new bloggerHi,<br />I've moved the blog into the new google template. Everyone logging in will be invited to move her/his account too. A nice aspect of the change is that the blog can now be searched easily!<br />Cheers, TomTom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1169557267803382372007-01-23T04:55:00.001-08:002007-01-23T05:01:07.813-08:00F lmerHi Marie,<br />do you really need to use the quasibinomial distribution? Is there evidence of overdispersion in the maximal model you fit?<br /><br />If not, then you can use the binomial distribution.<br />That allows you to use likelihood-ratio tests for the fixed and random effects.<br />using <br />anova(model,modelsimplified,test="Chisq")<br /><br />If there is clear overdispersion, then I guess the best thing to do would be resampling according to the model, if that works for the quasibinomial "distribution". You need to use mcmcsamp() for that.<br />Otherwise, it will become hand-work to get some of those hotly debated F-ratios out....<br />Cheers, TomTom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1169547716286570672007-01-23T02:13:00.000-08:002007-01-23T02:22:57.420-08:00lmer P-valuehi everybody<br />Do you know the simplest and correct way to produce P-values of lmer and of the ANOVA (F test) on the lmer model in R?<br />My lmer model is of this form:<br />model<-lmer(pref1~brood.size*logtarsus*cohort+(1birth.nest/foster.brood),family=quasibinomial,method='Laplace')<br />summary(model)<br />anova(model,test="F")<br />I read an extensive eplanation by Douglas Bates on R help archive named "[R] lmer, p-values and all that". However, nothing is mentionned regarding Laplace. So i do not know what to do.<br />Could someone give me advice?<br />many thanks<br />MarieMarie Holveckhttp://www.blogger.com/profile/16185407690012580518noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1168950443701485892007-01-16T04:25:00.000-08:002007-01-16T04:27:41.623-08:00Longitudinal DataFinally, the participants in the course <a href="http://www.rug.nl/biologie/onderzoek/onderzoekScholen/functionalEcology/phdCourses/multivariateStatistics">Longitudinal and Incomplete Data</a> have been invited. Transcripts for R of the SAS code presented during the course will be added to the weblog soon.<br />Cheers, TomTom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1168946697558998072007-01-16T03:21:00.000-08:002007-01-16T03:24:57.570-08:00Loglinear modelsHi Nicola,<br />the best manual is <a href="http://www.statslab.cam.ac.uk/~pat/Splusdiscrete2.pdf">http://www.statslab.cam.ac.uk/~pat/Splusdiscrete2.pdf</a>, based on the book by Agresti.<br />You can basically use glm with a poisson link for log-linear models, I guess it is also discussed in the <a href="http://www.stats.ox.ac.uk/pub/MASS4/">Venables and Ripley book</a>.<br />Good luck, TomTom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1168945538186814972007-01-16T02:58:00.000-08:002007-01-16T03:05:38.206-08:00Loglin models and contingency tablesHi Tom,<br /><br />I'm trying to use loglinear models to investigate contingency tables with frequencies, but I can't figure out how to get the standardized residuals, odds ratios etc. Also the help function doesn't mention them, and the R-manuals I have don't even mention loglin models at all.<br />Do I have to calculate them by hand or is there perhaps a good manual on loglin models?<br />Thanks in advance!<br />Nicola<br /><br />ps I have a 3-way table with one binary dependent factor and 2 explanatory factors with 3 levels each.nicolahttp://www.blogger.com/profile/02845233819066972346noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1167996658796850212007-01-05T03:20:00.000-08:002007-01-05T04:28:45.006-08:00lmer, random effect with four levelsDear Amir,<br />the "place" random effect in your dataset has only four levels, and they are nested within plants. The error which is returned states that the variance covariance matrix is singular, meaning that variation in intercept, slope and their covariance cannot be fitted. If you fit the random effect as (1|place), you will read that the random effect variance is effectively zero. <br />It does not make sense to include place both as fixed and as random, as in you second to fourth models. Since the number of levels of the random effects is very low, you could even just model them as fixed, using glm. The idea of random effects is that they are a representative sample from a large population, which is hard to defend with just two levels sampled per group!<br />More importantly, the distribution of your data does not look poisson at al. There are many zero observations. I would convert the counts to presence/absence data and analyse those as bernoulli probabilities. That works to some extent:<br /><br />b<-ifelse(Predators.Allspecies.>0,1,0)<br />model<-lmer(b~julian.date+offset(log.time)+(1|place),family=binomial)<br /><br />I've added you models to your post, so other visitors to this blog can also read and respond to your question,<br />all the best, TomTom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1167927384649984342007-01-04T08:08:00.000-08:002007-01-05T03:32:54.786-08:00Mixed effects model<pre wrap="">Beste Tom,<br />Tijdens de laatste GLM cursus hebben we samen een model gemaakt voor de<br />analyse van een dataset. Ik probeer het model nu te gebruiken, maar krijg<br />steeds fout meldingen in R. Helaas heb ik de R file van de analyse die we<br />samen gedaan hebben op een verkeerde wijze opgeslagen, waardoor ik niet<br />meer terug kan halen wat ik fout doe. Kan je misschien kijken wat er<br />mis gaat?<br /><br />Ik stuur de dataset en de Tin R file met het model die we gebruikt hebben<br />als attachment naar je e-mail account. De dataset gaat over voorkomen van insecten in 4 plantages, twee van ieder plant soort (Guava en Eucalyptus), gedurende 2<br />jaar. De data is verzameld door bomen in de plantages te bemonsteren en<br />daarom nemen we de zoektijd per boom als weegfactor voor de datapunt. Ik<br />probeer erachter te komen of de plantsoort effect heeft op het voorkomen<br />van de verschillende insecten groepen.<br /><br />Alvast bedankt,<br />Amir.</pre><br /><br />Here's Amir's code:<br /><br />setwd("D:/Amir-18-9-06/field-data/All data/analysis")<br />library(lme4)<br />data<-read.table("pdata.txt",header=T)<br />attach(data)<br />model<-lmer(Predators.Allspecies.~plant*julian.date+offset(log.time)<br />+(1+julian.date|place),family=poisson,method="Laplace")<br />model<-lmer(Predators.Allspecies.~place*julian.date+offset(log.time)<br />+(1+julian.date|place),family=poisson,method="Laplace")<br />model<-lmer(Predators.Allspecies.~place*julian.date+offset(log.time)<br />+(1+julian.date|place),family=quasipoisson,method="Laplace")<br />model<-lmer(Predators.Allspecies.~place*julian.date+offset(log.time)<br />+(1+julian.date|place/tree.code),family=poisson,method="Laplace")Amirhttp://www.blogger.com/profile/11171389993862330545noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1165484062485416762006-12-07T00:54:00.000-08:002009-11-13T15:48:51.700-08:00multinomHi Tessa,<br />To compare the fit of two models, you can first fit each one separately, and compare them using the anova() function. It will do a likelihood ratio test.<br />I think you have to construct your response differently, and not use "replica" to group observations.<br />The response variable has to be a matrix, with a row for each replica, and the counts on each patch in the columns. The explanatory variables will have as many entries as there are lines in the matrix. I now see that your covariates are patch-specific, which is more involved than what multinom() allows. Multinom() only allows an offset with that structure.<br />Maybe you will have to fit your data using a log-linear poisson model then, or function vglm() from the "VGAM" library can also handle multinomial data. Agresti, Categorical Data Analysis has a lot of info on the analysis of multinomial data.<br />Scripts for that book are available at <a href="http://www-stat.stanford.edu/%7Eowen/courses/306a/Splusdiscrete2.pdf">http://www-stat.stanford.edu/~owen/courses/306a/Splusdiscrete2.pdf</a>.<br />Cheers, TomTom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1165338474592345122006-12-05T09:07:00.000-08:002006-12-05T09:07:55.650-08:00Statistics Forum Functional Ecology Research School<a href="http://stats-feco.blogspot.com/"> </a><p class="MsoNormal">Dear Tom,</p> <p class="MsoNormal">I have tried to use the function multinom on my data. My experiment consists of multiple replicas of batches of 15 mites that distribute themselves over 5 patches. Thus, the response variable is multinomial distributed. I included 3 more explanatory variables as well (idad, nceggs and ideggs, assuming for now that they are independent).<o:p><br /></o:p></p> <p class="MsoNormal">My data looks like this (but then I have some more replica’s)</p> <p class="MsoNormal"><span style="" lang="PT-BR">replica<span style=""> </span>patch<span style=""> </span>ncad<span style=""> </span>idad<span style=""> </span>ideggs<span style=""> </span>nceggs<br />aaaa<span style=""> </span>a<span style=""> </span>2<span style=""> </span>0<span style=""> </span>0<span style=""> </span>2<br />aaaa<span style=""> </span>b<span style=""> </span>10<span style=""> </span>5<span style=""> </span>1<span style=""> </span>17<br />aaaa<span style=""> </span>c<span style=""> </span>0<span style=""> </span>1<span style=""> </span>0<span style=""> </span>2<br />aaaa<span style=""> </span>d<span style=""> </span>0<span style=""> </span>2<span style=""> </span>12<span style=""> </span>1<br />aaaa<span style=""> </span>e<span style=""> </span>0<span style=""> </span>7<span style=""> </span>4<span style=""> </span>2<br />bbbb<span style=""> </span>a<span style=""> </span>1<span style=""> </span>3<span style=""> </span>5<span style=""> </span>1<br />bbbb<span style=""> </span>b<span style=""> </span>9<span style=""> </span>3<span style=""> </span>1<span style=""> </span>8<br />bbbb<span style=""> </span>c<span style=""> </span>0<span style=""> </span>0<span style=""> </span>1<span style=""> </span>4<br />bbbb<span style=""> </span>d<span style=""> </span>2<span style=""> </span>0<span style=""> </span>11<span style=""> </span>6<br />bbbb<span style=""> </span>e<span style=""> </span>2<span style=""> </span>9<span style=""> </span>4<span style=""> </span>2</span></p> <p class="MsoNormal"><o:p> </o:p>My question are</p> <p class="MsoNormal">1. whether I have used the function multinom correctly? I have tried to find documentary which I could understand, but I failed. For example Venables and Ripley spent less then a page on it.</p> <p class="MsoNormal">2 how should I define “replica”. Now I have it as a factor, but I wonder if that is correct.</p> <p class="MsoNormal">3. how to do modelsimplification on the multinom function? Since the summary is a matrix and the drop1 function doesn’t work for multinom, I don’t know how to get to the minimal model. Maybe the most straightforward thing to do is start with the most simple model and add terms and see if the AIC decreases?<o:p></o:p></p> <p class="MsoNormal"><o:p></o:p>I will sent you my data file by email.<o:p><br /></o:p></p> <p class="MsoNormal">This is what I did:<o:p><br /></o:p></p> <p class="MsoNormal">library(MASS)<br />library(nnet)<br />setwd("n:\\Tessa\\experimenten\\SpatialdistrIdNc\\R")<br />mites=read.table("NcIdtryout.txt", header=T)<br />attach(mites)<br />names(mites)<br />head(mites)<br />model=multinom(patch~replica*idad*nceggs*ideggs, weight=ncad)<br />summary(model)</p><p class="MsoNormal"><o:p></o:p>Thanks in advance,</p> <p class="MsoNormal">Tessa<o:p></o:p></p> <a href="http://stats-feco.blogspot.com/">Statistics Forum Functional Ecology Research School</a>Tessahttp://www.blogger.com/profile/17419978293900752737noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1164886648228316422006-11-30T03:27:00.000-08:002006-11-30T03:37:28.250-08:00planned comparisonsDag Tom,<br /><br />het idee is om echt planned comparisons te doen, waardoor de error rate niet hoeft te worden aangepast (zoals wel gebeurt bij de simultaneous CI's methode).<br /><br />Ik vond devolgende manier (op 't net) om planned comparisons te doen, wat denk je daarvan?;<br /><br />cmat <- rbind("HL v LH" = c(0,1, -1, 0),"HH v LH/HL"=c(2,-1,-1,0), "LL v LH/HL"=c(0,-1,-1,2))<br />library(gmodels)<br />fit <-glm(ei2~kruising , contrasts=list("kruising"=make.contrasts(cmat)))<br />summary(fit)<br /><br />groet<br />Nicolanicolahttp://www.blogger.com/profile/02845233819066972346noreply@blogger.com2tag:blogger.com,1999:blog-22446344.post-1164750637076013722006-11-28T13:47:00.000-08:002006-11-28T13:50:37.096-08:00planned comparisonsHi Nicola, <br />you told me it's data on selection experiments. Are the data rally three blocks, meaning crosses between selection lines started from three different base populations?<br /><br />If you write it like this:<br />model<-aov(egg~block+block:cross)<br />there there are no main effects for cross, and you only want to compare crosses within blocks. That means that you only want contrasts between different groups in the same block.<br /><br />To compare these groups, first look at the coefficients of the model, and what they represent.<br />Reformulating the model like this probably makes comparisons easier:<br />model<-aov(egg~block:cross)<br />because all estimates are means per level.<br /><br />With four crosses per block, are the coefficients of block 1, the first four coefficients in the model?<br />Say, you want to make simultaneous confidence intervals for all pairwise comparisons for this block:<br /><br /># contrast matrix<br /><br />cmatrix<-matrix(rep(0,6*4),nrow=6)<br />cmatrix[1,1]<-1;cmatrix[1,2]<--1;<br />cmatrix[2,1]<-1;cmatrix[2,3]<--1;<br />cmatrix[3,1]<-1;cmatrix[3,4]<--1;<br />cmatrix[4,2]<-1;cmatrix[4,3]<--1;<br />cmatrix[5,2]<-1;cmatrix[5,4]<--1;<br />cmatrix[6,3]<-1;cmatrix[6,4]<--1;<br /><br />cmatrix<br /><br /># when this matrix is multiplied with the vector of coefficients, <br /># every element of the result, is an estimated difference between two means.<br /><br /># library(multcomp)<br /><br /># psubset selects the parameters of the block you want to investigate.<br /><br />simint(model,psubset=1:4,cmatrix=cmatrix)<br /><br />You should probably make a contrast matrix for all 18 pairwise comparisons simultaneously.<br /><br />cmatrix<-matrix(rep(0,3*6*4*3),nrow=3*6)<br />cmatrix[1,1]<-1;cmatrix[1,2]<--1;<br />cmatrix[2,1]<-1;cmatrix[2,3]<--1;<br />cmatrix[3,1]<-1;cmatrix[3,4]<--1;<br />cmatrix[4,2]<-1;cmatrix[4,3]<--1;<br />cmatrix[5,2]<-1;cmatrix[5,4]<--1;<br />cmatrix[6,3]<-1;cmatrix[6,4]<--1;<br /><br />cmatrix[7,5]<-1;cmatrix[7,6]<--1;<br />cmatrix[8,5]<-1;cmatrix[8,7]<--1;<br />cmatrix[9,5]<-1;cmatrix[9,8]<--1;<br />cmatrix[10,6]<-1;cmatrix[10,7]<--1;<br />cmatrix[11,6]<-1;cmatrix[11,8]<--1;<br />cmatrix[12,7]<-1;cmatrix[12,8]<--1;<br /><br />cmatrix[13,9]<-1;cmatrix[13,10]<--1;<br />cmatrix[14,9]<-1;cmatrix[14,11]<--1;<br />cmatrix[15,9]<-1;cmatrix[15,12]<--1;<br />cmatrix[16,10]<-1;cmatrix[16,11]<--1;<br />cmatrix[17,10]<-1;cmatrix[17,12]<--1;<br />cmatrix[18,11]<-1;cmatrix[18,12]<--1;<br /><br />simint(model,psubset=1:12,cmatrix=cmatrix)Tom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1164203931624688552006-11-22T05:31:00.000-08:002006-11-22T05:58:51.680-08:00planned comparisons in a nested ANOVAHi Tom and others,<br /><br />How do you perform planned comparisons within a nested ANOVA? I want to look at the differences between levels of the nested factor.<br /><br />When trying to build the appropriate contrast coefficient, I couldn't make matrices which were built from a combination of 2 factors (so main factor:nested factor). Is this the right way forward at all? And if so, how? A bookreference would be fine too!<br />Or would it be easier to just do t-tests?<br />Also, in the booklet of the course it was mentioned you could "construct appropriate independent variables" for planned comparisons. Is this the same as constructing contrast coefficients?<br /><br />Thanks in advance,<br />cheers<br />Nicolanicolahttp://www.blogger.com/profile/02845233819066972346noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1163588983788618022006-11-15T03:08:00.000-08:002006-11-15T03:09:43.800-08:00MultinomialHi Tessa,<br /><br /><br />Answers:<br /><br />1) how do I define my distribution in a glm when I have multinomial data (5 patches)?<br /><br />Well, the glm() function cannot handle multinomial data, but function multinom() from the "nnet" library can.<br />Another package which fits a multinomial probit model is "MNP", there is also the "multinomrob" library. I think multinom() is easiest to use.<br /><br /><br />2) the same question but for a generalized linear mixed model (I would like to include time as a factor)?<br /><br />When you need to add a random effect, "MNP" is the only option in R as far as I can see. Things become complicated in that case. Research on multinomial random effects models is still ongoing.<br />Some references: Multinomial logit random effects models, Statistical Modelling (2001) (J. Hartzel, A. Agresti, and B. Caffo); The analysis of ordered categorical data: An overview and a survey of recent developments, invited discussion paper for the Spanish Statistical Journal, TEST (2005) (I. Liu and A. Agresti).<br /><br />3) In Crawley they talk about GAM’s (general additive models). What I understand is that for these models you don’t need to define a distribution because they use non parametric smoothers. Can these models be used for multinomial data?<br /><br />Generalized additive models use non-parametric smoothers to model a non-linear dependence of a response variable on a predictor variable.<br />The distribution of the response variable is among the list of usual suspects: binomial, poisson, normal, gamma.<br />It has been claimed that package "gbm" (R) or "MART" (S-Plus) can do this for multinomial data, but I haven't tried those yet.<br /><br />4) Spatial distributions are difficult. Package "spatstat" could contain the solution. I am not very familiar with analysing spatial patterns.Tom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1162996853295793922006-11-08T06:40:00.000-08:002006-11-08T06:40:53.406-08:00multinomial distribution<a href="http://stats-feco.blogspot.com/"> </a><p class="MsoNormal">Hoi Tom and others</p> <p class="MsoNormal">I have a couple of questions about how to analyse my data. The data are on spatial distributions in 2 predatory mite species. The mites were given a choice between five patches. I would like to test for differences in distribution when both species are present and compare it with a control where only one species is present.</p> <p class="MsoNormal">Questions:</p> <p class="MsoNormal" style="margin-left: 0.5in; text-indent: -0.25in;"><!--[if !supportLists]--><span style="">1)<span style="font-family: "Times New Roman"; font-style: normal; font-variant: normal; font-weight: normal; font-size: 7pt; line-height: normal; font-size-adjust: none; font-stretch: normal;"> </span></span><!--[endif]-->how do I define my distribution in a glm when I have multinomial data (5 patches)?</p> <p class="MsoNormal" style="margin-left: 0.5in; text-indent: -0.25in;"><!--[if !supportLists]--><span style="">2)<span style="font-family: "Times New Roman"; font-style: normal; font-variant: normal; font-weight: normal; font-size: 7pt; line-height: normal; font-size-adjust: none; font-stretch: normal;"> </span></span><!--[endif]-->the same question but for a generalized linear mixed model (I would like to include time as a factor)?</p> <p class="MsoNormal" style="margin-left: 0.5in; text-indent: -0.25in;"><!--[if !supportLists]--><span style="">3)<span style="font-family: "Times New Roman"; font-style: normal; font-variant: normal; font-weight: normal; font-size: 7pt; line-height: normal; font-size-adjust: none; font-stretch: normal;"> </span></span><!--[endif]-->In <st1:place>Crawley</st1:place> they talk about GAM’s (general additive models). What I understand is that for these models you don’t need to define a distribution because they use non parametric smoothers. Can these models be used for multinomial data?</p> <p class="MsoNormal" style="margin-left: 0.5in; text-indent: -0.25in;"><!--[if !supportLists]--><span style="">4)<span style="font-family: "Times New Roman"; font-style: normal; font-variant: normal; font-weight: normal; font-size: 7pt; line-height: normal; font-size-adjust: none; font-stretch: normal;"> </span></span><!--[endif]-->The mites lay eggs during the experiment. I would like to know how the spatial distribution of the eggs, the eggs of the second species and the second species itself influence the spatial distribution of one mite. I thought about ancova, but my variables are not independent, so I suppose the correct way is to do several correlations? Or maybe multivariate analysis? I am also interested in interactions, so I wondered if it is possible to take this information in my model</p> <p class="MsoNormal" style="margin-left: 0.25in;"><o:p> </o:p></p> <p class="MsoNormal" style="margin-left: 0.25in;">Thanks in advance!</p> <p class="MsoNormal" style="margin-left: 0.25in;">Tessa</p> <p class="MsoNormal" style="margin-left: 0.25in;"><o:p> </o:p></p>Tessahttp://www.blogger.com/profile/17419978293900752737noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1162985838043994482006-11-08T03:30:00.000-08:002006-11-08T03:37:18.073-08:00random effectsHi Fransesc,<br />you seem to have only a small number of levels for the plot effects.<br />The consequence will be that estimates of variance between these random effects are very imprecise, and that it is not very solid to state that these levels are a representative sample from a large population of plot levels. <br />I would fit all effects as fixed.<br />Continuous variables such as time cannot define levels random effects, but they can appear in random regressions. Anyway, if you want to fit higher-order models of time, all you have to do is make new variables, like<br />timesq<-time^2<br />and add them to the model.<br />In a model with interactions of these time variables with treatment (or other factors), different polynomials of time will be fit for each treatment group.<br />Hope that helps,<br />TomTom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1162984780590891732006-11-08T01:37:00.000-08:002006-11-08T03:19:40.783-08:00mixed effectsHi Tom, and the rest of the experts ;-)<br /><br />I've been struggling quite a bit to fit a mixed effects model to my data. It is a RCB (or CR) design, with repeated measurements (RM) and, as said, both fixed and random effects. Treatment(fixed) = 2, plot (random) = 3 and nested in treatment. also there's replicate (random) = 2 and nested in plot, but this is the lowest level of variability, so I won't include this term.<br />Now here come the uncertainties: Time...I think is fixed...and treatment is also nested in time, because in the course of time the effect of the treatment seems to be different. The graph shows a large hump, which makes me suspect a higher-order term (time^2 or even time^3) . My problem is that I can't seem to find a relatively clear and simple explanation on how to deal with this. Crawley certainly doesn't treat it as it does the glm's. I reckon it goes a bit too far for a glm...<br />I hope I explained it clear enough...<br /><br />greetz,<br /> FrancescFranceschttp://www.blogger.com/profile/09910996248952030363noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1162811549747956262006-11-06T03:08:00.000-08:002006-11-06T03:12:29.756-08:00variable lengthsHi Cas,<br />something has gone wrong with attaching the data, or calculating new variables from old ones. The error message you get means that different variables in your model have different lenghts (so there's a mismatch).<br />for example, for a model of the form<br />glm(y~x+z)<br />you can check the lengths of x,z and y by<br />length(y)<br />length(x)<br />length(z)<br />The numbers returned should all be the same for a correct model specification. That makes it easy to spot where things went wrong.<br />good luck, TomTom VDhttp://www.blogger.com/profile/06116887666734281615noreply@blogger.com0tag:blogger.com,1999:blog-22446344.post-1162810883333191782006-11-06T02:58:00.000-08:002006-11-06T03:01:23.343-08:00Hi Tom and others of the GLM-R course,<br /><br />I am trying to run a GLM, but get the error saying that "variable lenghts differ (found for 'density')".The variable ''density" simply is a number from 3 to 7. Do you know what I am doing wrong?<br /><br />Cheers,<br /><br />Cas EikenaarCas Eikenaarhttp://www.blogger.com/profile/07999462243910385660noreply@blogger.com0