Re: [R] syntax for interactions in lme
> Host (fixed) > Sire (random) > Dam nested within Sire (random) > Host * Sire (random) > Host * Dam within Sire (random) > > So without the interactions I have: > > hogmodel = lme(gain ~ host, random = ~1|sire/dam) > > If I understand correctly, that "sire/dam" term gives me both > Sire and Dam within Sire as random factors. OK, so now I want > to add the two interactions (listed above)... Correct, for the interactions write: random = ~ host | sire/dam an interaction between a fixed and a random term can be interpreted as a variability in the fixed term for the different levels in the random term. Thus the 1 (which stands for the intercept) is exchanged with the fixed effect(s) with which interactions are of interest. This is easy if all the hierarchical levels of a nested random effects go into an interaction it is a bit more complicated if not. Say you only want the interaction of host and dam but not sire: random = list (~ 1 | sire, ~ host | dam) Hope this helps (it can all be found in Pinheiro and Bates and on the help pages). Lorenz __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] anova with models from glmmPQL
> I try to compare some models obtained from glmmPQL. > > model1 <- > glmmPQL(y~red*yellow+I(red^2)+I(yellow^2)+densite8+I(densite8^ > 2)+freq8_4 > +I(freq8_4^2), random=~1|num, binomial); > model2 <- > glmmPQL(y~red*yellow+I(red^2)+I(yellow^2)+densite8+I(densite8^ > 2)+freq8_4 > , random=~1|num, binomial); > anova(model1, model2) You try to compare models that differ in their fixed parts. This is not possible with the default method 'REML'. This would only be possible if you fitted your models using method 'ML' (see Pinheiro & Bates, 2000). In addition, if I understood a remark by José Pinheiro during one of his courses correctly, the anova comparisons are not save with a distribution other than normal. Thus, one should rely on the function intervals () to see whether confidence intervals of the parameters overlap zero or not. Perhaps someone else can comment on this issue? Regards, Lorenz - Lorenz Gygax Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] interpretation output glmmPQL
> We study the effect of several variables on fruit set for 44 > individuals (plants). For each individual, we have the number > of fruits, the number of flowers and a value for each variable. > ... > - Glm does not take account of the correlation between the > flowers of a unique individual. So we would like to add a > random effect 'individual' but the model2 (here after) gives an > output similar to the one of model1 for estimated coefficients > and p-values. > ... > Does it mean that there is no individual effect or is my model > not good (number of groups (individuals)=number of observations, > is it possible?). If you have only one observation per indiviudal plant, how could there be dependence within the plant? This would only make sense if your observations were the individual flowers. Data on those could be correlated within plant and then a random term for the plant is meaningful. Cheers, Lorenz - Lorenz Gygax Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] The mathematics inside lme()
> Now I want to evaluate GroupCov as a covariate to Treat. I > can then start with either m1 or m2 as base, but what is most > correct when GroupCov has only one value for each Group? > > m3 <- lm(Yield ~ Treat + GroupCov + Treat:GroupCov) > > gives the same fixed effects as > > m4 <- lme(Yield ~ Treat + GroupCov + Treat:GroupCov, > random =~1| Group) > > but this time the prob.values for GroupCov is much stronger > in m3 than in m4. Needless to say, anova(m3,m4) tells that m4 > is a better *model* than m3. But is it better for my purpose? Well, I do not actually know what your purpose is ... ... but in my oppinion the second model is much better (and I am tempted to say the first one is wrong). The crucial point here is that there is only one value of GroupCov in each Group. Thus the number of replications that provide degrees of freedom for the effect of GroupCov is the number of groups. m4 adjusts for this fact, has a lower df for the GroupCov and thus a lower p-value. In m3, you model as if all your observations are independent for all variables. This is actually the case for none but may become more visible for GroupCov because this variable is constant for all units within group (and thus this value is certainly not independent). Cheers, Lorenz - Lorenz Gygax, Dr. sc. nat. Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] testing equality of variances across groups in lme?
> I am fitting a two-level mixed model which assumes equality of > variance in the lowest-level residuals across groups. > The call is: > > fit3<-lme(CLnNAR~CLnRGR,data=meta.analysis, > + na.action="na.omit",random=~1+CLnRGR|study.code) I assume that CLnRGR is a factor and thus the groups which might have different variance in residuals are given by the different levels of CLnRGR!? > ... > know that one can test a nested model in which the residuals are > given weights, such as: You are on the right track here, just use: fit4<-lme(CLnNAR~CLnRGR,data=meta.analysis, na.action="na.omit", random=~CLnRGR|study.code, weights= varIdent (form= ~ 1 | CLnRGR)) with anova (fit3, fit4) you can test whether these weights improve the fits statistically. Lorenz - Lorenz Gygax, Dr. sc. nat. Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] logistic regression 3D-plot
> I tried to create a 3D surface showing the interaction between two > continuous explanatory variables; the response variable is > binary (0/1). > > The model is: > > model<-glm(incidence~sun*trees,binomial) > > then I used "wireframe" to create a 3D plot: > > xyz<-expand.grid(sun=seq(30,180,1),trees=seq(0,4000,10)) > > xyz$incidence<-as.vector(predict(model,xyz)) xyz$incidence<-as.vector(predict(model,xyz, type= "response")) should work > wireframe(incidence~sun*trees,xyz,scales=list(arrows=FALSE)) Cheers, Lorenz __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Split-split plot ANOVA
> I have been going over and over the examples in MASS > and the Pinheiro and Bates example, but cannot get my model > to run correctly with either aov or lme. > > Could someone give me a hand with the correct model statement? It would help to see some of the things you have tried already ... > First a description of the design. We are studying > germination rates for > various species under a variety of treaments. This is a > blocked split-split > plot design. The levels and treatments are: > > Blocks: 1-6 > > Whole plot treatment: >Overstory: Yes or No > > Split plot treatments: >Caging (to protect against seed predators): Yes or No >Herbaceous competition (i.e., grass): Yes or No > > Split-split plot treatment: >Tree species: 7 kinds > > The response variable is Lag, which is a indication of when > the seeds first germinated. I would try somthing like lme (fixed= Lag ~ Caging + herbaceous + tree, data= your.data, random= ~ 1 | Overstory/split/splitsplit) Perhaps you want/need to add some interactions as well. Overstory, split and splitsplit would be factors with specific levels for each of the plots, split plots and split-split plots, respectively. Thus what I attempted here is to separate the variables of the hierarchical design of data gathering (which go into the random effects) and the treatments (which go into the fixed effects). The degrees of freedom for the fixed effects are automatically adjusted to the correct level in the hierarchy. Did you try that? What did not work out with it? > Lastly, I have unbalanced data since some treatment > combinations never had any germination. In principle, the REML estimates in lme are not effected by unbalanced data. BUT I do not think that the missing germinations by themselves lead to an unbalanced data set: I assume it is informative that in some treatment combinations there was no germination. Thus, your lag there is something close to infinity (or at least longer than you cared to wait ;-). Thus, I would argue you have to somehow include these data points as well, otherwise you can only make a very restricted statement of the kind: if there was germination, this depended on such and such. > Since the data are highly nonnormal, I hope to do a > permutations test on the F-values for each main effect and > interaction in order to get my p-values. As these are durations a log transformation of your response might be enough. Regards, Lorenz - Lorenz Gygax, Dr. sc. nat. Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] binomia data and mixed model
> The experimental design was suppose to consist of 4
> treatments replicated 3 time, Source 1 and applied at 10 cm
> and source 2 applied at 20 cm. During the construction of the
> treatmetns the depths vary considerably so i can't test all my
> samples based on 10 and 20 cm any more the depths are now
> considered random and not fixed.
It is not really clear what you measuring ... but do you know the depth even
if it is not exactly 10 and 20 cm? Then, perhaps you could use this variable
in a continous fashion but still as fixed? I do not really see how you can
treat such a variable as random.
> The data is very non-normal (lots of zeros) therefore the only way
> to analyze it is to convert to binomial data.
I doubt this is the only way but certainly a valid one.
> Does any one know what type of analysis I should use? I was told
> that a NLmixed model would work but also a GLIM mixed model was
> appropriate. Is there any info using these in R.
I do not know about those but you can conduct binomial mixed effects models
by either using glmmPQL in library 'MASS' or GLMM in library 'lme4'.
Also do read:
author = {Pinheiro, Jose C and Bates, Douglas M},
title ={Mixed-Effects Models in {S} and {S}-{P}{L}{U}{S}},
publisher ={Springer},
year = {2000},
address = {New York}
and
AUTHOR = {Venables, W N and Ripley, B D},
TITLE = {Modern Applied Statistics with {S}},
PUBLISHER = {Springer},
YEAR = {2002},
ADDRESS = {New York},
EDITION = {fourth}
Regards, Lorenz
-
Lorenz Gygax, Dr. sc. nat.
Centre for proper housing of ruminants and pigs
Swiss Federal Veterinary Office
agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Special paper for postscript
> When I generate a "special" paper postscript image larger than > "a4" or "letter" using R, I can only see one-page portion of all > image, of course. > What will be the simple solution for this? Is there any way I > can set the bounding box information on the image? Or any other > suggestions? I often use 'special' paper for graphs that are smaller than A4. I believe that the postscript files thus created do have the correct bounding box. If you import such a graph it shows the right boundaries. Alternatively if I look at such a graph in ghostview I can set the boundaries myself in the Menu Media->user defined according to the bounding boxes in the postscript file. Thus, I am not sure what your problem is ... Regards, Lorenz - Lorenz Gygax, Dr. sc. nat. Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] random effects in lme
It is unlikely that your request will be answered faster if you post it a second time in exactly the same way ... > Suppose I have a linear mixed-effects model (from the package > nlme) with nested random effects (see below); how would I present > the results from the random effects part in a publication? > > Specifically, I´d like to know: > (1) What is the total variance of the random effects at each level? > Random effects: > Formula: ~ 1 | plotcode >(Intercept) > StdDev: 0.04176364 > > Formula: ~ 1 | treatment %in% plotcode > (Intercept) Residual > StdDev: 0.08660458 0.00833387 What is wrong with an estimted StdDev on the level of plotcode of 0.0418 and on the level of treatment (which is actually the same as to say that this is a block of plants within plotcode that received the same treatment?) of 0.087? > (2) How can I test the significance of the variance components? Calculate a model with an without the component of interest and compare the models using anova(). Perhaps you should read Pinheiro & Bates (2000)? > (3) Is there something like an "r squared" for the whole > model which I can state? None is provided and I do not know how easily it could be calculated. But the use of the measure can be questioned. It is an absolute measure on how much of the variability in the data is explained. Imagine you had true replicates (i.e. several measurements under absolutely identical situations). Imagine further that these measurements do nevertheless show some variability. If this variability was substantial, your r-squared would be low even though your model might pick up all the structure that you can find in the data. Thus you can only get as good as 'natural' variability in your data which is not considered by the r-squared measure. Thus, (corrected) r-squared values might tell you something if you compare different models based on the same data (in a similar way as the AIC and BIC criteria) but not if you compare completely different data sets. Regards, Lorenz - Lorenz Gygax, Dr. sc. nat. Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] lme and varFunc()
> I am currently analyzing a dataset using lme(). The model I
> use has the following structure:
>
> model<-lme(response~Covariate+TreatmentA+TreatmentB,
>random=~1|Block/Plot,method="ML")
>
> When I plot the residuals against the fitted values, I see a clear
> positive trend (meaning that the variance increases with the mean).
>
> I tried to solve this issue using weights=varPower(), but it
> doesn´t change the residual plot at all.
You are aware that you need to use something like
weigths= varPower (form= fitted (.))
and the plot residuals using e.g.
scatter.smooth (fitted (model), resid (model, type= 'n'))
Maybe the latter should also be ok with the default pearson residuals, but I
am not sure.
Possibly a look into the following would help?
@Book{Pin:00a,
author = {Pinheiro, Jose C and Bates, Douglas M},
title ={Mixed-Effects Models in {S} and {S}-{P}{L}{U}{S}},
publisher ={Springer},
year = {2000},
address = {New York}
}
> How would you implement such a positive trend in the variance? I´ve
> tried glmmPQL (which works great with poisson errors), but
> using glmmPQL I can´t do model simplification.
Well, what error distribution do you have / do you expect?
Regards, Lorenz
-
Lorenz Gygax, Dr. sc. nat.
Centre for proper housing of ruminants and pigs
Swiss Federal Veterinary Office
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] 3d bar plot
> This graph ->
> http://www.math.hope.edu/~tanis/dallas/images/disth36.gif
> is an example I found at
> http://www.math.hope.edu/~tanis/dallas/disth1.html
> created by Maple.
>
> Does anybody know how to create something similar in R?
>
> I have a feeling it could be possible using scatterplot3d
> (perhaps with type=h, the fourth example in help('scatterplot3d')?),
> but I cannot figure it out.
Sorry to butt in with a more fundamental question. Is this really the kind
of graph we want to cultivate and support? In my oppinion, it is hardly ever
necessary to have a graph in 3D or even in higher dimensions (one certain
exception is if one tries to spin a higly dimensional dataset in search of
patterns as you can do in ggobi and there might certainly be others).
At least the graph presented in the example does - in my eyes - not warrant
a 3D plot. Why not just draw curves for each of the n's in a plot of 'A'
against 'row'? This would enable a reader to make straightforward
comparisons of the curves and allow to estimate the height of the 'columns'
along the 'A'-axis much more easily.
Only because we can easily create 3D graphs, I do not believe that we should
use them often. Only if a careful evaluation of alternatives was not
promising success I would resign myself to using 3D graphs.
Lorenz
-
Lorenz Gygax, Dr. sc. nat.
Centre for proper housing of ruminants and pigs
Swiss Federal Veterinary Office
agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Multiple comparisons following nlme
> I need to do multiple comparisons following nlme analysis (Compare > the effects of different treatments on a response measured > repeatedly over time; > fixed = response ~ treat*time). If you have an interaction it does not really make sense to conduct a multiple comparison because the difference in treatment depends on time, i.e. this presumed post-hoc test could only give you a correct result for one of your points in time. Why not conduct this analysis and then interpret the pattern based on the estimates of your parameters and/or on a graphical display of your data? If your interaction is non-significant and you drop it, you still have a mulitvariate problem and I have never come across a multiple comparison test for such multivariate problems (but perhaps someone else has a pointer). In your case it might be enough to carefully decide on how the set contrasts. Then, an additional issue would be what kind of multiple comparisons to conduct (for univarite anova's there are at least a dozen methods). You can always conduct several to see which of the comparisons are highly significant and which ones might not be so strong. But usually you do not learn more than what you get when you interpret your parameters and/or graphs of your data. ... and by the way, I guess your model is using lme (linear mixed effects model) in package nlme and not actually an nlme (non-linear mixed effects model) itself. Regards, Lorenz - Lorenz Gygax, Dr. sc. nat. Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office, agroscope FAT Tänikon __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] GLMM and crossed effects
> I am analysing data with a dependent variable of insect > counts, a fixed effect of site and two random effects, day, > which is the same set of 10 days for each site, and then > transect, which is nested within site (5 each). And what exactly are you interested in? Just the differences between the two sites? Then you would have sampled on several days just to better estimate what is going on at those sites? Or are you also interested in a time course? > I am trying to fit the cross classified model using GLMM in lme4. I have yet to get around to working with lme4 (by the way is there some documentation that describes the changes and advantages of lme4 in comparison to nlme?), but assuming that the syntax is the same as in lme () or glmmPQL (): Why don't you use: GLMM (count ~ site, data= dat3, random= ~ 1 | site/trans, family= poisson) This compares your counts between sites considering that transects are nested in sites and that there are several measurements on each transect (the days). They are the repeated measurements on the lowest level (within trans), so you do not need to specify them explicitly. Acutally, you might not even need the site/ in your random term as this variable ist constant for each transect and thus the degrees of freedom are adjusted for the fixed effect of site (but then each of the transect need its own code). Well on further thought, you might necessarily need to leave this part of the random term away if you want to conduct any statistical test, because otherwise you only have an N= 1 for each site ... thus, I guess you need to assume that your transects are independent measures of your two sites (which means that you can conduct some statistical analysis, but you results hold only for these two specific sites and can not necessarily be generalised to other sites). Thus the model would be: GLMM (count ~ site, data= dat3, random= ~ 1 | trans2, family= poisson) or if you are interested in a time course you might try (and this explicitely models that these are the same days): GLMM (count ~ site*day, data= dat3, random= ~ 1 | trans2, family= poisson) I would argue, that you are either not interested in the days, then these are just your repeated measurements and it does not matter that they were exactly on the same days for the different transects (and then they are just implicitly nested in trans2) OR you are interested in them and then I would include them as a fixed effect, which is crossed with transect, i.e. all transects were sampled on all days. By the way, it is not clear to me from your description how trans2, ts and ts2 differ logically. > there are als ts1 and ts2 dummy variables, as was necessary > in the old lme. what are these necessary for? (But I have to admit that I usually feed a standard data.frame to lme and glmmPQL) > GLMM(count~site,data=dat3,random=list(day=~1,trans=~1|site, > family=poisson) I do not know whether you can write such a thing at all. If this has not changed a lot from nlme to lme4 you would need to write: random= list (~1 | day, ~ 1 | site/trans) but that you would implicitly define that site is nested in day, i.e. it would be the same as writing random= ~ 1 | day/site/trans which is not what you want. > #This does... but also note the differences in the summary > and VarCorr variance components... Here, you loose me completely. It is not clear to me what you compare and where you perceive a problem. Cheers, Lorenz - Lorenz Gygax, Dr. sc. nat. Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] The hidden costs of GPL software?
> So, is this analysis correct: are there hidden costs for free > software like R in the time required to learn it? At least > currently, for the people I know (biologists, ecologists, > oceanographers, ...), this is perfectly true. This is even an > insurmountable barrier for many of them I know, and they > have given up (they come back to Statistica, Systat, or S-PLUS > using exclusively functions they can reach through menus/dialog > boxes). I guess you are right, in that the steep initial learning curve could be smoothed for beginners. On the other hand I do not see how a GUI for R could cover more than the bare essentials because the available functionality is so vast. We also have S-Plus at our research institution and even there, I see, that people who do not know about the underlying code have difficulties in using the GUI. I personally believe that it is more a question how one is used to do statistics. Click and drag is the norm. (And I guess it is usually also the norm of how people/scientists use other Software.) In my eyes, using code instead, means that one is able to repeat the steps of an evaluation easily and to document at the same time what has been done. Very soon evaluations (and data handling) can be done far more efficiently than with click and drag. All these advantages outweigh the initial costs by several orders of magnitude. Thus, in my opinion it is more a question of education such that people might realize how they can work efficiently and cleanly. Perhaps one could even say that such an approach is more scientific because, in principal, it can be easily communicated and reproduced. It is, of course, easy for me to make these statements, as in the meantime I have been using S (S-Plus and R) for - gosh - over 10 years. But I see in some projects that I supervise that people get started easily with a snippet of code that I provide and the insight of the usefulness of such a work approach is usually easily within reach. Lorenz - Lorenz Gygax, Dr. sc. nat. Tel: +41 (0)52 368 33 84 / [EMAIL PROTECTED] Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] grouping for lme with nested repeated measurements
> 1.The simple problem ist that i have different Samples , from > which i make repeated measurements (each sample is measured 6 > times) and i repeat this experiment over several Days, so i > get the lme grouping term "random=~1|Days/Sample". I would rather specify this as: random= ~ 1 | Sample/Days, but I am not quite sure how this affects the model. > 2. Now i am measuring with 2 different measuring Apparatus > the same Sample each 6 times, to see how big the difference > from the appratus is. > Because Apparatus is on the level of repeated measurements i > cant write Days/Sample/Apparatus. the lme function offers a > list() feature to design the grouping, but i didnt understand > this, if it is the solution to the problem. Why would you want to include the Apparatus in the random effect? I assume that you are interested in differences and thus, this is a fixed effect: lme (fixed= response ~ apparatus, data= XX, random= ~ 1 | Sample/Days) Lorenz - Lorenz Gygax, [EMAIL PROTECTED] Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
FW: [R] glmmPQL and random factors
I have just realised that I sent this to Per only. For those interested on the list: -Original Message- From: Gygax Lorenz FAT Sent: Tuesday, September 14, 2004 4:35 PM To: 'Per Toräng' Subject: RE: [R] glmmPQL and random factors Hi Per, > glmmPQL(Fruit.set~Treat1*Treat2+offset(log10(No.flowers)), > random=~1|Plot, family=poisson, data=...) > > Plot is supposed to be nested in (Treat1*Treat2). > Is this analysis suitable? As far as I understand the methods and with my experience using such analyses, I would say that the model is ok the way you specified it. glmmPQL (and the underlying lme) is so intelligent (a thousand thanks to the developpers!) as to recognise if the treatments are fixed per plot, i.e. only one level of the two treatments appears in each plot. The denominator degrees of freedom in the anova table are adjusted automatically. I.e. your denominator df should be the number of plots minus five, the number of dfs you need for the fixed effects (Treat1, Treat2, the interaction, the covariate and the one df you always loose from the total of observations). > Moreover, what is the meaning of typing > random=~1|Plot compared to random=~Treat1*Treat2|Plot? The first version means, that the intercept / overall mean can vary from plot to plot. I.e. each plot may have another mean due to the fact that it grows somewhere else in addition to the differing treatments. The second version tries to model a difference in reaction to treatment 1 and 2 for each of the plots (which does not make sense in your case as each plot is only subjected to one kind of treatment). In a crossed design, i.e. if you could have treated your plants individually and had all treatment combinations in each of the plots, the first version implies that all the plots react in the same consistent way to the treatments. I.e. that the general level of each plot may be different, but the differences due to treatment are the same in each plot, the reaction of the plots are shifted but have the same shape (this is the same as saying that you only consider main effects of treatment and plot). The second version allows to estimate the reactions for each plot, i.e. in addition to a general shift, the treatments may have (slightly) different effects in each plot. This is the same as saying that you consider interactions between your fixed and random effects. See also the terrific book by Pinheiro & Bates (Mixed Effects Modelling in S and S-Plus, Springer, 2000). Cheers, Lorenz - Lorenz Gygax Tel: +41 (0)52 368 33 84 / [EMAIL PROTECTED] Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] drawing on axes
> I would like to make certain portions of axis lines thicker (as Tufte > suggests). How can I draw on axes? I only need a couple of line > segments on the left and bottom one. How about: plot (1, 1, xlim= c (-10, 10), bty= 'n') axis (1, at= c (-10, -5), labels= FALSE, tick= T, lwd= 5, tck= 0) axis (1, at= c (0, 3), labels= FALSE, tick= T, lwd= 5, tck= 0) axis (2, at= c (0.8, 1), labels= FALSE, tick= T, lwd= 5, tck= 0) Cheers, Lorenz - Lorenz Gygax Tel: +41 (0)52 368 33 84 / [EMAIL PROTECTED] Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] mixed models question
Hi Chris, > I am trying to fit the following linear model to logged per capita > fecundity data (ie number of babies per female) for a mouse: > > RsNRlS <- glm(formula = ln.fecundity ~ summer.rainfall + N + > lagged.rainfall + season, ) > > I am using this relationship in a simulation model, and the current > statistical model I have fit is unsatisfactory. The problem is I get a > global estimate of variance (MSE), but I think it varies across subsets > of the data. Specifically, seasons when there is lots of reproduction > (e.g. fall) tend to have high variance, while seasons with little > reproduction (e.g. summer) have small amounts of variance. I am > looking for a method for estimating the coefficients in my linear > model, and estimating a separate error for subsets of the data (ie for > each of the 4 seasons). The end goal is to take this linear model back > into my simulation model to make predictions about fecundity, but with > separate variance terms for subsets of the data. Are you using glm because you need a specific distribution family (such like poisson)? If not, you could possibly use gls with the argument weights= varFixed (~ season) With that you estimate your parameters and at the same time you allow for (and estimate) the different variances for the season. If you need the poisson distribution, I am not quite sure what to do. Perhaps glm also accepts this weight argument or perhaps you need to work with a generalised procedure of lme (either from one of the new lme packages or from MASS). Regards, Lorenz - Lorenz Gygax, Dr. sc. nat. Tel: +41 (0)52 368 33 84 / [EMAIL PROTECTED] Center for proper housing of ruminants and pigs Swiss Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland Fax : +41 (0)52 365 11 90 / Tel: +41 (0)52 368 31 31 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] multiple error strata in aov
Hi Murray > I am trying to perform a model 3 ANOVA for a 2 factor (say factor A and > factor B) anova in which factor A is fixed and factor B is random. > ... > In addition, I have tried using lme to perform this function, > but again without much success. What did not work? And, did you read Pinheiro & Bates? Personally, I find it easier to work with lme and this should be an easy one. What about lme (fixed= Y ~ A, random= ~ 1 | B) or lme (fixed= Y ~ A, random= ~ A | B)? Regards, Lorenz - Lorenz Gygax, Dr. sc. nat. Tel: +41 (0)52 368 33 84 / [EMAIL PROTECTED] Center for proper housing of ruminants and pigs Swiss Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland Fax : +41 (0)52 365 11 90 / Tel: +41 (0)52 368 31 31 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] lme newbie question
Hi Christoph, > Am I right in this lme implementation, when I want to investigate the > influence of the age.group, and the two conditions on the rt: > > my.lme <- lme(rt ~ age.group + angles * hands, data = > my.data, random = ~ 1 |subject) > > then I think I would have to compare the model above with a more > elaborated one, including more interactions: > > my.lme2 <- lme(rt ~ age.group * angles * hands, data = > my.data, random = ~ 1 |subject) > > and comparing them by performing a likelhood-ratio test, yes? If you compare these two models, you test whether the interactions of age.group with angles and hands and the three-way interaction all together make for a significant improvement of the model. Is that what you want? Also note: if you do this, you need to use the method ML and not the default REML. Or you start with the second model, use anova (my.lme2) and reduce the model stepwise. You can also ask whether there is a subject to subject variability in variables other than the intercept (i.e. interactions between your random and the fixed variables) and e.g. try things like random = ~ age.group + angles * hands | subject > I think, if I would like to generalize the influence of the experimental > conditions on the rt I should define angles and hands as a random effect, > yes? I do not see exactly what you aim at, here. Possibly, the second part of my answer above is an answer to this as well? > thanks for a short feedback. It seems, repeated-measures > anova's aren't a trivial topic in R :) They never are. But, after having read most of Pinheiro and Bates' book 'Mixed effects modelling in S and S-PLUS' (Springer), it seems easier to me than ever, because they use a consistent, integrated and concise approach. Regards, Lorenz - Lorenz Gygax, Dr. sc. nat. Tel: +41 (0)52 368 33 84 / [EMAIL PROTECTED] Center for proper housing of ruminants and pigs Swiss Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland Fax : +41 (0)52 365 11 90 / Tel: +41 (0)52 368 31 31 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] dfs in lme
Dear R-mixed-effects-modelers,
I could not answer this questions with the book by Pinheiro & Bates and did
not find anything appropriate in the archives, either ...
We are preparing a short lecture on degrees of freedom and would like to
show lme's as an example as we often need to work with these. I have a
problem in understanding how many dfs are needed if random terms are used
for explanatory variables in addition to the intercept (if I have understood
correctly that ist the same as saying that interactions between random and
fixed effects are considered). I tried the following code:
library ('nlme')
options (contrasts= c ('contr.treatment', 'contr.poly'))
# create fake data
data.df <- data.frame (gruppe= rep (1:4, rep (20, 4)))
# create response variable
data.df$zv <- rnorm (80, 2)
# create potential explanatory variables
data.df$explan <- rnorm (80, 2)
data.df$treat <- as.factor (sample (1:3, 80, T))
data.df$treat1 <- as.factor (sample (1:4, 80, T))
data.df$treat2 <- as.factor (sample (1:5, 80, T))
data.df$treat3 <- as.factor (sample (1:6, 80, T))
# with each of the explanatory variables
withoutInt <- lme (zv ~ explan, data= data.df, random= ~1 | gruppe)
withInt <- lme (zv ~ explan, data= data.df, random= ~ explan | gruppe)
anova (withoutInt)
anova (withInt)
anova (withoutInt, withInt)
There are two main things that I wonder about:
(1) the two anova() commands on the single models have the same residual
degrees of freedom even though the model withInt has an additional
explanatory variable. Why are the residual dfs not reduced?
(2) In the model comparison, it becomes visible that the model with 'explan'
in the random effect does indeed use more dfs. But I cannot figure out where
that number of dfs comes from. I thought that basically for each of the
levels in the grouping variable additional parameters are estimated? Thus, I
would expect somethind like df(interaction)= df(explanatory
variable)*df(random effect), but what I find is:
explanatory variabledelta-dfs of the model comparison
(= dfs of the interaction of the explanatory
variable with the random effect 'gruppe',
which has 4 levels, 3 dfs)
continuous (1 df)2
3 levels (2 dfs) 5
4 levels (3 dfs) 9
5 levels (4 dfs)14
6 levels (5 dfs)20
Can anyone point me in the right direction on where and how to answer these
questions?
Many thanks and regards, Lorenz
-
Lorenz Gygax
Tel: +41 (0)52 368 33 84 / [EMAIL PROTECTED]
Center for proper housing of ruminants and pigs
Swiss Veterinary Office
agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland
Fax : +41 (0)52 365 11 90 / Tel: +41 (0)52 368 31 31
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] bug in cor (..., use= ...)?
Dear R users, I have not found anything on this in the archives. Does anyone know whehther the parameter use= is not functioning in cor or enlighten me what it is supposed to do? My R version is "R version 1.8.1, 2003-11-21" on Windows 2000. I am hoping to be able to update to 1.9.1 as soon as it has appeared (we are not allowed here to install software on our own and thus I am trying to be able to have the .1 versions installed ...). Test code: x <- 1:10 y <- 2:11 x [1] <- NA y [10] <- 12 cor (x, y, use= 'all.obs', method= 'kendall') cor (x, y, use= 'complete.obs', method= 'kendall') cor (x, y, use= 'pairwise.complete.obs', method= 'kendall') As I understand, the first one of this should result in an error which it does not. All the results are the same and seemingly treat the NA as if it was 0. Any ideas are appreciated. Thanks and regards, Lorenz - Lorenz Gygax, Dr. sc. nat. Tel: +41 (0)52 368 33 84 / [EMAIL PROTECTED] Tag der offenen Tür, 11./12. Juni 2004: http://www.fat.ch/2004 Center for proper housing of ruminants and pigs Swiss Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland Fax : +41 (0)52 365 11 90 / Tel: +41 (0)52 368 31 31 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] glmmPQL - predict (..., type= 'response')?
Dear R users, Is there something like predict (..., type= 'response') for glmmPQL objects or how would I get fitted values on the scale of the response variable for the binomial and the poisson family? Any pointers are appreciated. Thanks, Lorenz - Lorenz Gygax, Dr. sc. nat. Tel: +41 (0)52 368 33 84 / [EMAIL PROTECTED] Tag der offenen Tür, 11./12. Juni 2004: http://www.fat.ch/2004 Center for proper housing of ruminants and pigs Swiss Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland Fax : +41 (0)52 365 11 90 / Tel: +41 (0)52 368 31 31 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Repeated measures regression
Why not start with:
@Book{Pin:00a,
author = {Pinheiro, Jose C and Bates, Douglas M},
title ={Mixed-Effects Models in {S} and {S}-{P}{L}{U}{S}},
publisher ={Springer},
year = {2000},
address = {New York}
}
Regards, Lorenz
> -Original Message-
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] Behalf Of knussear
> Sent: Wednesday, May 05, 2004 5:59 AM
> To: [EMAIL PROTECTED]
> Subject: [R] Repeated measures regression
>
>
> Hi List,
>
> Just wondering if there is such a thing as repeated measures
> regression, and if so, can R do it?
>
> I have repeated measurements of 10 individuals over a 45 day period,
> and I would like to regress their daily activity time against a daily
> environmental temperature. If I do so using averages of
> activity time I
> find a significant negative correlation, but I worry that because I
> have used the same 10 individuals for each daily mean that the daily
> averages are not independent.
>
> Can anyone help?
>
> Thanks
>
> Kenneth E. Nussear Phone 775 784-4565
> Biological Resources FAX 775 784-1369
> Research Center/315 [EMAIL PROTECTED]
> Reno, Nevada 89557http://www.brrc.unr.edu/~knussear/
>
> __
> [EMAIL PROTECTED] mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Mixed-effects model for nested design data
Dear Han, > I am using nlme for data from nested design. That is, "tows" are nested > within "trip", "trips" nested within "vessel", and "vessels" nested > within "season". I also have several covariates, say "tow_time", > "latitude" and "depth" > My model is >y = season + tow_time + latitude + depth + vessel(season) + > trip(season, vessel) + e > In SAS, the program would be > proc mixed NOCLPRINT NOITPRINT data=obtwl.x; > class vessel trip tow season depth; > model y = season depth latitude /solution; <--fixed effects > random vessel(season) trip(season vessel); > run; > My question is: How this nested mixed-effects model can be > fitted in R-> "nlme"? I do not know about SAS but I would guess that your model should be fitted as something like: lme (fixed= y ~ season + tow_time + latitude + depth, random= ~ 1 | season/vessel/trip) Maybe you should do some reading in the book by Pinheiro & Bates? They explain well how to set up models. Regards, Lorenz - Lorenz Gygax, Dr. sc. nat. Tel: +41 (0)52 368 33 84 / [EMAIL PROTECTED] Tag der offenen Tür, 11./12. Juni 2004: http://www.fat.ch/2004 Center for proper housing of ruminants and pigs Swiss Veterinary Office agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland Fax : +41 (0)52 365 11 90 / Tel: +41 (0)52 368 31 31 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] 3D histogram
Hi Tine,
> I have been trying to look for documentation on making a 3D
> histogram in R (so a histogram on which there are several distributions
> are plotted), but I cannot find anything on it. Is it possible to
> construct such histograms and should use a special package (f.e. R.Basic)?
In my opinion you should never use 3D graphs. Due to the perspective on
these graphs it is always difficult to compare different values of the
graphed data and usually it is also difficult to get a good estimate of the
actual values as the axes are not always aligned with the data of interst.
You might want to ask yourself whether your data is really as simple that it
can be shown in bars (or could you e.g. rather use boxlplots). If you want
to graph your data as in barplots, you could plot the values of your bars
using points in a 2D diagram and e.g. connect the points that would have
been in one row/column in the 3D diagram.
If you are interested in the graphical presentation of quantitative data, I
can only recommend:
@Book{Tuf:99,
author = {Tufte, Edward R},
title ={The Visual Display of Quantitative Information},
publisher ={Graphic Press},
year = {1999},
address = {Cheshire, Connecticut},
edition = {17th printing}
}
Hope this helps. Lorenz
-
Lorenz Gygax, Dr. sc. nat.
Tel: +41 (0)52 368 33 84 / [EMAIL PROTECTED]
Tag der offenen Tür, 11./12. Juni 2004: http://www.fat.ch/2004
Center for proper housing of ruminants and pigs
Swiss Veterinary Office
agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland
Fax : +41 (0)52 365 11 90 / Tel: +41 (0)52 368 31 31
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] modelling nested random effects with interactions in R
Hi Michael,
> I have a nested anova, with random factor "lakefac" nested within
> factor "fishfac" (fixed), with an additional fixed factor
> "Habfac". If I consider everything as fixed effects, it's addmittedly
> not the correct model, but I can at least get this to work:
> ...
> So now I try to run it using the lme4 package, treating
> lakefac as random;
I am not quite sure about lme4 but in limrary ('nlme') you would need to do
something like:
> lakernd <- lme(ltotinv ~ habfac * fishfac,
random = ~ 1 | fishfac/lakefac,
data= use)
with this model you have the fixed effect of habfac and fishfac and their
interaction, lme should get the df's of the model ok if you specify it that
way. In this way, the random term is only for the intercept (which to my
understanding is the same as saying there is no interaction between random
and fixed effects). If you want to include the interactions you need to do
something like this:
> lakernd.int <- lme(ltotinv ~ habfac * fishfac,
random = ~ habfac * fishfac | fishfac/lakefac,
data= use)
Thus, you specify random effects for the habfac, fishfac and their
interactions, which is the same as saying that there are interactions
between random and fixed effects.
Possibly you might need to do something like the following (because fishfac
is in the fixed and the random term, thus you specify interactions
seperately for both levels:
> lakernd.int <- lme(ltotinv ~ habfac * fishfac,
random = list (~ 1 | fishfac,
~ habfac * fishfac | lakefac)
data= use)
To test whether interactions are significant you need to compare models.
Thus the following would test whether all these included interactions
togethter lead to a statistically better model:
anova (lakernd, lakernd.int)
I hope this helps to get you on the right track.
Regards, Lorenz
-
Lorenz Gygax
Tel: +41 (0)52 368 33 84 / [EMAIL PROTECTED]
Center for proper housing of ruminants and pigs
Swiss Veterinary Office
agroscope FAT Tänikon, CH-8356 Ettenhausen / Switzerland
Fax : +41 (0)52 365 11 90 / Tel: +41 (0)52 368 31 31
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Hangup on save.image () / q ()
Dear all,
I am doing some data handling and tabualtions and then I am using glmmPQL to
fit a binomial hierarchical model (which I assign to a new object). If I do
a save.image () or a q ('yes') after this, I get the following warning
messages:
Warning messages:
1: namespaces may not be available when loading
2: names in persistent strings are currently ignored
Sometimes R seems to shut down properly, sometimes R hangs up during the
quit, in any case R cannot be restarted after this because the .RData seems
to be corrupt.
This is on R Version 1.7.0 (2003-04-16) on Windows and it happens both, if I
am using Emacs with ESS or the R-Gui.exe. .RData is about 200 KB before I do
the hierarchical model. After the hang up, the unusable file is about 1200
KB.
Any ideas what is happening or pointers where I should be looking (I tried
the R archives but did not seem to find the right key word)?
Many thanks, Lorenz
-
Lorenz Gygax
Tel: +41 52 368 33 84 / [EMAIL PROTECTED]
Center for proper housing of ruminants and pigs
Swiss Veterinary Office, FAT, CH-8356 Tänikon / Switzerland
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
