[R] The 'REP' term in AMMI{agricolae}
Dear all,
W2k, R 2.5.1
I am trying out the AMMI function in the agricolae package, to analyse the
dependence of environment for a certain cultivar. The function responds to
four basic variables:
ENV Environment
GEN Genotype
REP Replication
Y Response
My question concerns the 'REP' term. When I normally do an analysis of
variance, the replication would refer to repeated observations within the
each field trial. In the case I am analysing right now, I have five years
of observations for each 'ENV', in such a case it feels natural to use
year as the most important replication of the environment- especially as I
am interested in long term trends. This approach also seems to work
allright.
But the field trials are also structured in three main blocks, subdivided
into five 'lattice' blocks, a structure that is powerful in the analysis
of variance. (I use a random call in lme{nlme}).
Is it possible to use the block structure also in the AMMI analysis? If
that is possible, how should I code? I have tried to find out in the
documentation, but if it is stated there I do not understand it.
Thank you
/CG
--
CG Pettersson, PhD
Swedish University of Agricultural Sciences (SLU)
Dept. of Crop Production Ecology. Box 7043.
SE-750 07 Uppsala, Sweden
[EMAIL PROTECTED]
__
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
and provide commented, minimal, self-contained, reproducible code.
[R] Subsetting dataframes
Dear all!
W2k, R 2.5.1
I am working with an ongoing malting barley variety evaluation within
Sweden. The structure is 25 cultivars tested each year at four sites, in
field trials with three replicates and 'lattice' structure (the replicates
are divided into five sub blocks in a structured way). As we are normally
keeping around 15 varieties from each year to the next, and take in 10 new
for next year, we have tested totally 72 different varieties during five
years.
I store the data in a field trial database, and generate text tables with
the subset of data I want and import the frame to R. I take in all
cultivars in R and use 'subset' to select what I want to look at. Using
lme{nlme} works with no problems to get mean results over the years, but
as I now have a number of years I want to analyse the general site x
cultivar relation. I am testing AMMI{agricolae} for this and it seems to
work except for the subsetting. This is what happens:
If I do the subsetting like this:
x62_samvar - subset(x62_5, cn %in%
c(Astoria,Barke,Christina,Makof, Prestige,Publican,Quench))
A test run with AMMI seems to work in the first part:
AMMI(site, cn, rep, yield)
ANALYSIS AMMI: yield
Class level information
ENV: Hag Klb Bjt Ska
GEN: Astoria Prestige Makof Christina Publican Quench
REP: 1 2 3
Number of observations: 240
model Y: yield ~ ENV + REP%in%ENV + GEN + ENV:GEN
Analysis of Variance Table
Response: Y
DfSum Sq Mean Sq F valuePr(F)
ENV 3 120092418 40030806 90.0424 1.665e-06 ***
REP(ENV)8 3556620444578 0.5674 0.803923
GEN 5 21376142 4275228 5.4564 9.680e-05 ***
ENV:GEN15 28799807 1919987 2.4504 0.002555 **
Residuals 208 162973213783525
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Coeff var Mean yield
13.08629 6764.098
After this something goes wrong, as AMMI finds a cultivar name not
selected in the subsetting. (The plotting might go wrong anyhow, but I
haven´t got that far yet):
Error in model.frame.default(Terms, newdata, na.action = na.action, xlev =
object$xlevels) :
factor 'y' has new level(s) Arkadia
Looking at the dataframe using
edit(x62_samvar)
only shows the selected lines, but using levels() gives another answer as
levels(x62_samvar$cn)
gives back all 72 cultivar names used during the five years (starting with
Arcadia).
Where do I go wrong and how do I use subset in a proper way?
Thanks
/CG
--
CG Pettersson, PhD
Swedish University of Agricultural Sciences (SLU)
Dept. of Crop Production Ecology. Box 7043.
SE-750 07 Uppsala, Sweden
[EMAIL PROTECTED]
__
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
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Subsetting dataframes
Thanks a lot. But an ignorant R user, like me, needed the code example from Jim Holtman posted outside the list earlier today to understand that: x62_samvar$cn - x62_samvar$cn[,drop=TRUE] was the way to code. Thank you both! /CG On Thu, July 19, 2007 3:01 pm, Uwe Ligges said: CG Pettersson wrote: Dear all! W2k, R 2.5.1 I am working with an ongoing malting barley variety evaluation within Sweden. The structure is 25 cultivars tested each year at four sites, in /snip Where do I go wrong and how do I use subset in a proper way? So you have to drop the levels you are excluding. Example: x - factor(letters[1:4]) x x[1:2] x[1:2, drop=TRUE] Uwe Ligges Thanks /CG -- CG Pettersson, PhD Swedish University of Agricultural Sciences (SLU) Dept. of Crop Production Ecology. Box 7043. SE-750 07 Uppsala, Sweden [EMAIL PROTECTED] __ 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 and provide commented, minimal, self-contained, reproducible code.
[R] Why do I get a linebreak in the legend?
Dear all, W2k, R2.4.0 I want to place a legend in a regression plot, stating the adjusted R-square value. After some struggle with the coding I am nearly there, but only nearly. The best try so far is: legend(topleft, expression(paste(R[adj]^2), = 0.66)) This places the proper information in the legend, but I get a linebreak before the = 0.66 and I want the expression on one single line. All adjustments from this code I have tried so far either produces only half the expression or produces an error message. All the best and sorry for a trivial quastion /CG -- CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Crop Production Ekology. Box 7043. SE-750 07 Uppsala, Sweden [EMAIL PROTECTED] __ 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 and provide commented, minimal, self-contained, reproducible code.
[R] Plotting symbols with two positions?
Thanks a lot to Demitris for a prompt answer some minutes ago on another tread (see below). To avoid excess mails on the list, I move onto next question: I have another small plotting problem that confuses me. I want to plot results from a field trial series, using the numbers of the trials as symbols in the plot. pch = as.character(trial_no) works fine, but truncates the trial number to the first digit. As I have sixteen trials in the series I get into problems How do I squeeze in two positions as a symbol in a plot? All the best /CG On Thu, November 9, 2006 9:30 am, Dimitris Rizopoulos said: try the following: x - runif(100, -4, 4) y - 1 + 2 * x + rnorm(100, sd = 2) fit - lm(y ~ x) plot(x, y) abline(fit) legend(topleft, expression(paste(R[adj]^2, = 0.66))) ## or plot(x, y) abline(fit) legend(topleft, legend = substitute(R[adj]^2 == x, list(x = summary(fit)$adj.r.squared))) I hope it helps. Best, Dimitris Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm - Original Message - From: CG Pettersson [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Sent: Thursday, November 09, 2006 9:10 AM Subject: [R] Why do I get a linebreak in the legend? Dear all, W2k, R2.4.0 I want to place a legend in a regression plot, stating the adjusted R-square value. After some struggle with the coding I am nearly there, but only nearly. The best try so far is: legend(topleft, expression(paste(R[adj]^2), = 0.66)) This places the proper information in the legend, but I get a linebreak before the = 0.66 and I want the expression on one single line. All adjustments from this code I have tried so far either produces only half the expression or produces an error message. All the best and sorry for a trivial quastion /CG -- CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Crop Production Ekology. Box 7043. SE-750 07 Uppsala, Sweden [EMAIL PROTECTED] __ 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 and provide commented, minimal, self-contained, reproducible code. Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm -- CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Crop Production Ekology. Box 7043. SE-750 07 Uppsala, Sweden [EMAIL PROTECTED] __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] Plotting symbols with two positions?
Dear Dimitris, Thanks a lot, but I didn't really manage to apply it in my context, I first got an error message and then an ugly plot. I have also got the advice from Gabor Grothendieck to use letters instead of numbers, but I do prefer numbers as this is the normal way of referring to the individual trials. This is what I did, (and failed with): attach(DAT.nr) plot(jd.s,jd.h, type = no) Warning message: plot type 'no' will be truncated to first character in: plot.xy(xy, type, pch, lty, col, bg, cex, lwd, ...) This was not very hopeful, but as an empty frame was produced in the same time, I took a risk with this: text(jd.s, jd.h, pch=as.character(DAT.nr$t_no)) Which, to my surprise, produced a plot. The problem with it was _not_ that the numbers were truncated into the first (as could be expected) but that it seems like the line number, and not the trial number, were used. Not very nice as they number up to 576 and have 36 levels for each trial.. If I instead use the last expression, but as a plot: plot(jd.s, jd.h, pch=as.character(DAT.nr$t_no)) I am back where I begun with the correct numbers, but only the first position of it. What is happening? All the best /CG On Thu, November 9, 2006 10:19 am, Dimitris Rizopoulos said: try this: y - rnorm(16) plot(y, type = n) text(1:16, y, 1:16) Best, Dimitris - Original Message - From: CG Pettersson [EMAIL PROTECTED] To: Dimitris Rizopoulos [EMAIL PROTECTED] Cc: CG Pettersson [EMAIL PROTECTED]; r-help@stat.math.ethz.ch Sent: Thursday, November 09, 2006 10:09 AM Subject: Plotting symbols with two positions? Thanks a lot to Demitris for a prompt answer some minutes ago on another tread (see below). To avoid excess mails on the list, I move onto next question: I have another small plotting problem that confuses me. I want to plot results from a field trial series, using the numbers of the trials as symbols in the plot. pch = as.character(trial_no) works fine, but truncates the trial number to the first digit. As I have sixteen trials in the series I get into problems How do I squeeze in two positions as a symbol in a plot? All the best /CG On Thu, November 9, 2006 9:30 am, Dimitris Rizopoulos said: try the following: x - runif(100, -4, 4) y - 1 + 2 * x + rnorm(100, sd = 2) fit - lm(y ~ x) plot(x, y) abline(fit) legend(topleft, expression(paste(R[adj]^2, = 0.66))) ## or plot(x, y) abline(fit) legend(topleft, legend = substitute(R[adj]^2 == x, list(x = summary(fit)$adj.r.squared))) I hope it helps. Best, Dimitris Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm - Original Message - From: CG Pettersson [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Sent: Thursday, November 09, 2006 9:10 AM Subject: [R] Why do I get a linebreak in the legend? Dear all, W2k, R2.4.0 I want to place a legend in a regression plot, stating the adjusted R-square value. After some struggle with the coding I am nearly there, but only nearly. The best try so far is: legend(topleft, expression(paste(R[adj]^2), = 0.66)) This places the proper information in the legend, but I get a linebreak before the = 0.66 and I want the expression on one single line. All adjustments from this code I have tried so far either produces only half the expression or produces an error message. All the best and sorry for a trivial quastion /CG -- CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Crop Production Ekology. Box 7043. SE-750 07 Uppsala, Sweden [EMAIL PROTECTED] __ 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 and provide commented, minimal, self-contained, reproducible code. Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm -- CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Crop Production Ekology. Box 7043. SE-750 07 Uppsala, Sweden [EMAIL PROTECTED] Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm -- CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Crop Production Ekology. Box 7043. SE-750 07 Uppsala, Sweden [EMAIL PROTECTED] __ 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 and provide commented, minimal, self-contained, reproducible code.
[R] Model vs. Observed for a lme() regression fit using two variables
Dear all. R 2.3.1, W2k. I am working with a field trial series where, for the moment, I do regressions using more than one covariate to explain the protein levels in malting barley. To do this I use lme() and a mixed call, structured by both experiment (trial) and repetition in each experiment (block). Everything works fine, resulting in nice working linear models using two covariates. But how do I visualize this in an efficient and clear way? What I want is something like the standard output from all multivariate tools I have worked with (Observed vs. Predicted) with the least square line in the middle. It is naturally possible to plot each covariate separate, and also to use the 3d- sqatterplot in Rcmdr to plot both at the same time, but I want a plain 2d plot. Who has made a plotting method for this and where do I find it? Or am I missing something obvious here, that this plot is easy to achieve without any ready made methods? Cheers /CG -- CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences (SLU) Dept. of Crop Production Ecology. Box 7043. SE-750 07 UPPSALA, Sweden. +46 18 671428, +46 70 3306685 [EMAIL PROTECTED] __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] Model vs. Observed for a lme() regression fit using two variables
Hi Andrew, Thanks a lot, That would give me what I want. But using my own data and models resulted in this: plot(fitted(tcos31.c.cp, level=1), FCR.c$g.cp) Error in xy.coords(x, y, xlabel, ylabel, log) : 'x' and 'y' lengths differ This is quite correct, as there are some missing values in the covariate and I made the model using the 'na.action=na.omit' option. I know there is a way of using the model to fix this, but haven´t been able to get the code right during the afternoon. How do I code this and where should I have looked? Cheers /CG On Thu, September 7, 2006 12:03 pm, Andrew Robinson said: Hi CG, I think that the best pair of summary plots are 1) the fitted values without random effects against the observed response variable, and 2) fitted values with random effects against the observed response variable. The first plot gives a summary of the overall quality of the fixed effects of the model, the second gives a summary of the overall quality of the fixed effects and random effects of the model. eg fm1 - lme(distance ~ age, data = Orthodont) plot(fitted(fm1, level=0), Orthodont$distance) abline(0, 1, col=red) plot(fitted(fm1, level=1), Orthodont$distance) abline(0, 1, col=red) I hope that this helps. Andrew On Thu, Sep 07, 2006 at 11:35:40AM +0200, CG Pettersson wrote: Dear all. R 2.3.1, W2k. I am working with a field trial series where, for the moment, I do regressions using more than one covariate to explain the protein levels in malting barley. To do this I use lme() and a mixed call, structured by both experiment (trial) and repetition in each experiment (block). Everything works fine, resulting in nice working linear models using two covariates. But how do I visualize this in an efficient and clear way? What I want is something like the standard output from all multivariate tools I have worked with (Observed vs. Predicted) with the least square line in the middle. It is naturally possible to plot each covariate separate, and also to use the 3d- sqatterplot in Rcmdr to plot both at the same time, but I want a plain 2d plot. Who has made a plotting method for this and where do I find it? Or am I missing something obvious here, that this plot is easy to achieve without any ready made methods? Cheers /CG -- CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences (SLU) Dept. of Crop Production Ecology. Box 7043. SE-750 07 UPPSALA, Sweden. +46 18 671428, +46 70 3306685 [EMAIL PROTECTED] __ 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 and provide commented, minimal, self-contained, reproducible code. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 Email: [EMAIL PROTECTED] http://www.ms.unimelb.edu.au -- CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Crop Production Ekology. Box 7043. SE-750 07 Uppsala, Sweden [EMAIL PROTECTED] __ 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 and provide commented, minimal, self-contained, reproducible code.
[R] New behavior in estimable(), bug or feature?
Hello all! w2k, R-2.2.1. update.packages done today. I just started to work with a new dataset, using lme() (library nlme) and estimable() (library gmodels). I first wanted to establish the fixed effects for eight fertiliser treatments (variable treat) coded as A to H. Fitting and reducing a first model for grain yield went smoothly. When I wanted to look at the fixed effects with estimable() I got an error message claiming that I was using the wrong variable names, estimable wanted the variable names in the form usually given by fixed.effects(). I changed the names in my matrix to the requested form, but got the same error message. Then I changed to an old library with a variety trial material using a matrix I _know_ worked two months ago, I got the same type of error message. What is happening? Is there a new namecheck in estimable acting a little over-enthusiastic or what? Here are some output from R: trMat1 A 1 0 0 0 0 0 0 0 B 0 1 0 0 0 0 0 0 C 0 0 1 0 0 0 0 0 D 0 0 0 1 0 0 0 0 E 0 0 0 0 1 0 0 0 F 0 0 0 0 0 1 0 0 G 0 0 0 0 0 0 1 0 H 0 0 0 0 0 0 0 1 formula(m4.y) y.dm ~ treat - 1 fixed.effects(m4.y) treatA treatB treatC treatD treatE treatF treatG treatH 2065.267 4052.033 4571.479 4933.026 4680.980 5021.347 5063.306 5198.153 estimable(m4.y, trMat1) Error in FUN(newX[, i], ...) : Invalid parameter name(s): , , , , , , , Valid names are: treatA, treatB, treatC, treatD, treatE, treatF, treatG, treatH All the best /CG -- CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Crop Production Ekology. Box 7043. SE-750 07 Uppsala, Sweden [EMAIL PROTECTED] __ 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
[R] Failure when updating package nlme
Hello all, R2.2.1, W2k I posted a similar question last week after having problem with a group update (using update.packages()), without getting any answer. Now I have updated all packages separately, and found that the problem comes from nlme: Everything comes home, the md5 checksums are ok, but then the old version can´t be removed from my computer, leaving a nlme library with only the libs sublibrary left. As my R doesn´t start without nlme, I am stuck here. R doesn´t start with only the libs left and doesn´t start without nlme. What´s going on? It´s no big crisis as I have the nlme version from the R installation file, but I try to keep my installation as updated as possible, so it´s irritating. :) /CG -- CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Crop Production Ekology. Box 7043. SE-750 07 Uppsala, Sweden [EMAIL PROTECTED] __ 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
[R] Start problem after package update
Hello all! R2.2.1, W2k. When I used update.packages() today, R wanted to update VR, cluster and nlme packages. I did so using the Danish mirror. Everything worked as usual during update, but after closing R and trying to restart I got an error message saying that the .Rdata could not be restored. Re-installing from my downloaded .exe file cured this, but what happened and what package could have caused the problem? /CG -- CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Crop Production Ekology. Box 7043. SE-750 07 Uppsala, Sweden [EMAIL PROTECTED] __ 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
[R] using abline and a fitted 2nd degree formula
Hello all, R2.1.1, Wk2 I am doing some two-step plotting, first using plot() to illustrate the datapoints and then using abline() to place a trend line from a fitted model into the plot. Everything works well as long as the formula of the fitted model i of the type: m1 - lm(Dependent ~ Independent) then abline(m1) puts the proper straight line into the plot. But if I use: m2 - lm(Dependent ~ Independent + I(Independent^2)) abline(m2) produces a straight line, only from the first order term. Why, and what should I do about it? Cheers /CG -- CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Crop Production Ekology. Box 7043. SE-750 07 Uppsala, Sweden [EMAIL PROTECTED] __ 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
[R] Problems with abline adding regression line to a graph
Hello all,
R2.1.1, W2k
I try to make a plot of a simple regression model in this way:
with(njfA_bcd, {
+ plot(TC_OS.G31,Prot,cex = 2, col = red, xlab= TC/OS at GS32,
+ ylab=Grain crude protein (CP))
+ })
This part works well and produces the datapoints as red circles.
When I try to add a line, using a fitted linear model in a way
that works perfect with other variables in the same dateset the
following happens:
with(njfA_bcd, {
+abline(lm(predict(m1tc) ~ TC_OS.G31), lty = 1, col = red)
+ })
Error in model.frame(formula, rownames, variables, varnames, extras,
extranames, :
variable lengths differ
And this means?
There exists missing values for TC_OS.G31 in the dataset. From the
beginning
m1tc was a lm() object, which gave the same Error message. To try to fix the
problem I changed to lme() and used na.action=na.omit explicitely, but this
didn´t help.
Here is the summary of m1tc:
summary(m1tc)
Linear mixed-effects model fit by REML
Data: njfA_bcd
AIC BIClogLik
209.4914 219.0692 -100.7457
Random effects:
Formula: ~1 | Trial
(Intercept) Residual
StdDev:1.242184 0.6520464
Fixed effects: Prot ~ TC_OS.G31
Value Std.Error DF t-value p-value
(Intercept) 14.86209 0.957630 68 15.519662 0
TC_OS.G31 -24.22286 4.792801 68 -5.054008 0
Correlation:
(Intr)
TC_OS.G31 -0.935
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.68329774 -0.73751040 -0.05600477 0.68301243 2.21693174
Number of Observations: 83
Number of Groups: 14
What is happening and what shall I do about it?
Cheers
/CG
--
CG Pettersson, MSci, PhD Stud.
Swedish University of Agricultural Sciences (SLU)
Dept. of Crop Production Ecology. Box 7043.
SE-750 07 UPPSALA, Sweden.
+46 18 671428, +46 70 3306685
[EMAIL PROTECTED]
__
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
[R] Analyses of covariation with lme() or lm()
Hello all! I have a problem that calls for a better understanding, than mine, of how lme() uses the random part of the call. The dataset consists of eleven field trials (Trial) with three replicates (Block) and four fertiliser treatments (Treat). Analysing for example yield with lme() is easy: m1 - lme(Yield ~ Treat, data=data, random =~1| Trial/Block) giving estimates of Treat effects with good significances. If I compare m1 with the model without any random structure: m2 - lm(Yield ~ Treat, data=data), m1 is, naturally, much better than m2. So far so good. Now I have one (1) measure from each Trial, of soil factors weather and such, that I want to evaluate. Remember: only one value of the covariate for each Trial. The suggestion I have got from my local guru is to base this in m1 like: m3 - lme(Yield ~ Treat + Cov1 + Treat:Cov1, data=data, random =~1| Trial/Block) thus giving a model where the major random factor (Trial) is represented both as a (1) measure of Cov1 in the fixed part and by itself in the random part. Trying the simpler call: m4 - lm(Yield ~ Treat + Cov1 + Treat:Cov1, data=data) gives back basically the same fixed effects as m3, but with better significances for Cov1. Tested with anova(m3,m4) naturally gives the answer that m3 is better than m4. Ok, what about dealing with Trial in the fixed call? : m5 - lm(Yield ~ Trial + Treat + Cov1 + Treat:Cov1, data=data) lm() swallows this, but silently moves out Cov1 from the analysis, an action that feels very logical to me. My guru says that using the random call secures you from overestimating the p-values of the covariate. I fear that the risk is as big that you underestimate them with the same action. Working on a paper, I naturally want to be able to do some sort of discussion on the impact of covariates... ;-) What is the wise solution? Or, if this is trying to make other people do my homework, could anyone tell me where the homework is? (I´ve got both Pinhiero Bates and MASS as well as some others in the bookshelf.) Cheers /CG -- CG Pettersson MSci. PhD.Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Crop Production Ecology (VPE). http://www.slu.se/ [EMAIL PROTECTED] __ 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
[R] No ~ in JGR
R2.1.0 JGR 1.2 W2k Hello all! I´ve just installed JGR on my both R-equipped computers and am very pleased with the look and functionality. Except in one, very important, way. I can´t figure out how to get the ~ sign from the keyboard to the console. Copying it from old code works fine. Using the traditional GUI works as usual. I have a Swedish keyboard layout, where ~ shares key with ¨ and ^, all reguiring a space after the hit to produce the sign on the screen. The other signs from the key works, but AltGr + ~ followed by space just gives a blank. What do I do wrong? /CG -- CG Pettersson MSci. PhD.Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Ecology and Crop production sciences (EVP). http://www.slu.se/ [EMAIL PROTECTED] __ 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
[R] Correct effect plots from lme() objects
Hello all! R2.0.1, W2k I posted this question to the list last Sunday without getting any replies on the list. I got two off the list though, suggesting me to plot manually as a second step, from estimable() or intervals() objects respectively. As this was not really what I wished for, I take the risk to upset somebody with a trivial question, and re-post it (just a little edited). So, here it is again: I use lme() from nlme to make mean estimates from series of field experiments where not all treatments (like variety) are present in all experiments. (An old problem in variety evaluation). A typical call is: yield.lme - lme(Yield ~ Variety, data = data, na.action = na.omit, random = ~ 1 | Experiment/Block) This works well, even when observations are lacking. I have checked against the accepted method for doing this in Sweden, which is PROC MIXED in SAS, and the fitted fixed effects are more or less identical. I use estimable() from gmodels (gregmisc package) to extract estimates, standard errors and such. I use matrices with the variety names as row names, it works smooth. What I am unable to, as yet, is to make nice plots of the estimates for a given set of varieties. To use only the fixed call directly on the dataset works for many plooting functions, but produces the wrong graphs, as the structure is not used. The structure (the random call) has to be used, as there are NA:s in the dataset. Pinheiro Bates have a lot of graphics on lme objects, but they try to illustrate more sophisticated relations than my need. I´ve looked through gplots and the graphic parts of nlme without any hits. Probably, my difficulties are just due to my own lack of skill. Some standard plotting facility plotting directly from the lme object ought to work, but I don´t understand how. How should I do this? Cheers /CG -- CG Pettersson MSci. PhD.Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Ecology and Crop production sciences (EVP). http://www.slu.se/ [EMAIL PROTECTED] __ 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
[R] Plotting from lme() objects
Hello all! R2.0.1, W2k I use lme() from nlme to make mean estimates from series of field experiments where not all treatments (like variety) are present in all experiments. A typical call is: yield.lme - lme(Yield ~ Variety, data = data, na.action = na.omit, random = ~ 1 | Experiment/Block) This works well, even when observations are lacking. I have checked against the accepted method for doing this in Sweden, which is PROC MIXED in SAS, and the fitted fixed effects are more or less identical. I use estimable() from gmodels (gregmisc package) to extract estimates, standard errors and such. I use matrices with the variety names as row names, it works smooth. What I am unable to, as yet, is to make nice plots of the estimates for a given set of varieties. To use only the fixed call directly on the dataset works, but produces the wrong graph, as the structure is not used. Pinheiro Bates have a lot of graphics on lme objects, but they try to illustrate more sophisticated relations than my need. I´ve looked through gplots and the graphic parts of nlme without any hits. Probably, my difficulties are just due to my own lack of skill. Some standard plotting facility plotting directly from the lme object ought to work, but I don´t understand how. Another possiblility is to plot from the estimable() output. How should I do this? All suggestions are welcome. Cheers /CG CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences Dep. of Ecology and Crop Production. Box 7043 SE-750 07 Uppsala __ 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
[R] Postgraduate distance learning courses for MSc in Statistics
30 Januari, Tilo Blenk wrote: Does anybody has experience with postgraduate distance learning courses resulting in a MSc in Statistics? I am thinking about such a course to ... snip ... Sorry for posting something only distantly related to R but I hope there some people reading this mailing list who could help me. Moreover, R is one of the reasons why I got interested in statistics. Tilo Dear Tilo As it was surprisingly silent on this question, I can offer a tip for you, nearer than both England or USA: At KVL in Copenhagen, they give several courses with well developed distance methodology. They use both SAS and R, not totally parallell though.. Check out the following link: http://statmaster.sdu.dk/ Cheers /CG CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences Dep. of Ecology and Crop Production. Box 7043 SE-750 07 Uppsala Sweden __ 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] Conflicts using Rcmdr, nlme and lme4
Thanks a lot, both to Prof Bates and Prof Fox. I see at least some bits of the light now. One original question remaines though: Why do all these (22 of them) packages autoload in the first place? I don´t want them to. It must have something to do with the original loading from Rcmdr, as (for example) nlme does *not* autoload by itself no matter how many lme objects there are in a reopened workspace, this is good as far as I am concerned. It is especially strange as Rcmdr itself does not outoload. How do I supress this behavior? Thanks /CG 28 Jan, John Fox wrote: --- Dear CG, An addendum to my previous response: If you do decide to recompile the Rcmdr package to eliminate the dependency on nlme, as I suggested, you'll also have to edit the .onLoad() function in the package (which is in the file startup.R in the source package) to remove nlme from the vector of required packages. Regards, John John Fox Department of Sociology McMaster University Hamilton, Ontario Canada L8S 4M4 905-525-9140x23604 http://socserv.mcmaster.ca/jfox -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of CG Pettersson Sent: Friday, January 28, 2005 5:15 AM To: r-help@stat.math.ethz.ch Subject: [R] Conflicts using Rcmdr, nlme and lme4 Hello all! R2.0.1, W2k. All packages updated. I´m heavily dependant on using mixed models. Up til´now I have used lme() from nlme as I have been told to. Together with estimable() from gmodels it works smooth. I also often run Rcmdr, mostly for quick graphics. ...snip... CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences Dep. of Ecology and Crop Production. Box 7043 SE-750 07 Uppsala __ 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
[R] Conflicts using Rcmdr, nlme and lme4
Hello all! R2.0.1, W2k. All packages updated. I´m heavily dependant on using mixed models. Up til´now I have used lme() from nlme as I have been told to. Together with estimable() from gmodels it works smooth. I also often run Rcmdr, mostly for quick graphics. After using Rcmdr, on reopening the R workspace all help libraries for Rcmdr (22 !) loads, among them nlme, but not Rcmdr itself. Why? Now I saw on the list yesterday, that lmer() from lme4 offers better coding possibilities for crossed random factors, it feels natural to switch from lme() to lmer() if estimable() still works on the objects. So I installed lme4, but got a conflict with the (auto)loaded nlme: library(lme4) Loading required package: Matrix Loading required package: latticeExtra Error in fun(...) : Package lme4 conflicts with package nlme. To attach lme4 you must restart R without package nlme. Error: .onLoad failed in loadNamespace for 'lme4' Error in library(lme4) : package/namespace load failed for 'lme4' Several questions here: Why do all these packages autoload, and could I avoid this in some way? Why a conflict? Isn´t it possible to mask the conflicting parts of nlme when loading lme4? Why are the mixed models tools from the same author(s?) splitted between nlme and lme4? Are there any technical reasons not to include lmer() in nlme? Cheers /CG CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences Dep. of Ecology and Crop Production. Box 7043 SE-750 07 Uppsala __ 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
[R] Fine tuning scatterplot() (car package)
Hello all! System: R2.0.1, W2k. All packages apdated with update.packages(). I use the default graphic device for my system (no active choice). I´m trying to fine-tune scatterplot() graphs (package car), for publication. I want to control plotting symbols and colours to match plots from other, less flexible, systems. I also need to make all text and symbols bigger, to enable printing of the plots in smaller scale than the original outputs. My strategy has been to start in Rcmdr, using the command line produced there and rerun variants of the scatterplot() call directly in the console. After consulting MASS and CAR, the colour and characters were easy to fix by overrunning the defaults with col=c() and pch=c(). But the size of letters and figures beats me so far. I´ve tried cex and csi in the same positiones as col and pch on the line. cex doesn´t do anything as far as I can see, csi produces warning messages (which is some sort of effect!) but does nothing on the size. Trying to tune with these commands inside xlab and ylab calls (based on the defaults from the documentation) results in unspecified syntax error messages. My code in the console is this for the moment: scatterplot(Kernel.protein~TC.OS | Year, reg.line=FALSE, smooth=FALSE, labels=FALSE, boxplots=FALSE, span=0.5, by.groups=TRUE, col=c(,red,blue,black), pch=c(15,4,5), data=ecpa.f) Which produces nice plots, but with too small symbols, letters an numbers. Thanks! /CG CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences Dep. of Ecology and Crop Production. Box 7043 SE-750 07 Uppsala __ 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
[R] MIME decoding in Mozilla Thunderbird
Hello all! I recently switched mail program to Mozilla Thunderbird, running on W2k. Everything works fine, except the MIME-digests from this list. The decoding doesn´t work properly. I had contact with Martin Maechler some time ago, and he suggested a try on this list for ideas on how to do the decoding in Windows, even if it´s not a proper R-question. Outlook do the proper job on the MIME, but using that is to go a bit far! Or? /CG -- CG Pettersson MSci. PhD.Stud. Swedish University of Agricultural Sciences (SLU) Dep. of Ecology and Crop production sciences (EVP). http://www.slu.se/ [EMAIL PROTECTED] __ [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
[R] Subsetting with more than one criteria
Hello all! I have been using R for two years, but still have trivial problems. For the moment, I run R1.9.1 on a W2k computer. I work with a dataset from a variety evaluation project. Every year 25 varieties have been tested. The tested cultivars have been the same each year in all trials, but some have changed over the three year period. Now I want to do some calculations and plotting on subsets. The structure looks like this: YearADB Block Vcode Variety Yield Protein 2002SW0024 1 20226 Denise 584312.8 2002SW0024 1 9865Astoria 672911.4 2002SW0024 1 9622Barke 612112 2002SW0024 1 9604Cecilia 557912.7 2002SW0024 1 20223 Granta 559111.6 2002SW0024 1 20222 Class 559111.7 2002SW0024 1 9922Wiking 574412.5 2002SW0024 1 20103 Vortex 586310.6 . . And so on both down and sideways. Three years and four trials with three replicats each gives 900 lines. How do I use several criteria to subset this? Ast - subset(data, Variety == Astoria) works fine giving back the 36 lines where Astoria appears. But how do I pick two (or more) varieties? AND picking one (or two) years? Everything I have tried results in rubbish, like: AstBark - subset(data, Variety == c(Astoria,Barke)) which seems to work but gives back 19 and 18 lines of the varieties respectively. There exists 36 of each Thanks /CG CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences Dep. of Ecology and Crop Production. Box 7043 SE-750 07 Uppsala __ [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
[R] Data import, going from 8 to 550 columns
Hello all! I need to import a NIR dataset into R. It should be quite trivial, but I can´t make it work. (No problems with the text in the beginning, as # is recognised by read.table as the comment sign.) The thing I can´t get around is the CR that ends every line after column eight as the line in R should be 550 columns wide (including the JF-number). Every new line in R should begin with the JF2455 and so on. Naturally it is possible to re-shape the tables in Excel before import, but it is quite tedious and doesn´t feel right...! How do I make read.table to just go on reading on the next line when it comes to CR, and how do I make it use the double CR followed by a blank to begin the next line? The data-file(s) looks like this: #ID=Samples from soil scanning #SAMPLE_NUMBERS_PRESENT=Y #NX_VARIABLES=550 #NY_VARIABLES=0 #FIRST_WAVELENGTH=1300.00 #LAST_WAVELENGTH=2398.00 #WAVELENGTH_INCREMENT=2.00 JF2455 0.4367495 0.4365539 0.4363573 0.4361560 0.4359702 0.4357788 0.4355963 0.4354126 0.4352311 0.4350726 0.4349101 0.4347557 0.4346097 0.4344587 0.4343193 0.4341759 0.4340320 0.4338984 0.4337671 0.4336369 0.4335097 0.4333864 the original table is 8 columns wide, ended with a CR sixty four lines removed here 0.5015950 0.5020472 0.5026294 0.5033303 0.5041344 0.5049909 0.5059010 0.5067372 0.5075415 0.5082389 0.5089509 0.5095288 0.5101137 0.5106306 0.5111954 0.5116805 JF2456 0.3604568 0.3600681 0.3596676 0.3592694 0.3588919 0.3585098 0.3581379 0.3577725 0.3573992 0.3570563 0.3566975 0.3563588 0.3560365 0.3556931 0.3553730 0.3550543 0.3547286 0.3544230 0.3541073 0.3537982 0.3535004 0.3531921 0.3529077 0.3526271 0.3523493 0.3520919 0.3518271 0.3515673 0.3513192 0.3510693 0.3508208 0.3505693 ... and so on Thanks! /CG CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences Dep. of Ecology and Crop Production. Box 7043 SE-750 07 Uppsala __ [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] Is there a /ddfm=satterth for R?
Hello all! When you are working with a little more complicated models in SAS PROC MIXED, you often use the /ddfm=satterth call to make sure the df decomposition is done the best way possible. Running the same models in lme, without any special calls, results in warning messages about the df handling. Is anybody out there working with something like the /ddfm=satterth? It would be handy, or are there any reasons not to use it? /CG CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences Dep. of Ecology and Crop Production. Box 7043 SE-750 07 Uppsala __ [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 - problems with model
Thanks a lot for the answer! Now, I only have the last one left - How do I get round it? I knew about the missing cells in the design, but didn´t know how lme would react on them. In this case, I can remove the water:temp term, but how can I be sure that this is the right thing to do? Is the lm run without the random term enough for removing water:temp from the lme model, or do I have to do a PROC MIXED run with the random term to make that decision in a case like this? Is it possible (for me) to understand why MIXED accepts the design but not lme? They ought to get the same sort of problems, or have I missed something? /CG --- CG Pettersson [EMAIL PROTECTED] writes: Hello all! I´m working with some training datasets in a SAS-based course, trying to do the same things in lme that I do in PROC MIXED. Why don´t lme do an analysis on this dataset when I use the model water*temp? The trouble comes from the water:temp term, as it works with water+temp. The data are, indeed, assymetric but lm accepts the water:temp term giving results in the F-test near what PROC MIXED produces. MIXED accepts the model. The water:temp term can be removed from the model according to the F-test in SAS (and to the lm model without any random term). Doing so in both MIXED and lme gives reasonably similar results for both systems. What do the error message mean, and how can I get around this? Because of missing cells in the design xtabs(~water + temp, milk) temp water 100 110 120 140 1 3 3 3 0 2 3 0 3 3 3 3 3 0 3 the model matrix for the fixed effects is rank deficient. In lm the rank deficiency is detected and appropriate adjustments made coef(summary(lm(maill6 ~ water * temp, milk))) Estimate Std. Errort value Pr(|t|) (Intercept) 2.1767 0.1142339 19.0544730 2.218661e-13 water2 0.2833 0.1615511 1.7538308 9.647013e-02 water3 0.0533 0.1615511 0.3301329 7.451108e-01 temp110 0.1400 0.1615511 0.8665987 3.975669e-01 temp120 0.3133 0.1615511 1.9395305 6.827304e-02 temp140 0.2333 0.1615511 1.4443312 1.658280e-01 water3:temp110 -0.1867 0.2284678 -0.8170371 4.245898e-01 water2:temp120 0.0967 0.2284678 0.4231085 6.772282e-01 water2:temp140 0.2167 0.2284678 0.9483467 3.555125e-01 Notice that you would expect 6 degrees of freedom for the interaction term but only three coefficients are estimated. In lme it is much more difficult to compensate for such rank deficiencies because they could be systematic, like this, or they could be due to relative precision parameters approaching zero during the iterations. Because of this we just report the error (although admittedly we could be a bit more explicit about the nature of the problem - we are reporting the symptom that we detect, not the probable cause). The dataset: milk water temp rep maill4 maill6 maill8 taste4 taste6 taste8 1 1 100 1 2.90 2.13 2.39 10.1 10.09.6 2 1 100 2 2.19 2.20 2.27 11.09.3 11.0 3 1 100 3 2.13 2.20 2.41 10.17.09.6 4 1 110 1 2.13 2.34 2.41 11.0 10.59.8 5 1 110 2 2.32 2.27 2.25 11.0 11.3 11.2 6 1 110 3 2.13 2.34 2.429.4 10.79.0 7 1 120 1 2.00 2.49 2.71 11.1 11.2 11.4 8 1 120 2 2.41 2.49 2.46 11.6 11.79.6 9 1 120 3 2.22 2.49 2.73 10.7 10.3 10.2 10 2 100 1 2.13 2.41 2.49 11.1 10.8 11.2 11 2 100 2 2.49 2.34 2.53 11.1 11.29.2 12 2 100 3 2.80 2.63 3.338.39.77.8 13 2 120 1 2.38 2.85 2.06 11.9 11.2 11.2 14 2 120 2 2.61 2.70 2.70 11.7 10.8 11.0 15 2 120 3 2.77 3.06 3.25 10.99.09.4 16 2 140 1 2.56 2.84 3.10 10.7 11.29.8 17 2 140 2 2.63 2.61 2.81 10.8 11.0 11.6 18 2 140 3 2.99 3.28 3.759.29.69.6 19 3 100 1 2.60 2.24 2.32 10.88.4 10.8 20 3 100 2 2.06 2.11 2.20 11.0 11.2 11.8 21 3 100 3 1.98 2.34 2.80 10.3 10.2 10.6 22 3 110 1 1.91 2.06 2.29 11.0 11.49.4 23 3 110 2 1.98 1.98 2.15 10.0 11.8 10.6 24 3 110 3 1.98 2.51 2.819.39.2 10.2 25 3 140 1 2.27 2.42 2.72 10.8 11.6 12.0 26 3 140 2 2.27 2.20 2.41 11.2 11.0 11.4 27 3 140 3 2.20 2.77 3.06 10.5 10.2 10.0 The failing model: lme(maill6 ~ water * temp , random= ~1|rep, data = milk) Error in MEEM(object, conLin, control$niterEM) : Singularity in backsolve at level 0, block 1 The smaller (working) model: lme(maill6
[R] Problems loading dataset in Rcmdr
Hello all! I´ve been using Rcmdr for some time, as a quick way of producing graphics and basic statistics. I run R1.8.1, OS W2000. Two days ago the dataset loader stopped working. Normally, the button No active dataset is clickable to give you the opportunity to choose dataset to load in the Rcmdr context. Clicking on the button now produces this: Rcmdr Version 0.9-3 Error in parse(file, n, text, prompt) : parse error What is happening and what could I do to get out of it? (I have removed Rcmdr and downloaded it again from CRAN, as well as uninstalled R, reinstalling it from the copy of my download. Doesn´t help so far. Please help! /CG CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences Dep. of Ecology and Crop Production. Box 7043 SE-750 07 Uppsala __ [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] Black and white output from scatterplot (CAR package)
Hello! I want to produce outputs from scatterplot (CAR-package), that are printeble in black and white. The default output is very nice, but in color. I have found out how to make a greyscale instead. What I want is a good black and white output where the lines indicating regression for the groups, are possible to distinguish from eachother by traditional dotted lines (style) in black ink. Is this possible or do I have to do it by hand? /CG CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences Dep. of Ecology and Crop Production. Box 7043 SE-750 07 Uppsala __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Evaluating outer observations in an lme object.
Hello everybody! I´m working with a dataset from twelve fertilizer trials, where the technical fertilizer product and application method, but not the intensity of fertilization, is varied. (I´m using R1.7.1 and W2000.) The call: ejna1t4b.lme - lme( Yield ~ TrCode, data = ejna1t4, + random = ~ 1 | Trial/Block) works as far as I can understand well, the Block structure of the trials is used efficiently and everything looks nice according to plots of the object. Now I want to evaluate the influence of observations from the different experimental places (for example soil analyses or rainfall) - Could I do that without skipping the Trial/Block structure, or do I have to start from scratch again? The observed values will naturally only have one level for each Trial, so the term Trial/Block will host the effects of all observed (and unobserved) phenomena in each trial. Now I want to know where the effects come from. I´ve been looking for a text on this, both in MASS and Pinheiro Bates, without finding any. Any hints of where to look? Thanks /CG CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences Dep. of Ecology and Crop Production. Box 7043 SE-750 07 Uppsala __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Evaluating outer observations in an lme object.
Oh, the question might have been more precisely formulated I guess... update surely helps a lot in the practical work at the computer. But that is not my problem. The problem is where and how to introduce the variables in the command. I´ve tried things like: random = ~ P.AL | Trial/Block as the random call, looking for some sort of direct effect from the soil analysis P.AL. That results in a worse model than without the term, possibly becouse the effect of the soil already is, implicitlely, inside the term Trial. What I am looking for is a way of evaluating several factors that have just one observation for each Trial, still having the the nice Trial/Block structure present in the model. Or am I just stupid? /CG --- Might update help? spencer graves CG Pettersson wrote: Hello everybody! I´m working with a dataset from twelve fertilizer trials, where the technical fertilizer product and application method, but not the intensity of fertilization, is varied. (I´m using R1.7.1 and W2000.) The call: ejna1t4b.lme - lme( Yield ~ TrCode, data = ejna1t4, + random = ~ 1 | Trial/Block) works as far as I can understand well, the Block structure of the trials is used efficiently and everything looks nice according to plots of the object. Now I want to evaluate the influence of observations from the different experimental places (for example soil analyses or rainfall) - Could I do that without skipping the Trial/Block structure, or do I have to start from scratch again? The observed values will naturally only have one level for each Trial, so the term Trial/Block will host the effects of all observed (and unobserved) phenomena in each trial. Now I want to know where the effects come from. I´ve been looking for a text on this, both in MASS and Pinheiro Bates, without finding any. Any hints of where to look? Thanks /CG CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences Dep. of Ecology and Crop Production. Box 7043 SE-750 07 Uppsala __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help CG Pettersson, MSci, PhD Stud. Swedish University of Agricultural Sciences Dep. of Ecology and Crop Production. Box 7043 SE-750 07 Uppsala __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] na.action in model.tables and TukeyHSD
In 27/2, I got the following answer from Prof. Ripley: (The question is at the bottom) This ia already fixed in R-devel. The answer is the same: don't use na.omit implicitly: use it explicitly. I feel rather stupid for the moment, as I don´t understand an answer that looks very simple. What´s the code to do the trick using na.omit explicitly? (Preferably starting with my code in the question) I can´t get it to work, so my tries are not worth printing here... Thanks /CG On Wed, 26 Feb 2003, CG Pettersson wrote: Hello everybody! I use R 1.6.2 in Windows, and have a problem controlling the na.action. In a dataset with twelve trials, one of the trials lack any readings of the variable STS.SH (standing power at harvest) Fitting an aov() object with the call: led1t7sts.aov - aov(STS.SH ~ Trial/Block + Treatment + Treatment:Trial, data = led1t7, na.action=na.exclude) seems to work as it produces an object with 10 df for the factor Trial. But when I use model.tables or TukeyHSD on the object I get this: model.tables(led1t7sts.aov, means) Error in replications(paste(~, paste(names(tables), collapse = +)), : na.action must be a function I have tried to use na.action=na.exclude inside the model.tables call as well, without any bettering. I can naturally cope with the problem by taking the whole trial away from the dataset, but it doesn´t feel very sophisticated...;-) (Prof. Ripley answered a similar question from me two weeks ago. The answer was good but didn´t work as the reason of the error was the same as this time: a whole trial with only na:s in it). Thanks /CG CG Pettersson [EMAIL PROTECTED] CG Pettersson [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list http://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Updating R
Hi! The installing and updating functions for packages in R (install.package() , update.package()) really impresses me in behavior. For the moment I run R 1.5.1 and got warning messages after installing the multcomp package, as it is built under R 1.6.1. Is there any similar smooth ways of updating the R-version as there is for packages? I´ve looked in all places I can think of but can´t find any hints. Thanks /CG CG Pettersson [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list http://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Random or fixed factors
Hi! I work right now with my first real job in R. It´s a series of 12 fertilizer trials with three repetations. The treatments can be broken down to 2 x 2 factors (Product and Application). Experimenting both with lm() and aov() makes me a little confused as I don´t get any differences in the sum of squares. With aov() I have used a command of this type: aov(Yield ~ Product * Application + Error(Trial/Block)) With lm() I have used something like: lm(Yield ~ Trial/Block + Product * Application) The placing of the term Trial/Block in the lm() statement doesn´t seem to matter. How do I tell R that the Trial term should be random in the lm() statement? It´s quite obvious that the Error term in aov() is random, but in the lm()? Is that implied of the fact that it has a nested term inside? If so, and I don´t have plotwise data, how do I mark Trial as Random? /CG CG Pettersson [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list http://www.stat.math.ethz.ch/mailman/listinfo/r-help
