Re: [R] Font quality in base graphics

[Mark Difford]

Hi willemf,

Glad to hear that it helped.  Years ago (late-90s) I Linuxed, but have since
been forced into the Windows environment (where, however, I have the great
pleasure of being able to use MiKTeX and LyX, i.e. TeX/LaTeX). I therefore
can't help you further, except to say that I have never had a problem
controlling font sizes &c to my admittedly very demanding --- some people
say excessively demanding --- standards (and that's on Windows!). And I have
never had a problem with labels &c not being where they should be, or of the
size I want them to be, when I have built the graphic from "scratch." And
only very rarely have I encountered such problems when using canned graph
types.

In brief, what I am saying is that the problem almost certainly lies with
the way fonts &c are set up on your Linux box. Were this not the case, then
I can assure you that there would have many and varied sharply worded
statements on this list relating to the poor quality of R's graphs. And
there would have been just as many pointed, well-written rebukes, pointing
that  Yet there aren't. If you search the archives you will find that a
good many users migrated to R from other systems because of R's excellent
graphical subsystems. Look at the graphics in any of the many books now
published on using R, or that use R to elucidate problems Set your mind
at rest: look at your system setup, and the tools outside R that you are
using.

Hope it all works out. OpenOffice is now a very good suite of programs, but
if you want true quality of output then you really should be TeXing. Check
it out.

Bye, Mark.


willemf wrote:
> 
> Mark, your suggestion results in about 75% control over the plot. This is
> the best that I have managed to get it at, so thank you very much. In
> Linux you create a X11() device for screen output. Specifying identical
> device characteristics results in a fair degree of similarity between
> screen version and EPS version. However in this case, for instance, some
> labels along the X axis are omitted in the screen version and
> (thankfullly!) included in the Postscript version. Also, the relative
> sizes of caption font size and label font size are not identical in the
> two versions. I have learnt a few things in this exercise, so thanks you
> very much for the advice.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Font-quality-in-base-graphics-tp18465608p18483719.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Font quality in base graphics

[Mark Difford]

Hi willemf,

And perhaps I should have added (in case you are moving across systems) that
you should take a look at 

?embedFonts
http://cran.r-project.org/doc/Rnews/Rnews_2006-2.pdf

Ghostscript should be installed on your box, so there shouldn't be a
problem. Without additional information, it's not possible to help further.
Of course, you could send me the data and a script showing how you want it
plotted, and I would send you a PDF in return, showing you what R can do ;).

HTH, Mark.



Mark Difford wrote:
> 
> Hi willemf,
> 
> Glad to hear that it helped.  Years ago (late-90s) I Linuxed, but have
> since been forced into the Windows environment (where, however, I have the
> great pleasure of being able to use MiKTeX and LyX, i.e. TeX/LaTeX). I
> therefore can't help you further, except to say that I have never had a
> problem controlling font sizes &c to my admittedly very demanding --- some
> people say excessively demanding --- standards (and that's on Windows!).
> And I have never had a problem with labels &c not being where they should
> be, or of the size I want them to be, when I have built the graphic from
> "scratch." And only very rarely have I encountered such problems when
> using canned graph types.
> 
> In brief, what I am saying is that the problem almost certainly lies with
> the way fonts &c are set up on your Linux box. Were this not the case,
> then I can assure you that there would have many and varied sharply worded
> statements on this list relating to the poor quality of R's graphs. And
> there would have been just as many pointed, well-written rebukes, pointing
> that  Yet there aren't. If you search the archives you will find that
> a good many users migrated to R from other systems because of R's
> excellent graphical subsystems. Look at the graphics in any of the many
> books now published on using R, or that use R to elucidate problems
> Set your mind at rest: look at your system setup, and the tools outside R
> that you are using.
> 
> Hope it all works out. OpenOffice is now a very good suite of programs,
> but if you want true quality of output then you really should be TeXing.
> Check it out.
> 
> Bye, Mark.
> 
> 
> willemf wrote:
>> 
>> Mark, your suggestion results in about 75% control over the plot. This is
>> the best that I have managed to get it at, so thank you very much. In
>> Linux you create a X11() device for screen output. Specifying identical
>> device characteristics results in a fair degree of similarity between
>> screen version and EPS version. However in this case, for instance, some
>> labels along the X axis are omitted in the screen version and
>> (thankfullly!) included in the Postscript version. Also, the relative
>> sizes of caption font size and label font size are not identical in the
>> two versions. I have learnt a few things in this exercise, so thanks you
>> very much for the advice.
>> 
>> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Font-quality-in-base-graphics-tp18465608p18483965.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] help with bivariate density plot question

[Mark Difford]

Hi Daniela,

Spencer (? Graves) is not at home. Seriously, this is a list that many
people read and use. If you wish to elicit a response, then you would be
wise to give a better statement of what your difficulty is.

The function you enquire about is well documented with an example, see

##
library(GenKern)   ## load the package
?KernSur  ## get help on the function

You don't need to do anything special to get adaptive bandwidths, it's all
done for you (by the authors of the package). Just replace the x and y
values in the example with your values, and perhaps deal with any NAs in
your data set.

Should one moralize?: Well, it is generally true that you want help from
others...

HTH, Mark.


Daniela Carollo wrote:
> 
> Hi Spencer,
> 
> I have seen your name on the web site, and perhaps you can
> help me with my R problem.
> 
> I'm trying to use KernSur to put in evidence a substructure in a
> bidimensional plot. My problem is that, in order to get the density
> in the low density areas (in which the substructure is located) I should
> use different bandwidths. How I can do that?
> 
> Also, I think that the best choice for my case is to use the function
> "akerdmul" which perform the multivariate adaptive kernel density
> distribution. Are you familiar with this function?
> 
> Any help would be really appreciated.
> 
> Thank you very much!
> 
> Regards,
> 
> Daniela
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/help-with-bivariate-density-plot-question-tp18495958p18504638.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] calculate differences - strange outcome

[Mark Difford]

Hi Andreas,

It's because you are dealing with binary or floating point calculations, not
just a few apples and oranges, or an abacus (which, by the way, is an
excellent calculating device, and still widely used in some [sophisticated]
parts of the world).

http://en.wikipedia.org/wiki/Floating_point

HTH, Marks.


Kunzler, Andreas wrote:
> 
> Dear List,
> 
> I ran into some trouble by calculating differences. For me it is
> important that differences are either 0 or not.
> 
> So I don't understand the outcome of this calculation
> 
> 865.56-(782.86+0+63.85+18.85+0)
> [1] -1.136868e-13
> 
> I run R version 2.71 on WinXP
> 
> I could solve my problem by using
> round()
> but I would like to know the reason.
> Maybe someone can help me?
> 
> Thanx
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/calculate-differences---strange-outcome-tp18505010p18505498.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Newbie's question about lm

[Mark Difford]

Hi Ptit,

>> I would like to fit data with the following formula :
>> y=V*(1+alpha*(x-25))
>> where y and x are my data, V is a constant and alpha is the slope I'm
>> looking for.

Priorities first: lm() or ordinary least-square regression is a basically a
method for finding the best-fitting straight line through a set of points
(assuming that x is measured without error). So lm() determines the slope;
it's usually what you are trying to estimate.

Perhaps you are after abline:
?abline

You can feed this coefficients (i.e. your intercept and slope), and then
plot over an existing scatterplot.

##
plot(x,y)
abline(a=intercept, b=slope)  ## or: abline(coef=c(intercept, slope))

HTH, Mark.


Ptit_Bleu wrote:
> 
> Hello,
> 
> I would like to fit data with the following formula :
> y=V*(1+alpha*(x-25))
> where y and x are my data, V is a constant and alpha is the slope I'm
> looking for.
> 
> How to translate this into R-language ?
> At the moment, I only know : lm(y ~ x)
> 
> Sorry for such a basic question. I thought I could find the solution in a
> post but I have to confess that, up to know, I'm not able to understand
> the posts I read concerning lm (the level of the questions are too high
> for me and my english quite poor as well).
> 
> Have a nice evening,
> Ptit Bleu. 
> 

-- 
View this message in context: 
http://www.nabble.com/Newbie%27s-question-about-lm-tp18511262p18517239.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] R version 2.7.1 (warnings)

[Mark Difford]

Hi Jaap,

With all those packages loading it could take some time, unless it's a known
problem (?).  Why don't you do a vanilla start (add switch --vanilla to
startup) and do some simple core-related stuff. Then add packages
one-by-one...

Or: search through the source code of the packages for the error message...

HTH, Mark.


Van Wyk, Jaap wrote:
> 
> Hi
> I have installed the new version 2.7.1, and get the following messages. Is
> there some way I can "fix it" ?
> I use SciViews as well, so the problem may be with one of the libraries
> that is loaded at start-up:
>  
> Start-up:
> Loading required package: datasets
> Loading required package: utils
> Loading required package: grDevices
> Loading required package: graphics
> Loading required package: stats
> Loading required package: tcltk
> Loading Tcl/Tk interface ... done
> Loading required package: R2HTML
> Loading required package: svMisc
> Loading required package: svIO
> Loading required package: svViews
> During startup - Warning message:
> 'Sys.putenv' is deprecated.
> Use 'Sys.setenv' instead.
> See help("Deprecated") 
> 
> What should I do to get rid of the warning message?
>  
> Next, after loading some data in the usual way and making a simple plot, I
> get the following:
> Warning message:
> 'length(rescale)' differs between new and previous
>  ==> NOT changing 'rescale' 
> (This also happens at some other times, doing other things.)
>  
> Please help.
> Thank you.
>  
> Jacob L van Wyk
> Department of Statistics
> University of Johannesburg, APK
> Box 524
> Auckland Park 2006
> South Africa
> Office Tel: +27 11 559 3080
> Fax: +27 11 559 2832
>  
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/R-version-2.7.1-%28warnings%29-tp18542181p18542358.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] R version 2.7.1 (warnings)

[Mark Difford]

Hi Jaap,

Great stuff! As the old adage went, "Go well, go "

Bye, Mark.


Van Wyk, Jaap wrote:
> 
> Thanks, Mark, for the response.
> The problem is vith SciViews. It is not stable under the latest version of
> R.
> I found a solution by downloading the latest version of Tinn-R, which
> communicates with the latest version of R, and now I can carry on like
> normal (with both windows tiled horizontally).
>  
> Jacob
>  
>  
> Jacob L van Wyk
> Department of Statistics
> University of Johannesburg, APK
> Box 524
> Auckland Park 2006
> South Africa
> Office Tel: +27 11 559 3080
> Fax: +27 11 559 2832
>  
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/R-version-2.7.1-%28warnings%29-tp18542181p18552685.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Decoding subscripts/superscripts from CSVs

[Mark Difford]

Hi Nina,

> Subscripts, superscripts, and italics are encoded in carats, and I was
> wondering 
> if R can actually recognize those and print actual superscripts, etc.
> Here's an example:

I know that ladies are fond of diamonds (perhaps you mean carets?); though
it isn't quite clear what you want (though we do know what ...). Can R do
superscripts, subscripts and all the other plotmath stuff? Yes, very much
so; but it still isn't clear what your real question is.

See ?plotmath, in any event.

HTH, Mark.


naw3 wrote:
> 
> Hi,
> 
> I have a CSV file with various biological reactions. Subscripts,
> superscripts,
> and italics are encoded in carats, and I was wondering if R can actually
> recognize those and print actual superscripts, etc. Here's an example:
> 
>  S-adenosyl-L-methionine + rRNA  =  S-adenosyl-L-homocysteine +
> rRNA containing N6-methyladenine
> 
> Thanks,
> -Nina
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Decoding-subscripts-superscripts-from-CSVs-tp18598344p18598727.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Accessing individual variables from summary of lm

[Mark Difford]

Hi Ben,

Sorry (still a little out-of-tune), perhaps what you really need to know
about is ?"["

HTH, Mark.


Mark Difford wrote:
> 
> Hi Ben,
> 
>> If you wouldn't mind, how do I access the individual components inside
>> coefficients matrix?
> 
> What you want to know about is ?attributes
> 
> ##
> attributes(model)
> model$coefficients
> model$coefficients[1]
> model$coefficients[2:4]
> model$coefficients[c(1,5)]
> 
> HTH, Mark.
> 
> 
> ascentnet wrote:
>> 
>> Thanks for your help!  If you wouldn't mind, how do I access the
>> individual components inside coefficients matrix?
>> 
>> thanks,
>> Ben.
>> 
>> Paul Smith wrote:
>>> 
>>> On Tue, Jul 22, 2008 at 5:44 AM, ascentnet <[EMAIL PROTECTED]> wrote:
>>>> I am trying to access the R2, intercept, and slope out of the summary
>>>> from an
>>>> lm so that I can insert it into a database using RODBC.  How can I
>>>> access
>>>> these individually?
>>> 
>>> Try:
>>> 
>>> data(swiss)
>>> model <- lm(Infant.Mortality ~ .,data=swiss)
>>> summary(model)$r.square
>>> # The following returns a matrix with the coefficients
>>> summary(model)$coefficients
>>> 
>>> Paul
>>> 
>>> __
>>> R-help@r-project.org 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.
>>> 
>>> 
>> 
>> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Accessing-individual-variables-from-summary-of-lm-tp18582274p18604672.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Accessing individual variables from summary of lm

[Mark Difford]

Hi Ben,

> If you wouldn't mind, how do I access the individual components inside
> coefficients matrix?

What you want to know about is ?attributes

##
attributes(model)
model$coefficients
model$coefficients[1]
model$coefficients[2:4]
model$coefficients[c(1,5)]

HTH, Mark.


ascentnet wrote:
> 
> Thanks for your help!  If you wouldn't mind, how do I access the
> individual components inside coefficients matrix?
> 
> thanks,
> Ben.
> 
> Paul Smith wrote:
>> 
>> On Tue, Jul 22, 2008 at 5:44 AM, ascentnet <[EMAIL PROTECTED]> wrote:
>>> I am trying to access the R2, intercept, and slope out of the summary
>>> from an
>>> lm so that I can insert it into a database using RODBC.  How can I
>>> access
>>> these individually?
>> 
>> Try:
>> 
>> data(swiss)
>> model <- lm(Infant.Mortality ~ .,data=swiss)
>> summary(model)$r.square
>> # The following returns a matrix with the coefficients
>> summary(model)$coefficients
>> 
>> Paul
>> 
>> __
>> R-help@r-project.org 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.
>> 
>> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Accessing-individual-variables-from-summary-of-lm-tp18582274p18604576.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Coefficients of Logistic Regression from bootstrap - how to get them?

[Mark Difford]

Hi All,

It really comes down to a question of attitude: you either want to learn
something fundamental or core and so bootstrap yourself to a "better" place
(at least away from where you are), or you don't. As Marc said, Michal seems
to have erected a wall around his thinking.

I don't think it's fair to take pot shots at Frank for not wanting to
promote or further something he doesn't believe in. He's a regular
contributor to the list, who gives sound advice. He's also one of the few
experts on the list who is prepared to give statistical advice.

Regards, Mark.


Rolf Turner-3 wrote:
> 
> 
> On 23/07/2008, at 1:17 PM, Frank E Harrell Jr wrote:
> 
>> Michal Figurski wrote:
>>> Hmm...
>>> It sounds like ideology to me. I was asking for technical help. I  
>>> know what I want to do, just don't know how to do it in R. I'll go  
>>> back to SAS then. Thank you.
>>> -- 
>>> Michal J. Figurski
>>
>> You don't understand any of the theory and you are using techniques  
>> you don't understand and have provided no motivation for.  And you  
>> are the one who is frustrated with others.  Wow.
> 
>   Come off it guys.  It is indeed very frustrating when one asks ``How  
> can I do X''
>   and gets told ``Don't do X, do Y.''  It may well be the case that  
> doing X is
>   wrong-headed, misleading, and may cause the bridge to fall down, or  
> the world to
>   come to an end.  Fair enough to point this out --- but then why not  
> just tell
>   the poor beggar, who asked, how to do X?
> 
>   The only circumstance in which *not* telling the poor beggar how to  
> do X is
>   justified is that in which it takes considerable *work* to figure  
> out how to
>   do X.  In this case it is perfectly reasonable to say ``I think  
> doing X is
>   stupid so I am not going to waste my time figuring out for you how  
> to do it.''
> 
>   I don't know enough about the bootstrapping software (don't know  
> *anything*
>   about it actually) to know whether the foregoing circumstance  
> applies here.
>   But I suspect it doesn't.  And I suspect that you (Frank) could tell  
> Michal in
>   a few lines the answer to the question that he *asked* (as opposed,  
> possibly,
>   to the question that he should have asked).
> 
>   If it were my problem I'd just write my own bootstrapping function  
> to apply
>   to the problem in hand.  It can't be that hard ... just a for loop  
> and a
>   call to sample(...,replace=TRUE).
> 
>   If you can write macros in SAS then .
> 
>   cheers,
> 
>   Rolf
> 
> ##
> Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Coefficients-of-Logistic-Regression-from-bootstrap---how-to-get-them--tp18570684p18605881.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Histogram

[Mark Difford]

Hi Angelo,

Look carefully at package vcd; and at log-linear models (+ glm(...,
family=poisson)). For overdispersion there are more advanced methods.

HTH, Mark.


Angelo Scozzarella wrote:
> 
> Hi,
> 
> how can I treat data organised in classes and frequencies?
> 
> Ex.
> 
> class frequency
> 
> 20-23 9
> 23-25 7
> 26-28 5
> 29-31 5
> 32-34 3
> 
> 
> 
> Thanks
> 
> Angelo Scozzarella
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Histogram-tp18617432p18619628.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Just 2 more questions - for now!

[Mark Difford]

Hi Robin,

>> I ... can't get lm to work despite reading the help. I can get it to work
>> with a single 
>> explanatory variable, EG model <- lm(data$admissions~data$maxitemp)

I'll answer just the second of your questions. Advice: don't just read the
help file, look at the examples and run them; look at how the calls are
made. lm() and many other programs use (or can use) the ?formula interface,
so learn how to use it. For your "example," rather do

##
EGModel <- lm(admissions ~ maxitemp, data = data)

where data is your data set (to avoid confusion call it something else, e.g.
mydata). To add another variable/predictor, which should be in your data
frame (for simplicity) called data, simply do

EGModel <- lm(admissions ~ maxitemp+another.var, data = data)

If your maxitemp(s) occur at different sites, then you might be interested
in a nested model

EGModel1 <- lm(admissions ~ site/maxitemp - 1, data = data)
EGModel2 <- lm(admissions ~ site*maxitemp, data = data) ## alt.
parameterization
anova(EGModel1, EGModel2)
anova(EGModel, EGModel1, EGModel2, test="Chi")

This shows part of the power of using the formula interface.

HTH, Mark.


Williams, Robin wrote:
> 
> Hi all, 
> Thanks for the help with my previous post. 
> I have just two more questions for the minute. 
>   I think I said in a previous post that I like to use the terminal,
> i.e. run rterm.exe. On exiting the terminal, I am asked if I want to
> save the workspace. If I hit y (yes), the workspace is just saved as
> .rdata in my working directory, does anyone know how I can name it
> directly from the terminal? 
>   More importantly, I can't then open this file from the terminal.
> Obviously loading the file from windows explorer brings up the GUI.
> Anyone know the command I need? All I can think of doing is adding
> rterm.exe to my path and running it from the command prompt (adding the
> file as an argument to my command), but surely there is an easy way to
> do this from R? Of course I would like to have the terminal open and
> open and close various workspaces in one session, without wanting to
> restart R all the time. 
>  
>   Finally, I rather embarrasingly can't get lm to work despite reading
> the help. I can get it to work with a single explanatory variable, EG 
> model <- lm(data$admissions~data$maxitemp)
> but how to include a second explanatory variable (minitemp)? I keep
> getting errors. Surely I don't need to use
> c(data$maxitemp,data$minitemp) etc? 
>  
> All help greatly appreciated - I am getting there slowly!
>  
> 
> Robin Williams
> Met Office summer intern - Health Forecasting
> [EMAIL PROTECTED] 
> 
>  
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Just-2-more-questions---for-now%21-tp18629686p18630709.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] factor question

[Mark Difford]

Hi Edna,

Because I am "always" subsetting, I keep the following function handy

mydata[] <- lapply(mydata, function(x) if(is.factor(x)) x[,drop=T] else x)

This will strip out all factor levels that have been dropped by a previous
subsetting operation. For novice users of R (though I am not suggesting that
you are) it's priceless.

The original author, as far as I know, is Andy Liaw (author of the
randomForest package).

HTH, Mark.


Edna Bell wrote:
> 
> Hi!
> 
> Suppose I have a factor:
> 
>> fac1 <- factor(rep(c("dog","cat","tree"),c(2,3,4)))
>> fac1
> [1] dog  dog  cat  cat  cat  tree tree tree tree
> Levels: cat dog tree
>> length(fac1)
> [1] 9
> 
> Now I want to get rid of the dogs:
>> fac1 <- fac1[3:9]
>> fac1
> [1] cat  cat  cat  tree tree tree tree
> Levels: cat dog tree
>>
> How should I remove the dog from the level, please?
> 
> Thanks in advance,
> Edna Bell
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/factor-question-tp18638814p18639010.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] R-package install

[Mark Difford]

Hi Cindy,

>> Hi, I have a R package, which is a .tar.tar file. And it is said to be
>> source code for all 
>> platforms. ... I am wondering if I can use this package in Windows R.

If it is source code you would first need to compile it to binary format
before you can use it. Can you do that?

Which package is it?

Regards, Mark.


cindy Guo wrote:
> 
> Hi, I have a R package, which is a .tar.tar file. And it is said to be
> source code for all platforms. And the author said it should install on
> any
> system (but he didn't know about Windows). I am wondering if I can use
> this
> package in Windows R.
> Thank you,
> Cindy
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/R-package-install-tp18636993p18639516.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] extracting Pr>ltl from robcov/ols (Design)

[Mark Difford]

ldb (?),

>> I am trying to extract significance levels out of a robcov+ols call
>> However, I can't figure out how to get to the significance (Pr>ltl).
>> It is obviously calculating it because the call:

It's calculated in print.ols(). See extract below. To see the whole
function, do

print.ols

## Extract from print.ols in package Design
se <- sqrt(diag(x$var))
z <- x$coefficients/se
P <- 2 * (1 - pt(abs(z), rdf))
co <- cbind(x$coefficients, se, z, P)
dimnames(co) <- list(names(x$coefficients), c("Value", "Std. Error", 
"t", "Pr(>|t|)"))
print(co)

HTH, Mark.


ldb-5 wrote:
> 
> I am trying to extract significance levels out of a robcov+ols call.
> 
> For background: I am analysing data where multiple measurements(2 per
> topic) were taken from individuals(36) on their emotional reaction
> (dependent variable) to various topics (3 topics). Because I have
> several emotions and a rotation to do on the topics, I'd like to have
> the results pumped into a nice table.
> 
> answer<-robcov(ols(emotion ~ topic,x=TRUE,y=TRUE),individual)
> 
>>From the robcov help it warns me:
> 
>>Adjusted ols fits do not have the corrected standard errors printed with
print.ols. Use sqrt(diag(adjfit$var)) to get this, where adjfit is the
result of robcov.
> 
> So I can get to the standard error by:
> 
> answer.se<-sqrt(diag(answer$var))
> 
> I can get to the coefficients by:
> 
> answer.co<-coefficients(answer)
> 
> I can get to the t-values by:
> 
> answer.t<-coefficients(answer)/ sqrt(diag(answer$var))#t-value
> 
> However, I can't figure out how to get to the significance (Pr>ltl).
> It is obviously calculating it because the call:
> 
> answer
> 
> provides it, but I can't figure out where it is getting it from or how
> it is calculating.
> 
> thanks for any help.
> 
> ldb
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/extracting-Pr%3Eltl-from-robcov-ols-%28Design%29-tp18651539p18651964.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Beginning lm

[Mark Difford]

Hi Kevin,

>> Can anyone give me a short tutorial on the formula syntax? ... I am sorry
>> but I could not 
>> glean this information from the help page on lm.

You can give yourself a very good tutorial by reading ?formula and Chapter
12 of

file://localhost/C:/Program%20Files/R/R-2.7.1pat/doc/manual/R-intro.html

Mark.


rkevinburton wrote:
> 
> I have ussed lm to generate a basic line correlation:
> 
> fit = lm(hours.of.sleep ~ ToSleep)
> 
> Note: From "Bayesian Computation with R", Jim Albert, p. 7
> 
> I understand the simple y = mx + b line that this fits the data to. Now
> apparently I don't understand formulas. The documentation indicates that
> there is an implied "intercept" in the formula so now I want to try and
> fit the line to a second degree polynomial so I tried:
> 
> ft = lm(hours.of.sleep ~ ToSleep ^ 2 + ToSleep)
> 
> and I still seem to get results that indicate a slope intercept, y = mx +
> b, type of fit. Can anyone give me a short tutorial on the formula syntax?
> I would like to fit the data to 2nd and higher order polynomials, 1 / x,
> log(x), etc. I am sorry but I could not glean this information from the
> help page on lm.
> 
> Thank you.
> 
> Kevin
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Beginning-lm-tp18668114p18669143.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] A easy way to write formula

[Mark Difford]

Hi Jinsong and Thierry,

>> (x1 + x2 + x3) ^2 will give you the main effects and the interactions.

Although it wasn't specifically requested it is perhaps important to note
that (...)^2 doesn't expand to give _all_ interaction terms, only
interactions to the second order, so the interaction term x1:x2:x3 will be
missing, To get these, do

(x1 + x2 + x3)^3

Regards, Mark.


ONKELINX, Thierry wrote:
> 
> (x1 + x2 + x3) ^2 will give you the main effects and the interactions.
> So this will shorten your equation to 
> 
> ft <- lm(y ~ (x1 + x2 + x3) ^2 + I(x1 * x1)+ I(x2 * x2) + I(x3 * x3),
> mydata)
> 
> HTH,
> 
> Thierry
> 
> 
> 
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
> and Forest
> Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
> methodology and quality assurance
> Gaverstraat 4
> 9500 Geraardsbergen
> Belgium 
> tel. + 32 54/436 185
> [EMAIL PROTECTED] 
> www.inbo.be 
> 
> To call in the statistician after the experiment is done may be no more
> than asking him to perform a post-mortem examination: he may be able to
> say what the experiment died of.
> ~ Sir Ronald Aylmer Fisher
> 
> The plural of anecdote is not data.
> ~ Roger Brinner
> 
> The combination of some data and an aching desire for an answer does not
> ensure that a reasonable answer can be extracted from a given body of
> data.
> ~ John Tukey
> 
> -Oorspronkelijk bericht-
> Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
> Namens Jinsong Zhao
> Verzonden: zondag 27 juli 2008 9:39
> Aan: r-help@r-project.org
> Onderwerp: [R] A easy way to write formula
> 
> Hi
> 
> I have a data frame, including x1, x2, x3, and y. I use lm() to fit
> second-order linear model, like the following:
> 
> ft <- lm(y ~ x1 + x2 + x3 + I(x1 * x1) + I(x1 * x2) + I(x1 * x3) + I(x2
> * x2) + I(x2 * x3) + I(x3 * x3), mydata)
> 
> if the independent variable number is large, the formula will be very
> long. Is there a easy way to write formula like the above one? I have
> read the R introduction, however, I don't find a easy way.
> 
> Any hints will be appreciated. Thanks,
> 
> Jinsong
> 
> __
> R-help@r-project.org 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-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/A-easy-way-to-write-formula-tp18674075p18674632.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] graphing regression coefficients and standard errors

[Mark Difford]

Hi Murali,

>> I am interested in plotting my regression analysis results(regression
>> coefficients and 
>> standard errors obtained through OLS and Tobit models) in the form of
>> graphs.

plot(obj$lm) will give you a set of diagnostic plots. What you seem to be
after is ?termplot. Also look at John Fox's effects package and at Frank
Harrell's Design package (sub functions ols and plot.Design).

Regards, Mark.


Murali K wrote:
> 
> Hello,
> I am interested in plotting my regression analysis
> results(regression coefficients and standard errors obtained through OLS
> and
> Tobit models)
> in the form of graphs.What is the best way to accomplish this? Thanks.
> 
> 
> Murali Kuchibhotla
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/graphing-regression-coefficients-and-standard-errors-tp18660870p18685184.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] multiv

[Mark Difford]

Hi Ileana,

See this thread:

http://www.nabble.com/R-package-install-td18636993.html

HTH, Mark.


Somesan, Ileana wrote:
> 
>> Hello,
>> 
>> I want to install the package "multiv" which is not maintained any
>> more (found in the archive: multiv_1.1-6.tar.gz from 16 July 2003). I
>> have installed an older version of R 1.4.1 on my Windows Vista in
>> order to match the older package. 
>> 
>> As no binary is distributed, is it necessary to compile the package?
>> If so, how? 
>> 
>> Thanks a lot,
>> Ileana
>> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/multiv-tp18695040p18699787.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Help interpreting density().

[Mark Difford]

Hi Kevin,

>> The documentation indicates that the bw is essentially the sd.
>> > d <- density(rnorm(1000))

Not so. The documentation states that the following about "bw": "The kernels
are scaled such that this is the standard deviation of the smoothing
kernel...," which is a very different thing.

The default bandwidth used by density is ?bw.nrd0. Read that documentation
carefully and all might be clear.

HTH, Mark.


rkevinburton wrote:
> 
> I issue the following:
> 
>> d <- density(rnorm(1000))
>> d
> 
> and get:
> 
> Call:
> density.default(x = rnorm(1000))
> 
> Data: rnorm(1000) (1000 obs.);  Bandwidth 'bw' = 0.2235
> 
>x y
>  Min.   :-3.5157   Min.   :2.416e-05  
>  1st Qu.:-1.6892   1st Qu.:1.129e-02  
>  Median : 0.1373   Median :7.267e-02  
>  Mean   : 0.1373   Mean   :1.367e-01  
>  3rd Qu.: 1.9639   3rd Qu.:2.693e-01  
>  Max.   : 3.7904   Max.   :4.014e-01  
> 
> The documentation indicates that the bw is essentially the sd. Yet I have
> specified an sd of 1? How am I to interpret the ranges of the values? x
> ranges almost from -4 to +4 and y ranges from 0 to 0.4. The mean x is .1
> which isn't too awfully close to what I would expect (0.0). Then there is:
> 
>> d <- density(rpois(1000,0))
>> d
> 
> Call:
> density.default(x = rpois(1000, 0))
> 
> Data: rpois(1000, 0) (1000 obs.);   Bandwidth 'bw' = 0.2261
> 
>x y  
>  Min.   :-0.6782   Min.   :0.01979  
>  1st Qu.:-0.3391   1st Qu.:0.14073  
>  Median : 0.   Median :0.57178  
>  Mean   : 0.   Mean   :0.73454  
>  3rd Qu.: 0.3391   3rd Qu.:1.32830  
>  Max.   : 0.6782   Max.   :1.76436  
> 
> Here I am getting the mean that I expect from a Poisson distribuition but
> y ranges from 0 to 1.75. Again I am not sure what these numbers mean. How
> can I map the output to the standard distirbution description parameters?
> 
> Thank you.
> 
> Kevin
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Help-interpreting-density%28%29.-tp18704955p18706154.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Help interpreting density().

[Mark Difford]

Hi Kevin,

>> I still have my original question. How does the output relate to
>> estimating the parameters 
>> of a given density? I read that for a gausian kernal:

This isn't the place for such questions: you need to do some _basic_ reading
on the subject so that you begin to understand something about the method
you are messing about with. Basically (very basically) it's a smoothed out
histogram. And you will probably (?) know that a histogram is [still used]
to show you how a set of univariate data (random variable) is distributed.

Perhaps start with http://en.wikipedia.org/wiki/Kernel_density, then go
somewhere else. But since you have access to the web you really should have
found something like this yourself.

Regards, Mark.


rkevinburton wrote:
> 
> OK. Thank you for pointing out my mistake.
> 
> I still have my original question. How does the output relate to
> estimating the parameters of a given density? I read that for a gausian
> kernal:
> 
> bw.nrd0 implements a rule-of-thumb for choosing the bandwidth of a
> Gaussian kernel density estimator. It defaults to 0.9 times the minimum of
> the standard deviation and the interquartile range divided by 1.34 times
> the sample size to the negative one-fifth power (= Silverman's ‘rule of
> thumb’
> 
> But how does that relate to say a Poisson distribution or a two-parameter
> distribution like a normal, beta, or binomial distribution?
> 
> Thank you.
> 
> Kevin
> 
>  Mark Difford <[EMAIL PROTECTED]> wrote: 
>> 
>> Hi Kevin,
>> 
>> >> The documentation indicates that the bw is essentially the sd.
>> >> > d <- density(rnorm(1000))
>> 
>> Not so. The documentation states that the following about "bw": "The
>> kernels
>> are scaled such that this is the standard deviation of the smoothing
>> kernel...," which is a very different thing.
>> 
>> The default bandwidth used by density is ?bw.nrd0. Read that
>> documentation
>> carefully and all might be clear.
>> 
>> HTH, Mark.
>> 
>> 
>> rkevinburton wrote:
>> > 
>> > I issue the following:
>> > 
>> >> d <- density(rnorm(1000))
>> >> d
>> > 
>> > and get:
>> > 
>> > Call:
>> > density.default(x = rnorm(1000))
>> > 
>> > Data: rnorm(1000) (1000 obs.);  Bandwidth 'bw' = 0.2235
>> > 
>> >x y
>> >  Min.   :-3.5157   Min.   :2.416e-05  
>> >  1st Qu.:-1.6892   1st Qu.:1.129e-02  
>> >  Median : 0.1373   Median :7.267e-02  
>> >  Mean   : 0.1373   Mean   :1.367e-01  
>> >  3rd Qu.: 1.9639   3rd Qu.:2.693e-01  
>> >  Max.   : 3.7904   Max.   :4.014e-01  
>> > 
>> > The documentation indicates that the bw is essentially the sd. Yet I
>> have
>> > specified an sd of 1? How am I to interpret the ranges of the values? x
>> > ranges almost from -4 to +4 and y ranges from 0 to 0.4. The mean x is
>> .1
>> > which isn't too awfully close to what I would expect (0.0). Then there
>> is:
>> > 
>> >> d <- density(rpois(1000,0))
>> >> d
>> > 
>> > Call:
>> > density.default(x = rpois(1000, 0))
>> > 
>> > Data: rpois(1000, 0) (1000 obs.);   Bandwidth 'bw' = 0.2261
>> > 
>> >x y  
>> >  Min.   :-0.6782   Min.   :0.01979  
>> >  1st Qu.:-0.3391   1st Qu.:0.14073  
>> >  Median : 0.   Median :0.57178  
>> >  Mean   : 0.   Mean   :0.73454  
>> >  3rd Qu.: 0.3391   3rd Qu.:1.32830  
>> >  Max.   : 0.6782   Max.   :1.76436  
>> > 
>> > Here I am getting the mean that I expect from a Poisson distribuition
>> but
>> > y ranges from 0 to 1.75. Again I am not sure what these numbers mean.
>> How
>> > can I map the output to the standard distirbution description
>> parameters?
>> > 
>> > Thank you.
>> > 
>> > Kevin
>> > 
>> > __
>> > R-help@r-project.org 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.
>> > 
>> > 
>> 
>> -- 
>> View this message in context:
>> http://www.nabble.com/Help-interpreting-density%28%29.-tp18704955p18706154.html
>> Sent from the R

Re: [R] Help interpreting density().

[Mark Difford]

Hi Kevin,

Clicking on the link I sent gets me there (?), though things are pretty slow
at the moment. Perhaps try this related link, and from it get back to the
first one:

http://en.wikipedia.org/wiki/Density_estimation

You can also get to this via histogram, so search for that in Wiki, and
then...

HTH, Mark.


rkevinburton wrote:
> 
> Sorry I tried WikiPedia and only found:
> 
> Wikipedia does not have an article with this exact name.
> 
> I will try to find some other sources of information.
> 
> Kevin
> 
>  Mark Difford <[EMAIL PROTECTED]> wrote: 
>> 
> Hi Kevin,
> 
>>> I still have my original question. How does the output relate to
>>> estimating the parameters 
>>> of a given density? I read that for a gausian kernal:
> 
> This isn't the place for such questions: you need to do some _basic_
> reading
> on the subject so that you begin to understand something about the method
> you are messing about with. Basically (very basically) it's a smoothed out
> histogram. And you will probably (?) know that a histogram is [still used]
> to show you how a set of univariate data (random variable) is distributed.
> 
> Perhaps start with http://en.wikipedia.org/wiki/Kernel_density, then go
> somewhere else. But since you have access to the web you really should
> have
> found something like this yourself.
> 
> Regards, Mark.
> 
> 
> rkevinburton wrote:
>> 
>> OK. Thank you for pointing out my mistake.
>> 
>> I still have my original question. How does the output relate to
>> estimating the parameters of a given density? I read that for a gausian
>> kernal:
>> 
>> bw.nrd0 implements a rule-of-thumb for choosing the bandwidth of a
>> Gaussian kernel density estimator. It defaults to 0.9 times the minimum
>> of
>> the standard deviation and the interquartile range divided by 1.34 times
>> the sample size to the negative one-fifth power (= Silverman's ‘rule of
>> thumb’
>> 
>> But how does that relate to say a Poisson distribution or a two-parameter
>> distribution like a normal, beta, or binomial distribution?
>> 
>> Thank you.
>> 
>> Kevin
>> 
>>  Mark Difford <[EMAIL PROTECTED]> wrote: 
>>> 
>>> Hi Kevin,
>>> 
>>> >> The documentation indicates that the bw is essentially the sd.
>>> >> > d <- density(rnorm(1000))
>>> 
>>> Not so. The documentation states that the following about "bw": "The
>>> kernels
>>> are scaled such that this is the standard deviation of the smoothing
>>> kernel...," which is a very different thing.
>>> 
>>> The default bandwidth used by density is ?bw.nrd0. Read that
>>> documentation
>>> carefully and all might be clear.
>>> 
>>> HTH, Mark.
>>> 
>>> 
>>> rkevinburton wrote:
>>> > 
>>> > I issue the following:
>>> > 
>>> >> d <- density(rnorm(1000))
>>> >> d
>>> > 
>>> > and get:
>>> > 
>>> > Call:
>>> > density.default(x = rnorm(1000))
>>> > 
>>> > Data: rnorm(1000) (1000 obs.);  Bandwidth 'bw' = 0.2235
>>> > 
>>> >x y
>>> >  Min.   :-3.5157   Min.   :2.416e-05  
>>> >  1st Qu.:-1.6892   1st Qu.:1.129e-02  
>>> >  Median : 0.1373   Median :7.267e-02  
>>> >  Mean   : 0.1373   Mean   :1.367e-01  
>>> >  3rd Qu.: 1.9639   3rd Qu.:2.693e-01  
>>> >  Max.   : 3.7904   Max.   :4.014e-01  
>>> > 
>>> > The documentation indicates that the bw is essentially the sd. Yet I
>>> have
>>> > specified an sd of 1? How am I to interpret the ranges of the values?
>>> x
>>> > ranges almost from -4 to +4 and y ranges from 0 to 0.4. The mean x is
>>> .1
>>> > which isn't too awfully close to what I would expect (0.0). Then there
>>> is:
>>> > 
>>> >> d <- density(rpois(1000,0))
>>> >> d
>>> > 
>>> > Call:
>>> > density.default(x = rpois(1000, 0))
>>> > 
>>> > Data: rpois(1000, 0) (1000 obs.);   Bandwidth 'bw' = 0.2261
>>> > 
>>> >x y  
>>> >  Min.   :-0.6782   Min.   :0.01979  
>>> >  1st Qu.:-0.3391   1st Qu.:0.14073  
>>> >  Median : 0.   Median :0.57178  
>>> >  Mean   : 0.

Re: [R] Re peated Measure ANOVA-Old question

[Mark Difford]

Hi Chunhao,

>> I google the website and I found that there are three ways to perform  
>> repeated measure ANOVA: aov, lme and lmer.

It's also a good idea to search through the archives.

>> I use the example that is provided in the above link and I try
>> > tt<-aov(p.pa~group*time+Error(subject/time),data=P.PA)
>> > TukeyHSD(tt)
>> Error in UseMethod("TukeyHSD") : no applicable method for "TukeyHSD"

Because of ... + Error().

See the following thread, where this is answered:

http://www.nabble.com/Tukey-HSD-(or-other-post-hoc-tests)-following-repeated-measures-ANOVA-td17508294.html#a17559307

HTH, Mark.


ctu wrote:
> 
> Hi R users,
> I google the website and I found that there are three ways to perform  
> repeated measure ANOVA: aov, lme and lmer.
> http://www.mail-archive.com/[EMAIL PROTECTED]/msg58502.html
> But the questions are which one is good to use and how to do post-hoc
> test?
> 
> I use the example that is provided in the above link and I try
>> tt<-aov(p.pa~group*time+Error(subject/time),data=P.PA)
>> TukeyHSD(tt)
> Error in UseMethod("TukeyHSD") : no applicable method for "TukeyHSD"
> 
> Any suggestions?
> 
> Appreciate in advance~
> Chunhao
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Repeated-Measure-ANOVA-Old-question-tp18726146p18726722.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Re peated Measure ANOVA-Old question

[Mark Difford]

Hi Chunhao,

If you carefully read the posting that was referred to you will see that
lme() and not lmer() was used as an example (for using with the multcomp
package). lmer() was only mentioned as an aside... lmer() is S4 and doesn't
work with multcomp, which is S3.

Apropos of specifying random effects: an example was given in the posting
showing that aov(... +Error(...)) and lme(..., random =...) give the same
result. Use that as your guideline.

HTH, Mark.


ctu wrote:
> 
> Tank you (Anna & Mark) very much for this super information
> Let me confirm that IF I want to perform the RM-ANOVA
> I should use "lmer" and perform the post-hoc test by using "glht", right?
> 
> Because I am not so familiar with "lmer" so I have two more questions.
> Let me use that following example again,
> 1. Is the random effect statement correct? if no what will be?
> 2. What did I do wrong in glht?
> 
> 
> require(lme4);require(nlme)
> group <-rep(rep(1:2, c(5,5)), 3);time <-rep(1:3, rep(10,3));subject <-  
> rep(1:10,3)
> p.pa <- c(92, 44, 49, 52, 41, 34, 32, 65, 47, 58, 94, 82, 48, 60, 47,
> 46, 41, 73, 60, 69, 95, 53, 44, 66, 62, 46, 53, 73, 84, 79)
> P.PA <- data.frame(subject, group, time, p.pa)
> P.PA$group=factor(P.PA$group);P.PA$time=factor(P.PA$time)
> P.PA$subject=factor(P.PA$subject)
> 
> tlmer<-lmer (p.pa ~ time * group + (time*group | subject), data=P.PA )
> 
>> summary(glht(tlmer,linfct=mcp(group="Tukey")))
> Error in UseMethod("fixef") : no applicable method for "fixef"
> In addition: Warning message
> 
> Thank you very much
> Chunhao
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Repeated-Measure-ANOVA-Old-question-tp18726146p18727795.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Mixed effects model where nested factor is not the repeated across treatments lme???

[Mark Difford]

Hi Miki and Chunhao,

<< Rusers (Anna, and Mark {thank you guys}) provide me a vary valuable  
<< information.

Also see Gavin Simpson's posting earlier today: apparently multcomp does now
work with lmer objects (it's gone through phases of not working, then
working: it's still being developed). Beware, though, that random effects
are specified differently, so it's not as easy to recast an aov(... +
Error(...)) term structure as an equivalent random effect's structure.

HTH, Mark.


ctu wrote:
> 
> Hi Miki,
> I just got the same problem with you couple hours ago.
> Rusers (Anna, and Mark {thank you guys}) provide me a vary valuable  
> information.
> link to following address.
> 
> http://www.nabble.com/Tukey-HSD-(or-other-post-hoc-tests)-following-repeated-measures-ANOVA-td17508294.html#a17559307
> for the A vs. B, A vs. C
> You could install and download the multcomp package and perform the  
> post hoc test
> such as
> summary(glht(lmel,linfct=mcp(treatment="Tukey")))
> 
> hopefully it helps
> Chunhao
> 
> 
> Quoting M Ensbey <[EMAIL PROTECTED]>:
> 
>> Hi,
>>
>>
>>
>> I have searched the archives and can't quite confirm the answer to this.
>> I appreciate your time...
>>
>>
>>
>> I have 4 treatments (fixed) and I would like to know if there is a
>> significant difference in metal volume (metal) between the treatments.
>> The experiment has 5 blocks (random) in each treatment and no block is
>> repeated across treatments. Within each plot there are varying numbers
>> of replicates (random) (some plots have 4 individuals in them some have
>> 14 and a range in between). NOTE the plots in one treatment are not
>> replicated in the others.
>>
>>
>>
>> So I end up with a data.frame with 4 treatments repeated down one column
>> (treatment=A, B, C, D), 20 plots repeated down the next (block= 1 to 20)
>> and records for metal volume (metal- 124 of these)
>>
>> I have made treatment and block a factor. But haven't grouped them (do I
>> need to and how if so)
>>
>>
>>
>> The main question is in 3 parts:
>>
>>
>>
>> 1.   is this the correct formula to use for this situation:
>> lme1<-lme(metal~treatment,data=data,random=~1|block) (or is lme even the
>> right thing to use here?)
>>
>>
>>
>> I get:
>>
>>> summary(lme1)
>>
>> Linear mixed-effects model fit by REML
>>
>>  Data: data
>>
>>AIC  BIClogLik
>>
>>   365.8327 382.5576 -176.9163
>>
>>
>>
>> Random effects:
>>
>>  Formula: ~1 | block
>>
>> (Intercept)  Residual
>>
>> StdDev:   0.4306096 0.9450976
>>
>>
>>
>> Fixed effects: Cu ~ Treatment
>>
>>  Value Std.Error  DF   t-value p-value
>>
>> (Intercept)   5.587839 0.2632831 104 21.223688  0. ***
>>
>> TreatmentB -0.970384 0.3729675  16 -2.601792  0.0193 ***
>>
>> TreatmentC  -1.449250 0.3656351  16 -3.963651  0.0011 ***
>>
>> TreatmentD  -1.319564 0.3633837  16 -3.631323  0.0022 ***
>>
>>  Correlation:
>>
>>  (Intr) TrtmAN TrtmCH
>>
>> TreatmentB -0.706
>>
>> TreatmentC  -0.720  0.508
>>
>> TreatmentD  -0.725  0.511  0.522
>>
>>
>>
>> Standardized Within-Group Residuals:
>>
>> Min  Q1 Med  Q3 Max
>>
>> -2.85762206 -0.68568460 -0.09004478  0.56237152  3.20650288
>>
>>
>>
>> Number of Observations: 124
>>
>> Number of Groups: 20
>>
>>
>>
>> 2.   if so how can I get p values for comparisons between every
>> group... ie is A different from B, is A different from C, is A different
>> from D, is B different from C, is B different from D etc... is there a
>> way to get all of these instead of just "is A different from B, is A
>> different from C, is A different from D" which summary seems to give?
>> 3.   last of all what is the best way to print out all the residuals
>> for lme... I can get qqplot(lme1) is there a pre-programmed call for
>> multiple diagnostic plots like in some other functions...
>>
>>
>>
>>
>>
>> Thankyou so Much for your time
>>
>>
>>
>> It is much appreciated
>>
>> ;-)
>>
>>
>>
>> Miki
>>
>>
>>
>>
>>  [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org 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-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Mixed-effects-model-where-nested-factor-is-not-the-repeated-across-treatments-lmetp18732327p18738202.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http

Re: [R] Coefficients of Logistic Regression from bootstrap - how to get them?

[Mark Difford]

Hi Michal,

>> It is that you just don't estimate mean, or CI, or variance on PK profile
>> data! 
>> It is as if you were trying to estimate mean, CI and variance of a 
>> "Toccata_&_Fugue_in_D_minor.wav" file. What for? The point is in the
>> music! 
>> Would the mean or CI or variance tell you anything about that?

I could be way of beam, but it strikes me that one reason for your "question
marks" or "issues" is that at bottom you are trying to use an essentially
traditional tool to fit "loopy" data and you can't see how that can work
...? If so, take a look at the fda (functional data analysis) package. It's
based on the work of Ramsey & Silverman (2005 --- there is an earlier pub.)
and may be closer to what you want. Word of caution: the subject matter is
pretty technical. You can get a taste of what's on offer at

http://www.functionaldata.org/
http://ego.psych.mcgill.ca/misc/fda/

HTH, Mark.


Gustaf Rydevik wrote:
> 
> On Thu, Jul 31, 2008 at 4:30 PM, Michal Figurski
> <[EMAIL PROTECTED]> wrote:
>> Frank and all,
>>
>> The point you were looking for was in a page that was linked from the
>> referenced page - I apologize for confusion. Please take a look at the
>> two
>> last paragraphs here:
>> http://people.revoledu.com/kardi/tutorial/Bootstrap/examples.htm
>>
>> Though, possibly it's my ignorance, maybe it's yours, but you actually
>> missed the important point again. It is that you just don't estimate
>> mean,
>> or CI, or variance on PK profile data! It is as if you were trying to
>> estimate mean, CI and variance of a "Toccata_&_Fugue_in_D_minor.wav"
>> file.
>> What for? The point is in the music! Would the mean or CI or variance
>> tell
>> you anything about that? Besides, everybody knows the variance (or
>> variability?) is there and can estimate it without spending time on
>> calculations.
>> What I am trying to do is comparable to compressing a wave into mp3 - to
>> predict the wave using as few data points as possible. I have a bunch of
>> similar waves and I'm trying to find a common equation to predict them
>> all.
>> I am *not* looking for the variance of the mean!
>>
>> I could be wrong (though it seems less and less likely), but you keep
>> talking about the same irrelevant parameters (CI, variance) on and on.
>> Well,
>> yes - we are at a standstill, but not because of Davison & Hinkley's
>> book. I
>> can try reading it, though as I stated above, it is not even "remotely
>> related" to what I am trying to do. I'll skip it then - life is too
>> short.
>>
>> Nevertheless I thank you (all) for relevant criticism on the procedure
>> (in
>> the points where it was relevant). I plan to use this methodology
>> further,
>> and it was good to find out that it withstood your criticism. I will look
>> into the penalized methods, though.
>>
>> --
>> Michal J. Figurski
>>
> 
> I take it you mean the sentence:
> 
> " For example, in here, the statistical estimator is  the sample mean.
> Using bootstrap sampling, you can do beyond your statistical
> estimators. You can now get even the distribution of your estimator
> and the statistics (such as confidence interval, variance) of your
> estimator."
> 
> Again you are misinterpreting text. The phrase about "doing beyond
> your statistical estimators", is explained in the next sentence, where
> he says that using bootstrap gives you information about the mean
> *estimator* (and not more information about the population mean).
> And since you're not interested in this information, in your case
> bootstrap/resampling is not useful at all.
> 
> As another example of misinterpretation: In your email from  a week
> ago, it sounds like you believe that the authors of the original paper
> are trying to improve on a fixed model
> Figurski:
> "Regarding the "multiple stepwise regression" - according to the cited
> SPSS manual, there are 5 options to select from. I don't think they used
> 'stepwise selection' option, because their models were already
> pre-defined. Variables were pre-selected based on knowledge of
> pharmacokinetics of this drug and other factors. I think this part I
> understand pretty well."
> 
> This paragraph is wrong. Sorry, no way around it.
> 
> Quoting from the paper Pawinski etal:
> "  *__Twenty-six(!)* 1-, 2-, or 3-sample estimation
> models were fit (r2 
 0.341– 0.862) to a randomly
> selected subset of the profiles using linear regression
> and were used to estimate AUC0–12h for the profiles not
> included in the regression fit, comparing those estimates
> with the corresponding AUC0–12h values, calculated
> with the linear trapezoidal rule, including all 12
> timed MPA concentrations. The 3-sample models were
> constrained to include no samples past 2 h."
> (emph. mine)
> 
> They clearly state that they are choosing among 26 different models by
> using their bootstrap-like procedure, not improving on a single,
> predefined model.
> This procedure is statistically sound (more or less at least), and not
> controversial.
>

Re: [R] contour lines in windows device but neither in pdf nor in postscript

[Mark Difford]

Hi Patrizio,

>> # I can see contour lines in a window device but I can't see them in
>> files pdftry.pdf and pstry.ps

No problem with either format on my system, though I am using 2.7.1 patched.
It's not mentioned as a bug fix for the patched version so surely was
working in 2.7.1. Probably something local.

HTH, Mark. 

platform   i386-pc-mingw32
arch   i386   
os mingw32
system i386, mingw32  
status Patched
major  2  
minor  7.1
year   2008   
month  06 
day27 
svn rev46012  
language   R  
version.string R version 2.7.1 Patched (2008-06-27 r46012)


Patrizio Frederic wrote:
> 
> library(mvtnorm)
> x = seq(-4,4,length=201)
> xy= expand.grid(x,x)
> sigma = (diag(c(1,1))+1)/2
> d2= matrix(dmvnorm(xy,sigma=sigma),201)
> xsamp = rmvnorm(200,sigma=sigma)
> 
> contour(x,x,d2)
> points(xsamp,col=3,pch=16)
> 
> pdf("pdftry.pdf")
> contour(x,x,d2)
> points(xsamp,col=3,pch=16)
> dev.off()
> 
> postscript("pstry.ps")
> contour(x,x,d2)
> points(xsamp,col=3,pch=16)
> dev.off()
> 
> # I can see contour lines in a window device but I can't see them in
> files pdftry.pdf and pstry.ps
>> version
>_
> platform   i386-pc-mingw32
> arch   i386
> os mingw32
> system i386, mingw32
> status
> major  2
> minor  7.1
> year   2008
> month  06
> day23
> svn rev45970
> language   R
> version.string R version 2.7.1 (2008-06-23)
> 
> what's going wrong?
> 
> Thanks in advance,
> Regards,
> 
> Patrizio Frederic
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/contour-lines-in-windows-device-but-neither-in-pdf-nor-in-postscript-tp18771094p18771537.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] How to get the p-value from lmer on a longitudinal analysis

[Mark Difford]

Hi Ronaldo,

... lmer p-values

There are two packages that may help you with this and that might work with
the current implementation of lmer(). They are languageR and RLRsim.

HTH, Mark.


Bugzilla from [EMAIL PROTECTED] wrote:
> 
> Hi,
> 
> I have a modelo like this:
> 
> Yvar <- c(0, 0, 0, 0, 1, 0, 0, 0, 1, 2, 1, 1, 2, 3, 6, 6, 3, 3, 4)
> TIME <- 4:22
> ID <- rep("PlotA",19)
> m <- lmer(Yvar~TIME+(TIME|ID),family=poisson)
> anova(m)
> summary(m)
> 
> How to get the p-value for this case?
> 
> Thanks
> Ronaldo
> -- 
> Just because you're paranoid doesn't mean they AREN'T after you.
> --
>> Prof. Ronaldo Reis Júnior
> |  .''`. UNIMONTES/Depto. Biologia Geral/Lab. de Biologia Computacional
> | : :'  : Campus Universitário Prof. Darcy Ribeiro, Vila Mauricéia
> | `. `'` CP: 126, CEP: 39401-089, Montes Claros - MG - Brasil
> |   `- Fone: (38) 3229-8187 | [EMAIL PROTECTED] |
> [EMAIL PROTECTED]
> | http://www.ppgcb.unimontes.br/lbc | ICQ#: 5692561 | LinuxUser#: 205366
> --
> Favor NÃO ENVIAR arquivos do Word ou Powerpoint
> Prefira enviar em PDF, Texto, OpenOffice (ODF), HTML, or RTF.
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/How-to-get-the-p-value-from-lmer-on-a-longitudinal-analysis-tp18776696p18779630.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Lattice: Changing the names of conditional variables in strips to mathematical expressions

[Mark Difford]

Hi Andrew,

This does it as part of the call.  I have increased the height of the strip
and added italic for the second name only.

densityplot(~density|type,data=Query,plot.points="jitter",ref=TRUE,width="sj", 
panel=function(x, ...){ 
panel.grid(h=-1, v=-1) 
panel.densityplot(x, ...) },
par.strip.text=list(lines=1.5),## incr. line height
strip=strip.custom(factor.levels=c(expression((sqrt(Gower^{1}))),
   expression((sqrt(italic(Kulczynski)^{1}))

HTH, Mark.


Andrewjohnclose wrote:
> 
> Hi,
> 
> I know this has come up before, but I am having a hard time getting any of
> the solutions I have found to work!
> 
> I am trying to change the conditioning variable names in my xyplot from
> the default names as they appear from the data table - see following:
>  
> X=densityplot(~density|type,data=Query,plot.points="jitter",ref=TRUE,width="sj",
> panel=function(x, ...){
> panel.grid(h=-1, v=-1)
> panel.densityplot(x, ...)
> })
> 
> X
> 
> to 
> 
> dimnames(X)[[1]] <- expression((sqrt(Gower^{1})),(sqrt(Kulczynski))
> 
> X
> 
> The key issue is gaining the square root symbols and superscript position
> of the "1" ...
> 
> As an aside, it is possible to vary the font type within expression? Say
> setting the text  to italics but leaving the superscript "1" as ariel?
> 
> Suggestion would be greatly appreciated
> 
> Thank you
> 
> Andrew
> 
>  http://www.nabble.com/file/p18787745/Query.csv Query.csv 
> 

-- 
View this message in context: 
http://www.nabble.com/Lattice%3A-Changing-the-names-of-conditional-variables-in-strips-to-mathematical-expressions-tp18787745p18787895.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] How do I get an inset graph (i.e. graph within a graph)?

[Mark Difford]

Hi Arthur,

This can be done quite easily using the appropriate arguments listed under
?par; and there are other approaches. Ready-made functions exist in several
packages. I tend to use ?add.scatter from package ade4. It's a short
function, so it's easy to customize it, but it works well straight out of
the box.

HTH, Mark.


Arthur Roberts wrote:
> 
> Hi, all,
> 
> How do I get an inset graph (i.e. graph within a graph)?  Your input  
> is greatly appreciated.
> 
> Best wishes,
> Art
> University of Washington
> Department of Medicinal Chemistry
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/How-do-I-get-an-inset-graph-%28i.e.-graph-within-a-graph%29--tp18796563p18798809.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] subset inside a lattice plot using panel.lines

[Mark Difford]

Hi Michael,

>> Pulling my hair out here trying to get something very simple to work. ...

I can't quite see what you are trying to do [and I am not sure that you
clearly state it], but you could make things easier and simpler by (1)
creating a factor to identify your groups of rows more cleanly and (2) by
using the groups= argument in lattice.

## Something along these lines
copyDat <- originalDat
copyDat$newFac <- gl(2, 387, 774, labels=c("A","B"))

xyplot(response ~ predictor|maybe.condition, data=copyDat, groups=newFac)

I hope this gives you some ideas.

Mark.


Michael Hopkins wrote:
> 
> 
> 
> Hi R people [duplicate - sorry, just posted HTML by mistake]
> 
> Pulling my hair out here trying to get something very simple to work.   
> Have a data frame of 774 rows and want to plot first and second half  
> on same axes with different colours.  A variable is present call 'row'  
> created like this and checked to be OK:
> 
>   row <- seq( len=length( variable1 ) )
> 
> ...so I just want to plot the two subsets where row <= 387 and where  
> row > 387.  None of these work properly:
> 
> plots both over each other in correct colours
>   opt_plot2 <- function( var_string, c, units ) {
>   print( xyplot( get(var_string) ~ variable1 | variable2, panel =
>   function( x, y, ... ) {
>   panel.xyplot( x, y, type ="n", ... )
>   panel.grid( h=-1, v=-1, lty=3, lwd=1, col="grey" )
>   panel.lines( x, y, col = "red", subset = row <= 387 )
>   panel.lines( x, y, col = "dark green", subset = row > 387 )
>   },
>   }
> 
> plots both just in red
>   ... panel.lines( x[row <= 387], y[row <= 387], col = "red" )
>   panel.lines( x[row > 387], y[row > 387], col = "dark 
> green" )
> 
> first <- (row <= 387)
> second <- (row > 387)
> 
> plots both over each other in correct colours
>   ... panel.lines( x, y, col = "red", subset = first )
>   panel.lines( x, y, col = "dark green", subset = second )
> plots both just in red
>   ... panel.lines( x[first], y[first], col = "red" )
>   panel.lines( x[second], y[second], col = "dark green" )
> 
> 
> I'm feeling frustrated and a bit stupid but should this be so  
> difficult?  Any help or tips on what I am doing wrong greatly  
> appreciated.
> 
> TIA
> 
> Michael
> 
> __
> 
>  Hopkins Research  Touch the Future
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/subset-inside-a-lattice-plot-using-panel.lines-tp18813236p18813818.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Is there any way to make pretty tables in R to pdf?

[Mark Difford]

Hi Arthur,

>> I was wondering if there was a package that can make pretty R tables to
>> pdf.

You got through TeX/LateX, but PDF could be your terminus. Package Hmisc:

? summary.formula

and its various arguments and options. You can't get much better.

http://cran.za.r-project.org/doc/contrib/Harrell-statcomp-notes.pdf

HTH, Mark.


Arthur Roberts wrote:
> 
> Hi, all,
> 
> All your comments have been very useful.  I was wondering if there was  
> a package that can make pretty R tables to pdf.  I guess I could use  
> xtable, but I would like something a little more elegant.  Your input  
> is greatly appreciated.
> 
> Best wishes,
> Art
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Is-there-any-way-to-make-pretty-tables-in-R-to-pdf--tp18818187p18818675.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Is there any way to make pretty tables in R to pdf?

[Mark Difford]

Hi Arthur,

Sorry, sent you down the wrong track: this will help you to get there:

http://biostat.mc.vanderbilt.edu/twiki/pub/Main/StatReport/summary.pdf

Regards, Mark.


Arthur Roberts wrote:
> 
> Hi, all,
> 
> All your comments have been very useful.  I was wondering if there was  
> a package that can make pretty R tables to pdf.  I guess I could use  
> xtable, but I would like something a little more elegant.  Your input  
> is greatly appreciated.
> 
> Best wishes,
> Art
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Is-there-any-way-to-make-pretty-tables-in-R-to-pdf--tp18818187p18818837.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] re quest for fine panel axis controls in lattice

[Mark Difford]

Hi David,

>> Specifically, within each panel, I want to set the limits for x and y
>> equal to each other since it is paired data (using the max value of the
>> two).

In addition to the code Chuck Cleland sent you, you may want to "square"
things up by adding the argument: aspect = "iso" before the closing bracket

xyplot(..., aspect="iso")

Regards, Mark.


DavidChosid wrote:
> 
> I'm trying to use fine axis controls in lattice for each panel.
> Specifically, within each panel, I want to set the limits for x and y
> equal to each other since it is paired data (using the max value of the
> two).  Of course, I have no problems setting the limits for the entire
> plot but I am having trouble setting them for each specific panel.
> Could someone please provide me some guidance?  Thanks in advance.
> 
>  
> 
> David Chosid
> 
> Massachusetts Division of Marine Fisheries
> 
> 1213 Purchase St., 3rd Floor
> 
> New Bedford, MA 02740
> 
> 508-990-2860 x140
> 
> email: [EMAIL PROTECTED]
> 
>  
> 
>  
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/request-for-fine-panel-axis-controls-in-lattice-tp18830204p18831668.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Correlation dichotomous factor, continous (numerical) and ordered factor

[Mark Difford]

Hi Birgitle,

>> ... my variables are dichotomous factors, continuous (numerical) and
>> ordered factors. ...
>> Now I am confused what I should use to calculate the correlation using
>> all my variables
>> and how I could do that in R.

Professor Fox's package polycor will do this for you in a very nice way.

Regards, Mark.


Birgitle wrote:
> 
> Hello R-User!
> 
> I appologise in advance if this should also go into statistics but I am
> presently puzzled.
> I have a data.frame (about 300 rows and about 80 variables) and my
> variables are dichotomous factors, continuous (numerical) and ordered
> factors.
> 
> I would like to calculate the linear correlation between every pair of my
> variables, because I would like to perform a logistic regression (glm())
> without the correlation between variables.
> 
> I thought I could use for the continous (numerical) and ordered factor a
> spearman correlation that is using the ranks.
> 
> But I thought also that I have to use a contingency table for the
> dichotomous factors. 
> 
> I read also that it is possible to use a point-biserial correlation to
> calculate the correlation between dichotomous and continuous variables.
> 
> Now I am confused what I should use to calculate the correlation using all
> my variables and how I could do that in R.
> Is it possible with cor(), rcorr(), cormat() or other R-functions using
> one of the available correlation-coefficients.
> 
> I would be very happy if somebody could enlighten my darkness.
> 
> Many thanks in advance.
> 
> B.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Correlation-dichotomous-factor%2C-continous-%28numerical%29-and-ordered-factor-tp18852158p18852399.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] defining the distance between axis label and axis

[Mark Difford]

Hi Jörg,

>> I haven't found anything in par()...

No? Well don't bet your bottom $ on it (almost never in R). ?par (sub mgp).

Mark.


Jörg Groß wrote:
> 
> Hi,
> 
> How can I make the distance between an axis-label and the axis bigger?
> 
> I haven't found anything in par()...
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/defining-the-distance-between-axis-label-and-axis-tp18858861p18858935.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Problems using hetcor (polycor)

[Mark Difford]

Hi Birgitle,

You need to get this right if someone is going to spend their time helping
you. Your code doesn't work: You have specified more columns in colClasses
than you have in the provided data set.

> TestPart<-read.table("TestPart.txt", header=TRUE,row.names=1,
> na.strings="NA" ,colClasses = Classe81)
Error in scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings, 
: 
  line 1 did not have 82 elements

> TestPart<-read.table("TestPart.txt", header=TRUE,row.names=1,
> na.strings="NA")
Warning message:
In scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings,  :
  number of items read is not a multiple of the number of columns

> length(names(TestPart))
[1] 72

> length(Classe81)
[1] 82

That will never work.

HTH, Mark.


Birgitle wrote:
> 
> Sorry if this post should be long but I tried to give you a piece of my
> data to reproduce my error message using hetcor:
> 
> Fehler in result$rho : $ operator is invalid for atomic vectors
> Zusätzlich: Warning messages:
> 1: In polychor(x, y, ML = ML, std.err = std.err) :
>   1 row with zero marginal removed
> 2: In polychor(x, y, ML = ML, std.err = std.err) :
>   the table has fewer than 2 rows
> 
> 
> Error in result$rho : $ operator is invalid for atomic vectors
> Additional: Warning message:
> 1: In polychor(x, y, ML = ML, std.err = std.err) :
>   1 row with zero marginal removed
> 2: In polychor(x, y, ML = ML, std.err = std.err) :
>   the table has fewer than 2 rows
> 
> Use tab delimited data at the end of the post.
> Copy in Texteditor and save as TestPart.txt 
> 
> Then use the following code
> 
> library(methods)
> setClass("of")
> setAs("character", "of", function(from) as.ordered(from)) 
> 
> Classe81<-cclasses <- c(rep("factor", 64),  rep("numeric",6), rep ("of",
> 12))
> 
> TestPart<-read.table("TestPart.txt", header=TRUE,row.names=1,
> na.strings="NA" ,colClasses = Classe81)
> 
> str(TestPart)
> 
> library(polycor)
> 
> TestPart.hetcor<-hetcor(TestPart, use="complete.obs")
> 
> this will produce the above mentioned error message that I can not
> interprete.
> 
> Would be great if somebody could help me to solve the problem.
> 
> Thanks
> 
> B.
> 
> 
> 
> 
>   1   2   3   4   7   8   9   10  12  
> 13  14  15  16  17  18  19  21  22  23
>   25  27  28  29  30  31  33  34
> 3536  37  38  39  40  41  42  43  44  
> 45  46  47  48  49  50  51  52  53  54
>   55  56  58  59  60
> 6162  63  64  65  66  67  68  69  70  
> 71  72  73  74  75  76  77  78  79  80
> AX1   1   0   0   1   0   0   1   0   
> 0   0   0   0   0   0   1   0   1   0 
>   0   1   1   0   0   0   0   0   1   0   
> 0   0   1   0   0   0   0
> 0 1   0   0   0   0   0   0   0   0   
> 0   1   1   0   0   1   0   1   25  5 
>   9   1   8.5 2.5 3   5   2   2   3   3   
> 1   1   1   2   1
> 2
> BX1   1   0   0   1   0   0   1   NA  
> NA  NA  0   0   0   0   1   0   0   1 
>   0   0   1   0   NA  NA  NA  NA  NA  NA  
> NA
> NA0   0   0   1   0   NA  NA  NA  NA  
> NA  NA  NA  0   0   0   0   1   1   0 
>   0   0   1   1   NA  NA  6   1   3.25
> 2.25  5   5   2   2   3   3   1   1   1   
> 1   1   1
> CX1   1   0   0   1   0   0   1   1   
> 0   0   0   1   0   1   0   0   1   0 
>   1   0   0   0   0   1   1   0   0   0   
> 0   0   0   1   0   0   0
> 0 1   0   0   0   0   0   0   0   0   
> 1   0   0   0   0   1   0   0   15  3.5   
>   6   1   5.5 5.5 5   5   2   2   1   2   
> 1   1   1   1
> 2 2
> DX1   1   0   0   1   0   0   1   0   
> 0   0   0   0   0   0   1   0   0   1 
>   0   0   1   0   0   1   0   1   0   0   
> 0   1   0   0   0   0   1
> 1 0   0   0   0   0   0   0   0   0   
> 1   0   1   0   1   0   1   0   50  17.5  
>   7.5 2.5  

Re: [R] Problems using hetcor (polycor)

[Mark Difford]

Hi Birgitle,

It seems to be failing on those columns that have just a single "entry" (i.e
= 1, with the rest as 0; having just 1, an , and then 0s gets you
through). And there are other reasons for failure (in the call to get a
positive definite matrix).

The main problem lies in the calculation of standard errors for some
combinations. You can get a result for all columns by turning this off. You
can get standard errors for all the columns up to 60 by omitting columns 12,
23, and 44.  You need to work out the rest yourself by "columning" forwards
till you get a failure, then drop that column from the index. There probably
isn't enough information to calculate SEs for the columns that cause an
error.

## This works with default settings but without columns 12, 23 and 44
hetcor(TestPart[,c(1:11,13:22,24:43,45)], pd=T, std.err=T,
use="complete.obs")

## The first fails; the second works
hetcor(TestPart[,c(1:11,13:22,24:43,45:60)], pd=T, std.err=F)
hetcor(TestPart[,c(1:72)], pd=F, std.err=F)

HTH, Mark.


Birgitle wrote:
> 
> Sorry if this post should be long but I tried to give you a piece of my
> data to reproduce my error message using hetcor:
> 
> Fehler in result$rho : $ operator is invalid for atomic vectors
> Zusätzlich: Warning messages:
> 1: In polychor(x, y, ML = ML, std.err = std.err) :
>   1 row with zero marginal removed
> 2: In polychor(x, y, ML = ML, std.err = std.err) :
>   the table has fewer than 2 rows
> 
> 
> Error in result$rho : $ operator is invalid for atomic vectors
> Additional: Warning message:
> 1: In polychor(x, y, ML = ML, std.err = std.err) :
>   1 row with zero marginal removed
> 2: In polychor(x, y, ML = ML, std.err = std.err) :
>   the table has fewer than 2 rows
> 
> Use tab delimited data at the end of the post.
> Copy in Texteditor and save as TestPart.txt 
> 
> Then use the following code
> 
> library(methods)
> setClass("of")
> setAs("character", "of", function(from) as.ordered(from)) 
> 
> Classe81<-cclasses <- c(rep("factor", 64),  rep("numeric",6), rep ("of",
> 12))
> 
> TestPart<-read.table("TestPart.txt", header=TRUE,row.names=1,
> na.strings="NA" ,colClasses = Classe81)
> 
> str(TestPart)
> 
> library(polycor)
> 
> TestPart.hetcor<-hetcor(TestPart, use="complete.obs")
> 
> this will produce the above mentioned error message that I can not
> interprete.
> 
> Would be great if somebody could help me to solve the problem.
> 
> Thanks
> 
> B.
> 
> 
> 
> 
>   1   2   3   4   7   8   9   10  12  
> 13  14  15  16  17  18  19  21  22  23
>   25  27  28  29  30  31  33  34
> 3536  37  38  39  40  41  42  43  44  
> 45  46  47  48  49  50  51  52  53  54
>   55  56  58  59  60
> 6162  63  64  65  66  67  68  69  70  
> 71  72  73  74  75  76  77  78  79  80
> AX1   1   0   0   1   0   0   1   0   
> 0   0   0   0   0   0   1   0   1   0 
>   0   1   1   0   0   0   0   0   1   0   
> 0   0   1   0   0   0   0
> 0 1   0   0   0   0   0   0   0   0   
> 0   1   1   0   0   1   0   1   25  5 
>   9   1   8.5 2.5 3   5   2   2   3   3   
> 1   1   1   2   1
> 2
> BX1   1   0   0   1   0   0   1   NA  
> NA  NA  0   0   0   0   1   0   0   1 
>   0   0   1   0   NA  NA  NA  NA  NA  NA  
> NA
> NA0   0   0   1   0   NA  NA  NA  NA  
> NA  NA  NA  0   0   0   0   1   1   0 
>   0   0   1   1   NA  NA  6   1   3.25
> 2.25  5   5   2   2   3   3   1   1   1   
> 1   1   1
> CX1   1   0   0   1   0   0   1   1   
> 0   0   0   1   0   1   0   0   1   0 
>   1   0   0   0   0   1   1   0   0   0   
> 0   0   0   1   0   0   0
> 0 1   0   0   0   0   0   0   0   0   
> 1   0   0   0   0   1   0   0   15  3.5   
>   6   1   5.5 5.5 5   5   2   2   1   2   
> 1   1   1   1
> 2 2
> DX1   1   0   0   1   0   0   1   0   
> 0   0   0   0   0   0   1   0   0   1 
>   0   0   1   0   0   1   0 

Re: [R] Problems using hetcor (polycor)

[Mark Difford]

Hi Birgitle,

>> I am so sorry there was a little mistake in the code again.

This arrived _after_ I sent in my "fix" (really just an exploration of what
the problem might be). Perhaps you need to wrestle a bit more with the art
of coding ...;). I don't think it will change the work-around. You won't
have standard errors for everything, but you will get them for just about
everything.

Regards, Mark.


Birgitle wrote:
> 
> I am so sorry there was a little mistake in the code again.
> 
> Right code now:
> 
> library(methods)
> setClass("of")
> setAs("character", "of", function(from) as.ordered(from))
> 
> Classe72<-cclasses <- c(rep("factor", 55),  rep("numeric",6), rep ("of",
> 12))
> 
> TestPart<-read.table("TestPart.txt", header=TRUE,row.names=1,
> na.strings="NA" ,colClasses = Classe72)
> 
> library(polycor)
> 
> TestPart.hetcor<-hetcor(TestPart, use="complete.obs")
> 
> B.
> 
> 
> Birgitle wrote:
>> 
>> Thanks Mark and I am sorry that I forgot to adapt the Classe-vector.
>> 
>> This should work now
>> 
>> library(methods)
>> setClass("of")
>> setAs("character", "of", function(from) as.ordered(from))
>> 
>> Classe72<-cclasses <- c(rep("factor", 55),  rep("numeric",6), rep
>> ("factor", 12))
>> 
>> TestPart<-read.table("TestPart.txt", header=TRUE,row.names=1,
>> na.strings="NA" ,colClasses = Classe72)
>> 
>> library(polycor)
>> 
>> TestPart.hetcor<-hetcor(TestPart, use="complete.obs")
>> 
>> Mark Difford wrote:
>>> 
>>> Hi Birgitle,
>>> 
>>> You need to get this right if someone is going to spend their time
>>> helping you. Your code doesn't work: You have specified more columns in
>>> colClasses than you have in the provided data set.
>>> 
>>>> TestPart<-read.table("TestPart.txt", header=TRUE,row.names=1,
>>>> na.strings="NA" ,colClasses = Classe81)
>>> Error in scan(file, what, nmax, sep, dec, quote, skip, nlines,
>>> na.strings,  : 
>>>   line 1 did not have 82 elements
>>> 
>>>> TestPart<-read.table("TestPart.txt", header=TRUE,row.names=1,
>>>> na.strings="NA")
>>> Warning message:
>>> In scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings,  :
>>>   number of items read is not a multiple of the number of columns
>>> 
>>>> length(names(TestPart))
>>> [1] 72
>>> 
>>>> length(Classe81)
>>> [1] 82
>>> 
>>> That will never work.
>>> 
>>> HTH, Mark.
>>> 
>>> 
>>> Birgitle wrote:
>>>> 
>>>> Sorry if this post should be long but I tried to give you a piece of my
>>>> data to reproduce my error message using hetcor:
>>>> 
>>>> Fehler in result$rho : $ operator is invalid for atomic vectors
>>>> Zusätzlich: Warning messages:
>>>> 1: In polychor(x, y, ML = ML, std.err = std.err) :
>>>>   1 row with zero marginal removed
>>>> 2: In polychor(x, y, ML = ML, std.err = std.err) :
>>>>   the table has fewer than 2 rows
>>>> 
>>>> 
>>>> Error in result$rho : $ operator is invalid for atomic vectors
>>>> Additional: Warning message:
>>>> 1: In polychor(x, y, ML = ML, std.err = std.err) :
>>>>   1 row with zero marginal removed
>>>> 2: In polychor(x, y, ML = ML, std.err = std.err) :
>>>>   the table has fewer than 2 rows
>>>> 
>>>> Use tab delimited data at the end of the post.
>>>> Copy in Texteditor and save as TestPart.txt 
>>>> 
>>>> Then use the following code
>>>> 
>>>> library(methods)
>>>> setClass("of")
>>>> setAs("character", "of", function(from) as.ordered(from)) 
>>>> 
>>>> Classe81<-cclasses <- c(rep("factor", 64),  rep("numeric",6), rep
>>>> ("of", 12))
>>>> 
>>>> TestPart<-read.table("TestPart.txt", header=TRUE,row.names=1,
>>>> na.strings="NA" ,colClasses = Classe81)
>>>> 
>>>> str(TestPart)
>>>> 
>>>> library(polycor)
>>>> 
>>>> TestPart.hetcor<-hetcor(TestPart, use="complete.obs

Re: [R] Problems using hetcor (polycor)

[Mark Difford]

Hi Birgitle,

>> It seems than, that it is not possible to use all variables without
>> somehow 
>> imputing missing values.

It depends on what you are after. You can use the full data set if you set
std.err=F and pd=F. Then exclude the columns that cause it to falter and
redo with SEs turned on. You have the correlations; all you lack are SEs for
ca. 5 columns.

And I haven't tried your revised classification; that's for you to do.

Regards, Mark.


Birgitle wrote:
> 
> Many Thanks Mark for your answer.
> 
> It seems than, that it is not possible to use all variables without
> somehow imputing missing values.
> 
> But I will try which variables I can finally use.
> 
> Many thanks again.
> 
> B.
> 
> 
> Mark Difford wrote:
>> 
>> Hi Birgitle,
>> 
>> It seems to be failing on those columns that have just a single "entry"
>> (i.e = 1, with the rest as 0; having just 1, an , and then 0s gets
>> you through). And there are other reasons for failure (in the call to get
>> a positive definite matrix).
>> 
>> The main problem lies in the calculation of standard errors for some
>> combinations. You can get a result for all columns by turning this off.
>> You can get standard errors for all the columns up to 60 by omitting
>> columns 12, 23, and 44.  You need to work out the rest yourself by
>> "columning" forwards till you get a failure, then drop that column from
>> the index. There probably isn't enough information to calculate SEs for
>> the columns that cause an error.
>> 
>> ## This works with default settings but without columns 12, 23 and 44
>> hetcor(TestPart[,c(1:11,13:22,24:43,45)], pd=T, std.err=T,
>> use="complete.obs")
>> 
>> ## The first fails; the second works
>> hetcor(TestPart[,c(1:11,13:22,24:43,45:60)], pd=T, std.err=F)
>> hetcor(TestPart[,c(1:72)], pd=F, std.err=F)
>> 
>> HTH, Mark.
>> 
>> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Problems-using-hetcor-%28polycor%29-tp18867343p18870999.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Where is the archive? Calling a function within a function.

[Mark Difford]

Hi Kevin,

>> Where is the archive?

Start with this:

?RSiteSearch

HTH, Mark.


rkevinburton wrote:
> 
> I seem to remember this topic coming up before so I decided to look at the
> archive and realized that I didn't know where it was. Is there a
> searchable archive for this list? Thank you.
> 
> My question is calling a function from within a function. I have
> 
> smerge <- function(d1,d2) {
> temp <- merge(d1,d2,all=TRUE,by.x="DayOfYear",by.y="DayOfYear")
> return (xmerge(temp))
> }
> 
> But I get an error message indicating that xmerge could not be found when
> it is defined just above.
> 
> Thank you.
> 
> Kevin
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Where-is-the-archive--Calling-a-function-within-a-function.-tp18874953p18875170.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Tick marks that correspond with bars on barplot

[Mark Difford]

Hi Megan,

>> I would like to have an X-axis where the labels for the years line up
>> after every two bars
>> in the plot (there is one bar for hardwood, and another for softwood).

It isn't clear to me from your description what you really want (I found no
attachment)? What you seem to be trying to get is a tickmark/label at the
centre/midpoint of each yearly grouping. If so, then look at the examples
that the author of barplot() has provided, and _in particular_ read what the
author of the method says under "Value". You use the colMeans of the object
returned by barplot() to position labels.

Regards, Mark.


Megan J Bellamy wrote:
> 
> Hello all,
> 
> I have created a barplot that shows change in hardwood/softwood density
> from 1965 to 2005 in 5 year periods (1965,1970, etc). I would like to have
> an X-axis where the labels for the years line up after every two bars in
> the plot (there is one bar for hardwood, and another for softwood). Below
> is my script:
> 
> density<-read.table("F:\\Megan\\Vtest.csv", header=TRUE, sep=",")
> attach (density)
> barplot(DENSITY,YEAR, col=c("blue", "green", "green", "blue", "blue",
> "green", "blue", "green", "green", "blue", "green", "blue", "blue",
> "green" , "green", "blue", "blue", "green"))
> legend(1,85,c("Softwood", "Hardwood"), fill=c("blue", "green"), bty="n")
> title(ylab="Density of Softwood and Hardwood by percent (%)")
> title(xlab="Year of measurement")
> title(main="Change in Softwood and Hardwood Density between 1965-2005 for
> PSP #80")
> axis(2, at=NULL)
> 
> I have tried using the tck, axis.lty functions with no luck and have also
> tried
> axis(1, at=YEAR) # but only the first year (1965) comes up. 
> 
> Attached is the csv file. Any help with this would be greatly appreciated.
> Many thanks in advance,
> 
> Megan
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Tick-marks-that-correspond-with-bars-on-barplot-tp18890762p18892131.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] [lme4]Coef output with binomial lmer

[Mark Difford]

Hi Tom,

>> 1|ass%in%pop%in%fam

This is "non-standard," but as you have found, it works. The correct
translation is in fact

1|fam/pop/ass

and not 1|ass/pop/fam as suggested by Harold Doran. Dropping %,
ass%in%pop%in%fam reads [means] as: nest ass in pop [= pop/ass], and then
nest this in fam == fam/pop/ass

HTH, Mark.


T.C. Cameron wrote:
> 
> Dear R users  
>  
> I have built the following model
>  
> m1<-lmer(y~harn+foodn+(1|ass%in%pop%in%fam),family = "quasibinomial")
>  
> where y<-cbind(alive,dead)
>  
> where harn and foodn are categorical factors and the random effect is a
> nested term to represent experimental structure
> e.g.  Day/Block/Replicate
> ass= 5 level factor, pop= 2 populations per treatment factor in each
> assay, 7 reps per population
>  
> The model can be family = quasibinomial or binomial
>  
> My complete lack of understanding is in retrieving the coefficients for
> the fixed effects to back-transform the effects of my factors on
> proportional survival
>  
> I get the following output:
>> coef(m1)
> $`ass %in% pop %in% fam`
>   (Intercept)  harn1 harn2   foodn2
> FALSE   1.0322375 -0.1939521 0.0310434 0.810084
> TRUE0.5997679 -0.1939521 0.0310434 0.810084
>  
> Where FALSE and TRUE refer to some attribute of the random effect 
>  
> My hunch is that it refers to the Coefficients with (=TRUE) and without
> (=FALSE) the random effects?
>  
> Any help appreciated
>  
> 
> 
> 
> Dr Tom C Cameron
> Genetics, Ecology and Evolution
> IICB, University of Leeds
> Leeds, UK
> Office: +44 (0)113 343 2837
> Lab:+44 (0)113 343 2854
> Fax:+44 (0)113 343 2835
> 
> 
> Email: [EMAIL PROTECTED]
> Webpage: click here
>  
> 
>  
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/-lme4-Coef-output-with-binomial-lmer-tp18894407p18904468.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Simple lme/lmer random effects questions

[Mark Difford]

Hi Brandon,

>> ...is it sufficient to leave the values as they are or should I generate
>> unique names for all 
>> combinations of sleeve number and temperature, using something like
>> > data$sleeve.in.temp <- factor(with(data, temp:sleeve)[drop=TRUE])

You might be luckier posting this on

https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

You need to make a unique identifier for each of your five sets of sleeves
(or you need to nest sleeve in temp) if you wish to take account of that
structure (temp/sleeve) as a random effect. I could be mistaken, but I think
there are subtle differences between ~ 1|temp/sleeve and ~
1|data$sleeve.in.temp designs.

>> confused on how to actually set up the random effects term for the 
>> models. Given my experimental setup, using the lme syntax...

The difference between (1) random = ~ 1|sleeve, ... and (2) random = ~
1+temp|sleeve is that the first will give a random intercept for each level
of sleeve, but the slope is fixed to be the same as their is no random term
for slope. In the second specification there is a random term for slope, viz
temp, so both intercepts and slopes can --- and probably will --- vary.

## To get some idea of what's being done, look at the following example. 

## random intercepts, parallel slopes 
mod1 <- lme(distance ~ age, Orthodont, random = ~ 1 | Subject)  
## random intercepts, separate slopes 
mod2 <- lme(distance ~ age, Orthodont, random = ~ age | Subject)  

plot(augPred(mod1, primary=~age))  ## parallel slopes 
plot(augPred(mod2, primary=~age))  ## separate slopes

HTH, Mark.


Brandon Invergo wrote:
> 
> Hello,
> 
> I have two very rudimentary questions regarding the random effects terms 
> in the lme and lmer functions. I apologize if this also partly strays 
> into a general statistics question, but I'm a bit new to this all. So 
> hopefully it'll be a quick problem to sort out...
> 
> Here is my experimental setup: I raised butterflies in 5 different 
> testing chambers all set to different temperatures. Within the testing 
> chambers, the butterflies were held in 10 different sleeves, which were 
> rotated daily to compensate for microenvironmental effects. I measured 
> several traits of the butterflies and I am constructing models for each 
> trait (unfortunately, multivariate analysis isn't possible). In my 
> models, sex and temperature are fixed factors and the sleeve is a random 
> effect. Most of the response variables are normally distributed, but 
> there is one with a Gamma distribution (time until an event) and another 
> with poisson distribution (counts), so some models use lme while others 
> use lmer. I would like to determine if, despite the daily rotation, 
> there are still random effects from the individual sleeves. My two 
> questions (assuming I haven't already made grave errors in my 
> description of the setup) are:
> 
> 1) In my data file, the "sleeve" variable is just marked with a number 1 
> through 10; the temperature is noted in a different column, so the 50 
> sleeves do not have unique names, but rather there are 5 instances of 
> each of the 10 sleeve numbers. If sleeve is to be properly included in 
> the models as a random effect, is it sufficient to leave the values as 
> they are or should I generate unique names for all combinations of 
> sleeve number and temperature, using something like
>  > data$sleeve.in.temp <- factor(with(data, temp:sleeve)[drop=TRUE])
> 
> 
> 2) (this is the one that strays more into standard statistics territory, 
> sorry) I'm a bit confused on how to actually set up the random effects 
> term for the models. Given my experimental setup, using the lme syntax, 
> should it be:
>  > model <- lme(response ~ sex*temp, random=~temp|sleeve, data)
> or
>  > model <- lme(response ~ sex*temp, random=~1|sleeve, data)
> or something else? I've searched and searched, but everything I find 
> online seems to be significantly more advanced than what I'm doing, 
> leaving me even more confused than when I started!
> 
> 
> Thank you very much for your help!! I want to be sure I do this analysis 
> right
> Cheers,
> -brandon
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Simple-lme-lmer-random-effects-questions-tp18933698p18939376.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] ylab with an exponent

[Mark Difford]

>> what is the problem?

A solution is:

plot(1,2, ylab=expression(paste("insects ", m^2)))

The problem is very much more difficult to determine.


stephen sefick wrote:
> 
> plot(1,2, ylab= paste("insects", expression(m^2), sep=" "))
> 
> I get insects m^2
> I would like m to the 2
> 
> what is the problem?
> 
> -- 
> Let's not spend our time and resources thinking about things that are
> so little or so large that all they really do for us is puff us up and
> make us feel like gods. We are mammals, and have not exhausted the
> annoying little problems of being mammals.
> 
>   -K. Mullis
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/ylab-with-an-exponent-tp19003927p19004197.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] before-after control-impact analysis with R

[Mark Difford]

Hi Nikolaos,

>> My question again is: Why can't I reproduce the results? When I try a 
>> simple anova without any random factors:

Lack of a "right" result probably has to do with the type of analysis of
variance that is being done. The default in R is to use so-called Type I
tests, for good reason. SAS, I think, still uses a type of test that many
authorities consider to be meaningless, as a general, first-level, test.

There is a long, long thread on this, going back to about approx April/May
1999, when a suitable Ripplyism was coined, which I believe found its way
into the fortunes package. But

RSiteSearch("type III")

should do for a start. Also see

?drop1
library(car)
?Anova

HTH, Mark.


Nikolaos Lampadariou wrote:
> 
> Hello everybody,
> 
> In am trying to analyse a BACI experiment and I really want to do it 
> with R (which I find really exciting).  So, before moving on I though it 
> would be a good idea to repeat some known experiments which are quite 
> similar to my own. I tried to reproduce 2 published examples but without 
> much success. The first one in particular is a published dataset 
> analysed with SAS by McDonald et al. (2000). They also provide the 
> original data as well as the SAS code. I don't know much about SAS and i 
> really want to stick to R. So here follow 3 questions based on these 2 
> papers.
> 
> 
> Paper 1
> (McDonald et al. 2000. Analysis of count data from before-after 
> control-impact studies. J. Agric. Biol. Envir Stat., 5:262-279)
> 
> Data were collected from 1984, 1989, 1990, 1991, 1993, and 1996. Samples 
> from 1984 are considered as before an impact and the remaining years as 
>   after the impact. Each year, 96 transects were sampled (36 in the 
> oiled area and 60 in the non-oiled area; "0" is for oiled and "1" for 
> non-oiled). The authors compare 3 different ways of analysing the data 
> including glmm. The data can be reproduced with the following commands 
> (density is fake numbers but I can provide the original data since I've 
> typed them in anyway):
> 
>  >x<-c(rep("0",36),rep("1",60))
>  >oiled<-rep(x,6)
>  >year<-c(rep(1984,96), rep(1989,96), rep(1990,96), rep(1991,96), 
> rep(1993,96), rep(1996,96))
>  >density<-runif(576, min=0, max=10)
>  >transect<-c(rep(1:96,6))
>  >oil<-data.frame("oiled"=oiled, "transect"=transect, "year"=year, 
> "density"=density)
> 
> Question 1:
> I can reproduce the results of the repeated measures anova with:
>  >oil.aov1<-aov(density~factor(year)*factor(oiled)+Error(factor(transect))
> 
> But why is the following command not working?
>  >oil.aov2<-aov(density~oiled*year + Error(oiled/transect), data=oil)
> 
> After reading the R-help archive, as well as Chambers and Hasties 
> (Statistical Models in S) and Pinheiro's and Bates (Mixed effects models 
> in S and S-plus) I would expect that the correct model is the oil.aov2. 
> As you might see from the data, transect is nested within oiled and I 
> still get the wrong results when I try +Error(transect) or 
> +Error(oiled:transect) and many other variants.
> 
> Question 2 (on the same paper):
> The authors conclude that it is better to analyse the data with a 
> Generalized Linear (Mixed) Model Technique. I tried lme and after 
> reading Douglas Bates article (May 2005, vol. 5/1) also lmer as follows:
>  >oil.lme<-lme(density~year*oiled, random=~1|oiled/transect)
> or
>  >harvest.lmer<-lmer(abund~H*BA+BA/year+BA/year:H+site:BA+(1|H/site))
> but again no luck. I will get always some error messages or some 
> interactions missing.
> 
> 
> Paper 2
> (Keough & Quinn, 2000. Legislative vs. practical protection of an 
> intertidal shoreline in southeastern Australia. Ecol. Appl. 10: 871-881)
> 
> Data were collected from 1989, 1990, 1991, 1993, 1994, 1995, 1996, 
> 1997). the 1989-1991 are the before impact years and the other 4 years 
> are the after the impact years. Eight sites were samples (2 protected 
> and 6 impacted). In this dataset, site is nested within harvest (H) and 
> year is nested within before-after (BA). Also, site is considered as 
> random by the authors. The data (fake again) can be produced with the 
> following commands:
> 
>  >site<-c(rep(c("A1","A2", "RR1", "RR2", "WT1", "WT2", "WT3", "WT4"),8))
>  >H<-c(rep(c("exp", "exp", "prot", "pro", "exp", "exp", "exp", "exp"), 8))
>  >year<-c(rep(1989,8), rep(1990,8), rep(1991,8), rep(1993,8), 
> rep(1994,8), rep(1995,8), rep(1996,8), rep(1997,8))
>  >BA<-c(rep("bef",24), rep("after",40))
>  >abund<-runif(64, min=0, max=10)
>  >harvest<-data.frame(abund, BA, H, site, year)
> 
> Question 3.
> The authors use a complex anova design and here is part of their anova 
> table which shows the design and the df.
> 
> source.of.var  df
> Harvesting(H) 1, 6
> before-after(BA)  1, 6
> H x BA1, 6
> Site(H)   6, 31
> Year(BA)  6, 31
> Site x BA 6, 31
> Year x H  6, 31
> Resid.31
> 
> 
> My question again is: Why can't I reproduce the results? When I 

Re: [R] before-after control-impact analysis with R

[Mark Difford]

Hi ...

Sorry, an "e" was erroneously "elided" from Ripley...



Mark Difford wrote:
> 
> Hi Nikolaos,
> 
>>> My question again is: Why can't I reproduce the results? When I try a 
>>> simple anova without any random factors:
> 
> Lack of a "right" result probably has to do with the type of analysis of
> variance that is being done. The default in R is to use so-called Type I
> tests, for good reason. SAS, I think, still uses a type of test that many
> authorities consider to be meaningless, as a general, first-level, test.
> 
> There is a long, long thread on this, going back to about approx April/May
> 1999, when a suitable Ripplyism was coined, which I believe found its way
> into the fortunes package. But
> 
> RSiteSearch("type III")
> 
> should do for a start. Also see
> 
> ?drop1
> library(car)
> ?Anova
> 
> HTH, Mark.
> 
> 
> Nikolaos Lampadariou wrote:
>> 
>> Hello everybody,
>> 
>> In am trying to analyse a BACI experiment and I really want to do it 
>> with R (which I find really exciting).  So, before moving on I though it 
>> would be a good idea to repeat some known experiments which are quite 
>> similar to my own. I tried to reproduce 2 published examples but without 
>> much success. The first one in particular is a published dataset 
>> analysed with SAS by McDonald et al. (2000). They also provide the 
>> original data as well as the SAS code. I don't know much about SAS and i 
>> really want to stick to R. So here follow 3 questions based on these 2 
>> papers.
>> 
>> 
>> Paper 1
>> (McDonald et al. 2000. Analysis of count data from before-after 
>> control-impact studies. J. Agric. Biol. Envir Stat., 5:262-279)
>> 
>> Data were collected from 1984, 1989, 1990, 1991, 1993, and 1996. Samples 
>> from 1984 are considered as before an impact and the remaining years as 
>>   after the impact. Each year, 96 transects were sampled (36 in the 
>> oiled area and 60 in the non-oiled area; "0" is for oiled and "1" for 
>> non-oiled). The authors compare 3 different ways of analysing the data 
>> including glmm. The data can be reproduced with the following commands 
>> (density is fake numbers but I can provide the original data since I've 
>> typed them in anyway):
>> 
>>  >x<-c(rep("0",36),rep("1",60))
>>  >oiled<-rep(x,6)
>>  >year<-c(rep(1984,96), rep(1989,96), rep(1990,96), rep(1991,96), 
>> rep(1993,96), rep(1996,96))
>>  >density<-runif(576, min=0, max=10)
>>  >transect<-c(rep(1:96,6))
>>  >oil<-data.frame("oiled"=oiled, "transect"=transect, "year"=year, 
>> "density"=density)
>> 
>> Question 1:
>> I can reproduce the results of the repeated measures anova with:
>> 
>> >oil.aov1<-aov(density~factor(year)*factor(oiled)+Error(factor(transect))
>> 
>> But why is the following command not working?
>>  >oil.aov2<-aov(density~oiled*year + Error(oiled/transect), data=oil)
>> 
>> After reading the R-help archive, as well as Chambers and Hasties 
>> (Statistical Models in S) and Pinheiro's and Bates (Mixed effects models 
>> in S and S-plus) I would expect that the correct model is the oil.aov2. 
>> As you might see from the data, transect is nested within oiled and I 
>> still get the wrong results when I try +Error(transect) or 
>> +Error(oiled:transect) and many other variants.
>> 
>> Question 2 (on the same paper):
>> The authors conclude that it is better to analyse the data with a 
>> Generalized Linear (Mixed) Model Technique. I tried lme and after 
>> reading Douglas Bates article (May 2005, vol. 5/1) also lmer as follows:
>>  >oil.lme<-lme(density~year*oiled, random=~1|oiled/transect)
>> or
>>  >harvest.lmer<-lmer(abund~H*BA+BA/year+BA/year:H+site:BA+(1|H/site))
>> but again no luck. I will get always some error messages or some 
>> interactions missing.
>> 
>> 
>> Paper 2
>> (Keough & Quinn, 2000. Legislative vs. practical protection of an 
>> intertidal shoreline in southeastern Australia. Ecol. Appl. 10: 871-881)
>> 
>> Data were collected from 1989, 1990, 1991, 1993, 1994, 1995, 1996, 
>> 1997). the 1989-1991 are the before impact years and the other 4 years 
>> are the after the impact years. Eight sites were samples (2 protected 
>> and 6 impacted). In this dataset, site is nested within harvest (H) and 
>> year is nested within before-after (BA). Also, site is considered

Re: [R] Re lational Operators in Legend

[Mark Difford]

Hi Lorenzo,

>> ...but I would like to write that 5<=k<=15.

This is one way to do what you want

plot(1,1)
legend("topright", expression(paste(R[g]~k^{1/d[f]^{small}}~5<=k, {}<=15)))

HTH, Mark.


Lorenzo Isella wrote:
> 
> Dear All,
> I am sure that what I am asking can be solved by less than a
> one-liner, but I am having a hard time to get this right.
> Inside a legend environment, I do not have any problem if I issue the
> command:
> expression(""*R[g]*"~"* k^{1/d[f]^{small}}*",  "*k<= {} *"15")
> but I would like to write that 5<=k<=15. It has to be trivial, but so
> far attempts like:
> expression(""*R[g]*"~"* k^{1/d[f]^{small}}*", for  5 "*<={}*" "*k<= {}
> *"15")
> always return some error message and no improvement with all the
> variations I tried out.
> Does anyone know what I am doing wrong?
> Many Thanks
> 
> Lorenzo
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Relational-Operators-in-Legend-tp19032571p19033251.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Re lational Operators in Legend

[Mark Difford]

Hi Lorenzo,

I may (?) have left something out. It isn't clear what "~" is supposed to
mean; perhaps it is just a spacer, or perhaps you meant the following:

plot(1,1)
legend("topright", expression(paste(R[g] %~~% k^{1/d[f]^{small}},~5<=k,
{}<=15)))

HTH, Mark.


Mark Difford wrote:
> 
> Hi Lorenzo,
> 
>>> ...but I would like to write that 5<=k<=15.
> 
> This is one way to do what you want
> 
> plot(1,1)
> legend("topright", expression(paste(R[g]~k^{1/d[f]^{small}}~5<=k,
> {}<=15)))
> 
> HTH, Mark.
> 
> 
> Lorenzo Isella wrote:
>> 
>> Dear All,
>> I am sure that what I am asking can be solved by less than a
>> one-liner, but I am having a hard time to get this right.
>> Inside a legend environment, I do not have any problem if I issue the
>> command:
>> expression(""*R[g]*"~"* k^{1/d[f]^{small}}*",  "*k<= {} *"15")
>> but I would like to write that 5<=k<=15. It has to be trivial, but so
>> far attempts like:
>> expression(""*R[g]*"~"* k^{1/d[f]^{small}}*", for  5 "*<={}*" "*k<= {}
>> *"15")
>> always return some error message and no improvement with all the
>> variations I tried out.
>> Does anyone know what I am doing wrong?
>> Many Thanks
>> 
>> Lorenzo
>> 
>> __
>> R-help@r-project.org 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.
>> 
>> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Relational-Operators-in-Legend-tp19032571p19033652.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] How to solve this problem ?

[Mark Difford]

Hi Daren,

>> Small progress, ...

m4 <- list(m1=m1, m2=m2, m3=m3)
boxplot(m4)

It's always a good idea to have a look at your data first (assuming you
haven't). This shows that the reliable instrument is m2.

HTH, Mark.


Daren Tan wrote:
> 
> 
> Small progress, I am relying on levene test to check for equality of
> variances. Is my understanding correct, the larger the p-value, the more
> likely the variances are the same ?
> 
>> trt
>  [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4
> Levels: 1 2 3 4
>> levene.test(rep(rnorm(5), 4), trt, option="median")
> 
> Modified Robust Brown-Forsythe Levene-type test based on the
> absolute
> deviations from the median
> 
> data:  rep(rnorm(5), 4)
> Test Statistic = 0, p-value = 1
> 
> 
> 
>> From: [EMAIL PROTECTED]
>> To: [EMAIL PROTECTED]
>> Subject: How to solve this problem ?
>> Date: Wed, 20 Aug 2008 13:33:23 +
>>
>>
>> I have disabled html text editing mode, thanks to Prof. Ripley for the
>> kind reminder.
>>
>> Given three geological devices that takes 5 readings at 4 environmental
>> conditions (A to D). What will be the proper approach to select the most
>> reliable device ?
>>
>> m1 <- c(73,52,49,53,83,43,58,94,53,62,75,66,41,72,70,75,57,59,85,84)
>> m2 <- c(31,38,30,35,36,26,27,38,22,31,24,35,36,31,38,33,32,28,33,30)
>> m3 <- c(65,57,36,40,36,30,40,34,37,40,33,33,37,29,37,37,30,33,40,35)
>>
>> names(m1) <- rep(LETTERS[1:4], each=5)
>> names(m2) <- rep(LETTERS[1:4], each=5)
>> names(m3) <- rep(LETTERS[1:4], each=5)
>>
>> Before writing this email, I have tried to compare the sd for each device
>> at each condition, but ran into obstacle on how to formulate the null
>> hypothesis. Alternative solution tried was ANOVA, I am unsure whether it
>> can help, as it compares the differences in means of each group.
>>
>> Thanks
>>
>> _
>> Easily edit your photos like a pro with Photo Gallery.
>> http://get.live.com/photogallery/overview
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/How-to-solve-this-problem---tp19069553p19084239.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Test of Homogeneity of Variances

[Mark Difford]

Have you read the documentation to either of the functions you are using?

?bartlett.test

"Performs Bartlett's test of the null that the variances in each of the
groups (samples) are the same."

This explicitly tells you what is being tested, i.e. the null tested is that
var1 = var2.

?rnorm

Generates random deviates, so the answer [i.e. p-value] will (almost
certainly) differ each time. Only by setting a seed for the random number
generator to use will rnorm() generate the same number-set/distribution and
therefore the same p-value.

?set.seed

##
set.seed(7)
bartlett.test(rep(rnorm(5),times=4), rep(1:5, each=4))
set.seed(7)
bartlett.test(rep(rnorm(5),times=4), rep(1:5, each=4))
set(23)
bartlett.test(rep(rnorm(5),times=4), rep(1:5, each=4))
bartlett.test(rep(rnorm(5),times=4), rep(1:5, each=4))

HTH, Mark.


Daren Tan wrote:
> 
> 
> I am testing whether the sample variances are equal. When p-value < 0.05
> (alpha), should accept null hypothesis (sample variances are equal) or
> reject it ?
> 
> 
> The two new examples with each having same sample variances also puzzle
> me. Why are the p-values different ?
> 
> bartlett.test(rep(rnorm(5),times=4), rep(1:5, each=4))
> 
> 
> Bartlett test of homogeneity of variances
> 
> 
> data:  rep(rnorm(5), times = 4) and rep(1:5, each = 4)
> Bartlett's K-squared = 0.8681, df = 4, p-value = 0.929
> 
> 
> bartlett.test(rep(rnorm(5),times=4), rep(1:5, each=4))
> 
> 
> Bartlett test of homogeneity of variances
> 
> 
> data:  rep(rnorm(5), times = 4) and rep(1:5, each = 4)
> Bartlett's K-squared = 3.5599, df = 4, p-value = 0.4688
> 
> 
>> From: [EMAIL PROTECTED]
>> To: [EMAIL PROTECTED]; [EMAIL PROTECTED]
>> Date: Fri, 22 Aug 2008 11:25:36 -0400
>> Subject: RE: [R] Test of Homogeneity of Variances
>>
>> What are your hypotheses? Once you state what they are, interpretation
>> should be straightforward.
>>
>>
>>
>> -Original Message-
>> From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
>> On Behalf Of Daren Tan
>> Sent: Friday, August 22, 2008 11:18 AM
>> To: [EMAIL PROTECTED]
>> Subject: [R] Test of Homogeneity of Variances
>>
>>
>> I am testing the homogeneity of variances via bartlett.test and
>> fligner.test. Using the following example, how should I interpret the
>> p-value in order to accept or reject the null hypothesis ?
>>
>> set.seed(5)
>> x <- rnorm(20)
>> bartlett.test(x, rep(1:5, each=4))
>>
>>
>> Bartlett test of homogeneity of variances
>>
>>
>> data: x and rep(1:5, each = 4)
>> Bartlett's K-squared = 1.7709, df = 4, p-value = 0.7778
>>
>> fligner.test(x, rep(1:5, each=4))
>>
>> Fligner-Killeen test of homogeneity of variances
>>
>>
>> data: x and rep(1:5, each = 4)
>> Fligner-Killeen:med chi-squared = 1.0819, df = 4, p-value = 0.8971
>>
>> __
>> R-help@r-project.org 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.
>> This email message, including any attachments, is for ...{{dropped:6}}
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Test-of-Homogeneity-of-Variances-tp19109383p19112613.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] aov, lme, multcomp

[Mark Difford]

Hi Richard,

>> The tests give different Fs and ps. I know this comes up every once in a 
>> while on R-help so I did my homework. I see from these two threads:

This is not so, or it is not necessarily so. The error structure of your two
models is quite different, and this is (one reason) why the F- and p-values
are different.

For instance, try the following comparison:

## Example 
require(MASS) ## for oats data set 
require(nlme) ## for lme() 
require(multcomp)  ## for multiple comparison stuff 

Aov.mod <- aov(Y ~ N + V + Error(B/V), data = oats) 
Lme.mod <- lme(Y ~ N + V, random = ~1 | B/V, data = oats) 

summary(Aov.mod) 
anova(Lme.mod) 


See:
http://www.nabble.com/Tukey-HSD-(or-other-post-hoc-tests)-following-repeated-measures-ANOVA-td17508294.html#a17553029

The example itself is from MASS (Venables & Ripley).

HTH, Mark.


Richard D. Morey wrote:
> 
> I am doing an analysis and would like to use lme() and the multcomp 
> package to do multiple comparisons. My design is a within subjects 
> design with three crossed fixed factors (every participant sees every 
> combination of three fixed factors A,B,C). Of course, I can use aov() to 
> analyze this with an error term (leaving out the obvious bits):
> 
> y ~ A*B*C+Error(Subject/(A*B*C))
> 
> I'd also like to use lme(), and so I use
> 
> y ~ A*B*C, random= ~1|Subject
> 
> The tests give different Fs and ps. I know this comes up every once in a 
> while on R-help so I did my homework. I see from these two threads:
> 
> http://www.biostat.wustl.edu/archives/html/s-news/2002-05/msg00095.html
> http://134.148.236.121/R/help/06/08/32763.html
> 
> that this is the expected behavior because of the way grouping works 
> with lme(). My questions are:
> 
> 1. is this the correct random argument to lmer:
> 
>   anova(lme(Acc~A*B*C,random=list(Sub=pdBlocked(list(
>   pdIdent(~1),
>   pdIdent(~A-1),
>   pdIdent(~B-1),
>   pdIdent(~C-1,data=data))
> 
> 2. How much do the multiple comparisons depend on the random statement?
> 
> 3. I'm also playing with lmer:
> 
> Acc~A*B*C+(1|Sub)
> 
> Is this the correct lmer call for the crossed factors? If not, can you 
> point me towards the right one?
> 
> 4. I'm not too concerned with getting "correct" Fs from the analyses 
> (well, except for aov, where it is easy), I just want to make sure that 
> I am fitting the same model to the data with all approaches, so that 
> when I look at parameter estimates I know they are meaningful. Are the 
> multiple comparisons I'll get out of lme and lmer meaningful with fully 
> crossed factors, given that they are both "tuned" for nested factors?
> 
> Thanks in advance.
> 
> -- 
> Richard D. Morey
> Assistant Professor
> Psychometrics and Statistics
> Rijksuniversiteit Groningen / University of Groningen
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/aov%2C-lme%2C-multcomp-tp19144362p19145003.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Two envelopes problem

[Mark Difford]

...

To pick up on what Mark has said: it strikes me that this is related to the
simplex, where the bounded nature of the vector space means that normal
arithmetical operations (i.e. Euclidean) don't work---that is, they can be
used, but the results are wrong. Covariances and correlations for instance,
are wrong, something that Pearson noted more than a century ago.

Taking logs of ratios of numbers (say a number divdided by geometric mean of
the other numbers, as in Aitchison's set of transformations) in this space
transfers everything to Euclidean space, so "squaring" the problem. This is
why taking logs fixes things ??

Mark.



statquant wrote:
> 
> Duncan: I think I see what you're saying but the strange thing is that if
> you use the utility function log(x) rather than x, then the expected
> values
> are equal. Somehow, if you are correct and I think you are, then taking
> the
> log , "fixes" the distribution of x which is kind of odd to me. I'm sorry
> to
> belabor this non R related discussion and I won't say anything more about
> it
> but I worked/talked  on this with someone for about a month a few years
> ago
> and we gave up so it's interesting for me to see this again.
> 
>Mark
> 
> -Original Message-
> From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
> On
> Behalf Of Duncan Murdoch
> Sent: Tuesday, August 26, 2008 8:15 AM
> To: Jim Lemon
> Cc: r-help@r-project.org; Mario
> Subject: Re: [R] Two envelopes problem
> 
> On 26/08/2008 7:54 AM, Jim Lemon wrote:
>> Hi again,
>> Oops, I meant the expected value of the swap is:
>> 
>> 5*0.5 + 20*0.5 = 12.5
>> 
>> Too late, must get to bed.
> 
> But that is still wrong.  You want a conditional expectation, 
> conditional on the observed value (10 in this case).  The answer depends 
> on the distribution of the amount X, where the envelopes contain X and 
> 2X.  For example, if you knew that X was at most 5, you would know you 
> had just observed 2X, and switching would be  a bad idea.
> 
> The paradox arises because people want to put a nonsensical Unif(0, 
> infinity) distribution on X.  The Wikipedia article points out that it 
> can also arise in cases where the distribution on X has infinite mean: 
> a mathematically valid but still nonsensical possibility.
> 
> Duncan Murdoch
> 
> __
> R-help@r-project.org 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-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Two-envelopes-problem-tp19150357p19163195.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] random error with lme for repeated measures anova

[Mark Difford]

Hi Jean-Pierre,

A general comment is that I think you need to think more carefully about
what you are trying to get out of your analysis. The random effects
structure you are aiming for could be stretching your data a little thin.

It might be a good idea to read through the archives of the
R-sig-mixed-models mailing list

https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

This should give you a better idea of some of the issues involved.

If you can't source the source (P&B) then there are other documents that
might help. You could begin by locationing the directory where nlme package
is installed and looking at the scripts in the scripts subdirectory. These
are from P&B.

Baayen has the following draft document on his website. Though DRAFT'ed all
over the place, it is well worth a read, even if you are not interested in
linguistic analysis. I think it has now been published by Cambridge UP.

http://www.ualberta.ca/~baayen/publications/baayenCUPstats.pdf

Paul Bliese's document (Multilevel Modeling in R) also has some useful
sections (find it under the contributed documents section on CRAN).

HTH, Mark.


Jean-Pierre Bresciani wrote:
> 
> Hi,
> 
> what is the appropriate syntax to get the random error correct when
> performing repeated measures anova with 'lme'.
> 
> let's say i have 3 independent variables, with 'aov', i would write
> something like: aov(dep_var~(indep_var1*indep_var2*indep_var3) +
> Error(subject/(indep_var1*indep_var2*indep_var3)).
> 
> With 'lme' however, i can't find the right formula. i tried things like:
> lme(dep_var~(indep_var1*indep_var2*indep_var3), random = ~1|subject) or
> nesting my independent variables in 'subject', but those are obviously
> wrong with my design.
> 
> i'm quite clueless (and i haven't found any convincing piece of
> information about how to correctly use 'lme' or 'lmer'). So, any advice
> along that line is more than welcome.
> 
> JP 
> 

-- 
View this message in context: 
http://www.nabble.com/random-error-with-lme-for-repeated-measures-anova-tp19178239p19180027.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] convert princomp output to equation for plane?

[Mark Difford]

Hi Bill,

>> Since x, y,and z all have measurement errors attached, the proper way 
>> to do the fit is with principal components analysis, and to use the 
>> first component (called loadings in princomp output).

The easiest way for you to do this is to use the pcr [principal component
regression] function in the pls package. Be aware that unless you fit all
components you will be carrying out a form of penalized regression. A small
example follows (assumes that you have installed the pls package):

##
lm.mod <- lm(Ozone ~ Solar.R + Wind + Month, data=airquality)
pc.mod <- pcr(Ozone ~ Solar.R + Wind + Month, data=airquality)

lm.mod
coef(pc.mod, intercept = TRUE)
coef(pc.mod, ncomp=1, intercept = TRUE)
coef(pc.mod, ncomp=3, intercept = TRUE)

Regards, Mark.


William Simpson-2 wrote:
> 
> I want to fit something like:
> z = b0 + b1*x + b2*y
> 
> Since x, y,and z all have measurement errors attached, the proper way
> to do the fit is with principal components analysis, and to use the
> first component (called loadings in princomp output).
> 
> My dumb question is: how do I convert the princomp output to equation
> coefficients in the format above?
> 
> I guess another dumb question would be: how about getting the standard
> deviations of b0, b1, b2?
> 
> Thanks very much for any help.
> 
> Bill
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/convert-princomp-output-to-equation-for-plane--tp19182643p19197360.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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 II regression - how do I do it?

[Mark Difford]

Hi Danilo,

>> I need to do a model II linear regression, but I could not find out how!!

The smatr package does so-called model II (major axis) regression.

Regards, Mark.


Danilo Muniz wrote:
> 
> I need to do a model II linear regression, but I could not find out how!!
> 
> I tryed to use the lm function, but I did not discovered how to specify
> the
> model (type I or type II) to the function... could you help me?
> 
> -- 
> Danilo Muniz
> [Gruingas Abdiel]
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/model-II-regression---how-do-I-do-it--tp19224427p19225640.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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 II regression - how do I do it?

[Mark Difford]

Hi Danilo,

>> I need to do a model II linear regression, but I could not find out how!!

The smatr package does so-called model II (major axis) regression.

Regards, Mark.


Danilo Muniz wrote:
> 
> I need to do a model II linear regression, but I could not find out how!!
> 
> I tryed to use the lm function, but I did not discovered how to specify
> the
> model (type I or type II) to the function... could you help me?
> 
> -- 
> Danilo Muniz
> [Gruingas Abdiel]
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/model-II-regression---how-do-I-do-it--tp19224427p19225663.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Density estimates in modelling framework

[Mark Difford]

Hi Hadley,

There is also locfit, which is very highly regarded by some authorities
(e.g. Hastie, Tibs, and Friedman).

Cheers, Mark.


hadley wrote:
> 
> Hi all,
> 
> Do any packages implement density estimation in a modelling framework?
>  I want to be able to do something like:
> 
> dmodel <- density(~ a + b, data = mydata)
> predict(dmodel, newdata)
> 
> This isn't how sm or KernSmooth or base density estimation works.  Are
> there other packages that do density estimation?  Or is there some
> reason that this is a bad idea.
> 
> Hadley
> 
> -- 
> http://had.co.nz/
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Density-estimates-in-modelling-framework-tp19225727p19226399.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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 II regression - how do I do it?

[Mark Difford]

Hi Dylan,

>> While this topic is fresh, are there any compelling reasons to use Model
>> II 
>> regression?

The fact that it is the type of regression used in principal component
analysis makes it a compelling method. Compelling reason? It is used to take
account of measurement errors in both y _and_ x. Usually
practioners/analysts would sidestep this by putting what is measured with
least error on the x-axis and regress y on that. For instance, in doing
calibration curves in nutrient analysis where the analyte is measured using
a spec. Then the spec reading goes on x. I can't tell you how often I have
seen analysts put the (usually) inaccurately determined analyte on x and the
spec reading on y.

HTH, Mark.


Dylan Beaudette-2 wrote:
> 
> On Friday 29 August 2008, Mark Difford wrote:
>> Hi Danilo,
>>
>> >> I need to do a model II linear regression, but I could not find out
>> >> how!!
>>
>> The smatr package does so-called model II (major axis) regression.
>>
>> Regards, Mark.
> 
> While this topic is fresh, are there any compelling reasons to use Model
> II 
> regression?
> 
> Cheers,
> 
> Dylan
> 
> -- 
> Dylan Beaudette
> Soil Resource Laboratory
> http://casoilresource.lawr.ucdavis.edu/
> University of California at Davis
> 530.754.7341
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/model-II-regression---how-do-I-do-it--tp19224427p19226644.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] non-parametric Anova and tukeyHSD

[Mark Difford]

Hi Stephen,

See packages:

coin
nparcomp
npmc

There is also kruskalmc() in package pgirmess

Regards, Mark.


stephen sefick wrote:
> 
> I have insect data from twelve sites and like most environmental data
> it is non-normal mostly.  I would like to preform an anova and a means
> seperation like tukey's HSD in a nonparametric sense (on some sort of
> central tendency measure - median?).  I am searching around at this
> time on the internet.  Any suggestions, books, etc. would be greatly
> appreciated.
> 
> -- 
> Stephen Sefick
> Research Scientist
> Southeastern Natural Sciences Academy
> 
> Let's not spend our time and resources thinking about things that are
> so little or so large that all they really do for us is puff us up and
> make us feel like gods. We are mammals, and have not exhausted the
> annoying little problems of being mammals.
> 
>   -K. Mullis
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/non-parametric-Anova-and-tukeyHSD-tp19227579p19227986.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] non-parametric Anova and tukeyHSD

[Mark Difford]

That slipped right away from me before I could say  I meant to add this
link to a useful thread from Torsten Hothorn.

?RSiteSearch would probably have got you there.

http://tolstoy.newcastle.edu.au/R/help/05/06/5829.html


Mark Difford wrote:
> 
> Hi Stephen,
> 
> See packages:
> 
> coin
> nparcomp
> npmc
> 
> There is also kruskalmc() in package pgirmess
> 
> Regards, Mark.
> 
> 
> stephen sefick wrote:
>> 
>> I have insect data from twelve sites and like most environmental data
>> it is non-normal mostly.  I would like to preform an anova and a means
>> seperation like tukey's HSD in a nonparametric sense (on some sort of
>> central tendency measure - median?).  I am searching around at this
>> time on the internet.  Any suggestions, books, etc. would be greatly
>> appreciated.
>> 
>> -- 
>> Stephen Sefick
>> Research Scientist
>> Southeastern Natural Sciences Academy
>> 
>> Let's not spend our time and resources thinking about things that are
>> so little or so large that all they really do for us is puff us up and
>> make us feel like gods. We are mammals, and have not exhausted the
>> annoying little problems of being mammals.
>> 
>>  -K. Mullis
>> 
>> __
>> R-help@r-project.org 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.
>> 
>> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/non-parametric-Anova-and-tukeyHSD-tp19227579p19228111.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] ANCOVA/glm missing/ignored interaction combinations

[Mark Difford]

Hi Lara,

>> And I cant for the life of me work out why category one (semio1) is being
>> ignored, missing 
>> etc.

Nothing is being ignored Lara --- but you are ignoring the fact that your
factors have been coded using the default contrasts in R, viz so-called
treatment or Dunnett contrasts. That is, your intercept is semio1 and the
other coefficients are calculated off it, hence the name treatment
contrasts.

In general this type of coding is recommended for GLM models (but this can
raise deep divisions, as some don't consider them to be proper contrasts
because they are not orthogonal). Perhaps you should read up on how dummy
variables are coded.

To really understand this you will have to move away from ?contrasts as a
source. It might also help you to do the following on a simpler model:

?dummy.coef
dummy.coef(saved.glm/lm.model)

This shows how it is worked out.

HTH, Mark.


lara harrup (IAH-P) wrote:
> 
> Hi
> 
> I am using R version 2.7.2. on a windows XP OS and have a question
> concerning an analysis of covariance with count data I am trying to do,
> I will give details of a scaled down version of the analysis (as I have
> more covariates and need to take account of over-dispersion etc etc) but
> as I am sure it is only a simple problem but I just can't see how to fix
> it.
> 
> I have a data set with count data as the response (total) and two
> continuous covariates and one categorical explanatory variable (semio).
> When I run the following lines, after loading the data and assigning
> 'semio' as a factor, taking into account that I want to consider two way
> interactions:
> 
>> model.a<-glm(total~(temperature+humidity+semio)^2,family=poisson)
>> summary(model.a)
> 
> I get the output below. But not all interactions are listed there are 4
> semio categories, 1,2,3 and 4 but only 2,3,and 4 are listed in the
> summary (semio2,semio3,semio4). And I cant for the life of me work out
> why category one (semio1) is being ignored, missing etc.
> 
> Any help or suggestions would be most appreciated. Thanks in advance
> 
> Lara
> [EMAIL PROTECTED]
> 
> Call:
> glm(formula = total ~ (temperature + humidity + semio)^2, family =
> poisson)
> 
> Deviance Residuals: 
> Min   1Q   Median   3Q  Max  
> -22.212   -5.132   -2.4843.200   18.793  
> 
> Coefficients:
>   Estimate Std. Error z value Pr(>|z|)
> (Intercept)  23.848754   2.621291   9.098  < 2e-16 ***
> temperature  -1.038284   0.150465  -6.901 5.18e-12 ***
> humidity -0.264912   0.033928  -7.808 5.81e-15 ***
> semio2   22.852664   1.291806  17.690  < 2e-16 ***
> semio33.699901   1.349007   2.743   0.0061 ** 
> semio4   -1.851163   1.585997  -1.167   0.2431
> temperature:humidity  0.012983   0.001983   6.545 5.94e-11 ***
> temperature:semio2   -0.870430   0.037602 -23.149  < 2e-16 ***
> temperature:semio3   -0.060846   0.038677  -1.573   0.1157
> temperature:semio40.055288   0.046137   1.198   0.2308
> humidity:semio2  -0.114718   0.013369  -8.581  < 2e-16 ***
> humidity:semio3  -0.031692   0.013794  -2.298   0.0216 *  
> humidity:semio4   0.008425   0.016020   0.526   0.5990
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> 
> (Dispersion parameter for poisson family taken to be 1)
> 
> Null deviance: 10423.5  on 47  degrees of freedom Residual deviance:
> 2902.2  on 35  degrees of freedom
> AIC: 3086.4
> 
> Number of Fisher Scoring iterations: 7
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/ANCOVA-glm-missing-ignored-interaction-combinations-tp19285259p19286859.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] ANCOVA/glm missing/ignored interaction combinations

[Mark Difford]

And perhaps I should also have added: fit your model without an intercept and
look at your coefficients. You should be able to work it out from there
quite easily. Anyway, you now have the main pieces.

Regards, Mark.


Mark Difford wrote:
> 
> Hi Lara,
> 
>>> And I cant for the life of me work out why category one (semio1) is
>>> being ignored, missing 
>>> etc.
> 
> Nothing is being ignored Lara --- but you are ignoring the fact that your
> factors have been coded using the default contrasts in R, viz so-called
> treatment or Dunnett contrasts. That is, your intercept is semio1 and the
> other coefficients are calculated off it, hence the name treatment
> contrasts.
> 
> In general this type of coding is recommended for GLM models (but this can
> raise deep divisions, as some don't consider them to be proper contrasts
> because they are not orthogonal). Perhaps you should read up on how dummy
> variables are coded.
> 
> To really understand this you will have to move away from ?contrasts as a
> source. It might also help you to do the following on a simpler model:
> 
> ?dummy.coef
> dummy.coef(saved.glm/lm.model)
> 
> This shows how it is worked out.
> 
> HTH, Mark.
> 
> 
> lara harrup (IAH-P) wrote:
>> 
>> Hi
>> 
>> I am using R version 2.7.2. on a windows XP OS and have a question
>> concerning an analysis of covariance with count data I am trying to do,
>> I will give details of a scaled down version of the analysis (as I have
>> more covariates and need to take account of over-dispersion etc etc) but
>> as I am sure it is only a simple problem but I just can't see how to fix
>> it.
>> 
>> I have a data set with count data as the response (total) and two
>> continuous covariates and one categorical explanatory variable (semio).
>> When I run the following lines, after loading the data and assigning
>> 'semio' as a factor, taking into account that I want to consider two way
>> interactions:
>> 
>>> model.a<-glm(total~(temperature+humidity+semio)^2,family=poisson)
>>> summary(model.a)
>> 
>> I get the output below. But not all interactions are listed there are 4
>> semio categories, 1,2,3 and 4 but only 2,3,and 4 are listed in the
>> summary (semio2,semio3,semio4). And I cant for the life of me work out
>> why category one (semio1) is being ignored, missing etc.
>> 
>> Any help or suggestions would be most appreciated. Thanks in advance
>> 
>> Lara
>> [EMAIL PROTECTED]
>> 
>> Call:
>> glm(formula = total ~ (temperature + humidity + semio)^2, family =
>> poisson)
>> 
>> Deviance Residuals: 
>> Min   1Q   Median   3Q  Max  
>> -22.212   -5.132   -2.4843.200   18.793  
>> 
>> Coefficients:
>>   Estimate Std. Error z value Pr(>|z|)
>> (Intercept)  23.848754   2.621291   9.098  < 2e-16 ***
>> temperature  -1.038284   0.150465  -6.901 5.18e-12 ***
>> humidity -0.264912   0.033928  -7.808 5.81e-15 ***
>> semio2   22.852664   1.291806  17.690  < 2e-16 ***
>> semio33.699901   1.349007   2.743   0.0061 ** 
>> semio4   -1.851163   1.585997  -1.167   0.2431
>> temperature:humidity  0.012983   0.001983   6.545 5.94e-11 ***
>> temperature:semio2   -0.870430   0.037602 -23.149  < 2e-16 ***
>> temperature:semio3   -0.060846   0.038677  -1.573   0.1157
>> temperature:semio40.055288   0.046137   1.198   0.2308
>> humidity:semio2  -0.114718   0.013369  -8.581  < 2e-16 ***
>> humidity:semio3  -0.031692   0.013794  -2.298   0.0216 *  
>> humidity:semio4   0.008425   0.016020   0.526   0.5990
>> ---
>> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
>> 
>> (Dispersion parameter for poisson family taken to be 1)
>> 
>> Null deviance: 10423.5  on 47  degrees of freedom Residual deviance:
>> 2902.2  on 35  degrees of freedom
>> AIC: 3086.4
>> 
>> Number of Fisher Scoring iterations: 7
>> 
>> __
>> R-help@r-project.org 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.
>> 
>> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/ANCOVA-glm-missing-ignored-interaction-combinations-tp19285259p19286930.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Help understanding lrm function of Design library

[Mark Difford]

Hi Noah,

>> Could someone point me to an online resource where I could learn 
>> more?  (I'm big on trying to teach myself.) [about lrm funcction and the
>> Design library]

Go for the hardcopy, but you could look at Google Books if you are pressed.
There you will find a good preview of the text. Also look at Alzolla and
Harrell's excellent "An Introduction to S and the Hmisc and Design
Libraries" PDF, which you will find on CRAN under "Contributed
Documentation."

Regards, Mark.


Noah Silverman-3 wrote:
> 
> Hi,
> 
> I'm developing an experiment with logistic regression.
> 
> I've come across the lrm function in the Design  library.
> 
> While I understand and can use the basic functionality, there are a ton 
> of options that go beyond my knowledge.
> 
> I've carefully read the help page for lrm, but don't understand many of 
> the arguments and optional return values.
> (penalty, penalty.matrix, strata.penalty, var.penalty, weights)
> 
> Additionally, there are several items returned in the model that I don't 
> understand.
> (Max Derriv, Model L.R, C, Dxy, Gamma, Tau-a, R2, Briar, Wald Z)
> 
> 1) Could someone point me to an online resource where I could learn 
> more?  (I'm big on trying to teach myself.)
> 
> 2) Or, alternately, does anybody want to explain a few things for me?
> 
> Thanks!
> 
> --
> Noah
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Help-understanding-lrm-function-of-Design-library-tp25022388p25023794.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Contrasts within ANOVA frame (Repost)

[Mark Difford]

Hi Jun,

>> I have three levels for a factor names "StdLot" and want to make three
>> comparisons, 1 vs 2, 1 vs 3 and 2 vs 3.

With only three levels to your factor, the contrast matrix you are
specifying is over-parametrized (i.e. over-specified): it has 3 rows and 3
columns.

## Look at the default behaviour
contrasts(d3$StdLot) <- NULL
contrasts(d3$StdLot)

contrasts(d3$StdLot) <- contr.helmert
contrasts(d3$StdLot)

I think what you really want are linear contrasts, which you will find in
package car (?linear.hypothesis), in gmodels (?estimable) [based on Frank
Harrell's approach; so look at his package Design], and package multcomp,
amongst others.

Regards, Mark.


Jun Shen-3 wrote:
> 
> Would like to try my luck to see if I can catch your eyes.
> 
> I was trying to do some contrasts within ANOVA. I searched the archive and
> found a clue posted by Steffen Katzner
> ( http://tolstoy.newcastle.edu.au/R/help/06/01/19385.html)
> 
> I have three levels for a factor names "StdLot" and want to make three
> comparisons, 1 vs 2, 1 vs 3 and 2 vs 3.
> 
> First,
> contrasts(d3$StdLot,3)<-matrix(c(1,-1, 0,0,1,-1,1,0,-1),3,3)   #d3 is the
> data set. set up the contrast matrix
> 
> Second,
> aov(Bkg~StdLot,na.rm=T,data=d3,contrasts=contrasts(d3$StdLot))->mod.aov
> #ANOVA,
> 
> Finally,
> summary(mod.aov,split=list(StdLot=list('1 vs 2'=1,'2 vs 3'=2,'1 vs 3'=3)))
> #comparison summary
> 
> Here is the final result I got, the third comparison is missing. Does
> anyone
> have any idea what is wrong here? If I change the order of the comparisons
> it's always the third one missing.  So I guess it's not due to the data.
> Appreciate any comment.
> 
>   Df Sum Sq Mean Sq F valuePr(>F)
> StdLot 2  1.905   0.953 10.3769 3.710e-05 ***
>   StdLot: 1 vs 2   1  0.223   0.223  2.42390.1200
>   StdLot: 2 vs 3   1  1.683   1.683 18.3299 2.162e-05 ***
>   StdLot: 1 vs 3   1
> Residuals601 55.173   0.092
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 2 observations deleted due to missingness
> 
> 
> -- 
> Jun Shen PhD
> PK/PD Scientist
> BioPharma Services
> Millipore Corporation
> 15 Research Park Dr.
> St Charles, MO 63304
> Direct: 636-720-1589
> 
>   [[alternative HTML version deleted]]
> 
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Contrasts-within-ANOVA-frame-%28Repost%29-tp25052361p25057955.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] applying summary() to an object created with ols()

[Mark Difford]

Hi Benjamin,

>> Does anyone know how I can set the *datadist()* and the *options()* such 
>> that I will get access to all coefficients?

## Do this before you fit your models, i.e. tell datadist &c what data set
you are using.
d <- datadist( subset(aa, Jahr>=1957 & Jahr<=1966) )
options( datadist = "d" )
d

You can specify the limits yourself when you call plot.Design or
summary.Design, but this automates the process.

Regards, Mark.


Benjamin Volland wrote:
> 
> Hello R-list,
> 
> I am trying to calculate a ridge regression using first the *lm.ridge()* 
> function from the MASS package and then applying the obtained Hoerl 
> Kennard Baldwin (HKB) estimator as a penalty scalar to the *ols()* 
> function provided by Frank Harrell in his Design package.
> It looks like this:
>  > rrk1<-lm.ridge(lnbcpc ~ lntex + lnbeerp + lnwinep + lntemp + pop, 
> subset(aa, Jahr>=1957 & Jahr<=1966))
>  > f <- ols(lnbcpc ~ lntex + lnbeerp + lnwinep + lntemp + pop, 
> subset(aa, Jahr>=1957 & Jahr<=1966), penalty = rrk$kHKB)
>  > f
> 
> which returns
>  >Linear Regression Model
>  >
>  >ols(formula = lnbcpc ~ lntex + lnbeerp + lnwinep + lntemp + pop,
>  >data = subset(aa, Jahr >= 1957 & Jahr <= 1966), penalty = rrk$kHKB)
>  >
>  > n Model L.R.   d.f. R2  Sigma
>  >10  38.59  8.814 0.98390.02796
>  >
>  >Residuals:
>  >1 2 3 4 5 6 
> 7 8 910
>  >-0.014653 -0.002787  0.017515 -0.018145 -0.008757 -0.008035  0.006066  
> 0.045826 -0.001244 -0.015786
>  >
>  >Coefficients:
>  >Value Std. Error   t Pr(>|t|)
>  >Intercept  1.5240 3.3034  0.4613   0.8496
>  >lntex  0.3722 0.2071  1.7975   0.6801
>  >lnbeerp0.9085 0.5760  1.5771   0.6964
>  >lnwinep   -0.1458 0.1874 -0.7781   0.7863
>  >lntemp-0.0772 0.1344 -0.5743   0.8240
>  >pop   -4.1889 1.9286 -2.1720   0.6571
>  >
>  >Adjusted R-Squared: 0.2227
> 
> All in all beautiful (leaving aside that the results suck). The problem 
> starts when I want to write the obtained coefficients (incl. Std. 
> Errors, t-, and p-values) into a matrix.
> Via the *f$coef* command I can only access the betas (1st column) and 
> using the *summary(f)* function I get
>  > summary(f)
>  >Fehler in summary.Design(f) :
>  >adjustment values not defined here or with datadist for lntex lnbeerp 
> lnwinep lntemp pop
> 
> Does anyone know how I can set the *datadist()* and the *options()* such 
> that I will get access to all coefficients?
> 
> I tried:
>  > options(datadist=NULL)
>  > f <- ols(lnbcpc ~ lntex + lnbeerp + lnwinep + lntemp + pop, 
> subset(aa, Jahr>=1957 & Jahr<=1966), penalty = rrk$kHKB)
>  > d <- datadist(f)
> but got:
>  > Fehler in sort.list(unique(y)) : 'x' must be atomic for 'sort.list'
>  > Have you called 'sort' on a list?
> 
> In the R documentation on ?ols() it states concerning the values 
> returned: "the same objects returned from |lm| (/unless |penalty| or 
> |penalty.matrix| are given/ - then an abbreviated list is returned since 
> |lm.pfit| is used as a fitter)..." Unfortunately no information seems to 
> be available on lm.pfit.
> Does anyone know why the using that function leads to an abbreviated 
> return list? Is there a trick to circumvent that?
> 
> Thanks
> Benjamin Volland
> 
> P.S. Currently using R-version 2.7.1 on a Windows PC.
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/applying-summary%28%29-to-an-object-created-with-ols%28%29-tp25083729p25091255.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] using loglog link in VGAM or creating loglog link for GLM

[Mark Difford]

Hi Kendra,

>> I am trying to figure out how to apply a loglog link to a binomial  
>> model (dichotomous response variable with far more zeros than ones).

If I were you I would look at ?zeroinfl in package pscl.

Regards, Mark.


Kendra Walker wrote:
> 
> 
> 
> I am trying to figure out how to apply a loglog link to a binomial  
> model (dichotomous response variable with far more zeros than ones).   
> I am aware that there are several relevant posts on this list, but I  
> am afraid I need a little more help.  The two suggested approaches  
> seem to be: 1) modify the make.link function in GLM, or 2) use the  
> loglog or cloglog functions in the VGAM package.  Below are my  
> questions for each.  Responses to either are much appreciated.
> 
> 1) Modifying the make.link function:
>   In his Sat, 24 Jan 2009 post, Re: [R] glm binomial loglog (NOT  
> cloglog) link, William Simpson suggested a loglog={} insertion that  
> seems reasonable.   Being new to R, however, I am going in circles  
> trying to figure out make this seemingly simple modification to the  
> code.   If anyone has the patience to step me through the process, or  
> refer me to the relevant information, I would be very grateful.
> 
> 2) Using the loglog function in VGAM:
> I tried fitting the model:  m<- vglm(Y~X,  
> family=binomialff(link="loglog"), data = d)  but I get the following  
> error? Error in lm.fit(X_vlm, z_vlm, ...) : NA/NaN/Inf in foreign  
> function call (arg 4) In addition: Warning message:In log(log(theta))  
> : NaNs produced?.
> If I run the same model using cloglog as the link, I get a result that  
> looks like the result I get using cloglog in GLM.  This is a bad fit  
> for my data, however, as I have many more zeros than ones.   The help  
> document for loglog states that NaNs are produced when theta is close  
> to 1 unless earg is used.  I am confused as to how to properly use the  
> earg parameter (and why I do not need it for cloglog despite having  
> many zeros), leading me to wonder whether the loglog link here is  
> really what I think it is (the compliment of cloglog).  Again, any  
> insights as to what I am missing would be appreciated.
> 
> Many thanks,
> 
> 
> Kendra Walker,
> School of Natural Resources and the Environment,
> University of Michigan
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/using-loglog-link-in-VGAM-or-creating-loglog-link-for-GLM-tp25091476p25096521.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Between-group variance from ANOVA

[Mark Difford]

Hi Emma,

>> 

R gives you the tools to work this out.

## Example
set.seed(7)
TDat <- data.frame(response = c(rnorm(100, 5, 2), rnorm(100, 20, 2)))
TDat$group <- gl(2, 100, labels=c("A","B"))
with(TDat, boxplot(split(response, group)))
summary(aov(response ~ group, data=TDat))

Regards, Mark.


emj83 wrote:
> 
> can anyone advise me please?
> 
> 
> emj83 wrote:
>> 
>> I have done some ANOVA tables for some data that I have, from this I can
>> read the within-group variance. can anyone tell me how i may find out the
>> between-group variance?
>> 
>> Thanks Emma
>> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Between-group-variance-from-ANOVA-tp24954045p25121532.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Between-group variance from ANOVA

[Mark Difford]

Hi Emma,

>> ...from this I can read the within-group variance. can anyone tell me how
>> i may find 
>> out the between-group variance?

But it's in the table, above the "within-group" variance. Remember that F is
the ratio of these two quantities, i.e. the mean of the group variances
divided by the mean of the within-group variances . I will work with my
example since you never set seed so your answers are different from mine
(which really does not help matters).

set.seed(7) 
TDat <- data.frame(response = c(rnorm(100, 5, 2), rnorm(100, 20, 2))) 
TDat$group <- gl(2, 100, labels=c("A","B"))
summary(aov(response ~ group, data=TDat))

11225.25/3.64
[1] 3083.86

There is some rounding error on the mean squares (i.e. mean variances) but F
is correct. Using estimates calculated by a different route we have:

11225.249057/3.639801
[1] 3084.028

Does this answer your question?

Regards, Mark.


emj83 wrote:
> 
> I have done this in R and this is the following ANOVA table I get:
> 
>> summary(aov(response ~ group, data=TDat))
>  Df  Sum Sq Mean Sq F valuePr(>F)
> group 1 11203.5 11203.5  2505.0 < 2.2e-16 ***
> Residuals   198   885.5 4.5
> 
> The model is response(i,j)= group(i)+ error(i,j),
> 
> we assume that group~N(0,P^2) and error~N(0,sigma^2)
> 
> I know that sigma^2 is equal to 4.5, how do I find out P^2? 
> 
> In the problem that I am trying to apply this to, I have more than 2
> groups. I was hoping there would be a function that helps you do this that
> I don't know about.
> 
> 
> Thanks for your help Emma
> 
> 
> 
> 
> Mark Difford wrote:
>> 
>> Hi Emma,
>> 
>>>> 
>> 
>> R gives you the tools to work this out.
>> 
>> ## Example
>> set.seed(7)
>> TDat <- data.frame(response = c(rnorm(100, 5, 2), rnorm(100, 20, 2)))
>> TDat$group <- gl(2, 100, labels=c("A","B"))
>> with(TDat, boxplot(split(response, group)))
>> summary(aov(response ~ group, data=TDat))
>> 
>> Regards, Mark.
>> 
>> 
>> emj83 wrote:
>>> 
>>> can anyone advise me please?
>>> 
>>> 
>>> emj83 wrote:
>>>> 
>>>> I have done some ANOVA tables for some data that I have, from this I
>>>> can read the within-group variance. can anyone tell me how i may find
>>>> out the between-group variance?
>>>> 
>>>> Thanks Emma
>>>> 
>>> 
>>> 
>> 
>> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Between-group-variance-from-ANOVA-tp24954045p25129942.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Between-group variance from ANOVA

[Mark Difford]

Hi Emma,

...

I forgot to add the tabular ouput, which doesn't help either:

T.sum <- summary(aov(response ~ group, data=TDat))
print(T.sum)

 Df  Sum Sq Mean Sq F valuePr(>F)
group 1 11225.2 11225.23084 < 2.2e-16 ***
Residuals   198   720.7 3.6  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

unlist(T.sum)

unlist(T.sum)[5]/unlist(T.sum)[6]
Mean Sq1 
3084.028

Regards, Mark.


Mark Difford wrote:
> 
> Hi Emma,
> 
>>> ...from this I can read the within-group variance. can anyone tell me
>>> how i may find 
>>> out the between-group variance?
> 
> But it's in the table, above the "within-group" variance. Remember that F
> is the ratio of these two quantities, i.e. the mean of the group variances
> divided by the mean of the within-group variances . I will work with my
> example since you never set seed so your answers are different from mine
> (which really does not help matters).
> 
> set.seed(7) 
> TDat <- data.frame(response = c(rnorm(100, 5, 2), rnorm(100, 20, 2))) 
> TDat$group <- gl(2, 100, labels=c("A","B"))
> summary(aov(response ~ group, data=TDat))
> 
> 11225.25/3.64
> [1] 3083.86
> 
> There is some rounding error on the mean squares (i.e. mean variances) but
> F is correct. Using estimates calculated by a different route we have:
> 
> 11225.249057/3.639801
> [1] 3084.028
> 
> Does this answer your question?
> 
> Regards, Mark.
> 
> 
> emj83 wrote:
>> 
>> I have done this in R and this is the following ANOVA table I get:
>> 
>>> summary(aov(response ~ group, data=TDat))
>>  Df  Sum Sq Mean Sq F valuePr(>F)
>> group 1 11203.5 11203.5  2505.0 < 2.2e-16 ***
>> Residuals   198   885.5 4.5
>> 
>> The model is response(i,j)= group(i)+ error(i,j),
>> 
>> we assume that group~N(0,P^2) and error~N(0,sigma^2)
>> 
>> I know that sigma^2 is equal to 4.5, how do I find out P^2? 
>> 
>> In the problem that I am trying to apply this to, I have more than 2
>> groups. I was hoping there would be a function that helps you do this
>> that I don't know about.
>> 
>> 
>> Thanks for your help Emma
>> 
>> 
>> 
>> 
>> Mark Difford wrote:
>>> 
>>> Hi Emma,
>>> 
>>>>> 
>>> 
>>> R gives you the tools to work this out.
>>> 
>>> ## Example
>>> set.seed(7)
>>> TDat <- data.frame(response = c(rnorm(100, 5, 2), rnorm(100, 20, 2)))
>>> TDat$group <- gl(2, 100, labels=c("A","B"))
>>> with(TDat, boxplot(split(response, group)))
>>> summary(aov(response ~ group, data=TDat))
>>> 
>>> Regards, Mark.
>>> 
>>> 
>>> emj83 wrote:
>>>> 
>>>> can anyone advise me please?
>>>> 
>>>> 
>>>> emj83 wrote:
>>>>> 
>>>>> I have done some ANOVA tables for some data that I have, from this I
>>>>> can read the within-group variance. can anyone tell me how i may find
>>>>> out the between-group variance?
>>>>> 
>>>>> Thanks Emma
>>>>> 
>>>> 
>>>> 
>>> 
>>> 
>> 
>> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Between-group-variance-from-ANOVA-tp24954045p25130266.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] adding factor scores back to an incomplete dataset...

[Mark Difford]

Hi David,

>> I am doing a factor analysis on survey data with missing values. ... Is
>> there a way to subset
>> from my original data set that will work with factanal() and preserve the
>> original rows or that 
>> will allow me to append the factor scores back onto the original dataset
>> with the proper 
>> rows and NAs where there could be no data?

The "trick" is to call factanal with the formula interface, and to specify
na.action=na.exclude. This is documented under ?factanal &c and will put NAs
in the right place(s) to pad out scores.

## Dummy
T.fac <- factanal(~ SrvP+EdcT+EcnD+BdvP+EnvP, data=TTDat, factors=2,
na.action=na.exclude, scores="regression")

?na.exclude

Regards, Mark.


David G. Tully wrote:
> 
> I am sure there is a simple way to do the following, but i haven't been 
> able to find it. I am hoping a merciful soul on R-help could point me in 
> the right direction.
> 
> I am doing a factor analysis on survey data with missing values. to do 
> this, I run:
> 
> FA1<-factanal(na.omit(DATA), factors = X, rotation = 'oblimin', scores = 
> 'regression')
> 
> Now that I have my factors and factor scores, I want to add those scores 
> back to my original dataset so I can plot factor scores by demographics. 
> However, when I try to add the scores back to the original data frame, 
> the variables are of different lengths.
> 
> Is there a way to subset from my original data set that will work with 
> factanal() and preserve the original rows or that will allow me to 
> append the factor scores back onto the original dataset with the proper 
> rows and NAs where there could be no data?
> 
> Again, I apologize if I am missing something basic. I am a self taught R 
> user and couldn't find an answer to this question.
> 
> Thanks in advance,
>David
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/adding-factor-scores-back-to-an-incomplete-dataset...-tp25140959p25142328.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] adding factor scores back to an incomplete dataset...

[Mark Difford]

Hi David, Phil,

Phil Spector wrote:
>> David - 
>> Here's the easiest way I've been able to come up with.

Easiest? You are making unnecessary work for yourselves and seem not to
understand the purpose of ?naresid (i.e. na.action = na.exclude). Why not
take the simple route that I gave, which really is R's + factanal's route.
Using Phil's data as example:

##
dat = data.frame(matrix(rnorm(100),20,5)) 
dat[3,4] = NA 
dat[12,3] = NA

scrs <- factanal(~X1+X2+X3+X4+X5, data=dat,factors=2,scores='regression',
 na.action=na.exclude)$scores
TrueDat <- merge(dat,scrs,by=0,all.x=TRUE,sort=FALSE)
TrueDat

Regards, Mark.


David G. Tully wrote:
> 
> Thanks, Prof Spector. Your first solution works well for me.
> 
> Phil Spector wrote:
>> David -
>>Here's the easiest way I've been able to come up with. I'll provide 
>> some sample data to make things clearer (hint, hint):
>>
>>> dat = data.frame(matrix(rnorm(100),20,5))
>>> dat[3,4] = NA
>>> dat[12,3] = NA
>>> scrs = factanal(na.omit(dat),factors=2,scores='regression')$scores
>>> rownames(scrs) = rownames(na.omit(dat))
>>> newdat = merge(dat,scrs,by=0,all.x=TRUE,sort=FALSE)
>>
>> This will result in the observations with missing values being
>> at the end of the data frame.  If you want the original order
>> (assuming default row names), you could use
>>
>> newdat[order(as.numeric(newdat$Row.names)),]
>>
>> A somewhat more complicated approach is, in some sense, more direct:
>>
>>> dat$Factor1 = NA
>>> dat$Factor2 = NA
>>> dat[rownames(na.omit(dat[,-c(6,7)])),c('Factor1','Factor2')] = 
>> +
>> factanal(na.omit(dat[,-c(6,7)]),factors=2,scores='regression')$scores
>>
>> The order of the data is preserved.
>> - Phil Spector
>>  Statistical Computing Facility
>>  Department of Statistics
>>  UC Berkeley
>>  spec...@stat.berkeley.edu
>>
>>
>>
>>
>>
>>
>>
>>
>> On Tue, 25 Aug 2009, David G. Tully wrote:
>>
>>> I am sure there is a simple way to do the following, but i haven't 
>>> been able to find it. I am hoping a merciful soul on R-help could 
>>> point me in the right direction.
>>>
>>> I am doing a factor analysis on survey data with missing values. to 
>>> do this, I run:
>>>
>>> FA1<-factanal(na.omit(DATA), factors = X, rotation = 'oblimin', 
>>> scores = 'regression')
>>>
>>> Now that I have my factors and factor scores, I want to add those 
>>> scores back to my original dataset so I can plot factor scores by 
>>> demographics. However, when I try to add the scores back to the 
>>> original data frame, the variables are of different lengths.
>>>
>>> Is there a way to subset from my original data set that will work 
>>> with factanal() and preserve the original rows or that will allow me 
>>> to append the factor scores back onto the original dataset with the 
>>> proper rows and NAs where there could be no data?
>>>
>>> Again, I apologize if I am missing something basic. I am a self 
>>> taught R user and couldn't find an answer to this question.
>>>
>>> Thanks in advance,
>>>  David
>>>
>>> __
>>> R-help@r-project.org 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-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/adding-factor-scores-back-to-an-incomplete-dataset...-tp25140959p25204698.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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-help@r-project.org

[Mark Difford]

Hi John,

>> Has a test for bimodality been implemented in R?

You may find the code at the URL below useful. It was written by Jeremy
Tantrum (a PhD of Werner Stuetzle's). Amongst other things there is a
function to plot the unimodal and bimodal Gaussian smoothers closest to the
observed data. A dip-test statistic is also calculated.

Regards, Mark.

http://www.stat.washington.edu/wxs/Stat593-s03/Code/jeremy-unimodality.R



John Sansom wrote:
> 
> Has a test for bimodality been implemented in R?
> 
> Thanks, John
> 
> NIWA is the trading name of the National Institute of Water & Atmospheric
> Research Ltd.
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/test-for-bimodality-In-Reply-To%3D-tp25216164p25220627.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Problem with xtabs(), exclude=NULL, and counting NA's

[Mark Difford]

>> I must say that this is slightly odd behavior to require both
>> na.action= AND exclude=.  Does anyone know of a justification?

Not strange at all.

?options

na.action, sub head "Options set in package stats." You need to override the
default setting.


ws-7 wrote:
> 
>>> xtabs(~wkhp, x, exclude=NULL, na.action=na.pass)
>> wkhp
>>  20   30   40   45   60 
>>   1    1   10    1    3    4
> 
> Thanks!  I must say that this is slightly odd behavior to require both
> na.action= AND exclude=.  Does anyone know of a justification?
> Shouldn't it be changed?   Ah well, if R were internally
> consistent or corresponded to typical programming Unix practices, it
> just wouldn't feel the same ... 
> 
> Cheers!
> 
>> --
>> David Winsemius, MD
>> Heritage Laboratories
>> West Hartford, CT
>>
>>
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Problem-with-xtabs%28%29%2C-exclude%3DNULL%2C-and-counting-NA%27s-tp25304142p25305878.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] How can I appoint a small part of the whole data

[Mark Difford]

Yichih,

Answer 2 is "correct," because your indexing specification for 1 is wrong.
You also seem to have left out a comma.

##
mu1990$wage[mu1990$edu==2|mu1990$edu==3|mu1990$edu==4, ] ## like this
mu1990$wage[mu1990$edu%in%2:4, ]

You really could have worked this out for yourself by looking at the results
of your subsetting/indexing operation.

Mark.


Yichih Hsieh wrote:
> 
> Dear all,
> 
> I got another problem:
> 
> if education have five levels
> 
> edu=1
> edu=2
> edu=3
> edu=4
> edu=5
> 
> If I want to appoint y=edu2~4 in 1990
> which programs is correct?
> I tried this two programs, they both work, but two results is different.
> 
> 1.
> fig2b<-reldist(y=mu1990$wage[mu1990$edu==2|3|4],..)
> 
> 
> 2.
> fig2b<-reldist(y=mu1990$wage[mu1990$edu%in%2:4],..)
> 
> which one is correct?
> and why they have different results?
> 
> 
> All help high appreciated.
> 
> 
> best,
> Yichih
> 
> 2009/9/5 Yichih Hsieh 
> 
>>
>> Dear Petr,
>>
>> your suggestion is useful
>>
>> many thanks for your help !
>>
>>
>> best,
>> Yichih
>>
>> 2009/9/3 Petr PIKAL 
>>
>> Hi
>>>
>>> use any of suitable selection ways that are in R.
>>>
>>> E.g.
>>>
>>> data[data$gender==1, ]
>>>
>>> selects only female values
>>>
>>> data$wage[(data$gender==1)  & (data$race=1)] selects black female wages.
>>>
>>> and see also ?subset
>>>
>>> Regards
>>> Petr
>>>
>>> r-help-boun...@r-project.org napsal dne 03.09.2009 10:51:59:
>>>
>>> > Dear all,
>>> >
>>> > I have 1980~1990 eleven datas,
>>> > every year have three variables,
>>> > wage
>>> > gender(1=female, 2=male)
>>> > race(1=black, 2=white)
>>> >
>>> > My original commands is:
>>> >
>>> > fig2b<-reldist(y=mu1990$wage,yo=mu1980$wage,...)
>>> >
>>> > I have three questions:
>>> > 1. If I want to appoint y=women's wage in 1990
>>> > yo=women's wage in 1980
>>> > 2. If I want to appoint y=women's wage in 1990
>>> > yo=men's wage in 1990
>>> > 3. If I want to appoint y=black women's wage in 1990
>>> > yo=white women's wage in 1990
>>> >
>>> > How can I modify the commands?
>>> >
>>> > All help highly appreciated.
>>> >
>>> > Best,
>>> > Yichih
>>> >
>>> >
>>> > --
>>> > Yichih Hsieh
>>> >
>>> > e-mail : yichih.hs...@gmail.com
>>> >
>>> >[[alternative HTML version deleted]]
>>> >
>>> > __
>>> > R-help@r-project.org 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.
>>>
>>>
>>
>>
>> --
>> Yichih Hsieh
>>
>> e-mail : yichih.hs...@gmail.com
>>
> 
> 
> 
> -- 
> Yichih Hsieh
> 
> e-mail : yichih.hs...@gmail.com
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/How-can-I-appoint-a-small-part-of-the-whole-data-tp25272209p25308714.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Omnibus test for main effects in the face ofaninteraction containing the main effects.

[Mark Difford]

Hi John,

>> When Group is entered as a factor, and the factor has two levels, the 
>> ANOVA table gives a p value for each level of the factor. 

This does not (normally) happen so you are doing something strange. 

## From your first posting on this subject
fita<-lme(Post~Time+factor(Group)+factor(Group)*Time,
random=~1|SS,data=blah$alldata)

To begin with, what is blah$alldata? lme() asks for a data frame if the
formula interface is used. This looks like part of a list, or the column of
a data frame.

Have a look at the output below, from a syntactically correct model that has
essentially the same structure as yours. The data set should be loaded with
nlme so you can run it directly to see the result. Sex, with two levels, is
not split in the anova table.

##
str(Orthodont)
anova( lme(distance ~ age * Sex, data = Orthodont, random = ~ 1) )

Surely this is essentially what you are looking for? If Sex is not already a
factor, and it really is better to make it one when you build your data set,
then you can use as.factor as you have done, with the same result.

(Note: age * Sex expands to age + Sex + age:Sex, which equals the messy and
unnecessary age + Sex + age * Sex.)

Regards, Mark.


John Sorkin wrote:
> 
> Daniel,
> When Group is entered as a factor, and the factor has two levels, the
> ANOVA table gives a p value for each level of the factor. What I am
> looking for is the omnibus p value for the factor, i.e. the test that
> the factor (with all its levels) improves the prediction of the outcome.
> 
> You are correct that normally one could rely on the fact that the model 
> Post<-Time+as.factor(Group)+as.factor(Group)*Time
> contains the model 
> Post<-Time+as.factor(Group)
> and compare the two models using anova(model1,model2). However, my model
> is has a random effect, the comparison is not so easy. The REML
> comparions of nested random effects models is not valid when the fixed
> effects are not the same in the models, which is the essence of the
> problem in my case. 
> 
> In addition to the REML problem if one wants to perform an omnibus test
> for Group, one would want to compare nested models, one containing
> Group, and the other not containing group. This would suggest comparing
> Post<-Time+  as.factor(Group)*Time to
> Post<-Time+Group+as.factor(Group)*Time
> The quandry here is whether one should or not "allow" the first model as
> it is poorly specified - one term of the interaction,
> as.factor(Group)*Time, as.factor(Group) does not appear as a main effect
> - a no-no in model building. 
> John
> 
> 
> John Sorkin M.D., Ph.D.
> Chief, Biostatistics and Informatics
> Baltimore VA Medical Center GRECC,
> University of Maryland School of Medicine Claude D. Pepper OAIC,
> University of Maryland Clinical Nutrition Research Unit, and
> Baltimore VA Center Stroke of Excellence
> 
> University of Maryland School of Medicine
> Division of Gerontology
> Baltimore VA Medical Center
> 10 North Greene Street
> GRECC (BT/18/GR)
> Baltimore, MD 21201-1524
> 
> (Phone) 410-605-7119
> (Fax) 410-605-7913 (Please call phone number above prior to faxing)
> jsor...@grecc.umaryland.edu
 "Daniel Malter"  09/07/09 9:23 PM >>>
> John, your question is confusing. After reading it twice, I still cannot
> figure out what exactly you want to compare.
> 
> Your model "a" is the unrestricted model, and model "b" is a restricted
> version of model "a" (i.e., b is a hiearchically reduced version of a,
> or
> put differently, all coefficients of b are in a with a having additional
> coefficients). Thus, it is appropriate to compare the models (also
> called
> nested models).
> 
> Comparing c with a and d with a is also appropriate for the same reason.
> However, note that depedent on discipline, it may be highly
> unconventional
> to fit an interaction without all direct effects of the interacted
> variables
> (the reason for this being that you may get biased estimates).
> 
> What you might consider is:
> 1. Run an intercept only model
> 2. Run a model with group and time
> 3. Run a model with group, time, and the interaction
> 
> Then compare 2 to 1, and 3 to 2. This tells you whether including more
> variables (hierarchically) makes your model better.
> 
> HTH,
> Daniel
> 
> On a different note, if lme fits with "restricted maximum likelihood," I
> think I remember that you cannot compare them. You have to fit them with
> "maximum likelihood." I am pointing this out because lmer with
> restricted
> maximum likelihood by standard, so lme might too.
> 
> -
> cuncta stricte discussurus
> -
> 
> -Ursprüngliche Nachricht-
> Von: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
> Im
> Auftrag von John Sorkin
> Gesendet: Monday, September 07, 2009 4:00 PM
> An: r-help@r-project.org
> Betreff: [R] Omnibus test for main effects in the face of aninteraction
> containing the main effects.
> 
> R 2.9.1
> Windows XP
> 
> UPDATE,
> Even my first

Re: [R] Alternative to Scale Function?

[Mark Difford]

>> The scale function will return the mean and sd of the data.

By default. Read ?scale.

Mark.


Noah Silverman-3 wrote:
> 
> I think I just answered my own question.
> 
> The scale function will return the mean and sd of the data.
> 
> So the process is fairly simple.
> scale training data varaible
> note mean and sd from the scale
> then manually scale the test data using the mean and sd from the 
> training data.
> 
> That should make sure that a value is transformed the same regardless of 
> which data set it is in.
> 
> Do I have this correct, or can anybody contribute any more to the concept?
> 
> Thanks!
> 
> 
> --
> Noah
> 
> On 9/11/09 1:10 PM, Noah Silverman wrote:
>> Hi,
>>
>> Is there an alternative to the scale function where I can specify my 
>> own mean and standard deviation?
>>
>> I've come across an interesting issue where this would help.
>>
>> I'm training and testing on completely different sets of data.  The 
>> testing set is smaller than the training set.
>>
>> Using the standard scale function of R seems to introduce some error.  
>> Since it scales data WITHIN the set, it may scale the same number to 
>> different value since the range in the training and testing set may be 
>> different.
>>
>> My thought was to scale the larger training set of data, then use the 
>> mean and SD of the training data to scale the testing data according 
>> to the same parameters.  That way a number will transform to the same 
>> result regardless of whether it is in the training or testing set.
>>
>> I can't be the first one to have looked at this.  Does anyone know of 
>> a function in R or if there is a scale alternative where I can control 
>> the parameters?
>>
>> Thanks!
>>
>> -- 
>> Noah
>>
>> __
>> R-help@r-project.org 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-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Alternative-to-Scale-Function--tp25407625p25408289.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] could not find function "Varcov" after upgrade of R?

[Mark Difford]

Hi Zhu,

>> could not find function "Varcov" after upgrade of R?

Frank Harrell (author of Design) has noted in another thread that Hmisc has
changed... The problem is that functions like anova.Design call a function
in the _old_ Hmisc package called Varcov.default. In the new version of
Hmisc this is called something else (vcov.default).

The best way to fix this is to install the new (i.e. current) version of
Hmisc and Frank's replacement for his Design package, which is called rms
(Regression Modeling Strategies).

Regards, Mark.


zhu yao wrote:
> 
> I uses the Design library.
> 
> take this example:
> 
> library(Design)
> n <- 1000
> set.seed(731)
> age <- 50 + 12*rnorm(n)
> label(age) <- "Age"
> sex <- factor(sample(c('Male','Female'), n,
>   rep=TRUE, prob=c(.6, .4)))
> cens <- 15*runif(n)
> h <- .02*exp(.04*(age-50)+.8*(sex=='Female'))
> dt <- -log(runif(n))/h
> label(dt) <- 'Follow-up Time'
> e <- ifelse(dt <= cens,1,0)
> dt <- pmin(dt, cens)
> units(dt) <- "Year"
> dd <- datadist(age, sex)
> options(datadist='dd')
> Srv <- Surv(dt,e)
> 
> f <- cph(Srv ~ rcs(age,4) + sex, x=TRUE, y=TRUE)
> cox.zph(f, "rank") # tests of PH
> anova(f)
> # Error in anova.Design(f) : could not find function "Varcov"
> 
> 
> 
> Yao Zhu
> Department of Urology
> Fudan University Shanghai Cancer Center
> No. 270 Dongan Road, Shanghai, China
> 
> 
> 2009/9/12 Ronggui Huang 
> 
>> I cannot reproduce the problem you mentioned.
>>
>> >  ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
>> >  trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
>> >   group <- gl(2,10,20, labels=c("Ctl","Trt"))
>> >   weight <- c(ctl, trt)
>> >   anova(lm.D9 <- lm(weight ~ group))
>> > sessionInfo()
>> R version 2.9.2 (2009-08-24)
>> i386-pc-mingw32
>>
>> locale:
>> LC_COLLATE=Chinese (Simplified)_People's Republic of
>> China.936;LC_CTYPE=Chinese (Simplified)_People's Republic of
>> China.936;LC_MONETARY=Chinese (Simplified)_People's Republic of
>> China.936;LC_NUMERIC=C;LC_TIME=Chinese (Simplified)_People's Republic
>> of China.936
>>
>> attached base packages:
>> [1] stats graphics  grDevices utils datasets  methods   base
>>
>> 2009/9/12 zhu yao :
>> > After upgrading R to 2.9.2, I can't use the anova() fuction.
>> > It says "could not find function "Varcov" ".
>> > What's wrong with my computer? Help needed, thanks!
>> >
>> > Yao Zhu
>> > Department of Urology
>> > Fudan University Shanghai Cancer Center
>> > No. 270 Dongan Road, Shanghai, China
>> >
>> >[[alternative HTML version deleted]]
>> >
>> > __
>> > R-help@r-project.org 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.
>> >
>>
>>
>>
>> --
>> HUANG Ronggui, Wincent
>> Doctoral Candidate
>> Dept of Public and Social Administration
>> City University of Hong Kong
>> Home page: http://asrr.r-forge.r-project.org/rghuang.html
>>
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/could-not-find-function-%22Varcov%22-after-upgrade-of-R--tp25412881p25414257.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Triangular distribution for kernel regression

[Mark Difford]

Hi Brian,

>> I am trying to get fitted/estimated values using kernel regression and a 
>> triangular kernel.

Look at Loader's locfit package. You are likely to be pleasantly surprised.

Regards, Mark.


Bryan-65 wrote:
> 
> Hello,
> 
> I am trying to get fitted/estimated values using kernel regression and a
> triangular kernel.  I have found packages that easily fit values from a
> kernel regression (e.g. ksmooth) but do not have a triangular distribution
> option, and density estimators that have triangular distribution options
> that I can't seem to use to produce estimated values (e.g. density).  Any
> help is appreciated.
> 
> Bryan
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Triangular-distribution-for-kernel-regression-tp25416706p25417878.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] hel with factorial analysis with mixed data

[Mark Difford]

andreiabb wrote:

>> the message that I am getting is 

>> Error in AFDM (all_data_sub.AFDM, type=c(rep("s",1), rep("n",1), rep("n",
>> : 
>> unused arguments (s) (type=c("s", "n","n")) 

>> Can someone help me?

If you are in hel[l] then it is entirely your own fault. The error message
is clear and would have become even more clear had you bothered to
read/check the help file for AFDM for yourself. There is no argument "type=
" to AFDM.

##
?AFDM

I hope this will lighten your way.

Mark.


andreiabb wrote:
> 
> Dear forum
> 
> I am trying to use the package FactoMineR, since I have 2 class variables
> and 1 continuous variable I am trying to use AFDM, factorial analysis of
> mixed data.
> I have 3 variables, 2 qualitative and one quantitative:
> 
> the code I am using is 
> 
> all_data_sub.AFM<-all_data_sub [, c=("V4","comb","V2")]
> 
> res<-AFDM(all_data_AFDM,type=c(rep("s",1), rep("n",1), rep("n",1)), ncp=5,
> sup.var=3:3, graph=FALSE)
> 
> the message that I am getting is
> 
> Error in AFDM (all_data_sub.AFDM, type=c(rep("s",1), rep("n",1), rep("n",
> :
> unused arguments (s) (type=c("s", "n","n"))
> 
> Can someone help me?
> thanks
> 

-- 
View this message in context: 
http://www.nabble.com/hel-with-factorial-analysis-with-mixed-data-tp25453138p25453922.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] How to do a PCA with ordered variables?

[Mark Difford]

P. Branco wrote:

>> I have used the dudi.mix method from the ade4 package, but when I do the
>> $index it shows 
>> me that R has considered my variables as quantitative. 

>> What should I do?

You should make sure that they are encoded as ordered factors, which has
nothing to do with ade4's dudi.mix().

##
?ordered

Mark.



P.Branco wrote:
> 
> Hi,
> 
> I want to do a pca using a set of variables which are ordered. I have used
> the dudi.mix method from the ade4 package, but when I do the $index it
> shows me that R has considered my variables as quantitative.
> 
> What should I do?
> 

-- 
View this message in context: 
http://www.nabble.com/How-to-do-a-PCA-with-ordered-variables--tp25491950p25491969.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] The Claw Density and LOCFIT

[Mark Difford]

Hi Jaap,

>> Could anybody please direct me in finding an updated version of this
>> document, or help me
>> correct the code given in the file. The (out-of-date) code is as follows:

You are not helping yourself, or anyone else, by not including the error
messages you get when trying to execute "your" code. The matter is in fact
quite straightfowards: ev does not accept an evaluation structure called
grid. The current documentation for locfit does tell you this: ?locfit.raw
(sub ev). When I executed your code I got

>   fit1 <- locfit( ~ claw54, deg = 0, kern = "gauss", alpha = c(0, 0.315),
> ev = 
+ grid(100, ll = -3.5, ur = 2.7)) 
Error in grid(100, ll = -3.5, ur = 2.7) : 
  unused argument(s) (ll = -3.5, ur = 2.7)

This explicitly tells you that the problem lies with the call to ev =
grid(...), which should in fact be ev = lfgrid(...). The following works for
me and should do so for you.

fit1 <- locfit( ~ claw54, deg = 0, kern = "gauss", alpha = c(0, 0.315), ev = 
   lfgrid(100, ll = -3.5, ur = 2.7)) 

HTH,
Mark.

R version 2.10.0 Under development (unstable) (2009-06-20 r48806) 
i386-pc-mingw32 

locale:
[1] LC_COLLATE=English_South Africa.1252  LC_CTYPE=English_South Africa.1252   
[3] LC_MONETARY=English_South Africa.1252 LC_NUMERIC=C 
[5] LC_TIME=English_South Africa.1252

attached base packages:
[1] splines   stats graphics  grDevices utils datasets  methods  
base 

other attached packages:
[1] locfit_1.5-4lattice_0.17-25 akima_0.5-2 DPpackage_1.0-7
ade4_1.4-11 Design_2.2-0   
[7] survival_2.35-4 Hmisc_3.6-0

loaded via a namespace (and not attached):
[1] cluster_1.12.0 grid_2.10.0tools_2.10.0 


Van Wyk, Jaap wrote:
> 
> I am trying to reproduce Figure 10.5 of Loader's book: Local Regression
> and Likelihood. The code provided in the book does not seem to work.
> I have managed (a while ago) to get the accompanied R-code for the figures
> in the book (file called lffigs.R) from somewhere - cannot find it on the
> web anymore. The code in the .R script file does not work either.
> Could anybody please direct me in finding an updated version of this
> document, or help me correct the code given in the file. The (out-of-date)
> code is as follows:
> 
>data(claw54)
>   fit1 <- locfit( ~ claw54, deg = 0, kern = "gauss", alpha = c(0, 0.315),
> ev =
> grid(100, ll = -3.5, ur = 2.7))
>   fit2 <- locfit( ~ claw54, deg = 0, kern = "gauss", alpha = c(0, 0.985),
> ev =
> grid(100, ll = -3.5, ur = 2.7))
>   x <- seq(-3.5, 2.7, length.out = 200)
>   y <- dnorm(x, -1., 0.1) + dnorm(x, -0.5, 0.1) + dnorm(x, 0, 0.1) +
> dnorm(x,
> 0.5, 0.1) + dnorm(x, 1., 0.1)
>   y <- (y + 5 * dnorm(x))/10
>   plot(fit1, get.data = T, main = "h=0.315", ylim = c(0, max(y)))
>   lines(x, y, lty = 2)
>   plot(fit2, get.data = T, main = "h=0.985", ylim = c(0, max(y)))
>   lines(x, y, lty = 2)
> 
> THANKS FOR ANY ASSISTANCE.
> ps: This code differs from that in the book. I have tried both, without
> success. Even if I just use, for example,
> fit1 <- locfit( ~ claw54, deg = 0, kern = "gauss", alpha = c(0, 0.315))
> I do not get the same result.
> 
> Jacob L van Wyk, Dept. of Statistics, University of Johannesburg (APK),
> Box 524, Auckland Park, 2006.
> Office: +27 11 559 3080, Fax: +27 11 559 2499
> 
> 
> 
> This email and all contents are subject to the following disclaimer:
> 
> http://www.uj.ac.za/UJ_email_legal_disclaimer.htm
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/The-Claw-Density-and-LOCFIT-tp24218274p24220007.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Correlation Network Diagram?

[Mark Difford]

Hi Rory,

There are several. Have a look at the gR Task Views. There you will also
find a link to the statnet suite, where you will find links to a dedicated
set of jstatsoft articles.

Regards, Mark.


Rory Winston wrote:
> 
> Hi all
> 
> On page 39 of this paper [1] by Andrew Lo there is a very interesting  
> correlation network diagram (sorry I dont have a more convenient link to  
> the type of diagram I'm talking about). does anyone know of any package in  
> R that can generate these types of diagrams?
> 
> Cheers
> -- Rory
> 
> [1] http://web.mit.edu/alo/www/Papers/august07.pdf
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Correlation-Network-Diagram--tp24321104p24329465.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Differing Variable Length Inconsistencies in Random Effects/Regression Models

[Mark Difford]

Hi Aditi,

Parts of _your_ code for the solution offered by Jerome Goudet are wrong;
see my comments.

> famfit<-lmer(peg.no~1 + (1|family), na.action=na.omit, vcdf)   ## use:
> na.action=na.exclude
> resfam<-residuals(famfit) 
> for( i in 1:length(colms)) 
+ { 
+ print ("Marker", i) 
+ regfam<-abline(lm(resfam~i))## you need to use:
abline(lm(resfam~colms[,i]))
+ print(regfam) 
+ }

Corrected code:
famfit<-lmer(peg.no~1 + (1|family), na.action=na.exclude, vcdf)
resfam<-residuals(famfit) 
for( i in 1:length(colms)) 
{ 
print ("Marker", i) 
regfam<-abline(lm(resfam~colms[,i]))
}

This should work.

Regards, Mark.


A Singh wrote:
> 
> 
> Dear All,
> 
> I am quite new to R and am having a problem trying to run a linear model 
> with random effects/ a regression- with particular regard to my variable 
> lengths being different and the models refusing to compute any further.
> 
> The codes I have been using are as follows:
> 
> vc<-read.table("P:\\R\\Testvcomp10.txt",header=T)
>>> attach(vc)
>>
>> family<-factor(family)
>> colms<-(vc)[,4:13] ## this to assign the 10 columns containing marker
>> datato a new variable, as column names are themselves not in any
>> recognisable sequence
>>
>> vcdf<-data.frame(family,peg.no,ec.length,syll.length,colms)
>> library(lme4)
> 
>>> for (c in levels(family))
>> + {for (i in 1:length(colms))
>> +{ fit<-lmer(peg.no~1 + (1|c/i), vcdf)
>> +}
>> +summ<-summary(fit)
>> +av<-anova(fit)
>> +print(summ)
>> +print(av)
>> + }
>>
>> This gives me:
>>
>> Error in model.frame.default(data = vcdf, formula = peg.no ~ 1 + (1 +  :
>>  variable lengths differ (found for 'c')
> 
> 
> I had posted a similar message on the R mixed model list a few days ago, 
> with respect to my fundamental methods, and Jerome Goudet had kindly 
> referred me to an alternative approach using residuals obtained from a 
> random effects model in lmer(), and then doing regressions using those 
> [residuals being the dependent variable and my marker data columns the 
> independent variable].
> 
> The code for that is as follows:
> 
>  vc<-read.table("P:\\R\\Text 
> Files\\Testvcomp10.txt",header=T,sep="",dec=".",na.strings=NA,strip.white=T)
>> attach(vc)
>>
>> family<-factor(family)
>> colms<-(vc)[,4:13]
>>
>> names(vc)
>  [1] "male.parent"  "family"   "offspring.id" "P1L55""P1L73" 
> 
>  [6] "P1L74""P1L77""P1L91""P1L96""P1L98" 
> 
> [11] "P1L100"   "P1L114"   "P1L118"   "peg.no" 
> "ec.length"
> [16] "syll.length"
>>
>> vcdf<-data.frame(family, colms, peg.no, ec.length, syll.length)
>>
>> library(lme4)
> 
>> famfit<-lmer(peg.no~1 + (1|family), na.action=na.omit, vcdf)
>> resfam<-residuals(famfit)
>> for( i in 1:length(colms))
> + {
> + print ("Marker", i)
> + regfam<-abline(lm(resfam~i))
> + print(regfam)
> + }
> 
> This again gives me the error:
> 
> 
> [1] "Marker"
> Error in model.frame.default(formula = resfam ~ i, drop.unused.levels = 
> TRUE) :
>   variable lengths differ (found for 'i')
> 
> 
> My variables all have missing values somewhere or the other. The missing 
> values are not consistent for all individuals, i.e some individuals have 
> some values missing, others have others,
>  and as much as I have tried to use na.action=na.omit and na.rm=T, the 
> differing variable length problem is dogging me persistently..
> 
> I also tried to isolate the residuals, save them in a new variable (called 
> 'resfam' here), and tried to save that in the data.frame()->vcdf, that I 
> had created earlier.
> 
> The problem with that was that when the residuals were computed, lmer() 
> ignored missing data in 'peg.no' with respect to 'family', which is 
> obviously not the same data missing for say another variable E.g. 
> 'ec.length'- leading again to an inconsistency in variable lengths. 
> Data.frame would then not accept that addition at all to the previous set.
> This was fairly obvious right from the start, but I decided to try it 
> anyway. Didn't work.
> 
> I apologise if the solution to working with these different variable 
> lengths is obvious and I don't know it- but I don't know R that well at
> all.
> 
> My data files can be downloaded at the following location:
> 
>  (excel-
> .xlsx)
> 
> 
> (.txt file)
> 
> 
> Any pointers would be greatly appreciated, as this is holding me up loads.
> 
> Thanks a ton for your help,
> 
> Aditi
> 
> 
> 
> --
> A Singh
> aditi.si...@bristol.ac.uk
> School of Biological Sciences
> University of Bristol
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http:

Re: [R] Differing Variable Length Inconsistencies in Random Effects/Regression Models

[Mark Difford]

Perhaps I should have added the following: To see that it "works," run the
following:

famfit<-lmer(peg.no~1 + (1|family), na.action=na.exclude, vcdf)
resfam<-residuals(famfit) 
for( i in 1:length(colms)) 
{ 
print(coef(lm(resfam~colms[,i])))
}

Regards, Mark.


A Singh wrote:
> 
> 
> Dear All,
> 
> I am quite new to R and am having a problem trying to run a linear model 
> with random effects/ a regression- with particular regard to my variable 
> lengths being different and the models refusing to compute any further.
> 
> The codes I have been using are as follows:
> 
> vc<-read.table("P:\\R\\Testvcomp10.txt",header=T)
>>> attach(vc)
>>
>> family<-factor(family)
>> colms<-(vc)[,4:13] ## this to assign the 10 columns containing marker
>> datato a new variable, as column names are themselves not in any
>> recognisable sequence
>>
>> vcdf<-data.frame(family,peg.no,ec.length,syll.length,colms)
>> library(lme4)
> 
>>> for (c in levels(family))
>> + {for (i in 1:length(colms))
>> +{ fit<-lmer(peg.no~1 + (1|c/i), vcdf)
>> +}
>> +summ<-summary(fit)
>> +av<-anova(fit)
>> +print(summ)
>> +print(av)
>> + }
>>
>> This gives me:
>>
>> Error in model.frame.default(data = vcdf, formula = peg.no ~ 1 + (1 +  :
>>  variable lengths differ (found for 'c')
> 
> 
> I had posted a similar message on the R mixed model list a few days ago, 
> with respect to my fundamental methods, and Jerome Goudet had kindly 
> referred me to an alternative approach using residuals obtained from a 
> random effects model in lmer(), and then doing regressions using those 
> [residuals being the dependent variable and my marker data columns the 
> independent variable].
> 
> The code for that is as follows:
> 
>  vc<-read.table("P:\\R\\Text 
> Files\\Testvcomp10.txt",header=T,sep="",dec=".",na.strings=NA,strip.white=T)
>> attach(vc)
>>
>> family<-factor(family)
>> colms<-(vc)[,4:13]
>>
>> names(vc)
>  [1] "male.parent"  "family"   "offspring.id" "P1L55""P1L73" 
> 
>  [6] "P1L74""P1L77""P1L91""P1L96""P1L98" 
> 
> [11] "P1L100"   "P1L114"   "P1L118"   "peg.no" 
> "ec.length"
> [16] "syll.length"
>>
>> vcdf<-data.frame(family, colms, peg.no, ec.length, syll.length)
>>
>> library(lme4)
> 
>> famfit<-lmer(peg.no~1 + (1|family), na.action=na.omit, vcdf)
>> resfam<-residuals(famfit)
>> for( i in 1:length(colms))
> + {
> + print ("Marker", i)
> + regfam<-abline(lm(resfam~i))
> + print(regfam)
> + }
> 
> This again gives me the error:
> 
> 
> [1] "Marker"
> Error in model.frame.default(formula = resfam ~ i, drop.unused.levels = 
> TRUE) :
>   variable lengths differ (found for 'i')
> 
> 
> My variables all have missing values somewhere or the other. The missing 
> values are not consistent for all individuals, i.e some individuals have 
> some values missing, others have others,
>  and as much as I have tried to use na.action=na.omit and na.rm=T, the 
> differing variable length problem is dogging me persistently..
> 
> I also tried to isolate the residuals, save them in a new variable (called 
> 'resfam' here), and tried to save that in the data.frame()->vcdf, that I 
> had created earlier.
> 
> The problem with that was that when the residuals were computed, lmer() 
> ignored missing data in 'peg.no' with respect to 'family', which is 
> obviously not the same data missing for say another variable E.g. 
> 'ec.length'- leading again to an inconsistency in variable lengths. 
> Data.frame would then not accept that addition at all to the previous set.
> This was fairly obvious right from the start, but I decided to try it 
> anyway. Didn't work.
> 
> I apologise if the solution to working with these different variable 
> lengths is obvious and I don't know it- but I don't know R that well at
> all.
> 
> My data files can be downloaded at the following location:
> 
>  (excel-
> .xlsx)
> 
> 
> (.txt file)
> 
> 
> Any pointers would be greatly appreciated, as this is holding me up loads.
> 
> Thanks a ton for your help,
> 
> Aditi
> 
> 
> 
> --
> A Singh
> aditi.si...@bristol.ac.uk
> School of Biological Sciences
> University of Bristol
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Differing-Variable-Length-Inconsistencies-in-Random-Effects-Regression-Models-tp24502146p24506118.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/p

Re: [R] Split plot analysis problems

[Mark Difford]

Hi Jean-Paul,

>> ... since R is not able to extract residuals?

R can extract the residuals, but they are a "hidden" in models with an error
structure

##
str(aov(PH~Community*Mowing*Water + Error(Block)))
residuals(aov(PH~Community*Mowing*Water + Error(Block))$Block)
residuals(aov(PH~Community*Mowing*Water + Error(Block))$Within)

Regards, Mark.


Jean-Paul Maalouf wrote:
> 
> Hello,
> 
> I would be very grateful if someone could give me a hand with my split  
> plot design problems.
> 
> So here is my design :
> I am studying the crossed-effects of water (wet/dry) and mowing  
> (mowed/not-mowed = nm) on plant height (PH) within 2 types of plant  
> communities (Xerobromion and Mesobromion) :
> - Within each type of communities, I have localised 4 blocks
> - In each block, I have defined 4 plots in order to have the 4  
> possible treatments of both the water and mowing factors : nm/dry ;  
> mowed/dry ; mowed/wet ; nm/wet.
> 
> Here is my data table :
> 
> Community Block Mowing WaterPH
> 1   Mesob1  Mowed   Wet  7.40
> 2   Mesob1 nm   Wet 13.10
> 3   Mesob1  Mowed   Dry  5.55
> 4   Mesob1 nm   Dry 10.35
> 5   Mesob2 nm   Dry 10.70
> 6   Mesob2  Mowed   Dry  6.38
> 7   Mesob2 nm   Wet  9.75
> 8   Mesob2  Mowed   Wet  6.35
> 9   Mesob3 nm   Wet  9.60
> 10  Mesob3  Mowed   Dry  5.10
> 11  Mesob3 nm   Dry 10.05
> 12  Mesob3  Mowed   Wet  6.25
> 13  Mesob4 nm   Wet  9.00
> 14  Mesob4  Mowed   Wet  6.50
> 15  Mesob4 nm   Dry  7.75
> 16  Mesob4  Mowed   Dry  5.90
> 17  Xerob5 nm   Wet  7.69
> 18  Xerob5  Mowed   Wet  8.11
> 19  Xerob5 nm   Dry  3.98
> 20  Xerob5  Mowed   Dry  3.69
> 21  Xerob6 nm   Wet  5.24
> 22  Xerob6  Mowed   Wet  4.22
> 23  Xerob6 nm   Dry  6.55
> 24  Xerob6  Mowed   Dry  4.40
> 25  Xerob7  Mowed   Dry  3.79
> 26  Xerob7 nm   Dry  3.91
> 27  Xerob7 nm   Wet  9.00
> 28  Xerob7  Mowed   Wet  8.50
> 29  Xerob8  Mowed   Dry  3.33
> 30  Xerob8 nm   Wet  6.25
> 31  Xerob8  Mowed   Wet  8.00
> 32  Xerob8 nm   Dry  6.33
> 
> I actually have 2 questions :
> I wrote my model in two different ways, and there were differences in  
> P-Values according to the model written :
> 
> First model : summary(aov(PH~Community*Mowing*Water + Error(Block)))
> Error: Block
>Df Sum Sq Mean Sq F value   Pr(>F)
> Community  1 42.182  42.182  24.407 0.002603 **
> Residuals  6 10.370   1.728
> ---
> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
> 
> Error: Within
> Df Sum Sq Mean Sq F valuePr(>F)
> Mowing  1 40.007  40.007 21.1747 0.0002215 ***
> Water   1 23.120  23.120 12.2370 0.0025673 **
> Community:Mowing1 21.060  21.060 11.1467 0.0036554 **
> Community:Water 1  6.901   6.901  3.6524 0.0720478 .
> Mowing:Water1  1.611   1.611  0.8527 0.3680090
> Community:Mowing:Water  1  0.858   0.858  0.4542 0.5089331
> Residuals  18 34.008   1.889
> ---
> 
> - Second model (assuming that Mowing*Water are nested inside the Block  
> factor) :
> summary(aov(PH~Community*Mowing*Water + Error(Block/(Mowing*Water
> 
> Error: Block
>Df Sum Sq Mean Sq F value   Pr(>F)
> Community  1 42.182  42.182  24.407 0.002603 **
> Residuals  6 10.370   1.728
> ---
> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
> 
> Error: Block:Mowing
>   Df Sum Sq Mean Sq F valuePr(>F)
> Mowing1 40.007  40.007  37.791 0.0008489 ***
> Community:Mowing  1 21.060  21.060  19.893 0.0042820 **
> Residuals 6  6.352   1.059
> ---
> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
> 
> Error: Block:Water
>  Df  Sum Sq Mean Sq F value  Pr(>F)
> Water1 23.1200 23.1200  6.0725 0.04884 *
> Community:Water  1  6.9006  6.9006  1.8125 0.22685
> Residuals6 22.8439  3.8073
> ---
> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
> 
> Error: Block:Mowing:Water
> Df Sum Sq Mean Sq F value Pr(>F)
> Mowing:Water1 1.6110  1.6110  2.0085 0.2062
> Community:Mowing:Water  1 0.8581  0.8581  1.0697 0.3409
> Residuals   6 4.8126  0.8021
> 
> Both models give me interesting (but different!) results. Which one  
> would be the most appropriate?
> 
> Second question : How can I verify preliminary assumptions (normality  
> of residuals and variance homogeneity) in this kind of models?
> When I ask R to extract residuals, the answer is "NULL":
> 
>> residuals(aov(PH~Community*Mowing*Water + Error(Block/(Mowing*Water
> NULL
>> residuals(aov(PH~Community*Mowing*Water + Error(Block)))
> NULL
> 
> A huge thanks to the one who will rescue (or at least try to 

Re: [R] Split plot analysis problems

[Mark Difford]

Hi Jean-Paul,

>> However, I've tried both solutions on my model, and I got different
>> residuals :...
>> What could be the difference between the two?

There is no difference. You have made a mistake.

##
tt <- data.frame(read.csv(file="tt.csv", sep=""))  ## imports your data set
T.aov <- aov(PH~Community*Mowing*Water + Error(Block), data=tt)

summary(T.aov)  ## This matches your output

   Df Sum Sq Mean Sq F valuePr(>F)
Mowing  1 40.007  40.007 21.1747 0.0002215 ***
Water   1 23.120  23.120 12.2370 0.0025673 ** 
Community:Mowing1 21.060  21.060 11.1467 0.0036554 ** 
Community:Water 1  6.901   6.901  3.6524 0.0720478 .  
Mowing:Water1  1.611   1.611  0.8527 0.3680090
Community:Mowing:Water  1  0.858   0.858  0.4542 0.5089331
Residuals  18 34.008   1.889  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

##
length(residuals(T.aov <- aov(PH~Community*Mowing*Water + Error(Block),
data=tt)$Within))
[1] 24
length(proj(T.aov)[,"Residuals"])
[1] 24

## Details
residuals(T.aov <- aov(PH~Community*Mowing*Water + Error(Block),
data=tt)$Within)

   9   10   11   12   13  
14   15   16 
-1.107277057 -0.977277057 -0.007277057 -0.719777057 -1.021627882 
0.215872118 -1.621627882  0.508372118 
  17   18   19   20   21  
22   23   24 
 1.241311954  1.498811954 -0.616188046  0.483811954 -0.477827381
-1.660327381  2.684672619  1.924672619 
  25   26   27   28   29  
30   31   32 
-0.465198010 -1.735198010  1.502301990  0.839801990 -0.428515358
-0.751015358  0.836484642  1.181484642 

proj(T.aov)[,"Residuals"]
   9   10   11   12   13  
14   15   16   17   18   19  
20   21 
-1.107277057 -0.977277057 -0.007277057 -0.719777057 -1.021627882 
0.215872118 -1.621627882  0.508372118  1.241311954  1.498811954 -0.616188046 
0.483811954 -0.477827381 
  22   23   24   25   26  
27   28   29   30   31   32 
-1.660327381  2.684672619  1.924672619 -0.465198010 -1.735198010 
1.502301990  0.839801990 -0.428515358 -0.751015358  0.836484642  1.181484642 

Regards, Mark.


Jean-Paul Maalouf wrote:
> 
> Thanks Mark and Richard for your propositions on how to extract residuals.
> 
> However, I've tried both solutions on my model, and I got different  
> residuals :
> If we consider the within residuals :
> 
> Mark's solution gave me a vector of 24 residuals :
> 
>> as.vector(residuals(aov(PH~Community*Mowing*Water +
>> Error(Block))$Within))
>   [1] -1.107277057 -0.977277057 -0.007277057 -0.719777057 -1.021627882  
>   0.215872118 -1.621627882  0.508372118  1.241311954  1.498811954  
> -0.616188046
> [12]  0.483811954 -0.477827381 -1.660327381  2.684672619  1.924672619  
> -0.465198010 -1.735198010  1.502301990  0.839801990 -0.428515358  
> -0.751015358
> [23]  0.836484642  1.181484642
> 
> and Richard's solution gave me a vector 32 residuals :
> 
> do.call(data.frame,proj(aov(PH~Community*Water*Mowing +
> Error(Block->proj
> proj$Within.Residuals
>   [1] -0.216875  1.745625 -1.174375 -0.354375  0.800625  0.460625  
> -0.799375 -0.461875 -0.404375 -0.274375  0.695625 -0.016875 -0.541875   
> 0.695625 -1.141875
> [16]  0.988125  0.589375  0.846875 -1.268125 -0.168125 -1.095625  
> -2.278125  2.066875  1.306875 -0.500625 -1.770625  1.466875  0.804375  
> -0.638125 -0.960625
> [31]  0.626875  0.971875
> 
> What could be the difference between the two?
> 
> Regards
> 
> JPM
> 
> 
> Quoting "Richard M. Heiberger" :
> 
>> Jean-Paul Maalouf wrote:
>>> Do you have any idea on how can I verify preliminary assumptions in  
>>>  this model (normality of the residuals and variance homogeneity),   
>>> since R is not able to extract residuals?
>>
>> Of course, R extracts residuals.  Use the proj() function.  See ?proj
>> for the example
>> to get the projection of an aovlist object onto the terms of a linear
>> model
>>
>> proj(npk.aovE)
>>
>>
>> To get the projections into a simple data.frame, use
>> tmpE <- proj(npk.aovE)
>> tmpE.df <- do.call("data.frame", tmpE)
>> tmpE.df
>>
>> Mark Difford's solution effectively retrieved columns 3 and 10 from
>> tmpE.df
>>
>> Rich
> 
> 
> 
> -- 
> Jean-Paul Maalouf
> UMR 1202 BIOGECO
> Inra - Université Bordeaux 1
> Bâtiment B8, Avenue des Facultés
> 33405 Talence, France
> Tel : 05 40008772
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View th

Re: [R] how to correlate nominal variables?

[Mark Difford]

Hi Timo,

>> I need functions to calculate Yule's Y or Cramérs Index... Are such
>> functions existing?

Also look at assocstats() in package vcd.

Regards, Mark.


Timo Stolz wrote:
> 
> Dear R-Users,
> 
> I need functions to calculate Yule's Y or Cramérs Index, in order to
> correlate variables that are nominally scaled?
> 
> Am I wrong? Are such functions existing?
> 
> Sincerely,
> Timo
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/how-to-correlate-nominal-variables--tp18441195p24688304.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Adding picture to graph?

[Mark Difford]

Hi Rainer,

>> the question came up if it would be possible to add a picture
>> (saved on the HDD) to a graph (generated by plot()), which
>> we could not answer.

Yes. Look at package pixmap and, especially, at the examples sub s.logo() in
package ade4.

Regards, Mark.




Rainer M Krug-6 wrote:
> 
> Hi
> 
> while teaching R, the question came up if it would be possible to add
> a picture (saved on the HDD) to a graph (generated by plot()), which
> we could not answer.
> 
> 
> 
> It might easily kill a clean graph, but: is there a way of doing this,
> even one should not do it?
> 
> 
> On a similar line of thought: is it possibe to define own symbols so
> that they can be used in the plot function with pch=?
> 
> 
> Rainer
> 
> -- 
> Rainer M. Krug, Centre of Excellence for Invasion Biology,
> Stellenbosch University, South Africa
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Adding-picture-to-graph--tp24714724p24716913.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Selecting Bootstrap Method for Quantile Regression

[Mark Difford]

Hi Tom,

>> For example, if I want to use the "xy-pair bootstrap" how do I indicate
>> this in summary.rq?

The general approach is documented under summary.rq (sub se option 5).
Shorter route is boot.rq, where examples are given.

## ?boot.rq
y <- rnorm(50)
x <- matrix(rnorm(100),50)
fit <- rq(y~x,tau = .4)
summary(fit,se = "boot", bsmethod= "xy")

Regards, Mark.


Tom La Bone wrote:
> 
> The help page and vignette for summary.rq(quantreg) mention that there are
> three different bootstrap methods available for the se="bootstrap"
> argument, but I can't figure out how to select a particular method. For
> example, if I want to use the "xy-pair bootstrap" how do I indicate this
> in summary.rq?
> 
> Tom
> 

-- 
View this message in context: 
http://www.nabble.com/Selecting-Bootstrap-Method-for-Quantile-Regression-tp24737299p24744179.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] i'm so stuck with text file and contour plot

[Mark Difford]

Hannes,

>> been trying to read a text file that contains heading in the first line
>> in to R but cant.

You want the following:

##
TDat <- read.csv("small.txt", sep="\t")
TDat
str(TDat)

See ?read.csv

Regards, Mark.


hannesPretorius wrote:
> 
> Ok i feel pretty stupid.. been trying to read a text file that contains
> heading in the first line in to R but cant. all i need to do is make a
> contour plot for a friend but after weeks i feel like giving up.. i
> included the first few lines of the file.. any help will be great
> 
> Thanks
> 
> Hannes http://www.nabble.com/file/p24777697/small.txt small.txt 
> 

-- 
View this message in context: 
http://www.nabble.com/i%27m-so-stuck-with-text-file-and-contour-plot-tp24777697p24780416.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] i'm so stuck with text file and contour plot

[Mark Difford]

Hi David,

>> I think he may also need to add the header=TRUE argument:

No! The argument header= is not required in this case.

##
> TDat <- read.csv("small.txt", sep="\t") 
> str(TDat[,1:3])
'data.frame':   10 obs. of  3 variables:
 $ Placename: Factor w/ 10 levels "Aankoms","Aapieshoek",..: 1 2 3 4 5 6 7 8
9 10
 $ X_coord  : num  30.9 31.4 31.1 31.4 18.7 ...
 $ Y_coord  : num  -26.2 -27.4 -29 -29 -33.5 ...

?read.csv, sub header:
"If missing, the value is determined from the file format: header is set to
TRUE if and only if the first row contains one fewer field than the number
of columns."

regards, Mark.


David Winsemius wrote:
> 
> I think he may also need to add the header=TRUE argument:
> 
> tdat <- read.csv("http://www.nabble.com/file/p24777697/small.txt";,  
> header=TRUE, sep="\t")
> 
> Note: read.table with those arguments should have worked as well.
> 
> And then use names(tdat) <- c()
> 
> Perhaps along these lines:
> tdnames <- names(tdat)
> tdnames
> 
> #--don't paste--
> [1] "Placename"
> [2] "X_coord"
> [3] "Y_coord"
> [4] "Jan.to.Dec.2006.Stroke.Density.per.sq.km"
> [5] "Jan.to.Dec.2007.Stroke.Density.per.sq.km"
> [6] "Jan.to.Oct.2008.Stroke.Density.per.sq.km"
> [7] "Total.Strokes.per.sq.km.for.Jan.2006.to.Oct.2008"
> ###-
> 
>   names(tdat)[4:7] <- c("Strk.dens.2006", "Strk.dens.2007", "Strk.dens. 
> 2008", "cumStrk.2006_8")
> 
> # cannot use variable names that begin with numbers  
> without special efforts
> tdat   # now can be displayed more economically
> 
> -- 
> David
> 
> On Aug 2, 2009, at 2:10 PM, Mark Difford wrote:
> 
>>
>> Hannes,
>>
>>>> been trying to read a text file that contains heading in the first  
>>>> line
>>>> in to R but cant.
>>
>> You want the following:
>>
>> ##
>> TDat <- read.csv("small.txt", sep="\t")
>> TDat
>> str(TDat)
>>
>> See ?read.csv
>>
>> Regards, Mark.
>>
>>
>> hannesPretorius wrote:
>>>
>>> Ok i feel pretty stupid.. been trying to read a text file that  
>>> contains
>>> heading in the first line in to R but cant. all i need to do is  
>>> make a
>>> contour plot for a friend but after weeks i feel like giving up.. i
>>> included the first few lines of the file.. any help will be great
>>>
>>> Thanks
>>>
>>> Hannes http://www.nabble.com/file/p24777697/small.txt small.txt
>>>
>>
>> -- 
>> View this message in context:
>> http://www.nabble.com/i%27m-so-stuck-with-text-file-and-contour-plot-tp24777697p24780416.html
>> Sent from the R help mailing list archive at Nabble.com.
>>
>> __
>> R-help@r-project.org 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.
> 
> David Winsemius, MD
> Heritage Laboratories
> West Hartford, CT
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/i%27m-so-stuck-with-text-file-and-contour-plot-tp24777697p24786217.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] i'm so stuck with text file and contour plot

[Mark Difford]

And I meant to add, but somehow forgot, that the default for read.csv is
header=TRUE (which is different from read.table, where it is FALSE).

Regards, Mark.


Mark Difford wrote:
> 
> Hi David,
> 
>>> I think he may also need to add the header=TRUE argument:
> 
> No! The argument header= is not required in this case.
> 
> ##
>> TDat <- read.csv("small.txt", sep="\t") 
>> str(TDat[,1:3])
> 'data.frame':   10 obs. of  3 variables:
>  $ Placename: Factor w/ 10 levels "Aankoms","Aapieshoek",..: 1 2 3 4 5 6 7
> 8 9 10
>  $ X_coord  : num  30.9 31.4 31.1 31.4 18.7 ...
>  $ Y_coord  : num  -26.2 -27.4 -29 -29 -33.5 ...
> 
> ?read.csv, sub header:
> "If missing, the value is determined from the file format: header is set
> to TRUE if and only if the first row contains one fewer field than the
> number of columns."
> 
> regards, Mark.
> 
> 
> David Winsemius wrote:
>> 
>> I think he may also need to add the header=TRUE argument:
>> 
>> tdat <- read.csv("http://www.nabble.com/file/p24777697/small.txt";,  
>> header=TRUE, sep="\t")
>> 
>> Note: read.table with those arguments should have worked as well.
>> 
>> And then use names(tdat) <- c()
>> 
>> Perhaps along these lines:
>> tdnames <- names(tdat)
>> tdnames
>> 
>> #--don't paste--
>> [1] "Placename"
>> [2] "X_coord"
>> [3] "Y_coord"
>> [4] "Jan.to.Dec.2006.Stroke.Density.per.sq.km"
>> [5] "Jan.to.Dec.2007.Stroke.Density.per.sq.km"
>> [6] "Jan.to.Oct.2008.Stroke.Density.per.sq.km"
>> [7] "Total.Strokes.per.sq.km.for.Jan.2006.to.Oct.2008"
>> ###-
>> 
>>   names(tdat)[4:7] <- c("Strk.dens.2006", "Strk.dens.2007", "Strk.dens. 
>> 2008", "cumStrk.2006_8")
>> 
>> # cannot use variable names that begin with numbers  
>> without special efforts
>> tdat   # now can be displayed more economically
>> 
>> -- 
>> David
>> 
>> On Aug 2, 2009, at 2:10 PM, Mark Difford wrote:
>> 
>>>
>>> Hannes,
>>>
>>>>> been trying to read a text file that contains heading in the first  
>>>>> line
>>>>> in to R but cant.
>>>
>>> You want the following:
>>>
>>> ##
>>> TDat <- read.csv("small.txt", sep="\t")
>>> TDat
>>> str(TDat)
>>>
>>> See ?read.csv
>>>
>>> Regards, Mark.
>>>
>>>
>>> hannesPretorius wrote:
>>>>
>>>> Ok i feel pretty stupid.. been trying to read a text file that  
>>>> contains
>>>> heading in the first line in to R but cant. all i need to do is  
>>>> make a
>>>> contour plot for a friend but after weeks i feel like giving up.. i
>>>> included the first few lines of the file.. any help will be great
>>>>
>>>> Thanks
>>>>
>>>> Hannes http://www.nabble.com/file/p24777697/small.txt small.txt
>>>>
>>>
>>> -- 
>>> View this message in context:
>>> http://www.nabble.com/i%27m-so-stuck-with-text-file-and-contour-plot-tp24777697p24780416.html
>>> Sent from the R help mailing list archive at Nabble.com.
>>>
>>> __
>>> R-help@r-project.org 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.
>> 
>> David Winsemius, MD
>> Heritage Laboratories
>> West Hartford, CT
>> 
>> __
>> R-help@r-project.org 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.
>> 
>> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/i%27m-so-stuck-with-text-file-and-contour-plot-tp24777697p24786283.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] i'm so stuck with text file and contour plot

[Mark Difford]

Hannes,

>> When I read the entire text file in I get the following message

Then you have not followed the very simple instructions I gave you above,
which I repeat below. Or you have changed small.txt.

##
TDat <- read.csv("small.txt", sep="\t") 
TDat
str(TDat)

Mark.


hannesPretorius wrote:
> 
> When I read the entire text file in I get the following message
> 
>> x <- read.table('c:/small.txt', sep='\t', header=TRUE)
> Warning message:
> number of items read is not a multiple of the number of columns.
> 
> thanks.
> 
> 
> 
> 
> 
> 
> hannesPretorius wrote:
>> 
>> Ok i feel pretty stupid.. been trying to read a text file that contains
>> heading in the first line in to R but cant. all i need to do is make a
>> contour plot for a friend but after weeks i feel like giving up.. i
>> included the first few lines of the file.. any help will be great
>> 
>> Thanks
>> 
>> Hannes http://www.nabble.com/file/p24777697/small.txt small.txt 
>> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/i%27m-so-stuck-with-text-file-and-contour-plot-tp24777697p24824996.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Testing year effect in lm() ***failed first time, sending again

[Mark Difford]

Emmanuel,

>> somewhat incomplete  help pages : what in h*ll are valid arguments to
>> mcp() beyond "Tukey" ??? Curently, you'll have to dig in the source to
>> learn that...).

Not so: they are clearly stated in ?contrMat.

Regards, Mark.


Emmanuel Charpentier-3 wrote:
> 
> Le jeudi 30 juillet 2009 à 16:41 -0600, Mark Na a écrit :
>> Dear R-helpers,
>> 
>> I have a linear model with a year effect (year is coded as a factor),
>> i.e.
>> the parameter estimates for each level of my year variable have
>> significant
>> P values (see some output below) and I am interested in testing:
>> 
>> a) the overall effect of year;
>> b) the significance of each year vis-a-vis every other year (the model
>> output only tests each year against the baseline year).
> 
> install.packges("multcomp")
> help.start()
> 
> and use the vignettes ! They're good (and are a good complement to
> somewhat incomplete  help pages : what in h*ll are valid arguments to
> mcp() beyond "Tukey" ??? Curently, you'll have to dig in the source to
> learn that...).
> 
> HTH
> 
>   Emmanuel Charpentier
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Testing-year-effect-in-lm%28%29-***failed-first-time%2C-sending-again-tp24748526p24832337.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] PCA or CA

[Mark Difford]

Hi Paul,

>> I have a data set for which PCA based between group analysis (BGA) gives 
>> significant results but CA-BGA does not. 
>> I am having difficulty finding a reliable method for deciding which
>> ordination 
>> technique is most appropriate.

Reliability really comes down to you thinking about and properly defining
what _information_ you want to extract from your data set, which we know
nothing about. PCA and CA are fundamentally different. The classical use of
CA lies in the analysis of count-data (contingency tables), for which it
remains a brilliant method. It is also widely applied to analyzing normal n
x p matrices of the type usually analyzed by PCA. A doubly-centered PCA
would get you close to a CA of the normal n x p matrix (i.e. of an analysis
of the same matrix).

This is a biggish area, so grab your specs, and perhaps start with Jolliffe
(PCA) and Benzecri/Greenacre (CA).

Regards, Mark.


Paul Dennis-3 wrote:
> 
> 
> Dear all
> 
> I have a data set for which PCA based between group analysis (BGA) gives
> significant results but CA-BGA does not.
> 
> I am having difficulty finding a reliable method for deciding which
> ordination technique is most appropriate. 
> 
> I have been told to do a 1 table CA and if the 1st axis is>2 units go for
> CA if not then PCA.
> 
> Another approach is that described in the Canoco manual - perform DCA and
> then look at the length of the axes.  I used decorana in vegan and it
> gives axis lengths.  I assume that these are measured in SD units. Anyway
> the manual say if the axis length is <3 go for PCA,>4 use CA and if
> intermediate use either. 
> 
> Are either of these approaches good/valid/recommended or is there a better
> method?
> 
> Thanks
> 
> Paul  
> 
> _
> Get the best of MSN on your mobile
> http://clk.atdmt.com/UKM/go/147991039/direct/01/
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/PCA-or-CA-tp25668667p25670451.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] trouble with html() in Hmisc

[Mark Difford]

Hi Liviu,

>> > tmp <- latex(.object, cdec=c(2,2), title="") 
>> > class(tmp) 
>> [1] "latex" 
>> > html(tmp) 
>> /tmp/RtmprfPwzw/file7e72f7a7.tex:9: Warning: Command not found:
>> \tabularnewline 
>> Giving up command: \...@hevea@amper 
>> /tmp/RtmprfPwzw/file7e72f7a7.tex:11: Error while reading LaTeX: 
>> This array/tabular column has no specification

This has nothing to do with Hmisc or hevea. In the header/preample of all
LyX files you will, for instance, find the following:

## -
%% LyX 1.6.2 created this file.  For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[english]{article}
.
.
.
%% LyX specific LaTeX commands.
%% Because html converters don't know tabularnewline
\providecommand{\tabularnewline}{\\}

Regards, Mark.


Liviu Andronic wrote:
> 
> Dear all
> On my system html() conversion of a `latex()' object fails. Follows a
> dummy example:
>> require(Hmisc)
>> data(Angell)
>> .object <- cor(Angell[,1:2], use="complete.obs")
>> tmp <- latex(.object, cdec=c(2,2), title="")
>> class(tmp)
> [1] "latex"
>> html(tmp)
> /tmp/RtmprfPwzw/file7e72f7a7.tex:9: Warning: Command not found:
> \tabularnewline
> Giving up command: \...@hevea@amper
> /tmp/RtmprfPwzw/file7e72f7a7.tex:11: Error while reading LaTeX:
>   This array/tabular column has no specification
> Adios
> 
> 
> I have hevea-1.10 installed on Debian (according to the help page, it
> performs the conversion). Attached is teh produced .tex file. Is this
> a bug or would there be a way to work around this behaviour?
> Thank you
> Liviu
> 
> 
>> sessionInfo()
> R version 2.9.2 (2009-08-24)
> x86_64-pc-linux-gnu
> 
> locale:
> LC_CTYPE=en_GB.UTF-8;LC_NUMERIC=C;LC_TIME=en_GB.UTF-8;LC_COLLATE=en_GB.UTF-8;LC_MONETARY=C;LC_MESSAGES=en_GB.UTF-8;LC_PAPER=en_GB.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_GB.UTF-8;LC_IDENTIFICATION=C
> 
> attached base packages:
> [1] datasets  grDevices tcltk splines   graphics  utils stats
> [8] methods   base
> 
> other attached packages:
>  [1] fortunes_1.3-6   RcmdrPlugin.Export_0.3-0 Rcmdr_1.5-2
>  [4] car_1.2-15   hints_1.0.1-1boot_1.2-39
>  [7] relimp_1.0-1 xtable_1.5-5 Hmisc_3.7-0
> [10] survival_2.35-7
> 
> loaded via a namespace (and not attached):
> [1] cluster_1.12.0  grid_2.9.2  lattice_0.17-25 tools_2.9.2
>> system("uname -a")
> Linux debian-liv 2.6.30-1-amd64 #1 SMP Sat Aug 15 18:09:19 UTC 2009
> x86_64 GNU/Linux
> 
> 
> 
> 
> 
> -- 
> Do you know how to read?
> http://www.alienetworks.com/srtest.cfm
> Do you know how to write?
> http://garbl.home.comcast.net/~garbl/stylemanual/e.htm#e-mail
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/trouble-with-html%28%29-in-Hmisc-tp25721612p25728656.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] trouble with html() in Hmisc

[Mark Difford]

Liviu,

>> Does this mean that Hmisc produces LyX code (as opposed to LaTeX code)?

No. Sorry for confusing you: it means that html does not know what
tabularnewline means, it cannot interpret it. My reply showed where the
problem lies and how to fix it.

You need to add the following command to the preamble of the *.tex file:
\providecommand{\tabularnewline}{\\}

Regards, Mark.


Liviu Andronic wrote:
> 
> Hello
> 
> On 10/3/09, Mark Difford  wrote:
>> This has nothing to do with Hmisc or hevea.
>>
> Although I have LyX installed, I don't quite understand where LyX
> comes into play. The R code in the original e-mail takes a table-like
> object and transforms it into LaTeX; then html() calls hevea to conver
> .tex to .html and opens a browser to display the result.
> 
> 
>>  %% LyX specific LaTeX commands.
>>  %% Because html converters don't know tabularnewline
>>  \providecommand{\tabularnewline}{\\}
>>
> 
> If I print the `latex' object from the previous e-mail, I get the
> offending "tabularnewline".
>> data(Angell)
>> .object <- cor(Angell[,1:2], use="complete.obs")
>> tmp <- latex(.object, cdec=c(2,2), file="", title="")
> % latex.default(.object, cdec = c(2, 2), file = "", title = "")
> %
> \begin{table}[!tbp]
>  \begin{center}
>  \begin{tabular}{lrr}\hline\hline
> \multicolumn{1}{l}{}&\multicolumn{1}{c}{moral}&\multicolumn{1}{c}{hetero}\tabularnewline
> \hline
> moral&$ 1.00$&$-0.59$\tabularnewline
> hetero&$-0.59$&$ 1.00$\tabularnewline
> \hline
> \end{tabular}
> 
> \end{center}
> 
> \end{table}
> 
> 
> Does this mean that Hmisc produces LyX code (as opposed to LaTeX code)?
> Liviu
> 
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/trouble-with-html%28%29-in-Hmisc-tp25721612p25731240.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Trendline for a subset of data

[Mark Difford]

Hi Steve,

>> However, I am finding that ... the trendline ... continues to run beyond
>> this data segment 
>> and continues until it intersects the vertical axes at each side of the
>> plot.

Your "best" option is probably Prof. Fox's reg.line function in package car.

##
library(car)
?reg.line
reg.line

Regards, Mark.


smurray444 wrote:
> 
> 
> Dear all,
> 
> I am using abline(lm ...) to insert a linear trendline through a portion
> of my data (e.g. dataset[,36:45]). However, I am finding that whilst the
> trendline is correctly displayed and representative of the data portion
> I've chosen, the line continues to run beyond this data segment and
> continues until it intersects the vertical axes at each side of the plot.
> 
> How do I display the line so that it only runs between point 36 and 45 (as
> shown in the example above) as doesn't continue to display a line
> throughout the rest of the plot space?
> 
> Many thanks,
> 
> Steve
> 
> _
> View your other email accounts from your Hotmail inbox. Add them now.
> http://clk.atdmt.com/UKM/go/167688463/direct/01/
> __
> R-help@r-project.org 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Trendline-for-a-subset-of-data-tp25818425p25821972.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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.