Re: [R] Rcmdr scatter3d

[John Fox]

Dear chip,

When there is groups variable specified, scatter3d() colours the
regression surface, points, and residuals for each group according to
the colours specified in surface.col. (Setting surface=FALSE suppresses
the regression surfaces and residuals.)

You can save the rgl graph that scatter3d() produces as a bitmapped png
graphic [e.g., via Graphs -> 3D graph -> Save graph to file, which uses
rgl.snapshot()]. There is also an rgl.postscript() command, which
supports some vector-graphics formats, but I've been unable to use it
successfully.

I hope this helps,
 John

On Tue, 4 Sep 2007 13:43:44 -0700 (PDT)
 array chip <[EMAIL PROTECTED]> wrote:
> Hi, I am using the scatter3d function in Rcmdr to plot
> the first 3 principal components, I have a grouping
> variable of 2 groups, and tried to plot points with
> different colors, somehow I couldn't change the
> default colors of the 2 groups (blue and green)by
> using option points.col=c('red','blue'), what's the
> problem here?
> 
> scatter3d(all.pca$x[,2],all.pca$x[,3],all.pca$x[,1],
> surface=FALSE, residuals=TRUE, bg="white",
> axis.scales=F, grid=F, ellipsoid=F, xlab='PCA
> 2',ylab='PCA 3', zlab='PCA
> 1',sphere.size=1.5,groups=as.factor(c(rep(1,100),rep(2,50)))
> ,point.col=c('red','blue'))
> 
> I am also wondering if I can have a copy of the image
> in high resolution, just like copying a regular R plot
> in "Metafile"? 
> 
> Thanks
> 
> 
> 
>
>

> 
> that gives answers, not web links.
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Row-Echelon Form

[John Fox]

etric")
D <- rep(0, n)
L <- diag(n)
i <- 2:n
D[1] <- X[1, 1]
if (abs(D[1]) < tol) stop("matrix is numerically singular")
L[i, 1] <- X[i, 1]/D[1]
for (j in 2:(n - 1)){
k <- 1:(j - 1)
D[j] <- X[j, j] - sum((L[j, k]^2) * D[k])
if (abs(D[j]) < tol) stop("matrix is numerically singular")
i <- (j + 1):n
    L[i, j] <- (X[i, j] -
colSums(L[j, k] * t(L[i, k, drop=FALSE]) *
D[k]))/D[j]
}
k <- 1:(n - 1)
D[n] <- X[n, n] - sum((L[n, k]^2) * D[k])
if (abs(D[n]) < tol) stop("matrix is numerically singular")
L %*% diag(sqrt(D))
}




John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Peter Danenberg
> Sent: Monday, September 03, 2007 2:17 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Row-Echelon Form
> 
> I was looking for an R-package that would reduce matrices to 
> row-echelon form, but Google was not my friend; any leads?
> 
> If not, I wonder if the problem could be expressed in terms 
> of constraint satisfaction...
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] An issue with White's test in Anova

[John Fox]

Dear David,

You've found a bug in Anova() for linear models that was introduced some
time ago when the linear.hypothesis() function [which is used by Anova()]
was modified not to report sums of squares for "White-adjusted" tests. I
have to think about should be done in this case.

Sorry,
 John

--------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of David Kaplan
> Sent: Friday, August 24, 2007 3:55 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] An issue with White's test in Anova
> 
> Hi all,
> 
> I'm running White's test to correct for non-constant error 
> variance and I am using the Anova function in the package 
> CAR.  My command structure is
> 
>  > Anova(scireg3, white.adjust="hc3")
> 
> where scireg3 is an object from lm.
> 
> I get the message
> 
> "Error in SS[i] <- SS.term(names[i]) : nothing to replace with"
> 
> What does this mean and how do I fix it.
> 
> Thanks in advance.
> 
> David
> 
> 
> 
> 
> --
> ==
> =
> David Kaplan, Ph.D.
> Professor
> Department of Educational Psychology
> University of Wisconsin - Madison
> Educational Sciences, Room, 1061
> 1025 W. Johnson Street
> Madison, WI 53706
> 
> email: [EMAIL PROTECTED]
> homepage: 
> http://www.education.wisc.edu/edpsych/facstaff/kaplan/kaplan.htm
> Phone: 608-262-0836
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] perception of graphical data

[John Fox]

Dear Richard,

Though slightly dated, the following article is a nice summary of the
literature on graphical perception:

Lewandowsky, S & Spence, I. (1989) The perception of statistical
graphs. Sociological Methods and Research, 18, 200-242.

I hope this helps,
 John

On Fri, 24 Aug 2007 13:30:56 -0400
 "Yeh, Richard C" <[EMAIL PROTECTED]> wrote:
> Hello,
> 
> I apologize that this is off-topic.  I am seeking information on
> perception of graphical data, in an effort to improve the plots I
> produce.  Would anyone point me to literature reviews in this area?
>  (Or
> keywords to try on google?)  Is this located somewhere near cognitive
> science, psychology, human factors research?
> 
> For example, some specific questions I have are:
> 
> I recall as a child when I first saw a map where the areas of the
> containers (geographical states) were drawn as rectangles,
> proportional
> to a quantity other than land area.  Does anyone know of an algorithm
> for drawing such maps?  Would anyone know of a journal or reference
> where I can find studies on whether subjects reading these maps can
> accurately assess the meaning of the different areas, as [some of us]
> can assess different heights on a bar graph?  (What about areas in
> bar
> graphs with non-uniform widths?)
> 
> Scatter plots of microarray data often attempt to represent thousands
> or
> tens of thousands of points, but all I read from them are density and
> distribution --- the gene names cannot be shown.  At what point,
> would a
> sunflowerplot-like display or a smooth gradient be better?  When two
> data points drawn as 50% gray disks are small and tangent, are they
> perceptually equivalent to a single, 100% black disk?  Or a 50% gray
> disk with twice the area?  What problems are known about plotting
> with
> disks --- do viewers use the area or the diameter (or neither) to
> gauge
> weight?
> 
> 
> As you can tell, I'm a non-expert, mixing issues of data
> interpretation,
> visual perception, graphic representation.  Previously, I didn't have
> the flexibility of R's graphics, so I didn't need to think so much.
> I've read some of Edward S. Tufte's books, but found them more
> qualitative than quantitative.
> 
> Thanks!
> 
> Richard
> 
> 212-933-3305 / [EMAIL PROTECTED]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] SEM for categorical data

[John Fox]

Dear Upasna,

I apologize for responding so late, but I was out of town when you posted
your query to r-help.

The sem() function in the sem package can handle dichotomous and ordered
categorical observed variables via tetrachoric, polychoric, biserial, and
polyserial correlations computed by the polycor package. See in particular
the hetcor() function in that package, which can compute "heterogeneous"
correlation matrices. Standard errors for parameter estimates in the SEM can
then be computed by bootstrapping; see ?boot.sem in the sem package for an
example with a confirmatory factor-analysis model for ordinal observed
variables.

I hope that this helps,
 John

--- original message ---

Hi

I am looking for a structural equation modeling package in R which can be
used for categorical data. Is anyone aware of the existence of such a
package? Would appreciate any help on this.

Thank you
Upasna

-- 
-
Upasna Sharma
Research Scholar
Shailesh J. Mehta School of Management,
Indian Institute of Technology, Bombay
Powai, Mumbai - 400076, India

--------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Rcmdr window border lost

[John Fox]

Dear Andy and Peter,

My apologies for chiming in so late, but I was out of town and just caught
up with R-help today (and an impressive amount of material accumulates on
R-help in a couple of weeks!).

Although I develop the Rcmdr package on a Windows system, I test it under
Linux as well -- most recently on Ubuntu 7.04 with R 2.5.0. Curiously, the
windows for Rcmdr 1.3-5 look normal to me, as I just verified, so I'm unable
to reproduce the problem, and wonder whether anyone else has seen it.

As well, setting options(Rmcdr=list(suppress.X11.warnings=TRUE)) should
suppress the X11 warnings that appear on some systems, but I thought that
this problem disappeared with R 2.4.0. See ?Commander.

Regards,
 John

--- original messages --

Andy Weller wrote:
> OK, I tried completely removing and reinstalling R, but this has not 
> worked - I am still missing window borders for Rcmdr. I am certain that 
> everything is installed correctly and that all dependencies are met - 
> there must be something trivial I am missing?!
>
> Thanks in advance, Andy
>
> Andy Weller wrote:
>   
>> Dear all,
>>
>> I have recently lost my Rcmdr window borders (all my other programs have 
>> borders)! I am unsure of what I have done, although I have recently 
>> update.packages() in R... How can I reclaim them?
>>
>> I am using:
>> Ubuntu Linux (Feisty)
>> R version 2.5.1
>> R Commander Version 1.3-5
>>
>> 
This sort of behaviour is usually the fault of the window manager, not 
R/Rcmdr/tcltk. It's the WM's job to supply the various window 
decorations on a new window, so either it never got told that there was 
a window, or it somehow got into a confused state. Did you try 
restarting the WM (i.e., log out/in or reboot)? And which WM are we 
talking about?

Same combination works fine on Fedora 7, except for a load of messages 
saying

Warning: X11 protocol error: BadWindow (invalid Window parameter)


>> I have deleted the folder: /usr/local/lib/R/site-library/Rcmdr and 
>> reinstalled Rcmdr with: install.packages("Rcmdr", dep=TRUE)
>>
>> This has not solved my problem though.
>>
>> Maybe I need to reinstall something else as well?
>>
>> Thanks in advance, Andy
>> 


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] termplot with uniform y-limits

[John Fox]

Dear Bill,

Functions in the effects package are somewhat similar to termplot, and
allow you to specify the y-axis limits. For example, modify the last
line of the last example in example(all.effects) to 

plot(eff.pres, ask=FALSE, ylim=c(10, 70))

I hope this helps,
 John

--- original message ---

Does anyone have, or has anyone ever considered making, a version of
'termplot' that allows the user to specify that all plots should have
the same y-limits?

This seems a natural thing to ask for, as the plots share a y-scale.
 If
you don't have the same y-axes you can easily misread the comparative
contributions of the different components. 

Notes: the current version of termplot does not allow the user to
specify ylim.  I checked.

  the plot tools that come with mgcv do this by default.  Thanks
Simon.

--------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] SEM model fit

[John Fox]

prev, gam5,   0.8 
prev  -> enf,  gam6,   0.4

sem.enf.rq <- sem(ram = mdl.rq, S = hcor(dx),  N = nrow(dx), obs.v =
names(dx), raw = F, fixed = names(dx)[4:6], par.size = 's', maxiter =
1e3,
analytic = F, gradtol = 1e-10)  ##set raw to False
summary(obj = sem.enf.rq, dig = 3, conf = 0.9) 

Respectfully,

Frank Lawrence


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Inverse BoxCox transformation

[John Fox]

Dear Des,

The following should do the trick:

invBoxCox <- function(x, lambda)
if (lambda == 0) exp(x) else (lambda*x + 1)^(1/lambda)

I hope this helps,
 John

----
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Des Callaghan
> Sent: Monday, June 18, 2007 3:33 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Inverse BoxCox transformation
> 
> Hi,
>  
> I can't seem to find a function in R that will reverse a 
> BoxCox transformation. Can somebody help me locate one 
> please? Thanks in advance.
>  
> Best wishes,
> Des
>  
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] linear hypothesis test in gls model

[John Fox]

Dear Tine,

linear.hypothesis() currently has no method specifically for gls objects,
and so this usage invokes the default method. I'm not sure off-hand what's
appropriate for an F-test in this context (and indeed why the default test
is inappropriate). Can you describe the correct test or supply a reference?
I suspect that it shouldn't be hard to write a linear.hypothesis method for
gls objects that fixes up the result returned by linear.hypothesis.default. 

You might take a look at car:::linear.hypothesis.default to see that it does
-- the computations are pretty straightforward.

I hope this helps,
 John 

--------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Tine Huyghe
> Sent: Saturday, June 16, 2007 4:19 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] linear hypothesis test in gls model
> 
> Dear all,
> 
> For analysis of a longitudinal data set with fixed 
> measurement in time I built a gls model (nlme). For testing 
> hypotheses in this model I used the linear.hypothesis 
> function from the car package. A check with the results 
> obtained in SAS proc MIXED with a repeated statement revealed 
> an inconsistency in the results. The problem can be that the 
> linear.hypothesis function (1) only gives the asymptotic chi 
> square test and/or (2) only uses the residual error. Is there 
> another solution to testing linear hypotheses in a gls model?
> 
> Thanks in advance
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Refactor all factors in a data frame

[John Fox]

Dear Hilmar,

You could use something like

DF <- as.data.frame(lapply(DF, function (x) if (is.factor(x)) factor(x) else
x))

Where DF is the data frame.

I hope this helps,
 John

----
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Hilmar Berger
> Sent: Tuesday, June 05, 2007 8:20 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Refactor all factors in a data frame
> 
> Hi all,
> 
> Assume I have a data frame with numerical and factor 
> variables that I got through merging various other data 
> frames and subsetting the resulting data frame afterwards. 
> The number levels of the factors seem to be the same as in 
> the original data frames, probably because subset() calls 
> [.factor without drop = TRUE (that's what I gather from 
> scanning the mailing lists).
> 
> I wonder if there is a easy way to refactor all factors in 
> the data frame at once. I noted that fix(data_frame) does the 
> trick, however, this needs user interaction, which I'd like 
> to avoid. Subsequent write.table / read.table would be 
> another option but I'm not sure if R can guess the 
> factor/char/numeric-type correctly when reading the table.
> 
> So, is there any way in drop the unused factor levels from 
> *all* factors of a data frame without import/export ?
> 
> Thanks in advance,
> Hilmar
> 
> -- 
> 
> Hilmar Berger
> Studienkoordinator
> Institut für medizinische Informatik, Statistik und 
> Epidemiologie Universität Leipzig Härtelstr. 16-18
> D-04107 Leipzig
> 
> Tel. +49 341 97 16 101
> Fax. +49 341 97 16 109
> email: [EMAIL PROTECTED]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] [R-pkgs] Rcmdr 1.3-0 and RcmdrPlugins.TeachingDemos

[John Fox]

I'd like to announce a new version, 1.3-0, of the Rcmdr package. The Rcmdr
package provides a basic-statistics graphical user interface (GUI) to R. 

Beyond small changes and additions, this new version of the package makes
provision for "plug-ins" that permit extension of the Rcmdr GUI without
altering and rebuilding the Rcmdr source package or modifying the installed
package. An R Commander plug-in is an ordinary R package that (1) provides
extensions to the R Commander menus is a file named menus.txt located in the
package's etc directory; (2) provides call-back functions required by these
menus; and (3) in optional Log-Exceptions: and Models: fields in the
package's DESCRIPTION file, augments respectively the list of functions for
which printed output is suppressed and the list of model objects recognized
by the R Commander. The menus provided by a plug-in package are merged with
the standard Commander menus. 

Plug-in packages given in the R Commander plugins option (see ?Commander)
are automatically loaded when the Commander starts up. Plug-in packages may
also be loaded via the Commander "Tools -> Load Rcmdr plug-in(s)" menu; a
restart of the Commander is required to install the new menus. Finally,
loading a plug-in package when the Rcmdr is not loaded will load the Rcmdr
and activate the plug-in. 

An illustrative R Commander plug-in package, RcmdrPlugin.TeachingDemos
(providing a GUI to some of Greg Snow's TeachingDemos package), is now
available on CRAN. (I suggest using this naming convention -- RcmdrPlugin.*
-- so that plug-in packages will sort immediately below the Rcmdr package on
CRAN. This assumes, of course, that other people will be interested in
creating Rcmdr plugins!)

Because this is a new feature of the Rcmdr, feedback and suggestions would
be appreciated.

I'd like to acknowledge Richard Heiberger's suggestions for the design of
this plug-in facility.

John


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox

___
R-packages mailing list
[EMAIL PROTECTED]
https://stat.ethz.ch/mailman/listinfo/r-packages

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] GUI component Margin on tkcanvas, tkframe or tktoplevel

[John Fox]

Dear Hao,


> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Hao Liu
> Sent: Thursday, May 24, 2007 3:53 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] GUI component Margin on tkcanvas, tkframe or tktoplevel
> 
> Dear gurus:
> 
> I have a question on how to configure the margin layout on 
> tcl/tk GUI objects like tkcanvas, tkframe or tktoplevel.
> For example, if I want to leave a larger margin on the left 
> side of the GUI container, which one should I configure, the 
> toplevel or tkcanvas or tkframe?
> 
> top<-tktoplevel()
> canvas<-tkcanvas(top, relief=, borderwidth=...)

I've used the borderwidth argument (which, I believe, affects all of the
borders) in several places, including to tkframe(). You can try it out and
see whether you get the effect you want.

> 
> I also find the documentation for all the possible arguments 
> that could be passed to tktoplevel, tkframe, tkcanvas, tk.wm, 
> tkconfigure is not good. I wonder if there are good resource 
> for it... at the very least, I know I need to look at what 
> arguments that tck/tk takes... but I don't have time to go that far.
> 

I think that the tcltk package documentation is meant to be used alongside
Tcl/Tk documentation. The most useful source that I've found is Welch et
al., Practical Programming in Tcl and Tk.

I hope this helps,
 John

> Thanks
> Hao
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Selecting complementary colours

[John Fox]

Dear Thomas,

This seems simpler than the solution that I used, so I'll give it a
try.

Thanks,
 John

On Tue, 22 May 2007 09:01:01 -0700 (PDT)
 Thomas Lumley <[EMAIL PROTECTED]> wrote:
> On Mon, 21 May 2007, John Fox wrote:
> >
> > In retrospect, I didn't specify the problem clearly: What I want to
> be able
> > to do is to place text on a background of arbitrary (but known RGB)
> colour
> > so that the text is legible. I guess that this is better described
> as a
> > "contrasting" than a "complementary" colour.
> 
> Since luminance contrasts are necessary and sufficient for readable
> text, you could use white for dark colors and black for light colors.
> 
> Luminance is roughly proportional to  0.2*(R^2.4)+0.6*(G^2.4),
> suggesting something like
> 
> lightdark<-function (color)
> {
>  rgb <- col2rgb(color)/255
>  L <- c(0.2, 0.6, 0) %*% rgb
>  ifelse(L >= 0.2, "#60", "#A0")
> }
> 
> This uses a pale yellow for dark backgrounds and a dark blue for
> light backgrounds, and it seems to work reasonably well.
> 
>   -thomas


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Boostrap p-value in regression [indirectly related to R]

[John Fox]

Dear Wolfgang,

I agree that it's preferable to compute the two-sided p-value without
assuming symmetry. Another, equivalent, way of thinking about this is to use
t^2 for the two-sided test in place of t.

BTW, the formula used in my appendix (for the one-sided p-value) is from
Davison and Hinkley, I believe, and differs trivially from the one in Efron
and Tibshirani.

Regards,
 John


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of 
> Viechtbauer Wolfgang (STAT)
> Sent: Monday, May 21, 2007 10:41 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Boostrap p-value in regression [indirectly related to R]
> 
> Hello All,
> 
> Despite my preference for reporting confidence intervals, I 
> need to obtain a p-value for a hypothesis test in the context 
> of regression using bootstrapping. I have read John Fox's 
> chapter on bootstrapping regression models and have consulted 
> Efron & Tibshirani's "An Introduction to the Bootstrap" but I 
> just wanted to ask the experts here for some feedback to make 
> sure that I am not doing something wrong.
> 
> Let's take a simplified example where the model includes one 
> independent variable and the idea is to test H0: beta1 = 0 
> versus Ha: beta1 != 0.
> 
> 
> 
> ### generate some sample data
> 
> n  <- 50
> xi <- runif(n, min=1, max=5)
> yi <- 0 + 0.2 * xi + rnorm(n, mean=0, sd=1)
> 
> ### fit simple regression model
> 
> mod <- lm(yi ~ xi)
> summary(mod)
> b1  <- coef(mod)[2]
> t1  <- coef(mod)[2] / coef(summary(mod))[2,2]
> 
> ### 1000 bootstrap replications using (X,Y)-pair resampling
> 
> t1.star <- rep(NA,1000)
> 
> for (i in 1:1000) {
> 
>   ids<- sample(1:n, replace=TRUE)
>   newyi  <- yi[ids]
>   newxi  <- xi[ids]  
>   mod<- lm(newyi ~ newxi)
>   t1.star[i] <- ( coef(mod)[2] - b1) / coef(summary(mod))[2,2]
> 
> }
> 
> ### get bootstrap p-value
> 
> hist(t1.star, nclass=40)
> abline(v=t1, lwd=3)
> abline(v=-1*t1, lwd=3)
> 2 * mean( t1.star > abs(t1) )
> 
> 
> 
> As suggested in the chapter on bootstrapping regression 
> models by John Fox, the bootstrap p-value is 2 times the 
> proportion of bootstrap t-values (with b1 subtracted so that 
> we get the distribution under H0) larger than the absolute 
> value of the actual t-value observed in the data. 
> 
> Doesn't this assume that the bootstrap sampling distribution 
> is symmetric? And if yes, would it then not be more reasonable to
> calculate:
> 
> mean( abs(t1.star) > abs(t1) )
> 
> or in words: the number of bootstrap t-values that are more 
> extreme on either side of the bootstrap distribution than the 
> actual t-value observed?
> 
> Any suggestions or comments would be appreciated!
> 
> --
> Wolfgang Viechtbauer
>  Department of Methodology and Statistics  University of 
> Maastricht, The Netherlands  http://www.wvbauer.com
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Selecting complementary colours

[John Fox]

Dear Chuck,

This solution works reasonably well for me. Although it occasionally
produces an error, I'm able to trap that.

Thank you -- and to everyone else who responded.

John


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: Charles C. Berry [mailto:[EMAIL PROTECTED] 
> Sent: Monday, May 21, 2007 8:49 PM
> To: John Fox
> Cc: 'Deepayan Sarkar'; r-help@stat.math.ethz.ch
> Subject: Re: [R] Selecting complementary colours
> 
> On Mon, 21 May 2007, John Fox wrote:
> 
> > Dear Deepayan,
> >
> > I actually thought of the equivalent of this approach, but 
> it doesn't 
> > do quite what I want.
> >
> > In retrospect, I didn't specify the problem clearly: What I 
> want to be 
> > able to do is to place text on a background of arbitrary (but known 
> > RGB) colour so that the text is legible. I guess that this 
> is better 
> > described as a "contrasting" than a "complementary" colour.
> 
> John,
> 
> There may be no unique solution. (For gray, for example.)
> 
> I am not sure (in terms of color theory) that maximizing in 
> rgb space really is the right thing to do, but perhaps this 
> will help you:
> 
> > cval <- function(x,y) -sum((x-y)^2)
> > contrasting <- function(x) 
> > 
> optim(runif(3,0,255),cval,lower=0,upper=255,method="L-BFGS-B",y=x)$par
> > do.call(rgb,as.list(contrasting(col2rgb("gray"))/255))
> [1] "#00"
> > do.call(rgb,as.list(contrasting(col2rgb("gray"))/255))
> [1] "#FF"
> > do.call(rgb,as.list(contrasting(col2rgb("pink"))/255))
> [1] "#00FF00"
> 
> Regards,
> 
> Chuck
> 
> >
> 
> >
> > Your solution, for example breaks down for grays:
> >
> >> mycol <- "#88"
> >> do.call(rgb, as.list(1 - col2rgb(mycol) / 255))
> > [1] "#77"
> >
> > Thank you for the suggestion.
> >
> > John
> >
> > ----
> > John Fox, Professor
> > Department of Sociology
> > McMaster University
> > Hamilton, Ontario
> > Canada L8S 4M4
> > 905-525-9140x23604
> > http://socserv.mcmaster.ca/jfox
> > 
> >
> >> -Original Message-
> >> From: [EMAIL PROTECTED] 
> >> [mailto:[EMAIL PROTECTED] On Behalf Of Deepayan 
> >> Sarkar
> >> Sent: Monday, May 21, 2007 6:45 PM
> >> To: John Fox
> >> Cc: r-help@stat.math.ethz.ch
> >> Subject: Re: [R] Selecting complementary colours
> >>
> >> On 5/21/07, John Fox <[EMAIL PROTECTED]> wrote:
> >>> Dear r-helpers,
> >>>
> >>> I wonder whether, given the "#rrggbb" representation of a colour, 
> >>> there is a simple way to select the complementary colour,
> >> also expressed as a "#rrggbb"
> >>> string.
> >>>
> >>> Any suggestions would be appreciated.
> >>
> >> You want rgb2col. The following should work for any standard color
> >> specification:
> >>
> >>> mycol = "royalblue"
> >>> do.call(rgb, as.list(1 - col2rgb(mycol) / 255))
> >> [1] "#BE961E"
> >>
> >> -Deepayan
> >>
> >> __
> >> R-help@stat.math.ethz.ch mailing list 
> >> https://stat.ethz.ch/mailman/listinfo/r-help
> >> PLEASE do read the posting guide
> >> http://www.R-project.org/posting-guide.html
> >> and provide commented, minimal, self-contained, reproducible code.
> >>
> >
> > __
> > R-help@stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide 
> > http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
> 
> Charles C. Berry(858) 534-2098
>   Dept of 
> Family/Preventive Medicine
> E mailto:[EMAIL PROTECTED] UC San Diego
> http://biostat.ucsd.edu/~cberry/ La Jolla, San Diego 
> 92093-0901
> 
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Selecting complementary colours

[John Fox]

Dear Achim,

As I mentioned in my response to Deepayan's suggestion, I didn't specify the
original problem clearly: The object is to get contrasting colours, so that
when one is plotted over the other, the two will be readily distinguishable.
Your suggestions don't do this for neutral colours:

> x <- "#88"
> y_hcl <- as(hex2RGB(x), "polarLUV")
> [EMAIL PROTECTED], "H"] <- [EMAIL PROTECTED], "H"] + 180
> hex(y_hcl)
[1] "#88"
 
> y_hsv <- as(hex2RGB(x), "HSV")
> [EMAIL PROTECTED], "H"] <- [EMAIL PROTECTED], "H"] + 180
> hex(y_hsv)
[1] "#88"

Thank you for trying to help.

John


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
-------- 

> -Original Message-
> From: Achim Zeileis [mailto:[EMAIL PROTECTED] 
> Sent: Monday, May 21, 2007 7:07 PM
> To: John Fox
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] Selecting complementary colours
> 
> On Mon, 21 May 2007, John Fox wrote:
> 
> > Dear r-helpers,
> >
> > I wonder whether, given the "#rrggbb" representation of a colour, 
> > there is a simple way to select the complementary colour, 
> also expressed as a "#rrggbb"
> > string.
> 
> Is the complementary color uniquely defined? My understanding 
> is that you can take opposite colors on a color wheel, but 
> there are of course various color wheels available. With 
> "colorspace" you can experiment with this,
> e.g.:
>   x <- "#81A9D0"
>   y_hcl <- as(hex2RGB(x), "polarLUV")
>   [EMAIL PROTECTED], "H"] <- [EMAIL PROTECTED], "H"] + 180
>   y_hcl <- hex(y_hcl)
> which is a bit more balanced than
>   y_hsv <- as(hex2RGB(x), "HSV")
>   [EMAIL PROTECTED], "H"] <- [EMAIL PROTECTED], "H"] + 180
>   y_hsv <- hex(y_hsv)
> 
> hth,
> Z
> 
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Selecting complementary colours

[John Fox]

Dear Deepayan,

I actually thought of the equivalent of this approach, but it doesn't do
quite what I want.

In retrospect, I didn't specify the problem clearly: What I want to be able
to do is to place text on a background of arbitrary (but known RGB) colour
so that the text is legible. I guess that this is better described as a
"contrasting" than a "complementary" colour.

Your solution, for example breaks down for grays:

> mycol <- "#88"
> do.call(rgb, as.list(1 - col2rgb(mycol) / 255))
[1] "#77"

Thank you for the suggestion.

John


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Deepayan Sarkar
> Sent: Monday, May 21, 2007 6:45 PM
> To: John Fox
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] Selecting complementary colours
> 
> On 5/21/07, John Fox <[EMAIL PROTECTED]> wrote:
> > Dear r-helpers,
> >
> > I wonder whether, given the "#rrggbb" representation of a colour, 
> > there is a simple way to select the complementary colour, 
> also expressed as a "#rrggbb"
> > string.
> >
> > Any suggestions would be appreciated.
> 
> You want rgb2col. The following should work for any standard color
> specification:
> 
> > mycol = "royalblue"
> > do.call(rgb, as.list(1 - col2rgb(mycol) / 255))
> [1] "#BE961E"
> 
> -Deepayan
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Selecting complementary colours

[John Fox]

Dear r-helpers,

I wonder whether, given the "#rrggbb" representation of a colour, there is a
simple way to select the complementary colour, also expressed as a "#rrggbb"
string.

Any suggestions would be appreciated.

John

--------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Is it possible to pass a Tcl/Tk component as argument to afunction

[John Fox]

Dear Hao,

> -Original Message-
> From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] 
> Sent: Thursday, May 17, 2007 8:29 AM
> To: John Fox
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] Is it possible to pass a Tcl/Tk component as 
> argument to afunction
> 
> John:
> 
> Thanks for your reply, I spent some time on this and the 
> conclusion is it works:
> 
> top<- tktoplevel()
> mainFrame <- tkcanvas(top)
> 
> both top and mainFrame can be used as parameters to pass to 
> other function. The name, however, will conflict each other 
> if they are defined in the same environment, which means if 
> you have another top and mainFrame passed to another 
> function, the GUIs will get to the same container. To avaoid 
> this, better  use different name.
> 
> I wonder if there is a way for me to create an env and 
> eliminate an env dynamically, I will try to get some 
> information, but I definitely welcome some quick inputs...
> 

Yes, you can create and delete environments, though it's not clear to me why
you should have to do so:

> env <- new.env()
> env

> parent.env(env)

> rm(env)

Regards,
 John

> Thanks
> Hao
> 
> 
> John Fox wrote: 
> 
>   Dear Hao,
>   
>   You might take a look at how the Rcmdr package is 
> implemented with many
>   reusable elements. There is, for example, an 
> initializeDialog function.
>   
>   I hope this helps,
>John
>   
>   
>   John Fox, Professor
>   Department of Sociology
>   McMaster University
>   Hamilton, Ontario
>   Canada L8S 4M4
>   905-525-9140x23604
>   http://socserv.mcmaster.ca/jfox 
>    
>   
> 
> 
>   -Original Message-
>   From: [EMAIL PROTECTED] 
>   [mailto:[EMAIL PROTECTED] On 
> Behalf Of Hao Liu
>   Sent: Wednesday, May 16, 2007 8:58 AM
>   To: r-help@stat.math.ethz.ch
>   Subject: [R] Is it possible to pass a Tcl/Tk 
> component as 
>   argument to afunction
>   
>   hi! All:
>   
>   I wonder if someone has done this before...
>   
>   I am writing several functions that conduct statistical 
>   analysis using a GUI interface by Tcl/Tk, they 
> share many 
>   identical GUI components. What I am trying to 
> do now is to 
>   simplify the code by writing a GUI repository 
> for all the 
>   components they use, thus save effort for code 
> maintenance.
>   
>   Since they all use:
>   
>   mainFrame <- tkcanvas(top)
>   
>   --
>   
>   I wonder if I can write functions that take 
> mainFrame as an 
>   argument, and call those functions from other 
> place with 
>   initialized tkcanvas object. I did not see 
> example like this 
>   and from my *limited* experience with tcltk, I found it 
>   always need something to be initialized before 
> put to use, 
>   that makes me wonder if this idea will work... 
> if it does 
>   not, any work arounds? like using Macro?
>   
>   Thanks
>   Hao
>   
>   [[alternative HTML version deleted]]
>   
>   __
>   R-help@stat.math.ethz.ch mailing list
>   https://stat.ethz.ch/mailman/listinfo/r-help
>   PLEASE do read the posting guide 
>   http://www.R-project.org/posting-guide.html
>   and provide commented, minimal, self-contained, 
> reproducible code.
>   
>   
> 
>   
>   
> 
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Is it possible to pass a Tcl/Tk component as argument to afunction

[John Fox]

Dear Hao,

You might take a look at how the Rcmdr package is implemented with many
reusable elements. There is, for example, an initializeDialog function.

I hope this helps,
 John


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Hao Liu
> Sent: Wednesday, May 16, 2007 8:58 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Is it possible to pass a Tcl/Tk component as 
> argument to afunction
> 
> hi! All:
> 
> I wonder if someone has done this before...
> 
> I am writing several functions that conduct statistical 
> analysis using a GUI interface by Tcl/Tk, they share many 
> identical GUI components. What I am trying to do now is to 
> simplify the code by writing a GUI repository for all the 
> components they use, thus save effort for code maintenance.
> 
> Since they all use:
> 
> mainFrame <- tkcanvas(top)
> 
> --
> 
> I wonder if I can write functions that take mainFrame as an 
> argument, and call those functions from other place with 
> initialized tkcanvas object. I did not see example like this 
> and from my *limited* experience with tcltk, I found it 
> always need something to be initialized before put to use, 
> that makes me wonder if this idea will work... if it does 
> not, any work arounds? like using Macro?
> 
> Thanks
> Hao
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Some questions on repeated measures (M)ANOVA & mixed modelswith lme4

[John Fox]

Dear Marco,

You might also take a look at ?Anova (or ?Manova) in the car package; the
last examples are for a repeated-measures ANOVA using both MANOVA and
univariate approaches, the latter with GG and HF corrections.

I hope this helps,
 John


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Marco B
> Sent: Sunday, May 13, 2007 10:22 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Some questions on repeated measures (M)ANOVA & 
> mixed modelswith lme4
> 
> Dear R Masters,
> 
> I'm an anesthesiology resident trying to make his way through 
> basic statistics. Recently I have been confronted with 
> longitudinal data in a treatment vs. control analysis. My 
> dataframe is in the form of:
> 
> subj | group | baseline | time | outcome (long) or subj | 
> group | baseline | time1 |...| time6 | (wide)
> 
> The measured variable is a continuous one. The null 
> hypothesis in this analysis is that the Group factor does not 
> significantly influence the outcome variable. A secondary 
> null hypothesis is that the Group x Time interaction is not 
> significant, either. Visual of the group means indicates the 
> outcome measure decreases linearly (more or less) over time 
> from baseline values. The time==1...time==6 intervals are 
> equally-spaced and we have equal sample sizes for the groups.
> 
> I've done a little reading around and found (at least) four 
> possible approaches:
> 
> A. Linear mixed model using lme4 with random intercept and slope with
> lmer() or lme()
> 
> B. Repeated measures ANOVA using aov() with Error() 
> stratification (found in Baron & Li, 2006), something along 
> the lines of:
> aov(outcome ~ group * time + baseline + Error(subj+subj:time))
> 
> (from: http://cran.r-project.org/doc/contrib/Baron-rpsych.pdf, p. 41)
> 
> C. "Repeated measures" MANOVA as follows (using data in wide format):
> response <- cbind(time1,time2,time3,time4,time5,time6)
> mlmfit <- lm(response ~ group)
> mlmfit1 <- lm(response ~ 1)
> mlmfit0 <- lm(response ~ 0)
> # Test time*group effect
> anova.mlm(mlmfit, mlmfit1, X=~1, test="Spherical") # Test 
> overall group effect anova.mlm(mlmfit, mlmfit1, M=~1) # Test 
> overall time effect anova.mlm(mlmfit1, mlmfit0, X=~1, 
> test="Spherical")
> 
> (taken from http://tolstoy.newcastle.edu.au/R/help/05/11/15744.html)
> 
> Now, on with the questions:
> 
> 1. This is really a curiosity. I find lmer() easier to use 
> than lme(), but the former does not allow the user to model 
> the correlation structure of the data. I figure lmer() is 
> presently assuming no within-group correlation for the data, 
> which I guess is unlikely in my example. Is there a way to 
> compare directly (maybe in terms of
> log-likelihood?) similar models fitted in lme() and lmer()?
> 
> 2. Baron & Li suggest a painful (at least for me) procedure 
> to obtain Greenhouse-Geisser or Huyn-Feldt correction for the 
> ANOVA analysis they propose. Is there a package or function 
> which simplifies the procedure?
> 
> 3. I must admit that I don't understand solution C. I can 
> "hack" it to fit my model, and it seems to work, but I can't 
> seem to grasp the overall concept, especially regarding the 
> outer and/or inner projection matrices (M & X). Could anyone 
> point me to a basic explanation of the procedure?
> 
> 4. Provided the assumptions for ANOVA hold, or that 
> deviations from them are not horrible, am I correct in saying 
> that this procedure would be the most powerful one? How would 
> you choose solution A over solution B (or viceversa)?
> 
> My sincerest gratitude to anyone who will take the time to 
> answer my questions!
> 
> Best Regards,
> 
> Marco
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Weighted least squares

[John Fox]

Dear Hadley,

> -Original Message-
> From: hadley wickham [mailto:[EMAIL PROTECTED] 
> Sent: Wednesday, May 09, 2007 2:21 AM
> To: John Fox
> Cc: R-help@stat.math.ethz.ch
> Subject: Re: [R] Weighted least squares
> 
> Thanks John,
> 
> That's just the explanation I was looking for. I had hoped 
> that there would be a built in way of dealing with them with 
> R, but obviously not.
> 
> Given that explanation, it stills seems to me that the way R 
> calculates n is suboptimal, as demonstrated by my second example:
> 
> summary(lm(y ~ x, data=df, weights=rep(c(0,2), each=50))) 
> summary(lm(y ~ x, data=df, weights=rep(c(0.01,2), each=50)))
> 
> the weights are only very slightly different but the 
> estimates of residual standard error are quite different (20 
> vs 14 in my run)
> 

Observations with 0 weight are literally excluded, while those with very
small weight (relative to others) don't contribute much to the fit.
Consequently you get very similar coefficients but different numbers of
observations.

I hope this helps,
 John

> Hadley
> 
> On 5/8/07, John Fox <[EMAIL PROTECTED]> wrote:
> > Dear Hadley,
> >
> > I think that the problem is that the term "weights" has different 
> > meanings, which, although they are related, are not quite the same.
> >
> > The weights used by lm() are (inverse-)"variance weights," 
> reflecting 
> > the variances of the errors, with observations that have 
> low-variance 
> > errors therefore being accorded greater weight in the 
> resulting WLS regression.
> > What you have are sometimes called "case weights," and I'm 
> unaware of 
> > a general way of handling them in R, although you could 
> regenerate the 
> > unaggregated data. As you discovered, you get the same coefficients 
> > with case weights as with variance weights, but different 
> standard errors.
> > Finally, there are "sampling weights," which are inversely 
> > proportional to the probability of selection; these are 
> accommodated by the survey package.
> >
> > To complicate matters, this terminology isn't entirely standard.
> >
> > I hope this helps,
> >  John
> >
> > 
> > John Fox, Professor
> > Department of Sociology
> > McMaster University
> > Hamilton, Ontario
> > Canada L8S 4M4
> > 905-525-9140x23604
> > http://socserv.mcmaster.ca/jfox
> > 
> >
> > > -Original Message-
> > > From: [EMAIL PROTECTED] 
> > > [mailto:[EMAIL PROTECTED] On Behalf Of hadley 
> > > wickham
> > > Sent: Tuesday, May 08, 2007 5:09 AM
> > > To: R Help
> > > Subject: [R] Weighted least squares
> > >
> > > Dear all,
> > >
> > > I'm struggling with weighted least squares, where 
> something that I 
> > > had assumed to be true appears not to be the case.
> > > Take the following data set as an example:
> > >
> > > df <- data.frame(x = runif(100, 0, 100)) df$y <- df$x + 1 + 
> > > rnorm(100, sd=15)
> > >
> > > I had expected that:
> > >
> > > summary(lm(y ~ x, data=df, weights=rep(2, 100))) 
> summary(lm(y ~ x, 
> > > data=rbind(df,df)))
> > >
> > > would be equivalent, but they are not.  I suspect the 
> difference is 
> > > how the degrees of freedom is calculated - I had expected 
> it to be 
> > > sum(weights), but seems to be sum(weights > 0).  This seems 
> > > unintuitive to me:
> > >
> > > summary(lm(y ~ x, data=df, weights=rep(c(0,2), each=50))) 
> > > summary(lm(y ~ x, data=df, weights=rep(c(0.01,2), each=50)))
> > >
> > > What am I missing?  And what is the usual way to do a linear 
> > > regression when you have aggregated data?
> > >
> > > Thanks,
> > >
> > > Hadley
> > >
> > > __
> > > R-help@stat.math.ethz.ch mailing list 
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide
> > > http://www.R-project.org/posting-guide.html
> > > and provide commented, minimal, self-contained, reproducible code.
> > >
> >
> >
> >
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Weighted least squares

[John Fox]

Dear Adai,

> -Original Message-
> From: Adaikalavan Ramasamy [mailto:[EMAIL PROTECTED] 
> Sent: Tuesday, May 08, 2007 8:38 PM
> To: S Ellison
> Cc: [EMAIL PROTECTED]; [EMAIL PROTECTED]; R-help@stat.math.ethz.ch
> Subject: Re: [R] Weighted least squares
> 
> http://en.wikipedia.org/wiki/Weighted_least_squares gives a 
> formulaic description of what you have said.
> 
> I believe the original poster has converted something like this
> 
>   y x
>   0   1.1
>   0   2.2
>   0   2.2
>   0   2.2
>   1   3.3
>   1   3.3
>   2   4.4
>  ...
> 
> into something like the following
> 
>   y x freq
>   0   1.11
>   0   2.23
>   1   3.32
>   2   4.41
>  ...
> 
> Now, the variance of means of each row in table above is ZERO 
> because the individual elements that comprise each row are 
> identical. Therefore your method of using inverse-variance 
> will not work here.
> 
> Then is it valid then to use lm( y ~ x, weights=freq ) ?

No, because the weights argument gives inverse-variance weights not case
weights.

Regards,
 John

> 
> Regards, Adai
> 
> 
> 
> S Ellison wrote:
> > Hadley,
> > 
> > You asked
> >> .. what is the usual way to do a linear regression when you have 
> >> aggregated data?
> > 
> > Least squares generally uses inverse variance weighting. 
> For aggregated data fitted as mean values, you just need the 
> variances for the _means_. 
> > 
> > So if you have individual means x_i and sd's s_i that arise 
> from aggregated data with n_i observations in group i, the 
> natural weighting is by inverse squared standard error of the 
> mean. The appropriate weight for x_i would then be 
> n_i/(s_i^2). In R, that's n/(s^2), as n and s would be 
> vectors with the same length as x. If all the groups had the 
> same variance, or nearly so, s is a scalar; if they have the 
> same number of observations, n is a scalar. 
> > 
> > Of course, if they have the same variance and same number 
> of observations, they all have the same weight and you 
> needn't weight them at all: see previous posting!
> > 
> > Steve E
> > 
> > 
> > 
> > ***
> > This email and any attachments are confidential. Any use, 
> > co...{{dropped}}
> > 
> > __
> > R-help@stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide 
> > http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> > 
> > 
> > 
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Weighted least squares

[John Fox]

Dear Hadley,

I think that the problem is that the term "weights" has different meanings,
which, although they are related, are not quite the same. 

The weights used by lm() are (inverse-)"variance weights," reflecting the
variances of the errors, with observations that have low-variance errors
therefore being accorded greater weight in the resulting WLS regression.
What you have are sometimes called "case weights," and I'm unaware of a
general way of handling them in R, although you could regenerate the
unaggregated data. As you discovered, you get the same coefficients with
case weights as with variance weights, but different standard errors.
Finally, there are "sampling weights," which are inversely proportional to
the probability of selection; these are accommodated by the survey package. 

To complicate matters, this terminology isn't entirely standard.

I hope this helps,
 John


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of hadley wickham
> Sent: Tuesday, May 08, 2007 5:09 AM
> To: R Help
> Subject: [R] Weighted least squares
> 
> Dear all,
> 
> I'm struggling with weighted least squares, where something 
> that I had assumed to be true appears not to be the case.  
> Take the following data set as an example:
> 
> df <- data.frame(x = runif(100, 0, 100)) df$y <- df$x + 1 + 
> rnorm(100, sd=15)
> 
> I had expected that:
> 
> summary(lm(y ~ x, data=df, weights=rep(2, 100))) summary(lm(y 
> ~ x, data=rbind(df,df)))
> 
> would be equivalent, but they are not.  I suspect the 
> difference is how the degrees of freedom is calculated - I 
> had expected it to be sum(weights), but seems to be 
> sum(weights > 0).  This seems unintuitive to me:
> 
> summary(lm(y ~ x, data=df, weights=rep(c(0,2), each=50))) 
> summary(lm(y ~ x, data=df, weights=rep(c(0.01,2), each=50)))
> 
> What am I missing?  And what is the usual way to do a linear 
> regression when you have aggregated data?
> 
> Thanks,
> 
> Hadley
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Library & Package for Tobit regression

[John Fox]

Dear Sattar,

You can use the survreg function in the survival package, which is part of
the standard R distribution. [BTW, help.search("Tobit") would have led you
to that.] There are other possibilities as well -- e.g., the tobit function
in the VGAM package.

I hope this helps,
 John

--------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Abdus Sattar
> Sent: Thursday, May 03, 2007 10:09 PM
> To: R-help@stat.math.ethz.ch
> Subject: [R] Library & Package for Tobit regression
> 
> Hello R-Users:
>  
> I am want to use tobit regression for left censored 
> panel/longitudinal data. Could you please provide me the name 
> of "library" and/or "package" that will give me option of 
> fitting tobit regression model for longitudinal data? 
>  
> Thank you. 
>  
> Sattar
> 
> __
> 
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] R News, volume 7, issue 1 is now available --error in AMMI article

[John Fox]

Dear Yuandan,

The function definition given in the article doesn't produce an error for
me.

Regards,
 John


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: Lux Zhang [mailto:[EMAIL PROTECTED] 
> Sent: Thursday, April 26, 2007 7:56 AM
> To: John Fox
> Cc: R-help
> Subject: Re: [R] R News, volume 7,issue 1 is now available 
> --error in AMMI article
> 
> 
> 
> On 26/04/07, Lux Zhang <[EMAIL PROTECTED]> wrote:
> 
> 
> 
>   On 26/04/07, John Fox <[EMAIL PROTECTED] > wrote:
> 
>   Dear Yuandan,
>   
>   My attention was drawn by your claim of an 
> "error in AMMI article."
>   
>   As you say, the code for the AMMI function is 
> given directly in the article.
>   If the argument biplot is equal to 1, then a 
> biplot is drawn by the 
>   function, as is apparent from the section of 
> code labelled "## 5 - Biplots."
>   
>   Why is this an error?
> 
> 
> 
>   when loading this AMMI function, at the line "if 
> (biplot == 1) { ", R seems treating the 'biplot' as a 
> subject, the biplot function from stats package, inseatd of 
> treating it as argument for the AMMI function. 
>   
>   
>   here is the error messenge when I load it
>   
>   > source ("AMMI.R")
>   Error in parse(file, n = -1, NULL, "?") : syntax error at
>   51: ( bplot == 1 ) {
>   52: plot(1, type = 
> 
> 
> Sorry, I had another look, it seems something to do with plot 
> (1, type = 'n' ... [ i copy this code from the pdf file]
> 
> after change it to 
> 
> plot (1, type = "n" ... as below 
> 
>   if ( biplot == 1 ) {
> plot(1, type = "n", xlim = range(c(envir.mean, 
> var.mean)), ylim = range(c(E[,1], G[,1])), xlab = "Yield",
> 
> lt was loaded. 
> 
> 
> 
>  
> 
> 
>   
> 
>   Regards,
>   John
>   
>    
>   John Fox, Professor
>   Department of Sociology
>   McMaster University
>   Hamilton, Ontario
>   Canada L8S 4M4
>   905-525-9140x23604
>   http://socserv.mcmaster.ca/jfox 
> <http://socserv.mcmaster.ca/jfox> 
>    
>   
>   > -Original Message-
>   > From: [EMAIL PROTECTED]
>   > [mailto: [EMAIL PROTECTED] ] 
> On Behalf Of Lux Zhang
>   > Sent: Thursday, April 26, 2007 1:38 AM
>   > To: R-help; [EMAIL PROTECTED] 
>   > Subject: Re: [R] R News, volume 7,issue 1 is 
> now available 
>   > --error in AMMI article
>   >
>   > Hi,
>   >
>   > In this newsletter (Vol 7, 1),the article on 
> AMMI by Onofri
>   > and Ciriofolo presented a AMMI function.  One 
> of arguments
>   > for this function AMMI (Page 
>   > 17) is biplot. There is a biplot fucntion 
> from  {stats}
>   > package.  I guess they are not the same. 
> Could the authors
>   > clarify that?
>   >
>   > Thanks,
>   >
>   > Yuandan
>   >
>   >   [[alternative HTML version deleted]] 
>   >
>   > __
>   > R-help@stat.math.ethz.ch mailing list 
>   > https://stat.ethz.ch/mailman/listinfo/r-help 
>   > PLEASE do read the posting guide
>   > http://www.R-project.org/posting-guide.html 
>   > and provide commented, minimal, 
> self-contained, reproducible code. 
>   >
>   
>   
>   
> 
> 
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] R News, volume 7, issue 1 is now available --error in AMMI article

[John Fox]

Dear Yuandan,

My attention was drawn by your claim of an "error in AMMI article."

As you say, the code for the AMMI function is given directly in the article.
If the argument biplot is equal to 1, then a biplot is drawn by the
function, as is apparent from the section of code labelled "## 5 - Biplots."

Why is this an error?

Regards,
 John

--------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Lux Zhang
> Sent: Thursday, April 26, 2007 1:38 AM
> To: R-help; [EMAIL PROTECTED]
> Subject: Re: [R] R News, volume 7,issue 1 is now available 
> --error in AMMI article
> 
> Hi,
> 
> In this newsletter (Vol 7, 1),the article on AMMI by Onofri 
> and Ciriofolo presented a AMMI function.  One of arguments 
> for this function AMMI (Page
> 17) is biplot. There is a biplot fucntion from  {stats} 
> package.  I guess they are not the same. Could the authors 
> clarify that?
> 
> Thanks,
> 
> Yuandan
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] intersect more than two sets

[John Fox]

Dear Peter, Hadley, et al.,

Just to bring it back around, here's a recursive Fold():

> a <- letters[1:4]
> b <- letters[2:5]
> c <- letters[3:6]
> d <- letters[4:7]
> e <- letters[5:8]

> Fold <- function(f, x, y, ...){
+ if (missing(...)) f(x, y)
+ else f(x, Fold(f, y, ...))
+ }

> Fold(intersect, a, b)
[1] "b" "c" "d"
> Fold(intersect, a, b, c)
[1] "c" "d"
> Fold(intersect, a, b, c, d)
[1] "d"
> Fold(intersect, a, b, c, d, e)
character(0)
> 
> do.call(Fold, list(intersect, a, b, c, d))
[1] "d"

Regards,
 John


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Peter Dalgaard
> Sent: Tuesday, April 24, 2007 4:03 PM
> To: hadley wickham
> Cc: R Help; Weiwei Shi
> Subject: Re: [R] intersect more than two sets
> 
> hadley wickham wrote:
> > On 4/24/07, Weiwei Shi <[EMAIL PROTECTED]> wrote:
> >   
> >> assume t2 is a list of size 11 and each element is a 
> vector of characters.
> >>
> >> the following codes can get what I wanted but I assume 
> there might be 
> >> a one-line code for that:
> >>
> >> t3 <- t2[[1]]
> >> for ( i in 2:11){
> >> t3 <- intersect(t2[[i]], t3)
> >> }
> >>
> >> or there is no such "apply"?
> >> 
> >
> > The operation you want is called a fold 
> > 
> (http://en.wikipedia.org/wiki/Fold_%28higher-order_function%29), and 
> > if it was available in R, you'd be able to do:
> >
> > fold(t2, intersect)
> >
> > Unfortunately, it's not, but you could implement it as follows:
> >
> > fold <- function(x, fun) {
> > if (length(x) == 1) return(fun(x))
> > 
> > accumulator <- fun(x[[1]], x[[2]])
> > if (length(x) == 2) return(accumulator)
> >
> > for(i in 3:length(x)) {
> > accumulator <- fun(accumulator, x[[i]])
> > }
> > accumulator
> > }
> >
> > a <- list(c(1,3,5), c(1,3), c(1, 2, 5, 6)) fold(a, intersect)
> >
> >   
> 
> It's come up before. Gabor G posted this rather more succinct version:
> 
>  > Fold <- function(f, x, L) (for(e in L) x <- f(x, e))  > 
> Fold(intersect,a[[1]],a[-1]) [1] 1
> 
> or maybe prettier:
> 
>  > (E <- Fold(union, NULL, a))
> [1] 1 3 5 2 6
>  > Fold(intersect, E, a)
> [1] 1
> 
> 
> > Which is just a trivial generalisation of your code above
> >
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] intersect more than two sets

[John Fox]

Dear Weiwei Shi,

How about using recursion?

> intersection <- function(x, y, ...){
+ if (missing(...)) intersect(x, y)
+ else intersect(x, intersection(y, ...))
+ }
 
> a <- letters[1:4]
> b <- letters[2:5]
> c <- letters[3:6]
> d <- letters[4:7]
> e <- letters[5:8]
 
> intersection(a, b)
[1] "b" "c" "d"
> intersection(a, b, c)
[1] "c" "d"
> intersection(a, b, c, d)
[1] "d"
> intersection(a, b, c, d, e)
character(0)
 
> do.call(intersection, list(a, b, c))
[1] "c" "d"

I hope this helps,
 John


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Weiwei Shi
> Sent: Tuesday, April 24, 2007 2:59 PM
> To: R Help
> Subject: Re: [R] intersect more than two sets
> 
> assume t2 is a list of size 11 and each element is a vector 
> of characters.
> 
> the following codes can get what I wanted but I assume there 
> might be a one-line code for that:
> 
> t3 <- t2[[1]]
> for ( i in 2:11){
>   t3 <- intersect(t2[[i]], t3)
> }
> 
> or there is no such "apply"?
> 
> On 4/24/07, Weiwei Shi <[EMAIL PROTECTED]> wrote:
> > Hi,
> > I searched the archives and did not find a good solution to that.
> >
> > assume I have 10 sets and I want to have the common 
> character elements of them.
> >
> > how could i do that?
> >
> > --
> > Weiwei Shi, Ph.D
> > Research Scientist
> > GeneGO, Inc.
> >
> > "Did you always know?"
> > "No, I did not. But I believed..."
> > ---Matrix III
> >
> 
> 
> --
> Weiwei Shi, Ph.D
> Research Scientist
> GeneGO, Inc.
> 
> "Did you always know?"
> "No, I did not. But I believed..."
> ---Matrix III
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Fitting multinomial response in structural equation

[John Fox]

Dear adschai,

I'm not sure that I entirely follow what you want to do, but if the response
really is qualitative I don't see the sense in transforming it into a single
quantitative variable. Nor does the strategy of generating 24 dichotomous,
separately modelled responses make sense, since these are correlated. In
some circumstances, however, one can resolve a polytomous response into a
set of *nested* dichotomies, which are then independent of one another.
Finally, I wouldn't as a general matter recommend fitting any statistical
model to a 24-category response.

I suspect that you'd do well to find someone with whom you can about your
research problem.

Regards,
 John

--------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of 
> [EMAIL PROTECTED]
> Sent: Saturday, April 21, 2007 4:32 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Fitting multinomial response in structural equation
> 
> Hi - I am confronting a situation where I have a set of 
> structural equation and one or two of my responses are 
> multinomial. I understand that sem would not deal with the 
> unordered response. So I am thinking of the following two ways:
> 
> 1. Expanding my response to a new set of binary variables 
> corresponding to each label of my multinomial response. Then 
> use each of these as a separate response in my model. 
> However, since I have about 24 labels in this single 
> variable, it will be very expensive to do this way.
> 2. I am thinking of transforming this variable into a 
> continous-valued variable. I am thinking of using the 
> observed count to transform this variable using the probit 
> function. Then my new variable is just a step-wise function. 
> The trouble that I am struggling with is that this response 
> variable will also serve as a predictor in another equation 
> in my structural model. The interpretation of this equation 
> is not so straightforward for me. The coefficient of this 
> variable is no longer reading 'a unit change in this variable 
> holding everything else fixed corresponds to the x unit 
> change of the response'. All I can read from this method is 
> that when I change from one label to another, it means p 
> amount change in my step-wise-function predictor variable and 
> it corresponds to x unit change of the response holding 
> everything fixed.
> 
> The main purpose here for myself to post my question here is 
> to obtain your insight especially with respect to using sem 
> with the two approaches above. I would like to ensure that my 
> approaches make sense within the context of sem. Any 
> comments/opinions would be really appreciated. Thank you so 
> much in advance.
> 
> - adschai
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Suggestions for statistical computing course

[John Fox]

Dear Duncan and Giovanni,

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Duncan Murdoch
> Sent: Friday, April 20, 2007 10:13 AM
> To: Giovanni Petris
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] Suggestions for statistical computing course
> 
> On 4/20/2007 9:34 AM, Giovanni Petris wrote:

. . .

> 
> There are a couple of shareware/freeware editors (WinEDT, 
> Tinn-R) that have hooks to R.  WinEDT also has support for 
> TeX/LaTeX; if that's important to you, it might be worth the 
> cost/effort to install.  I'm less familiar with Tinn-R, but I 
> believe it's free, whereas WinEDT is not.
> 

Tinn-R is indeed free and also has some support for LaTeX. Information is
available at .

Regards,
 John

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to convert the lower triangle of a matrix to a symmetricmatrix

[John Fox]

Dear Ranjan,

If the elements are ordered by rows, then the following should do the trick:

X <- diag(p)
X[upper.tri(X, diag=TRUE)] <- elements
X <- X + t(X) - diag(diag(X))

If they are ordered by columns, substitute lower.tri() for upper.tri().

I hope this helps,
 John

----
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Ranjan Maitra
> Sent: Thursday, April 19, 2007 9:28 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] how to convert the lower triangle of a matrix to 
> a symmetricmatrix
> 
> Hi,
> 
> I have a vector of p*(p+1)/2 elements, essentially the lower 
> triangle of a symmetric matrix. I was wondering if there is 
> an easy way to make it fill a symmetric matrix. I have to do 
> it several times, hence some efficient approach would be very useful.
> 
> Many thanks and best wishes,
> Ranjan
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] general question about plotting multiple regression results

[John Fox]

Dear Thomas and Simon,

On Thu, 19 Apr 2007 07:37:28 -0700 (PDT)
 Thomas Lumley <[EMAIL PROTECTED]> wrote:
> On Thu, 19 Apr 2007, Simon Pickett wrote:
> 
> > Hi all,
> >
> > I have been bumbling around with r for years now and still havent
> come up
> > with a solution for plotting reliable graphs of relationships from
> a
> > linear regression.
> 
> termplot() does this for a range of regression models (without
> interaction 
> terms). The "effects" package does it better for linear regression
> models.
> 
>   -thomas
>

The effects package also works for generalized linear models (which, I
suppose, are arguably linear regression models).

Regards,
 John
 
> 
> > Here is an example illustrating my problem
> >
> > 1.I do a linear regression as follows
> >
> >
>
summary(lm(n.day13~n.day1+ffemale.yell+fmale.yell+fmale.chroma,data=surv))
> >
> > which gives some nice sig. results
> >
> > Coefficients:
> > Estimate Std. Error t value Pr(>|t|)
> > (Intercept)  -0.739170.43742  -1.690 0.093069 .
> > n.day11.004600.05369  18.711  < 2e-16 ***
> > ffemale.yell  0.224190.06251   3.586 0.000449 ***
> > fmale.yell0.258740.06925   3.736 0.000262 ***
> > fmale.chroma  0.235250.11633   2.022 0.044868 *
> >
> > 2. I want to plot the effect of "ffemale.yell", "fmale.yell" and
> > "fmale.chroma" on my response variable.
> >
> > So, I either plot the raw values (which is fine when there is a
> very
> > strong relationship) but what if I want to plot the effects from
> the
> > model?
> >
> > In this case I would usually plot the fitted values values against
> the raw
> > values of x... Is this the right approach?
> >
> >
>
fit<-fitted(lm(n.day13~n.day1+ffemale.yell+fmale.yell+fmale.chroma,data=fsurv1))
> >
> > plot(fit~ffemale.yell)
> >
> > #make a dummy variable across the range of x
> > x<-seq(from=min(fsurv1$ffemale.yell),to=max(fsurv1$ffemale.yell),
> length=100)
> >
> > #get the coefficients and draw the line
> > co<-coef(lm(fit~ffemale.yell,data=fsurv1))
> > y<-(co[2]*x)+co[1]
> > lines(x,y, lwd=2)
> >
> > This often does the trick but for some reason, especially when my
> model
> > has many terms in it or when one of the independent variables is
> only
> > significant when the other independent variables are in the
> equation, it
> > gives me strange lines.
> >
> > Please can someone show me the light?
> >
> > Thanks in advance,
> >
> > Simon.
> >
> >
> >
> >
> >
> >
> > Simon Pickett
> > PhD student
> > Centre For Ecology and Conservation
> > Tremough Campus
> > University of Exeter in Cornwall
> > TR109EZ
> > Tel 01326371852
> >
> > __
> > R-help@stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
> 
> Thomas Lumley Assoc. Professor, Biostatistics
> [EMAIL PROTECTED] University of Washington, Seattle
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Fit sem: problem of 'Error in solve.default(C[ind, ind]) : Lapack routine dgesv: system is exactly singular.'

[John Fox]

Dear adschai,

This model is almost surely underidentified, but it's hard to figure that
out for sure because of its very odd structure: you have the three latent
variables all influencing each other mutually and all with exactly the same
observed indicators, but the structural disturbances are all specified to be
uncorrelated. Additionally, there are observed exogenous variables that
affect the latent endogenous variables, but these exogenous variables are
also specified to be uncorrelated. It's hard for me to imagine that you
really intended this, and even if the model is identified, I seriously doubt
that you can fit it to the data. Finally, you should verify that the input
covariance matrix is positive-definite.

I think that the issues of model specification go well beyond how to use the
software, and I strongly suggest that you try to find someone local to talk
to about your research.

John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of 
> [EMAIL PROTECTED]
> Sent: Sunday, April 15, 2007 2:16 PM
> To: [EMAIL PROTECTED]
> Subject: [R] Fit sem: problem of 'Error in 
> solve.default(C[ind, ind]) : Lapack routine dgesv: system is 
> exactly singular.'
> 
> Hi - I would need help another time. I would the model and 
> wire the correlation matrix into the sem() method. I got the 
> system is exactly singular? Is my model underidentified? Here 
> I have my model below along with the result from traceback() 
> funciton. Any help would be really appreciated. Thank you again!
> 
> Average_Daily_Traffic -> Lat_Env, 
> gam1_1,   NA
> Average_Daily_Truck_Traffic -> Lat_Env,   
> gam1_2,   NA
> BridgeAges -> Lat_Env,
> gam1_3,   NA
> 
> Design_Load -> Lat_Specs, 
> gam2_1,   NA
> Type_of_Service_On_Bridge -> Lat_Specs,   
> gam2_2,   NA
> Railroad_Beneath -> Lat_Specs,
> gam2_3,   NA
> BridgeAges -> Lat_Specs,  
> gam2_4,   NA
> Highway_Beneath -> Lat_Specs, 
> gam2_5,   NA
> 
> Main_Structure_Material -> Lat_Design,
> gam3_1,   NA
> Main_Structure_Design -> Lat_Design,  
> gam3_2,   NA
> Length_of_Maximum_Span -> Lat_Design, 
> gam3_4,   NA
> Number_of_Spans_In_Main_Unit -> Lat_Design,   
> gam3_5,   NA
> BridgeAges -> Lat_Design, 
> gam3_6,   NA
> 
> Lat_Env -> Lat_Specs, 
> beta2_1,  NA
> Lat_Env -> Lat_Design,
> beta3_1,  NA
> Lat_Specs -> Lat_Env, 
> beta1_2,  NA
> Lat_Specs -> Lat_Design,  
> beta3_2,  NA
> Lat_Design -> Lat_Env,
> beta1_3,  NA
> Lat_Design -> Lat_Specs,  
> beta1_2,  NA
> 
> Lat_Env -> Operating_Rating,  NA, 
>   1
> Lat_Env -> Deck_Cond_Rating,  
> lamy2_1,  NA
> Lat_Env -> Superstructure_Cond_Rating,
> lamy3_1,  NA
> Lat_Env -> Substructure_Cond_Rating,  
> lamy4_1,  NA
> Lat_Specs -> Operating_Rating,NA, 
>   1
> Lat_Specs -> Deck_Cond_Rating,
> lamy2_2,  NA
> Lat_Specs -> Superstructure_Cond_Rating,  
> lamy3_2,  NA
> Lat_Specs -> Substructure_Cond_Rating,
> lamy4_2,  NA
> Lat_Design -> Operating_Rating,   NA, 
>   1
> Lat_Design -> Deck_Cond_Rating,   
> lamy2_3,  NA
> Lat_Design -> Superstructure_Cond_Rating, 
> lamy3_3,  NA
> Lat_Design -> Substructure_Cond_Rating,   
> lamy4_3,  NA
> 
> Lat_Env <-> Lat_Specs,
> psi2_1,   NA
> Lat_Env <-> Lat_Design,   
> psi3_1,   NA
> Lat_Specs <-> Lat_Design,

Re: [R] Fit sem model with intercept

[John Fox]

Dear adschai,

You needn't look too far, since the last example in ?sem is for a model with
an intercept. One would use the raw-moment matrix as input to sem(), either
entered directly or calculated with the raw.moments() function in the sem
package. The row/column of the raw-moment matrix is given a name just like
the other columns. You could use the name "1"; in the example, the name is
"UNIT".

As you say, however, you're using polychoric and polyserial correlations as
input. Since the origin and scale of the latent continuous variables
underlying the ordinal variables are entirely arbitrary, I can't imagine
what the purpose of a model with an intercept would be, but it's possible
that I'm missing something. If you think that this makes some sense, then
you could convert the correlations to raw moments by using the means and
standard deviations of the observed variables along with the means and
standard deviations that you assign to the latent variables derived from the
ordinal variables (the latter on what basis I can't imagine, but I suppose
you could fix them to 0s and 1s).

Finally, if the sem model that you show is meant to be a complete
specification, I notice that it includes no covariance components; moreover,
if this is the complete structural part of the model, then I think it is
underidentified, and the two parts of the model (those involving eta1 and
eta2) appear entirely separate.

I hope this helps,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of 
> [EMAIL PROTECTED]
> Sent: Sunday, April 15, 2007 4:28 AM
> To: [EMAIL PROTECTED]
> Subject: [R] Fit sem model with intercept
> 
> Hi - I am trying to fit sem model with intercepts. Here is 
> what I have in my model.
> 
> Exogeneous vars: x1 (continous), x2 (ordinal), x3 (ordinal), 
> x4(continuous) Endogeneous vars: y1 (continuous), y2 
> (ordinal), y3 (ordinal)
> 
> SEM model:
> x1 -> eta1; x2 -> eta1; x3 -> eta2; x4 -> eta2; eta1 -> 
> y1, eta1 -> y2, eta2 -> y2, eta2 -> y3 
> 
> However, in these arrow models, I don't know how to add 
> intercept onto it. I am trying to find an example code using 
> sem package on how to incorporate intercept but cannot find 
> any documents on the web. Or we can simply add something like 
> this? '1 ->  eta1'? This is my first question. 
> 
> Also, note that since my y2 and y3 are ordinal, I used the 
> 'hetcor' to calculate correlation of observed variables. 
> However, from the document, I would need to use the 
> covariance matrix rather then the correlation. And I need 
> additional column for 1. I am not sure how this matrix should 
> look like and how I can obtain this? If there is any example 
> you could point me to, I would really appreciate. Thank you.
> 
> - adschai
> 
>   [[alternative HTML version deleted]]
> 
> __
> [EMAIL PROTECTED] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Windows Vista issues

[John Fox]

Dear Brian,

> -Original Message-
> From: Prof Brian Ripley [mailto:[EMAIL PROTECTED] 
> Sent: Saturday, April 14, 2007 9:44 AM
> To: John Fox
> Cc: [EMAIL PROTECTED]
> Subject: RE: [R] Windows Vista issues
> 
> On Sat, 14 Apr 2007, John Fox wrote:
> 
> > Dear Brian,
> >
> > I don't yet have a Vista machine, but it occurs to me that some of 
> > these problems might be avoided if R installed by default into c:\R 
> > rather than c:\Program Files\R. Is that the case?
> 
> If C:\R is not a system area, the first para of 2) applies: 
> in fact I tried installing into E:\R as an unprivileged user. 
>  That's what I would recommend a sole user of R to do (and 
> have long done myself). 

Yes, I do that too (to c:\R, in my case). I don't know whether Vista
considers c:\R a system area; if not, making c:\R the default would avoid
problems for sole users.

> However, if you are a sysadmin 
> installing for many users of the machine, you will want to 
> install into a system area (or else a standard user is likely 
> to be able to uninstall R, for example).
> 

Assuming that c:\R is not a system area, the question is what's the best
default location for installing R. I suspect that most people who use R on
Windows are using it on sole-user machines, and are not used to
differentiating between administrator and other accounts. Vista will make
them attend to this distinction, and they may well get in the habit of
installing software as an administrator. This would lead to the kinds of
problems that have been surfacing: Difficulty installing packages as an
unprivileged user after doing a default install of R as an administrator.

Sysadmins, on the other hand, tend to be more sophisticated, and could
install R to c:\Program Files\R even if the default location is c:\R.

> My memory is that long ago the default for the R installer 
> was c:\R, and user pressure made us change it.  

Yes, I recall that. I agree that the change made sense at the time.

> BTW, the 
> default is only c:\Program Files\R for an administrator 
> account: for unprivileged users it is under c:\Documents and 
> Settings\ for XP and c:\Users\ for Vista.
> We have some control in the installer script of these things, 
> but have chosen to take the defaults which give the most 
> familiar experience to end users.
> 

I think that the general principle should be to choose defaults to cause the
least difficulty, especially to unsophisticated users. Perhaps the
implications of that have changed with Vista. As I said, I don't currently
have a Vista machine, so none of this is based on experience.

Regards,
 John

> The file access reporting issues were occurring in ordinary 
> user space.
> 
> Brian
> 
> 
> > Thank you, by the way, for pursuing these issues.
> >
> > John
> >
> > 
> > John Fox
> > Department of Sociology
> > McMaster University
> > Hamilton, Ontario
> > Canada L8S 4M4
> > 905-525-9140x23604
> > http://socserv.mcmaster.ca/jfox
> > 
> >
> >> -Original Message-
> >> From: [EMAIL PROTECTED] 
> >> [mailto:[EMAIL PROTECTED] On Behalf Of Prof Brian 
> >> Ripley
> >> Sent: Saturday, April 14, 2007 1:57 AM
> >> To: [EMAIL PROTECTED]
> >> Subject: [R] Windows Vista issues
> >>
> >> It seemed FUD [*] has been prevailing here and elsewhere on Vista 
> >> security features.
> >>
> >> I asked our sysadmins to set up a Vista box for me on which I have 
> >> access to all levels of accounts.  Many of the issues I found were 
> >> covered by earlier answers and all in the upcoming rw-FAQ 
> (currently 
> >> available at http://www.stats.ox.ac.uk/pub/R/rw-FAQ.html 
> and in the 
> >> 2.5.0
> >> pre-releases) but a quick reprise may help.
> >>
> >> Most of my testing was of 2.5.0 beta, but I did some quick 
> tests of 
> >> 2.4.1.
> >>
> >>
> >> 1) The R installer and uninstaller are from an 'unidentified 
> >> publisher'
> >> and you may have to agree that you trust them.  This is a 
> problem of 
> >> the Inno Setup installer kit we use.  An ultra-cautious sysadmin 
> >> could configure Vista to stop you installing via such a program.
> >>
> >>
> >> 2) Permission problems:
> >>
> >> If you install R as an ordinary user (into your own file
> >> space) you should see no permissions problems.  (There would have 
> >> been problems, including under XP, with some recent daily binary 
> >> builds as the inst

Re: [R] Windows Vista issues

[John Fox]

Dear Brian,

I don't yet have a Vista machine, but it occurs to me that some of these
problems might be avoided if R installed by default into c:\R rather than
c:\Program Files\R. Is that the case?

Thank you, by the way, for pursuing these issues.

John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Prof 
> Brian Ripley
> Sent: Saturday, April 14, 2007 1:57 AM
> To: [EMAIL PROTECTED]
> Subject: [R] Windows Vista issues
> 
> It seemed FUD [*] has been prevailing here and elsewhere on 
> Vista security features.
> 
> I asked our sysadmins to set up a Vista box for me on which I 
> have access to all levels of accounts.  Many of the issues I 
> found were covered by earlier answers and all in the upcoming 
> rw-FAQ (currently available at 
> http://www.stats.ox.ac.uk/pub/R/rw-FAQ.html and in the 2.5.0 
> pre-releases) but a quick reprise may help.
> 
> Most of my testing was of 2.5.0 beta, but I did some quick 
> tests of 2.4.1.
> 
> 
> 1) The R installer and uninstaller are from an 'unidentified 
> publisher' 
> and you may have to agree that you trust them.  This is a 
> problem of the Inno Setup installer kit we use.  An 
> ultra-cautious sysadmin could configure Vista to stop you 
> installing via such a program.
> 
> 
> 2) Permission problems:
> 
> If you install R as an ordinary user (into your own file 
> space) you should see no permissions problems.  (There would 
> have been problems, including under XP, with some recent 
> daily binary builds as the installer kit had changed one of 
> its defaults to disallow non-administrator installs, but 
> these have been fixed.)
> 
> I also encountered no problems installing R under the 
> Administrator account (normally hidden) and installing 
> packages under the same account.
> 
> Things are more complicated if you use an account which is in 
> the local administrator group (but is not Administrator 
> itself).  Such accounts are no longer (by default) equivalent 
> to Administrator, and run programs as ordinary user accounts. 
> They need to 'Run as Administrator' to do things in the 
> system area such as C:\Program Files.  You will be asked if 
> you want to run as administrator if you try to install 
> software such as R, but you will not be asked if you try to 
> install packages in the main R library (since asking is 
> something that applies to a program, not part of a particular 
> session).  One simple solution is to elevate your credentials 
> when running an R session to install packages in the same way 
> that you needed to when installing R.  (Unix and MacOS X 
> users will recognize a somewhat automated reincarnation of 'sudo'.)
> 
> It looks like the best practice will be to change the (full) 
> ownership of the R installation to the account used to 
> install it, something which would be standard practice in the 
> Unix world.  Also, we are encouraging people as from 2.5.0 to 
> install packages into a site or personal library where these 
> permission issues should not arise (except when updating 
> recommended packages).
> 
> 
> 3) The most worrying problem is that Vista is reporting quite 
> incorrectly file permissions through the POSIX interfaces 
> used by file.info() and file.access(), and furthermore 
> allowed me as a standard user to create directories in areas 
> over which it says I do not have write permission. We will 
> look further into possible solutions, but it seems the Win32 
> API functions are giving the same answers.
> 
> Problems with 'access' (the C call underlying file.access()) 
> mean that the MinGW compilers do not currently run on Vista 
> without a lot of hoop-jumping.
> 
> 
> [*] http://en.wikipedia.org/wiki/Fear,_uncertainty_and_doubt
> 
> -- 
> Brian D. Ripley,  [EMAIL PROTECTED]
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford, Tel:  +44 1865 272861 (self)
> 1 South Parks Road, +44 1865 272866 (PA)
> Oxford OX1 3TG, UKFax:  +44 1865 272595
> 
> __
> [EMAIL PROTECTED] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] sem is not "taking" the model

[John Fox]

Dear Dimitri,

This works find for me. It's not possible to know what you did wrong
without more information, but one guess (probably wrong) is that you
 didn't enter a blank line to exit from specify.model(), which simply
uses scan() to read the model:

> library(sem)
> model1 <- specify.model()
1: NECESSITY -> n1, NA, 1
2: NECESSITY -> n2, lam_nec_2, NA
3: NECESSITY -> n3, lam_nec_3, NA
4: NECESSITY -> n4, lam_nec_4, NA
5: FRUGALITY -> f1, NA, 1
6: FRUGALITY -> f2, lam_frug_2, NA
7: FRUGALITY -> f3, lam_frug_3, NA
8: TIME -> t1, NA, 1
9: n1 <-> n1, theta_n1, NA
10: n2 <-> n2, theta_n2, NA
11: n3 <-> n3, theta_n3, NA
12: n4 <-> n4, theta_n4, NA
13: f1 <-> f1, theta_f1, NA
14: f2 <-> f2, theta_f2, NA
15: f3 <-> f3, theta_f3, NA
16: t1 <-> t1, NA, 0.414
17: NECESSITY <-> NECESSITY, phi_NN, NA
18: FRUGALITY <-> FRUGALITY, phi_FF, NA
19: TIME <-> TIME, phi_TT, NA
20: NECESSITY <-> TIME, phi_NT, NA
21: NECESSITY <-> FRUGALITY, phi_NF, NA
22: FRUGALITY <-> TIME, phi_FT, NA
23: 
Read 22 records
> 

Regards,
 John

On Wed, 11 Apr 2007 11:41:25 -0700 (PDT)
 John Smith <[EMAIL PROTECTED]> wrote:
> A strange problem with sem:
> 
> I downloaded the sem library and then, I specified my simple
> measurement model (below). I highlighted it and ran it. It ran, but
> it did NOT tell me "22 lines read". And nothing works after that - it
> looks like it runs, but it does not produce anything...
> Did I make a mistake somewhere in the model? (notice, TIME has only 1
> indicator - t1, and I fixed t1's error variance at 0.414.)
> Thank you!
> 
> model1 <- specify.model()
> NECESSITY  -> n1,  NA,  1
> NECESSITY -> n2,  lam_nec_2, NA 
> NECESSITY -> n3,  lam_nec_3, NA
> NECESSITY -> n4,  lam_nec_4, NA
> FRUGALITY  -> f1,  NA,  1
> FRUGALITY  -> f2,  lam_frug_2, NA
> FRUGALITY  -> f3,  lam_frug_3, NA
> TIME   -> t1,  NA,  1
> n1  <-> n1,  theta_n1, NA
> n2  <-> n2,  theta_n2, NA
> n3  <-> n3,  theta_n3, NA
> n4  <-> n4,  theta_n4, NA
> f1  <-> f1,  theta_f1, NA
> f2  <-> f2,  theta_f2, NA
> f3  <-> f3,  theta_f3, NA
> t1  <-> t1,  NA,  0.414
> NECESSITY <-> NECESSITY, phi_NN, NA
> FRUGALITY <-> FRUGALITY, phi_FF, NA
> TIME  <-> TIME,  phi_TT, NA
> NECESSITY <-> TIME,  phi_NT, NA
> NECESSITY <-> FRUGALITY, phi_NF, NA
> FRUGALITY <-> TIME,  phi_FT, NA
> 
> 
>
>
____
> Finding fabulous fares is fun.
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] creating a path diagram in sem

[John Fox]

Dear Dmitri,

path.diagram() takes a fitted "sem" object as its first argument, not a
model-specification object:

> path.diagram(sem.anxiety1,minrank='a1,a2,a3,d1,d2,d3,f1,f2,f3',
maxrank='ANXIETY,DEPRESS,FEAR')

digraph "sem.anxiety1" {
  rankdir=LR;
  size="8,8";
  node [fontname="Helvetica" fontsize=14 shape=box];
  edge [fontname="Helvetica" fontsize=10];
  center=1;
  "ANXIETY" [shape=ellipse]
  "DEPRESS" [shape=ellipse]
  "FEAR" [shape=ellipse]
  "ANXIETY" -> "a1" [label=""];
  "ANXIETY" -> "a2" [label="lam_anx_2"];
  "ANXIETY" -> "a3" [label="lam_anx_3"];
  "DEPRESS" -> "d1" [label=""];
  "DEPRESS" -> "d2" [label="lam_dep_2"];
  "DEPRESS" -> "d3" [label="lam_dep_3"];
  "FEAR" -> "f1" [label=""];
  "FEAR" -> "f2" [label="lam_fear_2"];
  "FEAR" -> "f3" [label="lam_fear_3"];
}

This isn't terribly useful, by the way, unless you have the graphviz
software (reference on ?path.diagram) to draw the graph.

I hope this helps,
 John


On Wed, 11 Apr 2007 09:52:37 -0700 (PDT)
 John Smith <[EMAIL PROTECTED]> wrote:
> Hello,
> I finally run my measurement model in sem - successfully. Now, I am
> trying to print out the path diagram that is based on the results -
> but for some reason it's not working. Below is my script - but the
> problem is probably in my very last line:
> 
> # ANALYSIS OF ANXIETY, DEPRESSION, AND FEAR - LISREL P.31
> library(sem)
> # Creating the ANXIETY, DEPRESSION, AND FEAR intercorrelation matrix
> (KM)
> KM<-matrix(
> c(1,.8,.83,.2,.21,.19,.18,.18,.18,
> 0,1,.81,.22,.24,.18,.19,.19,.21,
> 0,0,1,.22,.19,.2,.2,.2,.22,
> 0,0,0,1,.84,.82,.22,.22,.21,
> 0,0,0,0,1,.84,.19,.18,.19,
> 0,0,0,0,0,1,.18,.18,.18,
> 0,0,0,0,0,0,1,.84,.82,
> 0,0,0,0,0,0,0,1,.81,
> 0,0,0,0,0,0,0,0,1), 9, 9)
> 
> # Creating the ANXIETY, DEPRESSION, AND FEAR SDs vector (SD)
> SD<-c(1.5, 2.4, 3.2, 2.3, 2.3, 2.6, 4.5, 4.7, 5.6)
> 
> # Calculating the Var-Covar matrix based on correlations and SDs
> 
> COVAR<-outer(SD, SD) * KM
> 
> # Providing variable names
>
rownames(COVAR)<-colnames(COVAR)<-c('a1','a2','a3','d1','d2','d3','f1','f2','f3')
>  
> # Specifying the measurement model 1
> model1 <- specify.model()
> ANXIETY -> a1, NA, 1
> ANXIETY -> a2, lam_anx_2, NA 
> ANXIETY -> a3, lam_anx_3, NA
> DEPRESS -> d1, NA, 1
> DEPRESS -> d2, lam_dep_2, NA
> DEPRESS -> d3, lam_dep_3, NA
> FEAR -> f1, NA, 1
> FEAR -> f2, lam_fear_2, NA
> FEAR -> f3, lam_fear_3, NA
> a1 <-> a1, theta_a1, NA
> a2 <-> a2, theta_a2, NA
> a3 <-> a3, theta_a3, NA
> d1 <-> d1, theta_d1, NA
> d2 <-> d2, theta_d2, NA
> d3 <-> d3, theta_d3, NA
> f1 <-> f1, theta_f1, NA
> f2 <-> f2, theta_f2, NA
> f3 <-> f3, theta_f3, NA
> ANXIETY <-> ANXIETY, phi_AA, NA
> DEPRESS <-> DEPRESS, phi_DD, NA
> FEAR <-> FEAR, phi_FF, NA
> ANXIETY <-> FEAR, phi_AF, NA
> ANXIETY <-> DEPRESS, phi_AD, NA
> DEPRESS <-> FEAR, phi_DF, NA
> 
> # Running the estimation using sem:
> sem.anxiety1<-sem(model1, COVAR, N=150, par.size="startvalues")
> 
> # Looking at the results:
> summary(sem.anxiety1)
> 
> # Calling modification indices:
> mod.indices(sem.anxiety1)
> summary(mod.indices(sem.anxiety1))
> 
> # Creating a path diagram
> path.diagram(model1,minrank='a1,a2,a3,d1,d2,d3,f1,f2,f3',
> maxrank='ANXIETY,DEPRESS,FEAR')
> 
> Thank you very much for your help!
> Dimitri
> 
> 
>
>

> TV dinner still cooling? 
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] sem vs. LISREL: sem fails

[John Fox]

Dear Dimitri,

Using the tools provided in the sem package [e.g., the debug argument to
sem()], I took a closer look at this convergence problem, which arises
because of poor start values selected for the variances of the latent
variables. It turns out that you can get a solution for your original model
fit to the covariance matrix by setting the argument par.size to
"startvalues" (as suggested in ?sem): 

sem.anxiety <- sem(model, COVAR, N=150, par.size="startvalues")

Regards,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: John Fox [mailto:[EMAIL PROTECTED] 
> Sent: Monday, April 09, 2007 9:46 PM
> To: 'John Smith'
> Cc: 'r-help@stat.math.ethz.ch'
> Subject: RE: [R] sem vs. LISREL: sem fails
> 
> Dear Dimitri,
> 
> You haven't done anything wrong: Your model is a 
> straightforward confirmatory factor analysis and it is 
> correctly specified. I suspect that sem() is picking poor 
> start values.
> 
> I can get a solution by specifying two alternative models 
> that are equivalent to yours:
> 
> (1) Rescaling the problem by using correlation-matrix input:
> 
> sem.anxiety <- sem(model, KM, N=150)
> 
> (2) Fixing the variances of the factors (latent variables) 
> rather than using reference indicators:
> 
> model.2 <- specify.model()
> ANXIETY -> a1,  lam_anx_1,  NA
> ANXIETY -> a2,  lam_anx_2,  NA
> ANXIETY -> a3,  lam_anx_3,  NA
> DEPRESS -> d1,  lam_dep_1,  NA
> DEPRESS -> d2,  lam_dep_2,  NA
> DEPRESS -> d3,  lam_dep_3,  NA
> FEAR -> f1, lam_fear_1, NA
> FEAR -> f2, lam_fear_2, NA
> FEAR -> f3, lam_fear_3, NA
> a1 <-> a1,  theta_a1,   NA
> a2 <-> a2,  theta_a2,   NA
> a3 <-> a3,  theta_a3,   NA
> d1 <-> d1,  theta_d1,   NA
> d2 <-> d2,  theta_d2,   NA
> d3 <-> d3,  theta_d3,   NA
> f1 <-> f1,  theta_f1,   NA
> f2 <-> f2,  theta_f2,   NA
> f3 <-> f3,  theta_f3,   NA
> ANXIETY <-> ANXIETY,NA, 1
> DEPRESS <-> DEPRESS,NA, 1
> FEAR <-> FEAR,  NA, 1
> ANXIETY <-> FEAR,   phi_AF, NA
> ANXIETY <-> DEPRESS,phi_AD, NA
> DEPRESS <-> FEAR,   phi_DF, NA
> 
> # Running the estimation using sem:
> sem.anxiety.2 <- sem(model.2, COVAR, N=150)
> 
> 
> A couple of small points unrelated to the problem you 
> experienced: (1) You didn't need to load the MASS package, 
> since you didn't appear to use anything in it; (2) comments 
> in R are prefixed by #, not !.
> 
> I hope this helps,
>  John
> 
> 
> John Fox
> Department of Sociology
> McMaster University
> Hamilton, Ontario
> Canada L8S 4M4
> 905-525-9140x23604
> http://socserv.mcmaster.ca/jfox
>  
> 
> > -Original Message-
> > From: [EMAIL PROTECTED] 
> > [mailto:[EMAIL PROTECTED] On Behalf Of John Smith
> > Sent: Monday, April 09, 2007 5:28 PM
> > To: r-help@stat.math.ethz.ch
> > Subject: [R] sem vs. LISREL: sem fails
> > 
> > I am new to R. 
> > I just tried to recreate in R (using sem package and the identical 
> > input data) a solution for a simple measurment model I have found 
> > before in LISREL. LISREL had no problems and converged in just 3 
> > iterations.
> > In sem, I got no solution, just the warning message:
> > 
> > "Could not compute QR decomposition of Hessian.
> > Optimization probably did not converge.
> >  in: sem.default(ram = ram, S = S, N = N, param.names = pars, 
> > var.names = vars, "
> > 
> > What does it mean? Maybe I am doing something wrong?
> > 
> > I have 3 latent factors (Anxiety, Depression, and Fear) - 
> each of them 
> > has 3 observed indicators (a1, a2, a3; d1, d2, d3, and f1, f2, f3) 
> > Below is my script in R:
> > 
> > ! ANALYSIS OF ANXIETY, DEPRESSION, AND FEAR - LISREL P.31
> > 
> > ! Creating the ANXIETY, DEPRESSION, AND FEAR 
> intercorrelation matrix 
> > (KM):
> > KM<-matrix(
> > c(1,.8,.83,.2,.21,.19,.18,.18,.18,
> > 0,1,.81,.22,.24,.18,.19,.19,.21,
> > 0,0,1,.22,.19,.2,.2,.2,.22,
> > 0,0,0,1,.84,.82,.22,.22,.21,
> > 0,0,0,0,1,.84,.19,.18,.19,
> > 0,0,0,0,0,1,.18,.18,.18,
> > 0,0,0,0,0,0,1,.84,.82,
> > 0,0,0,0,0,0,0,1,.81,
> > 0,0,0,0,0,0,0,0,1)

Re: [R] sem vs. LISREL: sem fails

[John Fox]

Dear Dimitri,

> -Original Message-
> From: John Smith [mailto:[EMAIL PROTECTED] 
> Sent: Tuesday, April 10, 2007 8:11 AM
> To: John Fox
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] sem vs. LISREL: sem fails
> 
> Dear Scott,

(John, actually.)

> thanks a lot - all your comments are extremely valuable!
> One quick follow-up question - the thing is, I was really 
> trying to compare sem and LISREL solutions, which is why I 
> was trying to specify exactly the same model. Hence, I'd 
> prefer to stick to the Var-Covar matrix, rather than 
> Correlation matrix. Also, I really like fixing the metric for 
> latent variables along the metric of one of the observed 
> variables - just for practical and didactical reasons 
> (although I'll try to run your suggestion as well - in both 
> sem and LISREL).

You should be able to do what you want to do; that you can't is a weakness
of sem(). When I have some time, I'll take another look at how sem()
computes start values.

> Any advice/guidance on picking the starting values for sem if 
> I still want to estimate my original model or run into a 
> similar problem in the future? I know how to do it (in terms 
> of entering something else but NA in the last column) - but 
> what starting values should I pick?

Of course, if I knew the answer to this question I'd implement it. What you
could do, however, is start with one of the solutions that I suggested and
then rescale to the metric that you want. I've thought about implementing
that internally in sem() -- rescaling to correlations and then translating
back at the end -- but have to figure out how to handle equality
constraints.

Regards,
 John

> Thanks a lot again for your help!
> Dimitri
> 
> 
> - Original Message 
> From: John Fox <[EMAIL PROTECTED]>
> To: John Smith <[EMAIL PROTECTED]>
> Cc: r-help@stat.math.ethz.ch
> Sent: Monday, April 9, 2007 9:46:03 PM
> Subject: RE: [R] sem vs. LISREL: sem fails
> 
> 
> Dear Dimitri,
> 
> You haven't done anything wrong: Your model is a 
> straightforward confirmatory factor analysis and it is 
> correctly specified. I suspect that
> sem() is picking poor start values.
> 
> I can get a solution by specifying two alternative models 
> that are equivalent to yours:
> 
> (1) Rescaling the problem by using correlation-matrix input:
> 
> sem.anxiety <- sem(model, KM, N=150)
> 
> (2) Fixing the variances of the factors (latent variables) 
> rather than using reference indicators:
> 
> model.2 <- specify.model()
> ANXIETY -> a1,  lam_anx_1,  NA
> ANXIETY -> a2,  lam_anx_2,  NA
> ANXIETY -> a3,  lam_anx_3,  NA
> DEPRESS -> d1,  lam_dep_1,  NA
> DEPRESS -> d2,  lam_dep_2,  NA
> DEPRESS -> d3,  lam_dep_3,  NA
> FEAR -> f1, lam_fear_1, NA
> FEAR -> f2, lam_fear_2, NA
> FEAR -> f3, lam_fear_3, NA
> a1 <-> a1,  theta_a1,   NA
> a2 <-> a2,  theta_a2,   NA
> a3 <-> a3,  theta_a3,   NA
> d1 <-> d1,  theta_d1,   NA
> d2 <-> d2,  theta_d2,   NA
> d3 <-> d3,  theta_d3,   NA
> f1 <-> f1,  theta_f1,   NA
> f2 <-> f2,  theta_f2,   NA
> f3 <-> f3,  theta_f3,   NA
> ANXIETY <-> ANXIETY,NA, 1
> DEPRESS <-> DEPRESS,NA, 1
> FEAR <-> FEAR,  NA, 1
> ANXIETY <-> FEAR,   phi_AF, NA
> ANXIETY <-> DEPRESS,phi_AD, NA
> DEPRESS <-> FEAR,   phi_DF, NA
> 
> # Running the estimation using sem:
> sem.anxiety.2 <- sem(model.2, COVAR, N=150)
> 
> 
> A couple of small points unrelated to the problem you 
> experienced: (1) You didn't need to load the MASS package, 
> since you didn't appear to use anything in it; (2) comments 
> in R are prefixed by #, not !.
> 
> I hope this helps,
> John
> 
> 
> John Fox
> Department of Sociology
> McMaster University
> Hamilton, Ontario
> Canada L8S 4M4
> 905-525-9140x23604
> http://socserv.mcmaster.ca/jfox
>  
> 
> > -Original Message-
> > From: [EMAIL PROTECTED] 
> > [mailto:[EMAIL PROTECTED] On Behalf Of John Smith
> > Sent: Monday, April 09, 2007 5:28 PM
> > To: r-help@stat.math.ethz.ch
> > Subject: [R] sem vs. LISREL: sem fails
> > 
> > I am new to R. 
> > I just tried to recreate in R (using sem package and the identical 
> > input data) a solution for a simple measurment model I have found 
> > before in LISREL. LISREL had no problems and converged in just 3 
> > iterations.
>

Re: [R] Dealing with large nominal predictor in sem package

[John Fox]

Dear adschai,

> -Original Message-
> From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] 
> Sent: Monday, April 09, 2007 8:30 PM
> To: John Fox; r-help@stat.math.ethz.ch
> Subject: Re: RE: [R] Dealing with large nominal predictor in 
> sem package
> 
> Hi John,
> 
> Additional two questions on this sem package:
> (1) The tsls is based on maximum likelihood or OLS?

Neither; it does two-stage least-squares (2SLS).

> (2) I am trying to find goodness of fit for the result of 
> tsls. Somehow, I don't see it in the documentation. Would you 
> please provide some examples? 

I'm not sure what you mean by goodness of fit. If you have in mind an
R^2-like measure; you can always use 1 -
error-variance/variance-of-endogenous-variable, but this is not guaranteed
to be positive.

> (3) If I would like to diagnostic of model selection, says 
> use AIC criteria, it is a bit unclear for me how I can apply 
> this on structural equation model as it is composed of 
> multiple equations rather than one. And is there any 
> functionality in sem that does this? 

Since there's no likelihood for 2SLS estimation, I don't see how you could
get an AIC. On the other hand, sem() fits by full-information
maximum-likelihood (FIML). It prints out the BIC; you could compute the AIC
if you liked.

John

> Any help would be really appreciated. Thank you.
> 
> - adschai
> 
> - Original Message -
> From: John Fox
> Date: Monday, April 9, 2007 8:04 am
> Subject: RE: [R] Dealing with large nominal predictor in sem package
> To: [EMAIL PROTECTED]
> Cc: r-help@stat.math.ethz.ch
> 
> > Dear adschai,
> > 
> > It's not possible to know from your description exactly what you're 
> > doing, but perhaps the following will help:
> > 
> > (1) I presume that your nominal variable is exogenous, 
> since otherwise 
> > it wouldn't be sensible to use 2SLS.
> > 
> > (2) You don't have to make your own dummy regressors for a nominal 
> > variable; just represent it in the model as a factor as you would, 
> > e.g., in lm().
> > 
> > (3) Do you have at least as many instrumental variables 
> (including the 
> > dummy
> > regressors) as there are structural coefficients to 
> estimate? If not, 
> > the structural equation is underidentified, which will produce the 
> > error that you've encountered.
> > 
> > I hope this helps,
> > John
> > 
> > 
> > John Fox
> > Department of Sociology
> > McMaster University
> > Hamilton, Ontario
> > Canada L8S 4M4
> > 905-525-9140x23604
> > http://socserv.mcmaster.ca/jfox
> > 
> > 
> > > -Original Message-
> > > From: [EMAIL PROTECTED] 
> > > [mailto:[EMAIL PROTECTED] On Behalf Of 
> > > [EMAIL PROTECTED]
> > > Sent: Sunday, April 08, 2007 11:07 PM
> > > To: r-help@stat.math.ethz.ch
> > > Subject: [R] Dealing with large nominal predictor in sem package
> > > 
> > > Hi,
> > > 
> > > I am using tsls function from sem package to estimate a 
> model which 
> > > includes large number of data. Among its predictors, it has a 
> > > nominal data which has about 10 possible values. So I expand this 
> > > parameter into 9-binary-value predictors with the coefficient of 
> > > base value equals 0. I also have another continuous predictor.
> > > 
> > > The problem is that, whenever I run the tsls, I will get 
> 'System is 
> > > computationally singular' error all the time. I'm 
> wondering if there 
> > > is anyway that I can overcome this problem? Please kindly 
> suggest. 
> > > Thank you so much in advance.
> > > 
> > > - adschai
> > > 
> > > [[alternative HTML version deleted]]
> > > 
> > > __
> > > R-help@stat.math.ethz.ch mailing list 
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide
> > > http://www.R-project.org/posting-guide.html
> > > and provide commented, minimal, self-contained, reproducible code.
> > > 
> > 
> > 
> > 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Dealing with large nominal predictor in sem package

[John Fox]

Dear adschai,

> -Original Message-
> From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] 
> Sent: Monday, April 09, 2007 7:33 PM
> To: John Fox
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: RE: [R] Dealing with large nominal predictor in 
> sem package
> 
> Hi John,
> 
> Thank you. I think (2) from your explanation hits the right 
> point. The reason is that when I made my own dummy variables 
> and my original nominal variable has 10 possible values, it 
> makes each each observed exogeneous variable vector of mine 
> has 9 zeros and 1 one value. And I have about 40 
> observations. So it will make the matrix almost zero.
> 

I'm afraid that I don't follow that, unless you're saying that some of the
levels of the factor have very few observations in them.

> One more question. If I have a nominal response, I guess the 
> tsls would no longer work. How can I go around with this? 

If the response is ordinal, then you can use sem() with
polyserial/polychoric correlations. Otherwise, the sem package won't handle
it.

> Says, I have 3 equations in my structure model whose 
> responses are continuous whereas another one has multinominal 
> response. Thank you so much.
> 

As I said, neither tsls() nor sem() will handle an unordered response.

John

> - adschai
> 
> - Original Message -
> From: John Fox
> Date: Monday, April 9, 2007 8:04 am
> Subject: RE: [R] Dealing with large nominal predictor in sem package
> To: [EMAIL PROTECTED]
> Cc: r-help@stat.math.ethz.ch
> 
> > Dear adschai,
> > 
> > It's not possible to know from your description exactly what you're 
> > doing, but perhaps the following will help:
> > 
> > (1) I presume that your nominal variable is exogenous, 
> since otherwise 
> > it wouldn't be sensible to use 2SLS.
> > 
> > (2) You don't have to make your own dummy regressors for a nominal 
> > variable; just represent it in the model as a factor as you would, 
> > e.g., in lm().
> > 
> > (3) Do you have at least as many instrumental variables 
> (including the 
> > dummy
> > regressors) as there are structural coefficients to 
> estimate? If not, 
> > the structural equation is underidentified, which will produce the 
> > error that you've encountered.
> > 
> > I hope this helps,
> > John
> > 
> > 
> > John Fox
> > Department of Sociology
> > McMaster University
> > Hamilton, Ontario
> > Canada L8S 4M4
> > 905-525-9140x23604
> > http://socserv.mcmaster.ca/jfox
> > 
> > 
> > > -Original Message-
> > > From: [EMAIL PROTECTED] 
> > > [mailto:[EMAIL PROTECTED] On Behalf Of 
> > > [EMAIL PROTECTED]
> > > Sent: Sunday, April 08, 2007 11:07 PM
> > > To: r-help@stat.math.ethz.ch
> > > Subject: [R] Dealing with large nominal predictor in sem package
> > > 
> > > Hi,
> > > 
> > > I am using tsls function from sem package to estimate a 
> model which 
> > > includes large number of data. Among its predictors, it has a 
> > > nominal data which has about 10 possible values. So I expand this 
> > > parameter into 9-binary-value predictors with the coefficient of 
> > > base value equals 0. I also have another continuous predictor.
> > > 
> > > The problem is that, whenever I run the tsls, I will get 
> 'System is 
> > > computationally singular' error all the time. I'm 
> wondering if there 
> > > is anyway that I can overcome this problem? Please kindly 
> suggest. 
> > > Thank you so much in advance.
> > > 
> > > - adschai
> > > 
> > > [[alternative HTML version deleted]]
> > > 
> > > __
> > > R-help@stat.math.ethz.ch mailing list 
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide
> > > http://www.R-project.org/posting-guide.html
> > > and provide commented, minimal, self-contained, reproducible code.
> > > 
> > 
> > 
> > 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] sem vs. LISREL: sem fails

[John Fox]

Dear Dimitri,

You haven't done anything wrong: Your model is a straightforward
confirmatory factor analysis and it is correctly specified. I suspect that
sem() is picking poor start values.

I can get a solution by specifying two alternative models that are
equivalent to yours:

(1) Rescaling the problem by using correlation-matrix input:

sem.anxiety <- sem(model, KM, N=150)

(2) Fixing the variances of the factors (latent variables) rather than using
reference indicators:

model.2 <- specify.model()
ANXIETY -> a1,  lam_anx_1,  NA
ANXIETY -> a2,  lam_anx_2,  NA
ANXIETY -> a3,  lam_anx_3,  NA
DEPRESS -> d1,  lam_dep_1,  NA
DEPRESS -> d2,  lam_dep_2,  NA
DEPRESS -> d3,  lam_dep_3,  NA
FEAR -> f1, lam_fear_1, NA
FEAR -> f2, lam_fear_2, NA
FEAR -> f3, lam_fear_3, NA
a1 <-> a1,  theta_a1,   NA
a2 <-> a2,  theta_a2,   NA
a3 <-> a3,  theta_a3,   NA
d1 <-> d1,  theta_d1,   NA
d2 <-> d2,  theta_d2,   NA
d3 <-> d3,  theta_d3,   NA
f1 <-> f1,  theta_f1,   NA
f2 <-> f2,  theta_f2,   NA
f3 <-> f3,  theta_f3,   NA
ANXIETY <-> ANXIETY,NA, 1
DEPRESS <-> DEPRESS,NA, 1
FEAR <-> FEAR,  NA, 1
ANXIETY <-> FEAR,   phi_AF, NA
ANXIETY <-> DEPRESS,phi_AD, NA
DEPRESS <-> FEAR,   phi_DF, NA

# Running the estimation using sem:
sem.anxiety.2 <- sem(model.2, COVAR, N=150)


A couple of small points unrelated to the problem you experienced: (1) You
didn't need to load the MASS package, since you didn't appear to use
anything in it; (2) comments in R are prefixed by #, not !.

I hope this helps,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of John Smith
> Sent: Monday, April 09, 2007 5:28 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] sem vs. LISREL: sem fails
> 
> I am new to R. 
> I just tried to recreate in R (using sem package and the 
> identical input data) a solution for a simple measurment 
> model I have found before in LISREL. LISREL had no problems 
> and converged in just 3 iterations. 
> In sem, I got no solution, just the warning message:
> 
> "Could not compute QR decomposition of Hessian.
> Optimization probably did not converge.
>  in: sem.default(ram = ram, S = S, N = N, param.names = pars, 
> var.names = vars, "
> 
> What does it mean? Maybe I am doing something wrong?
> 
> I have 3 latent factors (Anxiety, Depression, and Fear) - 
> each of them has 3 observed indicators (a1, a2, a3; d1, d2, 
> d3, and f1, f2, f3) Below is my script in R:
> 
> ! ANALYSIS OF ANXIETY, DEPRESSION, AND FEAR - LISREL P.31
> 
> ! Creating the ANXIETY, DEPRESSION, AND FEAR intercorrelation 
> matrix (KM):
> KM<-matrix(
> c(1,.8,.83,.2,.21,.19,.18,.18,.18,
> 0,1,.81,.22,.24,.18,.19,.19,.21,
> 0,0,1,.22,.19,.2,.2,.2,.22,
> 0,0,0,1,.84,.82,.22,.22,.21,
> 0,0,0,0,1,.84,.19,.18,.19,
> 0,0,0,0,0,1,.18,.18,.18,
> 0,0,0,0,0,0,1,.84,.82,
> 0,0,0,0,0,0,0,1,.81,
> 0,0,0,0,0,0,0,0,1), 9, 9)
> 
> ! Creating the ANXIETY, DEPRESSION, AND FEAR Standard 
> Deviations vector (SD):
> SD<-c(1.5, 2.4, 3.2, 2.3, 2.3, 2.6, 4.5, 4.7, 5.6)
> 
> ! Calculating the Var-Covar matrix based on correlations and SDs:
> library(MASS)
> COVAR<-outer(SD, SD) * KM
> 
> ! Creating variable names
> rownames(COVAR)<-colnames(COVAR)<-c('a1','a2','a3','d1','d2','
> d3','f1','f2','f3')
> 
> ! Specifying the measurement model to estimate:
> model<-specify.model()
> ANXIETY -> a1, NA, 1
> ANXIETY -> a2, lam_anx_2, NA 
> ANXIETY -> a3, lam_anx_3, NA
> DEPRESS -> d1, NA, 1
> DEPRESS -> d2, lam_dep_2, NA
> DEPRESS -> d3, lam_dep_3, NA
> FEAR -> f1, NA, 1
> FEAR -> f2, lam_fear_2, NA
> FEAR -> f3, lam_fear_3, NA
> a1 <-> a1,theta_a1, NA
> a2 <-> a2,theta_a2, NA
> a3 <-> a3,theta_a3, NA
> d1 <-> d1,   theta_d1, NA
> d2 <-> d2,   theta_d2, NA
> d3 <-> d3,   theta_d3, NA
> f1 <-> f1, theta_f1, NA
> f2 <-> f2, theta_f2, NA
> f3 <-> f3, theta_f3, NA
> ANXIETY <-> ANXIETY,   phi_AA, NA
> DEPRESS <-> DEPRESS, phi_DD, NA

Re: [R] Dealing with large nominal predictor in sem package

[John Fox]

Dear adschai,

It's not possible to know from your description exactly what you're doing,
but perhaps the following will help: 

(1) I presume that your nominal variable is exogenous, since otherwise it
wouldn't be sensible to use 2SLS. 

(2) You don't have to make your own dummy regressors for a nominal variable;
just represent it in the model as a factor as you would, e.g., in lm(). 

(3) Do you have at least as many instrumental variables (including the dummy
regressors) as there are structural coefficients to estimate? If not, the
structural equation is underidentified, which will produce the error that
you've encountered.

I hope this helps,
 John

----
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of 
> [EMAIL PROTECTED]
> Sent: Sunday, April 08, 2007 11:07 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Dealing with large nominal predictor in sem package
> 
> Hi,
> 
> I am using tsls function from sem package to estimate a model 
> which includes large number of data. Among its predictors, it 
> has a nominal data which has about 10 possible values. So I 
> expand this parameter into 9-binary-value predictors with the 
> coefficient of base value equals 0. I also have another 
> continuous predictor. 
> 
> The problem is that, whenever I run the tsls, I will get 
> 'System is computationally singular' error all the time. I'm 
> wondering if there is anyway that I can overcome this 
> problem? Please kindly suggest. Thank you so much in advance.
> 
> - adschai
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] MANOVA with repeated measurements

[John Fox]

Dear Mark,

Take a look at the Anova() function in the car package. There's an example
of a repeated-measures MANOVA in ?Anova (the last example). Also see
?anova.mlm.

I hope this helps,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Mark 
> E. Olszewski
> Sent: Thursday, April 05, 2007 5:03 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] MANOVA with repeated measurements
> 
> Hello,
> 
> I have a question regarding performing manova.  I have an 
> experiment where I want to measure 10 output variables for 3 
> different measurement methods.  Since each of the methods 
> requires some user interaction, I would also like to include 
> repeated measures for each of the output variables to include 
> intraobserver variability in the design.  
> 
> How can I perform such a repeated-measures manova using R?  I 
> think I have found how to perform the non-repeated version of 
> this experiment based on previous posts, but I am having 
> trouble generalizing to the repeated-measures version (maybe 
> I am just too tired today).
> 
> Thanks in advance,
> 
> Mark
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] about systemfit

[John Fox]

Dear Martin,

Note that the F-test provided by linear.hypothesis() in the car package
appears to agree with the Wald test provided by both linear.hypothesis() and
waldtest.systemfit(), while ftest.systemfit() gives a much smaller p-value
and an F that's even larger than the Wald chisquare. This suggests to me
that the first place to look for a problem is in ftest.systemfit().

Regards,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Martin Ivanov
> Sent: Thursday, April 05, 2007 1:02 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] about systemfit
> 
> Thank you very much for your responsiveness. Here are the 
> tests that show the same results, as they must:
> 
> linear.hypothesis(model=adfResulm,hypothesis.matrix=Rrestr10,t
> est="Chisq"):
>  Res.Df RSS  Df Sum of Sq  Chisq Pr(>Chisq)
> 1127  7.3782
> 2137  7.6848 -10   -0.3066 5.2769  0.872
> 
> waldtest.systemfit(object=adfResu,R.restr=Rrestr10):
> 
> Wald-test for linear parameter restrictions in equation systems
> Wald-statistic: 5.277
> degrees of freedom: 10
> p-value: 0.8719 
> 
> So the Chisq test with linear.hypothesis from car and the 
> waldtest.systemfit give the same output. But the  F test with 
> linear.hypothesis from car and ftest.systemfit give 
> absolutely different results, as I demonstrated in my first 
> post. I have no idea why.
> 
> Regards,
> Martin
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] about systemfit

[John Fox]

Dear Martin,

I'll address only part of your question, which is how to get the code for
linear.hypothesis() in the car package: linear.hypothesis() is an S3 generic
function with several methods:

> library(car)
> methods("linear.hypothesis")
[1] linear.hypothesis.default* linear.hypothesis.glm*
[3] linear.hypothesis.lm*  linear.hypothesis.mlm*

   Non-visible functions are asterisked

So the function has methods for glm, lm, and mlm objects, along with a
default method, which, as explained in ?linear.hypothesis, should work with
models for which coef() and vcov() methods exist.

To see linear.hypothesis.default(), which isn't exported from the car
namespace, you can, e.g., type

 car:::linear.hypothesis.default

I'm not sure why you get different results for the F-test but similar
results for the Wald test. Which F-test seems to square with the Wald test
(you don't show the Wald tests)?

I hope this helps,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Martin Ivanov
> Sent: Thursday, April 05, 2007 11:28 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] about systemfit
> 
> Hello. I am still a newbie in R. Excuse me if I am asking 
> something obvious. My efforts to get an answer through 
> browsing the mailing archives failed. I want to perform an 
> augmented Dickey-Fuller test and to obtain AIC and BIC and to 
> be able to impose some linear restrictions on the ADF 
> regression so as to decide the correct order of 
> autoregression. However I could find no obvious way to impose 
> linear restrictions or to obtain AIC from the the result of 
> ADF.test from uroot. That is why I turned to systemfit. I ran 
> the ADF regression with systemfit and obtained the same 
> coefficient estimates as through ADF.test (as it had to be). 
> Unfortunately I could not find how to extract AIC from the 
> result of systemfit, so I evaluated the ADF regression by lm. 
> So far so good. However the results of ftest.systemfit and 
> linear.hypothesis from the "car" package are very different, 
> while the results from waldtest.systemfit and 
> linear.hypothesis coincide. I have no explanation for this 
> issue a  nd I could not see the code of linear.hypothesis. 
> When I type "linear.hypothesis" I get:
> function (model, ...)
> UseMethod("linear.hypothesis")
> 
> 
> but when I type "ftest.systemfit", I do see the actual code. Why? 
> Anyway, here are the results in more detail:
> 
> eqns<-list(eq = y ~ trend+ x1 + x[,1] + x[,2] + x[,3] + x[,4] 
> + x[,5] + x[,6] + x[,7] + x[,8] + x[,9] + x[,10] + x[,11] + 
> x[,12] + x[,13]
> 
> Rrestr10<-matrix(0,10,16);Rrestr10[1,16]=Rrestr10[2,15]=Rrestr
10[3,14]=Rrestr10[4,13]=Rrestr10[5,12]=Rrestr10[6,11]=Rrestr10[7,10]=Rrestr1
0[8,9]=Rrestr10[9,8]=Rrestr10[10,7]=1
> 
> adfResc<-systemfit(method="OLS",eqns=eqns,R.restr=Rrestr10)
> 
> adfResu<-systemfit(method="OLS",eqns=eqns)
> 
> adfResulm<-lm(formula=eqns$eq)
> 
> ftest.systemfit( object=adfResu, R.restr=Rrestr10) :
> 
> F-test for linear parameter restrictions in equation systems
> F-statistic: 9.083
> degrees of freedom of the numerator: 10 degrees of freedom of 
> the denominator: 127
> p-value: 3.449e-11 
> 
> linear.hypothesis(model=adfResulm,hypothesis.matrix=Rrestr10,t
> est="F"):
> 
>  Res.Df RSS  Df Sum of Sq  F Pr(>F)
> 1127  7.3782
> 2137  7.6848 -10   -0.3066 0.5277  0.868
> 
> As I said, the results of
> the chisquare test with linear.hypothesis and the 
> waldtest.systemfit coincide.
> I have one more problem. This is the output of lrtest.systemfit:
> lrtest.systemfit(resultc=adfResc,resultu=adfResu)
> 
> Likelihood-Ratio-test for parameter restrictions in equation systems
> LR-statistic:  
> degrees of freedom: 
> p-value:  
> 
> Why do I get empty values?
> In summary, I need to understand why the two ftests give 
> different results; why lrtest.systemfit gives empty output; 
> is there some way to extract AIC and BIC from object of class 
> systemfit or from the result of ADF.test.
> 
> Excuse me if I am asking something too obvious, but I am 
> really at a loss.
> 
> Any suggestions on any of the above questions will be welcomed.
> 
> Regards,
> Martin
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read

Re: [R] p value for coefficients in multinomial model

[John Fox]

Dear Xingwang Ye,

As it says on the help page for Anova(), "Be very careful in formulating the
model for type-III tests, or the hypotheses tested will not make sense."

 It looks as if you fit the model with the default contrast coding for the
factors, contr.treatment. To get sensible "Type-III" sums of squares, you
have to use contrasts that are orthogonal in the row basis of the model
matrix, such as contr.sum. 

SAS takes a different approach from R, using a deficient-rank
parametrization of the model, and then in effect placing a restriction on
the parameters that is equivalent to contr.treatment; it calculates
sums-of-squares in a manner that's independent of the parametrization, while
Anova() does not.

I hope this helps,
 John

----
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Xingwang Ye
> Sent: Thursday, April 05, 2007 9:01 AM
> To: r-help
> Subject: [R] p value for coefficients in multinomial model
> 
> Dear all,
>
> 1)how can I easily get p value for the coefficients of 
> factors in a multinomial model?
> 2)why the p values for "type III" test with Anova are not 
> identical to that from SAS?
>  
> for example:  
> 
>  A,B and C are categorical variables,but the proportions of 
> each level in each categorical variables are not balance. Y 
> is a nominal variables(>=3 categories); 
>   
> To do multinomial logistic regression in SAS, I tried proc 
> logistic data=a;  
> class A B C Y/param = ref ref = first;  
> model Y=A B C/link=logit;
> run;  
> 
> In R, I tried this:  
> library(car) 
> library(nnet)
> b <- as.data.frame(lapply(a[,c("A","B","C","Y")], as.factor))
> mod.mul<- (Y~A+B+C,data=b)  
> summary(mod.mul,cor=F,Wald=T) 
> Anova(mod.mul,type="III")
> 
> I can get the identical coefficients and std.errors in R as 
> in SAS. how can I easily get the p value for each 
> coeefficients? However,the "LR Chisq" for each variables are 
> similar but not identical to that from SAS, why?  
> 
> 
> BTW,please forward your kind answers to my email. Thank you.  
>   
> 
> Best wishes
> yours, sincerely 
> Xingwang Ye
> PhD candidate 
> Research Group of Nutrition Related Cancers and Other Chronic 
> Diseases  
> Institute for Nutritional Sciences,  
> Shanghai Institutes of Biological Sciences, 
> Chinese Academy of Sciences 
> P.O.Box 32 
> 294 Taiyuan Road 
> Shanghai 200031 
> P.R.CHINA
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] rcmdr help

[John Fox]

Dear Daniel,

I'm not a Debian user, but your description of what you did makes me wonder
whether (1) you ever loaded the Rcmdr package, and (2) whether you installed
its dependencies. After starting up R, you can load (and start) the R
commander via the command library(Rcmdr). Prior to that, you should be able
to install the Rcmdr package along with its dependencies from within R via
the command install.packages("Rcmdr", dependencies=TRUE).

Finally, if you haven't already done so, you'd probably save yourself some
time and trouble by reading the "Introduction to R" manual that ships with
R.

I hope this helps,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Daniel Fodor
> Sent: Wednesday, April 04, 2007 3:24 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] rcmdr help
> 
> Dear All,
> 
> I use Debian etch and downloaded R from its repo. Although it 
> is running I would like to make the first steps using rcmdr. 
> I went through on wikis and howtos, installed rcmdr but it is 
> not running. Apparently it is properly installed, there is an 
> Rcmdr library in /usr/lib/R/site-library but I cannot start 
> it neither from command line nor by clicking on applications 
> -> debian 
> -> apps -> math -> GNU R. Clicking on it starts R without 
> GUI. Tcl 8.4 
> -> was
> already installed and there are no broken packages. I would 
> appreciate any kind of help in starting R with gui. Please, 
> excuse me if it is a trivial question, I am newbie to R.
> 
> Best wishes
> 
> Daniel
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] partial R

[John Fox]

Dear Pedram,

I think that you're confusing component+residual ("partial-residual") plots
with added-variable ("partial-regression") plots. The latter shows residuals
for Y and a particular X regressed on all the other X's, and the correlation
in an added-variable plot is therefore the partial correlation between a
particular X and Y "controlling for" all the other X's. This is exactly what
the partial.cor() function computes. If I follow correctly what you want,
you just have to square the partial correlations.

Regards,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: Pedram Rowhani [mailto:[EMAIL PROTECTED] 
> Sent: Monday, April 02, 2007 11:06 AM
> To: John Fox; 'Michael Kubovy'
> Cc: r-help@stat.math.ethz.ch
> Subject: RE: [R] partial R
> 
> Dear John and Michael,
> 
> sorry for the poor explanation of my question.
> 
> What I am looking for is the partial R-squares, to estimate 
> the proportion of unexplained variation of y that becomes 
> explained with the addition of variable x cr.plots does give 
> me the plots but I am looking for a way to get the actual 
> partial R squares, which would correspond to the R squares of 
> those plots
> 
> I hope that things are a bit more clear now
> 
> I guess I can calculate the SSE's and SSR's of the different 
> models myself and then use the formula for partial R square 
> but things would be just simpler if there was already a 
> function implemented in R
> 
> Thanks again
> Pedram
> 
> At 16:20 02/04/2007, John Fox wrote:
> 
> 
>   Dear Michael and Pedram,
>   
>   I'm afraid that Pedram's question is unclear, since the 
> subject line refers
>   to "partial R" (which might have been intended as 
> "partial r," or "partial
>   correlation"), while the message itself refers to 
> "partial regression
>   coefficients." 
>   
>   The latter are simply the coefficients returned by 
> lm(); there is a
>   partial.cor() function in the Rcmdr package that 
> computes partial
>   correlations. It's so simple, that I'll just reproduce it here:
>   
>   partial.cor <-
>   function (X, ...) 
>   {
>   R <- cor(X, ...)
>   RI <- solve(R)
>   D <- 1/sqrt(diag(RI))
>   R <- -RI * (D %o% D)
>   diag(R) <- 0
>   rownames(R) <- colnames(R) <- colnames(X)
>   R
>   }
>   
>   cr.plots() in the car package produces partial-residual plots
>   ("component+residual plots") but returns neither 
> partial correlations nor
>   partial-regression coefficients.
>   
>   Regards,
>John
>   
>   
>   John Fox
>   Department of Sociology
>   McMaster University
>   Hamilton, Ontario
>   Canada L8S 4M4
>   905-525-9140x23604
>   http://socserv.mcmaster.ca/jfox 
> <http://socserv.mcmaster.ca/jfox>  
>    
>   
>   > -Original Message-
>   > From: [EMAIL PROTECTED] 
>   > [ mailto:[EMAIL PROTECTED] 
> <mailto:[EMAIL PROTECTED]> ] On Behalf Of 
> Michael Kubovy
>   > Sent: Monday, April 02, 2007 9:16 AM
>   > To: Pedram Rowhani
>   > Cc: r-help@stat.math.ethz.ch
>   > Subject: Re: [R] partial R
>   > 
>   > 
>   > On Apr 2, 2007, at 5:49 AM, Pedram Rowhani wrote:
>   > 
>   > > i am wondering if there is a command in R that will 
> give me the 
>   > > partial regression coefficients
>   > 
>   > To answer your question, you could have started with 
>   > RSiteSearch('partial regression')
>   > 
>   > It's then likely that you would figured out that one 
> way to proceed is
>   > install.packages('car')
>   > ?cr.plots
>   > 
>   > (You may have to restart R to get the help on a 
> newly-installed
>   > package.)
>   > _
>   > Professor Michael Kubovy
>   > University of Virginia
>   > Department of Psychology
>   > USPS: P.O.Box 400400Charlottesville, VA 22904-4400
>   > Parcels:Room 102Gilmer Hall
>   >  McCormick RoadCharlottesville

Re: [R] partial R

[John Fox]

Dear Michael and Pedram,

I'm afraid that Pedram's question is unclear, since the subject line refers
to "partial R" (which might have been intended as "partial r," or "partial
correlation"), while the message itself refers to "partial regression
coefficients." 

The latter are simply the coefficients returned by lm(); there is a
partial.cor() function in the Rcmdr package that computes partial
correlations. It's so simple, that I'll just reproduce it here:

partial.cor <-
function (X, ...) 
{
R <- cor(X, ...)
RI <- solve(R)
D <- 1/sqrt(diag(RI))
R <- -RI * (D %o% D)
diag(R) <- 0
rownames(R) <- colnames(R) <- colnames(X)
R
}

cr.plots() in the car package produces partial-residual plots
("component+residual plots") but returns neither partial correlations nor
partial-regression coefficients.

Regards,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Michael Kubovy
> Sent: Monday, April 02, 2007 9:16 AM
> To: Pedram Rowhani
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] partial R
> 
> 
> On Apr 2, 2007, at 5:49 AM, Pedram Rowhani wrote:
> 
> > i am wondering if there is a command in R that will give me the 
> > partial regression coefficients
> 
> To answer your question, you could have started with 
> RSiteSearch('partial regression')
> 
> It's then likely that you would figured out that one way to proceed is
> install.packages('car')
> ?cr.plots
> 
> (You may have to restart R to get the help on a newly-installed
> package.)
> _
> Professor Michael Kubovy
> University of Virginia
> Department of Psychology
> USPS: P.O.Box 400400Charlottesville, VA 22904-4400
> Parcels:Room 102Gilmer Hall
>  McCormick RoadCharlottesville, VA 22903
> Office:B011+1-434-982-4729
> Lab:B019+1-434-982-4751
> Fax:+1-434-982-4766
> WWW:http://www.people.virginia.edu/~mk9y/
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Bonferroni p-value greater than 1

[John Fox]

Dear Horace,

The Bonferonni p-value is obtained from the "unadjusted" p-value by
multiplying the latter by the number of observations, and provides a
conservative (although usually quite accurate) outlier test. When the
adjusted p-value exceeds 1 you can take that as an indication that there are
no unusually large studentized residuals (and indeed that the largest
studentized residual is smaller than one would expect under the standard
linear-model assumptions). 

I hope this helps,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Horace Tso
> Sent: Wednesday, March 28, 2007 6:36 PM
> To: 'R R-help'
> Subject: [R] Bonferroni p-value greater than 1
> 
> Hi folks,
> 
> I use the outlier.test in package car to test a lm model and 
> the bonferroni p value returned is shown as NA. When the 
> object is typed it indicates the p value is greater than 1. 
> I'm not sure how to interpret it. 
> 
> Thanks in advance.
> 
> Horace W. Tso
> 
> 
> > outlier.test(mod)$test
> max|rstudent|df  unadjusted p  Bonferroni p
>2.04106376   18.0.05618628NA 
> 
> > outlier.test(mod)
> 
> max|rstudent| = 2.041064, degrees of freedom = 18,
> unadjusted p = 0.05618628, Bonferroni p > 1
> 
> Observation: 1 
> 
> The lm model looks fine to me,
> 
> > summary(mod)
> 
> Call:
> lm(formula = x ~ ind, na.action = na.fail)
> 
> Residuals:
> Min  1Q  Median  3Q Max 
> -1.2082 -0.5200  0.1309  0.5725  0.9593 
> 
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 59.845860.31900   187.6  < 2e-16 ***
> ind -0.167680.02541-6.6 2.57e-06 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> 
> Residual standard error: 0.705 on 19 degrees of freedom
> Multiple R-Squared: 0.6963, Adjusted R-squared: 0.6803 
> F-statistic: 43.56 on 1 and 19 DF,  p-value: 2.57
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] fitted probabilities in multinomial logistic regression are identical for each level

[John Fox]

Dear Bob,

If I understand correctly what you've done, the "newdata" that you're
using to get predicted values includes the three values of the response
variable, which are irrelevant to the predictions and cause each
prediction to be repeated three times.

I hope that this helps,
 John

On Tue, 27 Mar 2007 06:19:12 +1000
 Bob Green <[EMAIL PROTECTED]> wrote:
> 
> I was hoping for some advice regarding possible explanations for the 
> fitted probability values I obtained for a multinomial logistic 
> regression. The analysis aims to predict whether Capgras delusions 
> (present/absent) are associated with group (ABH, SV, homicide; values
> 
> = 1,2,3,), controlling for previous violence. What has me puzzled is 
> that for each combination the fitted probabilities are identical. I 
> haven't seen this in the worked examples I have come across and was 
> interested to know if this is a problem or what might be the cause 
> for this. I ran an analysis with another independent variable and 
> obtained a similar pattern.
> 
> Any assistance with this is appreciated
> 
> Bob Green
> 
>  > predictors <- expand.grid(group=1:3, in.acute.danger = c("y","n"),
> 
> violent.convictions = c("y","n"))
>  > p.fit <- predict(mod.multacute, predictors, type='probs')
>  > p.fit
> 1 2 3
> 1  0.4615070 0.3077061 0.2307869
> 2  0.4615070 0.3077061 0.2307869
> 3  0.4615070 0.3077061 0.2307869
> 4  0.7741997 0.1290310 0.0967693
> 5  0.7741997 0.1290310 0.0967693
> 6  0.7741997 0.1290310 0.0967693
> 7  0.4230927 0.3846055 0.1923017
> 8  0.4230927 0.3846055 0.1923017
> 9  0.4230927 0.3846055 0.1923017
> 10 0.7058783 0.1647063 0.1294153
> 11 0.7058783 0.1647063 0.1294153
> 12 0.7058783 0.1647063 0.1294153
> 
> 
>  > mod.multacute <- multinom(group ~ in.acute.danger * 
> violent.convictions, data = kc,  na.action = na.omit)
> # weights:  15 (8 variable)
> initial  value 170.284905
> iter  10 value 131.016050
> final  value 130.993722
> converged
>  > summary(mod.multacute, cor=F, Wald=T)
> Call:
> multinom(formula = group ~ in.acute.danger * violent.convictions,
>  data = kc, na.action = na.omit)
> 
> Coefficients:
>(Intercept) in.acute.dangery violent.convictionsy 
> in.acute.dangery:violent.convictionsy
> 2   -1.4552791.3599055   -0.3364982 
>   0.02651913
> 3   -1.6964160.9078901   -0.3830842 
>   0.47860722
> 
> Std. Errors:
>(Intercept) in.acute.dangery violent.convictionsy 
> in.acute.dangery:violent.convictionsy
> 2   0.29680820.52820770.6162498 
>0.9936493
> 3   0.32798380.63125690.6946869 
>1.1284891
> 
> Value/SE (Wald statistics):
>(Intercept) in.acute.dangery violent.convictionsy 
> in.acute.dangery:violent.convictionsy
> 2   -4.903094 2.574566   -0.5460419 
>   0.02668862
> 3   -5.172256 1.438226   -0.5514486 
>   0.42411327
> 
> Residual Deviance: 261.9874
> AIC: 277.9874
>  > Anova (mod.multacute)
> Anova Table (Type II tests)
> 
> Response: group
>  LR Chisq Df Pr(>Chisq)
> in.acute.danger  10.9335  2   0.004225 **
> violent.convictions   0.5957  2   0.742430
> in.acute.danger:violent.convictions   0.1895  2   0.909600
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Effect display of proportional odds model

[John Fox]

Dear Jan,

I've placed a copy of the published version of the paper at
<http://socserv.socsci.mcmaster.ca/jfox/Misc/polytomous-effect-displays/inde
x.html>, along with the R code for computing effects and their standard
errors and R code for the graphs in the paper. As you'll see, the functions
provided do not construct graphs automatically, as those in the effects
package do for linear and generalized linear models. Eventually, I'd like to
incorporate this material in the effects package.

Regards,
 John

----
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Jan Wijffels
> Sent: Friday, March 23, 2007 7:36 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Effect display of proportional odds model
> 
> Dear useRs,
> I very much like the effect display of the proportional odds 
> model on page 29 (Figure 8) of the following paper by John Fox:
> http://socserv.mcmaster.ca/jfox/Papers/logit-effect-displays.pdf
> It really gives a very concise overview of the model. I would 
> like to use it to illustrate the proportional odds mixed 
> models we fit here for a project on Diabetes but I can't seem 
> to reproduce the plot. Does anyone have code for the plot? 
> Maybe John Fox himself? I would appreciate it very much.
> Thanks,
> Jan
>  
> Jan Wijffels
> University Center for Statistics
> W. de Croylaan 54
> 3001 Heverlee
> Belgium
> tel: +32 (0)16 322784
> fax: +32 (0)16 322831
>  <http://www.kuleuven.be/ucs> http://www.kuleuven.be/ucs
>  
> 
> 
> Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Effect display of proportional odds model

[John Fox]

Dear Jan,

First, I inadvertently removed material on these displays from my web site
when the paper was published in Sociological Methodology 2006. I'll update
and repost the material some time in the next couple of days, including a
copy of the published paper (with a link on my home page). 

Second, the appendix to the paper and the originally posted examples didn't
include the code for Figure 8 (which is Figure 10 in the published version
of the paper). I'll add that.

Regards,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Jan Wijffels
> Sent: Friday, March 23, 2007 7:36 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Effect display of proportional odds model
> 
> Dear useRs,
> I very much like the effect display of the proportional odds 
> model on page 29 (Figure 8) of the following paper by John Fox:
> http://socserv.mcmaster.ca/jfox/Papers/logit-effect-displays.pdf
> It really gives a very concise overview of the model. I would 
> like to use it to illustrate the proportional odds mixed 
> models we fit here for a project on Diabetes but I can't seem 
> to reproduce the plot. Does anyone have code for the plot? 
> Maybe John Fox himself? I would appreciate it very much.
> Thanks,
> Jan
>  
> Jan Wijffels
> University Center for Statistics
> W. de Croylaan 54
> 3001 Heverlee
> Belgium
> tel: +32 (0)16 322784
> fax: +32 (0)16 322831
>  <http://www.kuleuven.be/ucs> http://www.kuleuven.be/ucs
>  
> 
> 
> Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to get "lsmeans"?

[John Fox]

Dear Frank,

The object is to focus on a high-order term of the model, holding other
terms "constant" at typical values; in the case of a factor, a "typical
value" is ambiguous, so an average is taken over the levels of the factor.
If the factor is, e.g., gender, one could produce an estimate of the average
fitted value for a population composed equally of men and women, or of a
population composed of men and women in proportion to their distribution in
the data. Otherwise, one would have to produce separate sets of fitted
values for men and women, with the number of such sets increasing as the
number of levels of the factors held constant increase. On the scale of the
linear predictor, these sets would differ only by constants.

Regards,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: Frank E Harrell Jr [mailto:[EMAIL PROTECTED] 
> Sent: Thursday, March 22, 2007 8:41 AM
> To: John Fox
> Cc: 'Prof Brian Ripley'; 'r-help'; 'Chuck Cleland'; 'Liaw, Andy'
> Subject: Re: [R] how to get "lsmeans"?
> 
> John Fox wrote:
> > Dear Brian et al.,
> > 
> > My apologies for chiming in late: It's been a busy day.
> > 
> > First some general comments on "least-squares means" and 
> "effect displays."
> ... much good stuff deleted ...
> 
> John - the one thing I didn't get from your post is a 
> motivation to do all that as opposed to easy-to-explain 
> predicted values.  I would appreciate your thoughts.
> 
> Thanks
> Frank
> 
> -- 
> Frank E Harrell Jr   Professor and Chair   School of Medicine
>   Department of Biostatistics   
> Vanderbilt University

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to get "lsmeans"?

[John Fox]

Dear Frank,

> -Original Message-
> From: Frank E Harrell Jr [mailto:[EMAIL PROTECTED] 
> Sent: Thursday, March 22, 2007 10:51 AM
> To: John Fox
> Cc: 'Prof Brian Ripley'; 'r-help'; 'Chuck Cleland'; 'Liaw, Andy'
> Subject: Re: [R] how to get "lsmeans"?
> 
> John Fox wrote:
> > Dear Frank,
> > 
> > The object is to focus on a high-order term of the model, holding 
> > other terms "constant" at typical values; in the case of a 
> factor, a 
> > "typical value" is ambiguous, so an average is taken over 
> the levels of the factor.
> > If the factor is, e.g., gender, one could produce an 
> estimate of the 
> > average fitted value for a population composed equally of men and 
> > women, or of a population composed of men and women in 
> proportion to 
> > their distribution in the data. Otherwise, one would have 
> to produce 
> > separate sets of fitted values for men and women, with the 
> number of 
> > such sets increasing as the number of levels of the factors held 
> > constant increase. On the scale of the linear predictor, 
> these sets would differ only by constants.
> > 
> > Regards,
> >  John
> 
> Makes sense.  I just set other factors to the mode, and if it 
> is important to see estimates for other categories, I give 
> estimates for each factor level.  If I want to uncondition on 
> a variable (not often) I take the proper weighted average of 
> predicted values.
> 

The mode seens arbitrary to me, but I don't think that there's a unique
right answer here. The weighted average (using sample proportions) is what
the effect() function does.

Regards,
 John

> Cheers
> Frank
> 
> > 
> > 
> > John Fox
> > Department of Sociology
> > McMaster University
> > Hamilton, Ontario
> > Canada L8S 4M4
> > 905-525-9140x23604
> > http://socserv.mcmaster.ca/jfox
> > --------
> > 
> >> -Original Message-
> >> From: Frank E Harrell Jr [mailto:[EMAIL PROTECTED]
> >> Sent: Thursday, March 22, 2007 8:41 AM
> >> To: John Fox
> >> Cc: 'Prof Brian Ripley'; 'r-help'; 'Chuck Cleland'; 'Liaw, Andy'
> >> Subject: Re: [R] how to get "lsmeans"?
> >>
> >> John Fox wrote:
> >>> Dear Brian et al.,
> >>>
> >>> My apologies for chiming in late: It's been a busy day.
> >>>
> >>> First some general comments on "least-squares means" and
> >> "effect displays."
> >> ... much good stuff deleted ...
> >>
> >> John - the one thing I didn't get from your post is a 
> motivation to 
> >> do all that as opposed to easy-to-explain predicted 
> values.  I would 
> >> appreciate your thoughts.
> >>
> >> Thanks
> >> Frank
> >>
> >> -- 
> >> Frank E Harrell Jr   Professor and Chair   School 
> of Medicine
> >>   Department of Biostatistics   
> >> Vanderbilt University
> > 
> > 
> 
> 
> -- 
> Frank E Harrell Jr   Professor and Chair   School of Medicine
>   Department of Biostatistics   
> Vanderbilt University

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to get "lsmeans"?

[John Fox]

Dear Bob,

I think that you make a valid point -- and I've included "Type-III" tests in
the Anova() function in the car package, for example, though not entirely
for consistency with SAS -- but my object in writing the effects package was
different. I wanted to provide a general approach to visualizing terms in
complex models with linear predictors. If I can with reasonable effort
provide a solution that's the same as "least-squares means" for models for
which "least-squares means" are defined then I'll do so at some point, but
duplicating what SAS does was not my goal.

Regards,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of 
> Muenchen, Robert A (Bob)
> Sent: Thursday, March 22, 2007 9:36 AM
> To: R-help@stat.math.ethz.ch
> Subject: Re: [R] how to get "lsmeans"?
> 
> Hi All,
> 
> Perhaps I'm stating the obvious, but to increase the use of R 
> in places where SAS & SPSS dominate, it's important to make 
> getting the same answers as easy as possible. That includes 
> things like lsmeans and type III sums of squares. I've read 
> lots of discussions here on sums of squares & I'm not 
> advocating type III use, just looking at it from a marketing 
> perspective. Too many people look for excuses to not change.
> The fewer excuses, the better.
> 
> Of course this is easy for me to say, as I'm not the one who 
> does the work! Much thanks to those who do.
> 
> Cheers,
> Bob
> 
> =
>   Bob Muenchen (pronounced Min'-chen), Manager
>   Statistical Consulting Center
>   U of TN Office of Information Technology
>   200 Stokely Management Center, Knoxville, TN 37996-0520
>   Voice: (865) 974-5230  
>   FAX:   (865) 974-4810
>   Email: [EMAIL PROTECTED]
>   Web:   http://oit.utk.edu/scc, 
>   News:  http://listserv.utk.edu/archives/statnews.html
> =
> 
> > -Original Message-
> > From: [EMAIL PROTECTED] [mailto:r-help- 
> > [EMAIL PROTECTED] On Behalf Of John Fox
> > Sent: Wednesday, March 21, 2007 8:59 PM
> > To: 'Prof Brian Ripley'
> > Cc: 'r-help'; 'Chuck Cleland'
> > Subject: Re: [R] how to get "lsmeans"?
> > 
> > Dear Brian et al.,
> > 
> > My apologies for chiming in late: It's been a busy day.
> > 
> > First some general comments on "least-squares means" and "effect 
> > displays."
> > The general idea behind the two is similar -- to examine 
> fitted values 
> > corresponding to a term in a model while holding other terms to
> typical
> > values -- but the implementation is not identical. There are also
> other
> > similar ideas floating around as well. My formulation is 
> more general 
> > in the sense that it applies to a wider variety of models, 
> both linear 
> > and otherwise.
> > 
> > "Least-squares means" (a horrible term, by the way: in a 
> 1980 paper in 
> > the American Statistician, Searle, Speed, and Milliken 
> suggested the 
> > more descriptive term "population marginal means") apply to factors 
> > and combinations of factors; covariates are set to mean 
> values and the 
> > levels of other factors are averaged over, in effect applying equal 
> > weight to each level. (This is from memory, so it's 
> possible that I'm 
> > not getting it quite right, but I believe that I am.) In my effect 
> > displays, each level of
> a
> > factor is weighted by its proportion in the data. In models 
> in which 
> > least-squares means can be computed, they should differ from the 
> > corresponding effect display by a constant (if there are different 
> > numbers of observations in the different levels of the factors that 
> > are held constant).
> > 
> > The obstacle to computing either least-squares means or effect
> displays
> > in R
> > via predict() is that predict() wants factors in the "new 
> data" to be 
> > set to particular levels. The effect() function in the 
> effects package 
> > bypasses
> > predict() and works directly with the model matrix, 
> averaging over the 
> > columns that pertain to a factor (and reconstructing 
> interactions as 
> > necessary). As mentioned, this has the effect of setting 
> 

Re: [R] how to get "lsmeans"?

[John Fox]

Dear Brian et al.,

My apologies for chiming in late: It's been a busy day.

First some general comments on "least-squares means" and "effect displays."
The general idea behind the two is similar -- to examine fitted values
corresponding to a term in a model while holding other terms to typical
values -- but the implementation is not identical. There are also other
similar ideas floating around as well. My formulation is more general in the
sense that it applies to a wider variety of models, both linear and
otherwise.

"Least-squares means" (a horrible term, by the way: in a 1980 paper in the
American Statistician, Searle, Speed, and Milliken suggested the more
descriptive term "population marginal means") apply to factors and
combinations of factors; covariates are set to mean values and the levels of
other factors are averaged over, in effect applying equal weight to each
level. (This is from memory, so it's possible that I'm not getting it quite
right, but I believe that I am.) In my effect displays, each level of a
factor is weighted by its proportion in the data. In models in which
least-squares means can be computed, they should differ from the
corresponding effect display by a constant (if there are different numbers
of observations in the different levels of the factors that are held
constant).

The obstacle to computing either least-squares means or effect displays in R
via predict() is that predict() wants factors in the "new data" to be set to
particular levels. The effect() function in the effects package bypasses
predict() and works directly with the model matrix, averaging over the
columns that pertain to a factor (and reconstructing interactions as
necessary). As mentioned, this has the effect of setting the factor to its
proportional distribution in the data. This approach also has the advantage
of being invariant with respect to the choice of contrasts for a factor.

The only convenient way that I can think of to implement least-squares means
in R would be to use deviation-coded regressors for a factor (that is,
contr.sum) and then to set the columns of the model matrix for the factor(s)
to be averaged over to 0. It may just be that I'm having a failure of
imagination and that there's a better way to proceed. I've not implemented
this solution because it is dependent upon the choice of contrasts and
because I don't see a general advantage to it, but since the issue has come
up several times now, maybe I should take a crack at it. Remember that I
want this to work more generally, not just for levels of factors, and not
just for linear models.

Brian is quite right in mentioning that he suggested some time ago that I
use critical values of t rather than of the standard normal distribution for
producing confidence intervals, and I agree that it makes sense to do so in
models in which the dispersion is estimated. My only excuse for not yet
doing this is that I want to undertake a more general revision of the
effects package, and haven't had time to do it. There are several changes
that I'd like to make to the package. For example, I have results for
multinomial and proportional odds logit models (described in a paper by me
and Bob Andersen in the 2006 issue of Sociological Methodology) that I want
to incorporate, and I'd like to improve the appearance of the default
graphs. But Brian's suggestion is very straightforward, and I guess that I
shouldn't wait to implement it; I'll do so very soon.

Regards,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Prof 
> Brian Ripley
> Sent: Wednesday, March 21, 2007 12:03 PM
> To: Chuck Cleland
> Cc: r-help
> Subject: Re: [R] how to get "lsmeans"?
> 
> On Wed, 21 Mar 2007, Chuck Cleland wrote:
> 
> > Liaw, Andy wrote:
> >> I verified the result from the following with output from JMP 6 on 
> >> the same data (don't have SAS: don't need it):
> >>
> >> set.seed(631)
> >> n <- 100
> >> dat <- data.frame(y=rnorm(n), A=factor(sample(1:2, n, 
> replace=TRUE)),
> >>   B=factor(sample(1:2, n, replace=TRUE)),
> >>   C=factor(sample(1:2, n, replace=TRUE)),
> >>   d=rnorm(n))
> >> fm <- lm(y ~ A + B + C + d, dat)
> >> ## Form a data frame of points to predict: all 
> combinations of the ## 
> >> three factors and the mean of the covariate.
> >> p <- data.frame(expand.grid(A=1:2, B=1:2, C=1:2)) p[] <- lapply(p, 
> >> factor) p <- cbind(p, d=mean

Re: [R] arrowhead styles

[John Fox]

Dear Hendrik,

Arrows() is the first entry turned up by RSiteSearch("arrows",
restrict="functions"), so its effort to hide isn't very successful.

Regards,
 John

----
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Hendrik Fuss
> Sent: Saturday, March 17, 2007 4:17 PM
> To: RHelp
> Subject: Re: [R] arrowhead styles
> 
> Mikkel Grum:
> >> I've been using the arrows() function in plots a lot, but I'm not 
> >> happy with the appearance of the arrow heads. I understand that 
> >> arrows
> >> () doesn't offer more sophisticated arrowhead shapes like e.g. a 
> >> filled triangle, possibly with choice of angle at the point. Does 
> >> anyone know an easy way to achieve this?
> > There's also an Arrows() function in the package IDPMisc.
> 
> The IDPMisc Arrows() are actually looking very good, and work 
> as drop- in replacement. Pity they're hiding in a library, in 
> which I wouldn't have looked for arrows.
> 
> many thanks
> Hendrik
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Rcmdir does not work in SciVews-R

[John Fox]

Dear Kihwang,

I believe that you'll have to use an older version of the Rcmdr package with
SciViews; you can download version 1.1-2 of the Rcmdr from the SciViews web
page <http://www.sciviews.org/SciViews-R/>. For more information, you might
contact the SciViews developers.

I hope this helps,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Kihwang Lee
> Sent: Tuesday, March 13, 2007 12:12 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Rcmdir does not work in SciVews-R
> 
> Hi,
> 
> I have recently installed version 2.4.1 of R and also 
> installed SciViews-R 0.8.9. It worked fine. Then I installed 
> the Rcmdr 1.2-7 from CRAN . However the 'R commander menu' in 
> the dockable panel does not work. If I start Rcmdr in R 
> without SciViews, it works perfectly.
> 
> Is there any way that I can use Rcmdr with SciViews?
> 
> Many thanks in advanve.
> 
> Kihwang
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] No fit statistics for some models using sem

[John Fox]

Dear David and Ista,

I haven't looked at this model carefully, but the fact that the df are
0 suggests that the model is just-identified and therefore necessarily
perfectly reproduces the covariances among the observed variables.
Removing a parameter would over-identify the model, making possible the
computation of the missing fit statistics.

Regards,
 John

On Wed, 7 Mar 2007 18:31:09 +
 "David Barron" <[EMAIL PROTECTED]> wrote:
> It's not the correlation as such that is the problem; it's because
> you
> only have 10 degrees of freedom "available" with four observed
> variables, and you are estimating 10 parameters, which is why you get
> a chi square of zero.  When you remove any one free parameter (such
> as
> the correlation), the model becomes identified.
> 
> On 07/03/07, Ista Zahn <[EMAIL PROTECTED]> wrote:
> > Hi,
> >
> > New to both R and SEM, so this may be a very simple question. I am
> > trying to run a very simple path analysis using the sem package.
> > There are 2 exogenous (FARSCH, LOCUS10) and 2 endogenous (T_ATTENT,
> > RMTEST) observed variables in the model.  The idea is that T_ATTENT
> > mediates the effect of FARSCH and LOCUS10 on RMTEST. The RAM
> > specification I used is
> >
> > FARSCH -> T_ATTENT, y1x1, NA
> > LOCUS10 -> T_ATTENT, y1x2, NA
> > FARSCH -> RMTEST10, y2x1, NA
> > LOCUS10 -> RMTEST10, y2x2, NA
> > T_ATTENT -> RMTEST10, y2y1, NA
> > FARSCH <-> FARSCH, x1x1, NA
> > LOCUS10 <-> LOCUS10, x2x2, NA
> > T_ATTENT <-> T_ATTENT, y1y1, NA
> > RMTEST10 <-> RMTEST10, y2y2, NA
> > LOCUS10 <-> FARSCH, x2x1, NA
> >
> > This model runs, but using the summary function does not return the
> > usual model fit statistics, only the following:
> >
> > Model Chisquare =  0   Df =  0 Pr(>Chisq) = NA
> >   Chisquare (null model) =  8526.8   Df =  6
> >   Goodness-of-fit index =  1
> >   BIC =  0
> >
> > If I omit the last line from the RAM specification(i.e., delete
> > "LOCUS10 <-> FARSCH, x2x1, NA"), I DO get all the usual statistics:
> >
> >   Model Chisquare =  1303.7   Df =  1 Pr(>Chisq) = 0
> >   Chisquare (null model) =  8526.8   Df =  6
> >   Goodness-of-fit index =  0.95864
> >   Adjusted goodness-of-fit index =  0.58639
> >   RMSEA index =  0.30029   90% CI: (NA, NA)
> >   Bentler-Bonnett NFI =  0.84711
> >   Tucker-Lewis NNFI =  0.082726
> >   Bentler CFI =  0.84712
> >   BIC =  1294.1
> >
> > My understanding is the you should always put in the correlation
> > between exogenous predictors, but when I do this I don't get fit
> > statistics. Can anyone help me understand what is happening here?
> >
> > Thank you,
> >
> > Ista
> >
> > __
> > R-help@stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
> 
> 
> -- 
> =
> David Barron
> Said Business School
> University of Oxford
> Park End Street
> Oxford OX1 1HP
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] sem: standardized covariance estimates

[John Fox]

Dear Mike,

I don't believe that I've provided a way to get these correlations
automatically, but you can, of course, just divide the covariance by the
product of the standard deviations.

I hope this helps,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Mike Allerhand
> Sent: Tuesday, March 06, 2007 11:14 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] sem: standardized covariance estimates
> 
> Dear all,
> 
> How do I get the standardized covariance (the correlation) 
> between two latent variables?
> 'standardized.coefficients' gives standardized path 
> coefficients, but not covariances.
> The covariance estimates are easily obtained from fit$coeff 
> or 'summary',  but EQS reports both the covariance and the 
> correlation, how can I get that?
> 
> best wishes,  Mike
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] SEM - standardized path coefficients?

[John Fox]

Dear Tim,

See ?standardized.coefficients (after loading the sem package).

Regards,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Tim Holland
> Sent: Wednesday, February 28, 2007 12:35 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] SEM - standardized path coefficients?
> 
> Hello -
> 
> Does anybody know how to get the SEM package in R to return 
> standardized path coefficients instead of unstandardized 
> ones?  Does this involve changing the covariance matrix, or 
> is there an argument in the SEM itself that can be changed?
> 
> Thank you,
> Tim
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] enhanced question / standardized coefficients

[John Fox]

Dear Simon,

In my opinion, standardized coefficients only offer the illusion of
comparison for quantitative explanatory variables, since there's no deep
reason that the standard deviation of one variable has the same meaning as
the standard deviation of another. Indeed, if the variables are in the same
units of measurement in the first place, permitting direct comparison of
unstandardized coefficients, then separate standardization of each X is like
using a rubber ruler.

That said, as you point out, it makes no sense to standardize the dummy
regressors for a factor, so you can just standardize the quantitative
variables (Y and X's) in the regression equation.

I hope that this helps,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Simon P. Kempf
> Sent: Wednesday, February 07, 2007 9:27 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] enhanced question / standardized coefficients
> 
> Hello,
> 
>  
> 
> I would like to repost the question of Joerg:
> 
>  
> 
> Hello everybody, 
> 
> a question that connect to the question of Frederik Karlsons 
> about 'how to stand. betas' 
> With the stand. betas i can compare the influence of the 
> different explaning variables. What do i with the betas of 
> factors? I can't use the solution of JohnFox, because there 
> is no sd of an factor. How can i compare the influence of the 
> factor with the influence of the numeric variables? 
> 
> I got the same problem. In my regression equation there are 
> several categorical variables and  I would like to compute 
> the standard coefficients. How can I do this?
> 
>  
> 
> Simon
> 
>  
> 
>  
> 
>  
> 
>  
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Reference to dataframe and contents

[John Fox]

Dear Rene,

There are several ways to do what you want, including get(s1)[,S2].

I hope this helps,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Rene Braeckman
> Sent: Sunday, February 04, 2007 3:43 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Reference to dataframe and contents
> 
> This is probably easy for experienced users but I could not 
> find a solution.
>  
> I have several R scripts that process several columns of a 
> dataframe (several dataframes and columns actually, but 
> simplified for my question).
> 
> References such as:
>  
> myDF$myCol
>  
> are all over. I like to automate this for other dataframes 
> and columns by defining a reference only once in the 
> beginning of the script.
> 
> One idea is to use strings:
> 
> s1 <- "myDF"
> S2 <- "myCol"
> 
> My question is how to construct the equivalent of myDF$myCol 
> that can be used as such. Or is there a better solution?
> 
> Thanks. The help and discussions on this forum are the best.
> Rene
> -
> Rene Braeckman, PhD
> Clinical Pharmacology & Pharmacometrics
> Irvine, California, USA
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] efficient code. how to reduce running time?

[John Fox]

Dear Brian,


> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Prof 
> Brian Ripley
> Sent: Monday, January 22, 2007 11:06 AM
> To: Charilaos Skiadas
> Cc: John Fox; r-help@stat.math.ethz.ch
> Subject: Re: [R] efficient code. how to reduce running time?
> 
> On Mon, 22 Jan 2007, Charilaos Skiadas wrote:
> 
> > On Jan 21, 2007, at 8:11 PM, John Fox wrote:
> >
> >> Dear Haris,
> >>
> >> Using lapply() et al. may produce cleaner code, but it won't 
> >> necessarily speed up a computation. For example:
> >>
> >>> X <- data.frame(matrix(rnorm(1000*1000), 1000, 1000)) y <- 
> >>> rnorm(1000)
> >>>
> >>> mods <- as.list(1:1000)
> >>> system.time(for (i in 1:1000) mods[[i]] <- lm(y ~ X[,i]))
> >> [1] 40.53  0.05 40.61NANA
> >>>
> >>> system.time(mods <- lapply(as.list(X), function(x) lm(y ~ x)))
> >> [1] 53.29  0.37 53.94NANA
> >>
> > Interesting, in my system the results are quite different:
> >
> > > system.time(for (i in 1:1000) mods[[i]] <- lm(y ~ X[,i]))
> > [1] 192.035  12.601 797.094   0.000   0.000
> > > system.time(mods <- lapply(as.list(X), function(x) lm(y ~ x)))
> > [1]  59.913   9.918 289.030   0.000   0.000
> >
> > Regular MacOSX install with ~760MB memory.
> 
> But MacOS X is infamous for having rather specific speed 
> problems with its malloc, and so gives different timing 
> results from all other platforms.
> We are promised a solution in MacOS 10.5.
> 

Thanks for the clarification.

> Both of your machines seem very slow compared to mine:
> 
> > system.time(for (i in 1:1000) mods[[i]] <- lm(y ~ X[,i]))
> user  system elapsed
>   11.011   0.250  11.311
> > system.time(mods <- lapply(as.list(X), function(x) lm(y ~ x)))
> user  system elapsed
>   13.463   0.260  13.812
> 
> and that on a 64-bit platform (AMD64 Linux FC5).
> 

As you can see from the specs (in a previous message), my system is quite
old, which probably accounts for at least part of the difference. The ratios
of the user times for my and your system aren't too different though:

> 53.29/40.53  # mine
[1] 1.314829

> 13.463/11.011  # yours
[1] 1.222686

Regards,
 John

> -- 
> Brian D. Ripley,  [EMAIL PROTECTED]
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford, Tel:  +44 1865 272861 (self)
> 1 South Parks Road, +44 1865 272866 (PA)
> Oxford OX1 3TG, UKFax:  +44 1865 272595
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] efficient code. how to reduce running time?

[John Fox]

Dear Haris,

My timings were on a 3 GHz Pentium 4 system with 1 GB of memory running Win
XP SP2 and R 2.4.1.

I'm no expert on these matters, and I wouldn't have been surprised by
qualitatively different results on different systems, but this difference is
larger than I would have expected. One thing that seems particularly
striking in your results is the large difference between elapsed time and
user CPU time, making me wonder what else was going on when you ran these
examples.

Regards,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: Charilaos Skiadas [mailto:[EMAIL PROTECTED] 
> Sent: Monday, January 22, 2007 10:00 AM
> To: John Fox
> Cc: r-help@stat.math.ethz.ch; 'miraceti'
> Subject: Re: [R] efficient code. how to reduce running time?
> 
> On Jan 21, 2007, at 8:11 PM, John Fox wrote:
> 
> > Dear Haris,
> >
> > Using lapply() et al. may produce cleaner code, but it won't 
> > necessarily speed up a computation. For example:
> >
> >> X <- data.frame(matrix(rnorm(1000*1000), 1000, 1000)) y <- 
> >> rnorm(1000)
> >>
> >> mods <- as.list(1:1000)
> >> system.time(for (i in 1:1000) mods[[i]] <- lm(y ~ X[,i]))
> > [1] 40.53  0.05 40.61NANA
> >>
> >> system.time(mods <- lapply(as.list(X), function(x) lm(y ~ x)))
> > [1] 53.29  0.37 53.94NANA
> >
> Interesting, in my system the results are quite different:
> 
>  > system.time(for (i in 1:1000) mods[[i]] <- lm(y ~ X[,i]))
> [1] 192.035  12.601 797.094   0.000   0.000
>  > system.time(mods <- lapply(as.list(X), function(x) lm(y ~ x)))
> [1]  59.913   9.918 289.030   0.000   0.000
> 
> Regular MacOSX install with ~760MB memory.
> 
> > In cases such as this, I don't even find the code using *apply() 
> > easier to read.
> >
> > Regards,
> >  John
> 
> Haris
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] efficient code. how to reduce running time?

[John Fox]

Dear Haris,

Using lapply() et al. may produce cleaner code, but it won't necessarily
speed up a computation. For example:

> X <- data.frame(matrix(rnorm(1000*1000), 1000, 1000))
> y <- rnorm(1000)
> 
> mods <- as.list(1:1000)
> system.time(for (i in 1:1000) mods[[i]] <- lm(y ~ X[,i]))
[1] 40.53  0.05 40.61NANA
> 
> system.time(mods <- lapply(as.list(X), function(x) lm(y ~ x)))
[1] 53.29  0.37 53.94NANA

In cases such as this, I don't even find the code using *apply() easier to
read.

Regards,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of 
> Charilaos Skiadas
> Sent: Sunday, January 21, 2007 7:01 PM
> To: miraceti
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] efficient code. how to reduce running time?
> 
> On Jan 21, 2007, at 5:55 PM, miraceti wrote:
> 
> > Thank you all for lookin at it.
> > I'll fix the code to preallocate the objects.
> > and I wonder if there is a way to call anova on all the 
> columns at the 
> > same time..
> > Right now I am calling (Y~V1, data) from V1 to V50 thru a loop.
> > I tried (Y~., data) but it gave me different values from 
> the results I 
> > get when I call them separately, So I can't help but call 
> them 25,000 
> > times...
> 
> Have you looked at lapply, sapply, apply and friends?
> 
> Haris
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] efficient code. how to reduce running time?

[John Fox]

Dear Mira,

I didn't work through your code in detail, but I did notice that you're
doing something that's very inefficient in R -- building up objects
element-by-element, e.g., by indexing beyond their current length. Instead,
you can preallocate the objects and simply replace elements. For example,
create the vector F in your perm.F() as F <- numeric(nperms) rather than as
an empty vector. (BTW, I'd probably not name a variable "F" since this is
usually a synonym for the logical value FALSE.) There's a similar problem in
your handling of maxF, which you build up column-by-column via cbind().
(BTW, is maxF really a matrix?) You also needlessly recompute max(F), never
appear to use MSSid, and end lines with unnecessary semicolons.

I hope that this helps,
 John

----
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of miraceti
> Sent: Sunday, January 21, 2007 12:38 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] efficient code. how to reduce running time?
> 
> Hi,
> I am new to R.
> and even though I've made my code to run and do what it needs to .
> It is taking forever and I can't use it like this.
> I was wondering if you could help me find ways to fix the 
> code to run faster.
> Here are my codes..
> the data set is a bunch of 0s and 1s in a data.frame.
> What I am doing is this.
> I pick a column and make up a new column Y with values 
> associated with that column I picked.
> then I remove the column.
> and see if any other columns in the data have a significant 
> association with the Y column I've generated.
> If you see anything obvious that takes a long time, any 
> suggestions would be appreciated.
> thanks a lot.
> 
> Mira
> 
> --
> 
> #sub function for finding index
> rfind <- function(x)seq(along=x)[x*(1-x)>MAF*(1-MAF)]
> 
> #sub function for permutation test
> perm.F = function(y,x,nperms,sites)
> {
>   maxF = c();
>   for (i in 1:nperms)
>   {
> F=numeric(S)   #create an empty vector to store the F-values
> newY=sample(y,length(y)) #permute the cancer types
> newX = cbind(x, newY);
> # anova for all sites
> for ( i in sites )
> {
>   a <- anova(lm(newY~factor(newX[,i])));
>   F[i] <- a$`F value`[1];
> }
> MSSid <- which (F == max(F)); # index of MSS (Most 
> Significant Site)
> maxF = cbind(maxF,max(F));
>   }
>   maxF;
> }
> 
> 
> # set the output file
> sink("/tmp/R.out.3932.100")
> # load the dataset
> snp = read.table(file("/tmp/msoutput.3932.100"))
> #print (snp);
> 
> # pi: desired proportion of variation due to QTN pi = 0.05; 
> print (paste("pi:", pi)); MAF = 0.05; print (paste("MAF:", 
> MAF)); # S: number of segregating sites S = length(snp[1,]); 
> # N: number of samples N = length(snp[,1]); Dips = 
> sample(1:N,N) DipA = Dips[1:50] DipB = Dips[51:100] disnp = 
> snp[DipA,]+snp[DipB,] snp = as.data.frame(disnp, 
> row.names=NULL); N = length(snp[,1]);
> 
> # get allele freq for all SNPs
> allele_f <- mean(snp[,1:S])/2;
> print (allele_f);
> sites = rfind(allele_f);
> print(sites);
> 
> # collapse sites that have perfect correlation newsites <- 
> sites; for (i in 1:(length(sites)-1)) {
>   for (j in (i+1):length(sites))
>   {
> test = (snp[sites[i]] == snp[sites[j]])
> if ( all(test) || all(!test) )
> {
>   print (paste("perfect correlation with", sites[i]));
>   print (paste("removing alleles", sites[j]));
>   newsites <- newsites[newsites!=sites[j]];
> }
>   }
> }
> sites <- newsites;
> print(sites);
> 
> # QTN: the site nearest right to S/4
> sitesid = floor(length(sites)/4);
> QTNid = sites[sitesid];
> QTN = snp[,QTNid];
> 
> print (paste("QTN:", names(snp[QTNid]))); print (QTN);
> 
> # remove QTN from sites
> sites <- sites [ sites != QTNid ];
> print(sites);
> print (paste("Number of usable SNPs:", length(sites)));
> 
> # p: allele frequency of QTN
> p0 = allele_f[QTNid];
> p = min(p0, 1-p0);
> print (paste("QTN_freq:", p));
> 
> # z: random normal deviate
> z = rnorm(N, mean = 0, sd = 1);
> # foreach sample give quantitative phenotype # each row is a 
> sample # phenotype value depends on QTN genotype, pi, p, and 
> z Y <- sqrt(10-(10*pi))*z + QTN*sqrt((10

Re: [R] extend summary.lm for hccm?

[John Fox]

Dear Achim,

> -Original Message-
> From: Achim Zeileis [mailto:[EMAIL PROTECTED] 
> Sent: Monday, December 25, 2006 8:38 AM
> To: John Fox
> Cc: 'Dirk Eddelbuettel'; 'ivo welch'; r-help@stat.math.ethz.ch
> Subject: Re: [R] extend summary.lm for hccm?
> 
> On Sun, 24 Dec 2006, John Fox wrote:
> 
> > Dear Dirk and Ivo,
> >
> > It's true that the sandwich package has more extensive 
> facilities for 
> > sandwich variance estimation than the hccm function in the car 
> > package, but I think that the thrust of Ivo's question is to get a 
> > model summary that automatically uses the adjusted standard errors 
> > rather than having to compute them and use them "manually."
> 
> I've written the function coeftest() in package "lmtest" for 
> this purpose.
> With this you can
>   coeftest(obj, vcov = sandwich)

Since coeftest() already exists in a package, I think that this is a
preferable solution.

> or you can put this into a specialized summary() function as 
> John suggested (where you probably would want the F statistic 
> to use the other vcov as well). 

Good point. Here's a modified version that also fixes up the F-test:

summaryHCCM.lm <- function(model, type=c("hc3", "hc0", "hc1", "hc2", "hc4"),

...){
if (!require(car)) stop("Required car package is missing.")
type <- match.arg(type)
V <- hccm(model, type=type)
sumry <- summary(model)
table <- coef(sumry)
table[,2] <- sqrt(diag(V))
table[,3] <- table[,1]/table[,2]
table[,4] <- 2*pt(abs(table[,3]), df.residual(model), lower.tail=FALSE)
sumry$coefficients <- table
p <- nrow(table)
hyp <- if (has.intercept(model)) cbind(0, diag(p - 1)) else diag(p)
sumry$fstatistic[1] <- linear.hypothesis(model, hyp,
white.adjust=type)[2,"F"]
print(sumry)
cat("Note: Heteroscedasticity-consistant standard errors using
adjustment",
type, "\n")
}

Regards,
 John

> See also function waldtest() 
> in "lmtest",
> linear.hypothesis() in "car" and
>   vignette("sandwich", package = "sandwich")
> 
> Although this works, it is still a nuisance to use a 
> different function and not summary() directly. In addition, 
> it would also be nice to plug in different vcovs into 
> confint() or predict() methods. Of course, one could write 
> different generics or overload the methods in certain packages.
> However, I guess that many practitioners want to use 
> different vcov estimators - especially in the social and 
> political scieneces, and econometrics etc. - so that this 
> might justify that the original "lm" (and
> "glm") methods are extended to allow for plugging in 
> different vcov matrices. Maybe we could try to convince 
> R-core to include somthing like this?
> 
> Best wishes,
> Z
> 
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] extend summary.lm for hccm?

[John Fox]

Dear Frank,

If I remember Freedman's recent paper correctly, he argues that sandwich
variance estimator, though problematic in general, is not problematic in the
case that White described -- an otherwise correctly specified linear model
with heteroscedasticity estimated by least-squares. 

Regards,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Frank 
> E Harrell Jr
> Sent: Sunday, December 24, 2006 10:29 AM
> To: John Fox
> Cc: r-help@stat.math.ethz.ch; 'ivo welch'; 'Dirk Eddelbuettel'
> Subject: Re: [R] extend summary.lm for hccm?
> 
> John Fox wrote:
> > Dear Dirk and Ivo,
> > 
> > It's true that the sandwich package has more extensive 
> facilities for 
> > sandwich variance estimation than the hccm function in the car 
> > package, but I think that the thrust of Ivo's question is to get a 
> > model summary that automatically uses the adjusted standard errors 
> > rather than having to compute them and use them "manually."  Here's 
> > one approach, which could be modified slightly if Ivo wants 
> "hc0" as 
> > the default; it could also be adapted to use the sandwich package.
> > 
> > I hope this helps,
> >  John
> 
> Another approach:
> 
> library(Design)  # also requires library(Hmisc) f <- ols(y ~ 
> x1 + x2, x=TRUE, y=TRUE)
> f <- robcov(f)   # sandwich; also allows clustering.  Also see bootcov
> anova(f) # all later steps use sandwich variance matrix
> summmary(f)
> contrast(f, list(x1=.5), list(x1=.2))
> 
> BUT note that sandwich covariance matrix estimators can have 
> poor mean squared error (a paper by Bill Gould in Stata 
> Technical Bulletin years ago related to logistic regression 
> showed an example with a 100-fold increase in the variance of 
> a variance estimate) and can give you the right estimate of 
> the wrong quantity (reference below).
> 
> Frank Harrell
> 
> @Article{free06so,
>author =   {Freedman, David A.},
>title ={On the so-called ``{Huber} sandwich 
> estimator'' and
> ``robust standard errors''},
>journal =  The American Statistician,
>year = 2006,
>volume =   60,
>pages ={299-302},
>annote =   {nice summary of derivation of sandwich
> estimators;questions why we should be interested in getting 
> the right variance of the wrong parameters when the model 
> doesn't fit} }
> 
> 
> > 
> > --- snip --
> > 
> > summaryHCCM <- function(model, ...) UseMethod("summaryHCCM")
> > 
> > summaryHCCM.lm <- function(model, type=c("hc3", "hc0", 
> "hc1", "hc2", 
> > "hc4"),
> > 
> > ...){
> > if (!require(car)) stop("Required car package is missing.")
> > type <- match.arg(type)
> > V <- hccm(model, type=type)
> > sumry <- summary(model)
> > table <- coef(sumry)
> > table[,2] <- sqrt(diag(V))
> > table[,3] <- table[,1]/table[,2]
> > table[,4] <- 2*pt(abs(table[,3]), df.residual(model), 
> lower.tail=FALSE)
> > sumry$coefficients <- table
> > print(sumry)
> > cat("Note: Heteroscedasticity-consistant standard errors using 
> > adjustment",
> > type, "\n")
> > }
> > 
> > 
> > John Fox
> > Department of Sociology
> > McMaster University
> > Hamilton, Ontario
> > Canada L8S 4M4
> > 905-525-9140x23604
> > http://socserv.mcmaster.ca/jfox
> > 
> > 
> >> -Original Message-
> >> From: [EMAIL PROTECTED] 
> >> [mailto:[EMAIL PROTECTED] On Behalf Of Dirk 
> >> Eddelbuettel
> >> Sent: Saturday, December 23, 2006 11:16 PM
> >> To: ivo welch
> >> Cc: r-help@stat.math.ethz.ch
> >> Subject: Re: [R] extend summary.lm for hccm?
> >>
> >>
> >> On 23 December 2006 at 20:46, ivo welch wrote:
> >> | I wonder whether it is possible to extend the summary
> >> method for the
> >> | lm function, so that it uses an option "hccm" (well, model
> >> "hc0").  In
> >> | my line of work, it is p

Re: [R] extend summary.lm for hccm?

[John Fox]

Dear Dirk and Ivo,

It's true that the sandwich package has more extensive facilities for
sandwich variance estimation than the hccm function in the car package, but
I think that the thrust of Ivo's question is to get a model summary that
automatically uses the adjusted standard errors rather than having to
compute them and use them "manually."  Here's one approach, which could be
modified slightly if Ivo wants "hc0" as the default; it could also be
adapted to use the sandwich package.

I hope this helps,
 John

--- snip --

summaryHCCM <- function(model, ...) UseMethod("summaryHCCM")

summaryHCCM.lm <- function(model, type=c("hc3", "hc0", "hc1", "hc2", "hc4"),

...){
if (!require(car)) stop("Required car package is missing.")
type <- match.arg(type)
V <- hccm(model, type=type)
sumry <- summary(model)
table <- coef(sumry)
table[,2] <- sqrt(diag(V))
table[,3] <- table[,1]/table[,2]
table[,4] <- 2*pt(abs(table[,3]), df.residual(model), lower.tail=FALSE)
sumry$coefficients <- table
print(sumry)
cat("Note: Heteroscedasticity-consistant standard errors using
adjustment",
type, "\n")
}


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Dirk 
> Eddelbuettel
> Sent: Saturday, December 23, 2006 11:16 PM
> To: ivo welch
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] extend summary.lm for hccm?
> 
> 
> On 23 December 2006 at 20:46, ivo welch wrote:
> | I wonder whether it is possible to extend the summary 
> method for the 
> | lm function, so that it uses an option "hccm" (well, model 
> "hc0").  In 
> | my line of work, it is pretty much required in reporting of 
> almost all 
> | linear regressions these days, which means that it would be 
> very nice 
> | not to have to manually library car, then sqrt the diagonal, and 
> | recompute T-stats; instead, I would love to get everything 
> in the same 
> | format as the current output---except errors heteroskedasticity 
> | adjusted.
> | 
> | easy or hard?
> 
> Did you consider the 'sandwich' package? A simple
> 
>   > install.packages("sandwich")
>   > library(sandwich)
>   > ?vcovHC
>   > ?vcovHAC
>   
> should get you there.
> 
> Dirk
> 
> --
> Hell, there are no rules here - we're trying to accomplish something. 
>   -- Thomas A. Edison
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] effect plot

[John Fox]

Dear Ambrogi,


> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of 
> Ambrogi Federico
> Sent: Tuesday, December 19, 2006 12:03 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] effect plot
> 
> Dear R users,
> 
> Is there a simple way to use the effect function 
> (library(effects)) with a GEE estimated model?

Not that I'm aware. 

Sorry,
 John

> 
>  
> 
> library(gee)
> 
> library(MASS)
> 
> library(effects)
> 
> attach(epil)
> 
>  
> 
> b = gee(y ~ lbase*trt + lage + V4, family=poisson, id=subject,
> corstr="exchangeable")
> 
>  
> 
> plot(effect("lbase*trt", b))
> 
> # Errore in effect("lbase*trt", b) : nessun metodo 
> applicabile per "effect"
> 
>  
> 
>  
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] unique columns of a matrix

[John Fox]

Dear Roman,

You can use unique(X, MARGIN=2). See ?unique for details.

I hope this helps,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Roman 
> Akhter Ahmed
> Sent: Monday, December 18, 2006 8:42 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] unique columns of a matrix
> 
> Dear all,
> 
> I have a matrix of repeating columns in R, for example a matrix X is
> 
>   [,1]   [,2]   [,3]   [,4]
> [1,]   1  1  1  1
> [2,]   1  1  2  2
> 
> I want to store unique columns of the matrix X in a new matrix Y. 
> Therefore, Y will be
> 
>   [,1]   [,2]  
> [1,]   1  1
> [2,]   1  2
> 
> It will be really appreciated if you can provide me some 
> function for this job.
> Thanks for your time and effort in advance, Roman
> 
> --
> --
> Roman Akhter Ahmed (Ph.D. Candidate)
> Department of Econometrics and Business Statistics Room 659, 
> Building 11 (East Wing), Clayton Campus Monash University, 
> Victoria 3800, Australia
> Ph.: +61 3 9905 8346 (W), +61 3 9543 1958 (R)
> Web: http://www.buseco.monash.edu.au/staff/profile.php?uid=rahmed
> --
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] How to load Rcmd without Commander() ?

[John Fox]

Dear Ronggui,

The Commander() function is called in the Rcmdr .onAttach() function, so
library(Rcmdr) automatically starts up the GUI. You could use
Rcmdr:::reliability(S), where S is the covariance matrix among the items,
without libraryI(Rcmdr). Perhaps somewhat better would be to make a copy of
the function in the global environment via reliability <-
Rcmdr:::reliability.

I hope this helps,
 John

----
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of ronggui
> Sent: Friday, December 15, 2006 2:54 AM
> To: R-Help
> Subject: [R] How to load Rcmd without Commander() ?
> 
> I would like to use the reliability() function in Rcmdr 
> package, but not the GUI, so I would to load Rcmd without 
> Commander() running automatically.
> 
> Is there any easy way to get it? Thanks.
> 
> 
> --
> Ronggui Huang
> Department of Sociology
> Fudan University, Shanghai, China
> 黄荣贵
> 复旦大学社会学系
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Small Rcmdr issue

[John Fox]

Dear Alan,

I've noticed that problem too, and have no idea why it happens. The function
that's executed when you select "Exit -> From Commander and R" is pretty
simple:

closeCommanderAndR <- function(){
response <- closeCommander()
if (response == "cancel") return()
cat("\n")
quit(save="no")
}

If anyone knows how I can fix the problem, I'd be happy to so do.

Regards,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Alan Jackson
> Sent: Sunday, December 10, 2006 11:21 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Small Rcmdr issue
> 
> Just installed Rcmdr on a Linux box.
> 
> When I exit from both Rcmdr and R together from the regular 
> Rcmdr exit menu, it leaves my xterm with "stty -echo", so 
> that nothing that is typed appears.
> If I exit only Rcmdr, and then exit R normally, everything is 
> okay (stty echo).
> 
> Other than that minor irritant, everything seems to work just fine.
> 
> Linux 2.6.16-gentoo-r3
> Athlon 64 X2 3800+
> R version 2.4.0 (2006-10-03)
> Rcmdr Version 1.2-5
> 
> --
> --
> -
> | Alan K. Jackson| To see a World in a Grain of 
> Sand  |
> | [EMAIL PROTECTED]  | And a Heaven in a Wild Flower, 
> |
> | www.ajackson.org   | Hold Infinity in the palm of 
> your hand |
> | Houston, Texas | And Eternity in an hour. - 
> Blake   |
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Multilevel Modeling in R

[John Fox]

Dear Jeff,

On Mon, 4 Dec 2006 13:11:32 -0500
 "Jeff Miller" <[EMAIL PROTECTED]> wrote:
> Wow, would someone please send pdf links like that for SEM?
>

R's SEM capabilities aren't as strong as its mixed-model capabilities,
but here are some links:

<http://socserv.socsci.mcmaster.ca/jfox/Misc/sem/SEM-paper.pdf> 

<http://socserv.socsci.mcmaster.ca/jfox/Books/Companion/appendix-sems.pdf>

<http://socserv.socsci.mcmaster.ca/jfox/Courses/Oxford-2006/index.html>

I hope this helps,
 John

> Thanks,
> Jeff
> 
> -Original Message-
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] On Behalf Of Chuck Cleland
> Sent: Monday, December 04, 2006 1:01 PM
> To: Matthew Bridgman
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] Multilevel Modeling in R
> 
> Matthew Bridgman wrote:
> > Can anyone recommend a good text or resource for learning how to do
>  
> > Multilevel modeling in R?
> 
>   Here are a few other resources in addition to Pinheiro & Bates
> (2000):
> 
> http://finzi.psych.upenn.edu/R/library/mlmRev/doc/MlmSoftRev.pdf
> 
> http://cran.r-project.org/doc/packages/multilevel.pdf
> 
> http://stat.ethz.ch/CRAN/doc/Rnews/Rnews_2003-3.pdf [Lockwood, Doran
> &
> McCaffrey]
> 
>
http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-mixed-models.pd
> f
> 
> > Thanks,
> >Matt
> > 
> > __
> > R-help@stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> 
> -- 
> Chuck Cleland, Ph.D.
> NDRI, Inc.
> 71 West 23rd Street, 8th floor
> New York, NY 10010
> tel: (212) 845-4495 (Tu, Th)
> tel: (732) 512-0171 (M, W, F)
> fax: (917) 438-0894
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Box Tidwell / Error Message / Error in parse(file, n, text, prompt) : syntax error in

[John Fox]

Dear Simon,

I think that the most likely answer is that there's an error in the command
that you wrote (which is being sourced from a file?), but it's hard to tell
from the information that you supply. It might help to do a traceback()
after the error. I doubt that the error depends on the data.

I hope this helps,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Simon P. Kempf
> Sent: Monday, December 04, 2006 2:53 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Box Tidwell / Error Message / Error in 
> parse(file, n, text,prompt) : syntax error in
> 
> Dear R-Users,
> 
>  
> 
> I used the box.tidwell () function of the car Package. 
> 
>  
> 
> So far everything is fine. However, if the number of dummy 
> variables in the part not to be transformed (other.x formula) 
> exceeds a certain level (around 70), I receive the following 
> error message:
> 
>  
> 
> Error in parse(file, n, text, prompt) : syntax error in 
> 
>  
> 
> What did I miss? And how can I solve this problem?
> 
>  
> 
> I read some of the messages in the mail archives about this 
> error message, but to be honest I am rather a novice in R, I 
> did not understand them.
> 
>  
> 
> Some background information:
> 
>  
> 
> -  The data set does not have any missing values
> 
> -  The data set contains more than 19000 observations.
> 
>  
> 
> Thanks in advance and if you need more information about the 
> data, please let me know,
> 
>  
> 
> Simon
> 
>  
> 
>  
> 
>  
> 
>  
> 
> Simon P. Kempf 
> 
> Dipl.-Kfm. MScRE Immobilienvkonom (ebs)
> 
> Wissenschaftlicher Assistent
> 
>  
> 
> B|ro:
> 
> IREBS Immobilienakademie
> 
> c/o ebs Immobilienakademie GmbH
> 
> Berliner Str. 26a
> 
> 13507 Berlin
> 
>  
> 
> Privat:
> 
> Dunckerstra_e 60
> 
> 10439 Berlin
> 
>  
> 
> Mobil: 0176 7002 6687
> 
> Email:  <mailto:[EMAIL PROTECTED]> [EMAIL PROTECTED]
> 
>  
> 
> 
>   [[alternative HTML version deleted]]
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] lmer and a response that is a proportion

[John Fox]

Dear Cameron,

Given your description, I thought that this might be the case. 

I'd first examine the distribution of the response variable to see what it
looks like. If the values don't push the boundaries of 0 and 1, and their
distribution is unimodal and reasonably symmetric, I'd consider analyzing
them directly using normally distributed errors. If the values do stack up
near 0, 1, or both, I'd consider a transformation, or perhaps a different
family (depending on the pattern); in particular, if they stack up near both
0 and 1, a logit or similar transformation could help. Finally, if you have
many values of 0, 1, or both, then a transformation isn't promising (and,
indeed, the logit wouldn't be defined for these values). In any event, I'd
check diagnostics after a preliminary fit.

I hope this helps,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: Cameron Gillies [mailto:[EMAIL PROTECTED] 
> Sent: Sunday, December 03, 2006 6:31 PM
> To: Prof Brian Ripley; John Fox
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] lmer and a response that is a proportion
> 
> Dear Brian and John,
> 
> Thanks for your insight.  I'll clarify a couple of things 
> incase it changes your advice.
> 
> My response is a ratio of two measures taken during a bird's 
> path, which varies from 0  to 1, so I cannot convert it 
> columns of the number of successes.  It has to be reported as 
> the proportion.  I could logit transform it to make it 
> normal, but I am trying to avoid that so I can analyze it directly.
> 
> The subjects are individual birds and I have a range of 
> sample sizes from each bird (from 8 to >200, average of about 
> 75 measurements/bird).
> 
> Thanks!
> Cam
> 
> 
> On 12/3/06 3:47 PM, "Prof Brian Ripley" <[EMAIL PROTECTED]> wrote:
> 
> > On Sun, 3 Dec 2006, John Fox wrote:
> > 
> >> Dear Cameron,
> >> 
> >>> -Original Message-
> >>> From: [EMAIL PROTECTED] 
> >>> [mailto:[EMAIL PROTECTED] On Behalf Of Cameron 
> >>> Gillies
> >>> Sent: Sunday, December 03, 2006 1:58 PM
> >>> To: r-help@stat.math.ethz.ch
> >>> Subject: [R] lmer and a response that is a proportion
> >>> 
> >>> Greetings all,
> >>> 
> >>> I am using lmer (lme4 package) to analyze data where the 
> response is 
> >>> a proportion (0 to 1).  It appears to work, but I am wondering if 
> >>> the analysis is treating the response appropriately - 
> i.e. can lmer 
> >>> do this?
> >>> 
> >> 
> >> As far as I know, you can specify the response as a proportion, in 
> >> which case the binomial counts would be given via the weights 
> >> argument -- at least that's how it's done in glm(). An alternative 
> >> that should be equivalent is to specify a two-column matrix with 
> >> counts of "successes" and "failures" as the response. 
> Simply giving 
> >> the proportion of successes without the counts wouldn't be 
> appropriate.
> >> 
> >>> I have used both family=binomial and quasibinomial - is one more 
> >>> appropriate when the response is a proportion?  The coefficient 
> >>> estimates are identical, but the standard errors are larger with 
> >>> family=binomial.
> >>> 
> >> 
> >> The difference is that in the binomial family the 
> dispersion is fixed 
> >> to 1, while in the quasibinomial family it is estimated as a free 
> >> parameter. If the standard errors are larger with family=binomial, 
> >> then that suggests that the data are underdispersed 
> (relative to the 
> >> binomial); if the difference is substantial -- the factor 
> is just the 
> >> square root of the estimated dispersion -- then the 
> binomial model is 
> >> probably not appropriate for the data.
> > 
> > John's last deduction is appropriate to a GLM, but not 
> necessarily to 
> > a GLMM. I don't have detailed experience with lmer for 
> binomial, but I 
> > do for various other fitting routines for GLMM.  Remember 
> there are at 
> > least two sources of randomness in a GLMM, and let us keep 
> it simple 
> > and have just a subject effect and a measurement error.  Then if 
> > over-dispersion is happening within subjects, forcing the binomial 
> > dispersion (at the measurement level) to 1 tends to increase the 
> > estimate of the subject-level variance component to 
> compensate, and in 
> > turn increase some of the standard errors.
> > 
> > (Please note the 'tends' in that para, as the details of 
> the design do 
> > matter.  For cognescenti, think about plot and sub-plot 
> treatments in 
> > a split-plot design.)
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] lmer and a response that is a proportion

[John Fox]

Dear Cameron,

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Cameron Gillies
> Sent: Sunday, December 03, 2006 1:58 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] lmer and a response that is a proportion
> 
> Greetings all,
> 
> I am using lmer (lme4 package) to analyze data where the 
> response is a proportion (0 to 1).  It appears to work, but I 
> am wondering if the analysis is treating the response 
> appropriately - i.e. can lmer do this?
>

As far as I know, you can specify the response as a proportion, in which
case the binomial counts would be given via the weights argument -- at least
that's how it's done in glm(). An alternative that should be equivalent is
to specify a two-column matrix with counts of "successes" and "failures" as
the response. Simply giving the proportion of successes without the counts
wouldn't be appropriate.
 
> I have used both family=binomial and quasibinomial - is one 
> more appropriate when the response is a proportion?  The 
> coefficient estimates are identical, but the standard errors 
> are larger with family=binomial.
>

The difference is that in the binomial family the dispersion is fixed to 1,
while in the quasibinomial family it is estimated as a free parameter. If
the standard errors are larger with family=binomial, then that suggests that
the data are underdispersed (relative to the binomial); if the difference is
substantial -- the factor is just the square root of the estimated
dispersion -- then the binomial model is probably not appropriate for the
data.

I hope this helps,
 John
 
> Thanks very much for any insight you may have!
> Cam
> 
> 
> Cam Gillies
> PhD Candidate
> Biological Sciences
> University of Alberta
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] error in hetcor function (polycor package)?

[John Fox]

Dear Brendan,

That's curious, because the use argument to hetcor() works fine for me (see
below). Is it possible that you tried to use this argument without
specifying a data frame as the first argument to hetcor? If so, please see
?hetcor. If not, it would help to have an example.

I hope this helps,
 John

 snip -

> set.seed(12345) # adapting the example in ?hetcor
> R <- matrix(0, 4, 4)
> R[upper.tri(R)] <- runif(6)
> diag(R) <- 1
> R <- cov2cor(t(R) %*% R)
> round(R, 4)  # population correlations
   [,1]   [,2]   [,3]   [,4]
[1,] 1. 0.5848 0.5718 0.6233
[2,] 0.5848 1. 0.7374 0.6249
[3,] 0.5718 0.7374 1. 0.5923
[4,] 0.6233 0.6249 0.5923 1.
> data <- rmvnorm(1000, rep(0, 4), R)
> round(cor(data), 4)   # sample correlations
   [,1]   [,2]   [,3]   [,4]
[1,] 1. 0.5933 0.5659 0.6088
[2,] 0.5933 1. 0.7334 0.6230
[3,] 0.5659 0.7334 1. 0.5802
[4,] 0.6088 0.6230 0.5802 1.
> x1 <- data[,1]
> x2 <- data[,2]
> y1 <- cut(data[,3], c(-Inf, .75, Inf))
> y2 <- cut(data[,4], c(-Inf, -1, .5, 1.5, Inf))
> data <- data.frame(x1, x2, y1, y2)
> data[1,1] <- NA
> hetcor(data) 

Two-Step Estimates

Correlations/Type of Correlation:
   x1  x2 y1 y2
x1  1 Pearson Polyserial Polyserial
x2 0.5932   1 Polyserial Polyserial
y1 0.5952  0.7409  1 Polychoric
y2  0.624  0.6316 0.5708  1

Standard Errors:
x1  x2 y1
x1   
x2 0.02053   
y1 0.03092 0.02296   
y2 0.02027 0.01995 0.0374

n = 999 

P-values for Tests of Bivariate Normality:
   x1 x2  y1
x1  
x2 0.4782   
y1 0.4023 0.8871
y2 0.1166 0.5077 0.05526

> hetcor(data, use = "complete.obs")

Two-Step Estimates

Correlations/Type of Correlation:
   x1  x2 y1 y2
x1  1 Pearson Polyserial Polyserial
x2 0.5932   1 Polyserial Polyserial
y1 0.5952  0.7409  1 Polychoric
y2  0.624  0.6316 0.5708  1

Standard Errors:
x1  x2 y1
x1   
x2 0.02053   
y1 0.03092 0.02296   
y2 0.02027 0.01995 0.0374

n = 999 

P-values for Tests of Bivariate Normality:
   x1 x2  y1
x1  
x2 0.4782   
y1 0.4023 0.8871
y2 0.1166 0.5077 0.05526

> hetcor(data, use = "pairwise.complete.obs")

Two-Step Estimates

Correlations/Type of Correlation:
   x1  x2 y1 y2
x1  1 Pearson Polyserial Polyserial
x2 0.5932   1 Polyserial Polyserial
y1 0.5952  0.7409  1 Polychoric
y2  0.624  0.6317 0.5711  1

Standard Errors/Numbers of Observations:
x1  x2  y1   y2
x1 999 999 999  999
x2 0.0205310001000 1000
y1 0.03092 0.022951000 1000
y2 0.02027 0.01994 0.03738 1000

P-values for Tests of Bivariate Normality:
   x1 x2  y1
x1  
x2 0.4952   
y1 0.4023  0.878
y2 0.1166 0.5255 0.05615


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of BRENDAN KLICK
> Sent: Friday, December 01, 2006 1:42 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] error in hetcor function (polycor package)?
> 
> I have been using the hetcor function in the polycor package. 
>  When I don't specify the use option everything runs 
> smoothly.  However, when I specify use either as 
> "pairwise.complete.obs" or "complete.obs" I get this error
>  
> Error in optim(rho, f, control = control, hessian = TRUE, method =
> "BFGS") : 
> non-finite value supplied by optim
> 
> Is this an error in the hetcor function or am I missing something. 
> Thanks for your help.
>  
> Brendan Klick
> Johns Hopkins University School of Medicine [EMAIL PROTECTED]
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Box Tidwell / Error Message

[John Fox]

Dear Simon,

It's hard to tell without the data and more information about the nature of
the variables, but I suspect that the program is running into numerical
difficulties because of a flat likelihood at the maximum. Is age2 by any
chance age^2? How is year (in either form) related to age?

I hope this helps,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Simon P. Kempf
> Sent: Friday, December 01, 2006 3:30 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Box Tidwell / Error Message
> 
> Dear R-Users,
> 
>  
> 
> I used the box.tidwell () function of the car Package. 
> 
>  
> 
> When I used the following formula:
> 
>  
> 
> semi.sub.in.mi1.boxtidwell_h<-box.tidwell(RENT_LG ~ 
> I(age+1)+I(age2+1)+X06A
> + I(X08B+1) + I(X22+1) + I(X24+1) + X31A, ~B_YEAR + C_X01 + C_X14 + 
> + C_X19 +
> C_X29A +C_X21 + C_X23 + D_X12 + D_X17 + D_X18 + D_X25 + D_X27 
> + D_X30 +
> D_X32 + D_X35, data = semi.sub.in.mi1)
> 
>  
> 
> everything is fine.
> 
>  
> 
> However, when I replaced the time dummy variable:
> 
>  
> 
> semi.sub.in.mi1.boxtidwell_h<-box.tidwell(RENT_LG ~ 
> I(age+1)+I(age2+1)+X06A
> + I(X08B+1) + I(X22+1) + I(X24+1) + X31A, ~B_HALF + C_X01 + C_X14 + 
> + C_X19 +
> C_X29A +C_X21 + C_X23 + D_X12 + D_X17 + D_X18 + D_X25 + D_X27 
> + D_X30 +
> D_X32 + D_X35, data = semi.sub.in.mi1)
> 
>  
> 
> I get the following error message:
> 
>  
> 
> Error in lm.fit(x, y, offset = offset, singular.ok = 
> singular.ok, ...) : 
> 
> NA/NaN/Inf in foreign function call (arg 1)
> 
>  
> 
> When I use the following formula (I deleted the I(age2+1) term):
> 
>  
> 
> semi.sub.in.mi1.boxtidwell_h<-box.tidwell(RENT_LG ~ I(age+1)+X06A +
> I(X08B+1) + I(X22+1) + I(X24+1) + X31A, ~B_HALF + C_X01 + 
> C_X14 + C_X19 + C_X29A +C_X21 + C_X23 + D_X12 + D_X17 + D_X18 
> + D_X25 + D_X27 + D_X30 +
> D_X32 + D_X35, data = semi.sub.in.mi1)
> 
>  
> 
> It works.
> 
>  
> 
> Some background information:
> 
>  
> 
> -  The formula with the predictors to be transformed 
> contains only
> variables which are >0.
> 
> -  The data set does not have any missing values
> 
> -  B_YEAR is a factor with 10 levels
> 
> -  B_HALF is a factor with 20 levels
> 
> -  The data set contains more than 19000 observations.
> 
>  
> 
> Now, I am bit confused. Why does the function works when I 
> use B_YEAR respecitvely why does it work with B_HALF when I 
> delete I(age2+1)
> 
>  
> 
> Thanks in advance,
> 
>  
> 
> Simon
> 
>  
> 
>  
> 
>  
> 
>  
> 
> Simon P. Kempf 
> 
> Dipl.-Kfm. MScRE Immobilienvkonom (ebs)
> 
> Wissenschaftlicher Assistent
> 
>  
> 
> B|ro:
> 
> IREBS Immobilienakademie
> 
> c/o ebs Immobilienakademie GmbH
> 
> Berliner Str. 26a
> 
> 13507 Berlin
> 
>  
> 
> Privat:
> 
> Dunckerstra_e 60
> 
> 10439 Berlin
> 
>  
> 
> Mobil: 0176 7002 6687
> 
> Email:  <mailto:[EMAIL PROTECTED]> [EMAIL PROTECTED]
> 
>  
> 
> 
>   [[alternative HTML version deleted]]
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] String question

[John Fox]

Dear Carmen,

length(grep("[^01]", bits)) == 0

should do the trick, returning TRUE if the string contains only 0's and 1's.

I hope this helps,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Carmen Meier
> Sent: Wednesday, November 29, 2006 9:06 AM
> To: R-help
> Subject: [R] String question
> 
> Hi to all
> I would to determinate whether bits is a binary code and I 
> would to find out the which bit is set to 1
> 
> bits <-"00110110"
> I found to detect whether there are only numbers
> all.digits(bits)
> but is there any function to detect whether there are only 0 
> and 1 in the string
> 
> And how could I get the f.e the third "bit" from the right hand side
> 
> With regards Carmen
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Random effects ANOVA

[John Fox]

Dear Luis,

I believe that lmer() doesn't print out standard errors or Wald tests for
variance-covariance components because these tests aren't reliable. You can
omit each random effect from the model in turn, and use anova() to perform a
likelihood-ratio test for the corresponding variance and covariance
components, comparing each null model with the full model.

I hope this helps,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Luis Salasar
> Sent: Friday, November 24, 2006 1:57 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Random effects ANOVA
> 
> Hi,
> 
> I used the command lmer, in package lme4, to fit a random 
> effects ANOVA. But i didn't get the p-values of significance 
> tests of variance components. Does anyone know how  to do it? Thanks,
> 
> Luis Ernesto
> 
>   
> -
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to create this design matrix?

[John Fox]

Dear Michael,

If I follow what you want to do, which is not altogether certain, what
you're describing as the "design matrix" is intended to constitute the
left-hand side of the model. That is, your d1,1, ..., d1,12 are the 12
response variables. If so, I'm not sure why you think that my original
suggestion won't give you the answer, but I guess that at this point I
should leave to someone else who understands more clearly what you want to
sort it out.

Regards,
 John

----
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: Michael [mailto:[EMAIL PROTECTED] 
> Sent: Wednesday, November 15, 2006 9:49 PM
> To: John Fox
> Subject: Re: [R] how to create this design matrix?
> 
> How about I make a design matrix as follows:
> 
> 1  d1,1  0   0 0   0 0   0   0
> 00   1   d1,2  0   0 0   0   0
> 00   0   0 1   d1,3  0   0   0
> ...
> ...
> 
> 000  0 0   0  ...  ... 1  d1,12 
> 
> 
> The above matrix will work for the 1st row of Y data;
> 
> d1, 1 means the 1st column of the 1st row of data; d1, 12 
> means the 12th column of the 1st row of data.
> 
> 
> But now since I have N samples(rows) of Y data; how do I 
> conceptually do a 3D design matrix? 
> 
> 
> On 11/15/06, John Fox <[EMAIL PROTECTED]> wrote:
> 
>   Dear Michael,
>   
>   > -Original Message-
>   > From: Michael [mailto:[EMAIL PROTECTED]
>   > Sent: Wednesday, November 15, 2006 2:40 PM
>   > To: John Fox
>   > Cc: R-help@stat.math.ethz.ch
>   > Subject: Re: [R] how to create this design matrix?
>   >
>   > There are 12 response variables, columns 1 to 12 are response
>   > variables, i.e., these are y's, they all regress to the 13th
>   > column, which is the predictor, i.e. the X.
>   >
>   
>   Right.
>   
>   > Let's take column 1, call this Y1, and there are n rows(n
>   > samples) of it,
>   >
>   > I need Y1= b0_1 + b1_1* X + epsilon, where X is the 
> 13th column
>   >
>   > Similarly, for column 1 to column 12, we do the above,
>   >
>   > Y12= b0_12 + b1_12 * X + epsilon, where Y12 is the 
> 12th column, 
>   >
>   > they all have different b0's and b1's.
>   
>   Right.
>   
>   > Totally there are 24 b0's and b1's.
>   >
>   
>   Yes.
>   
>   > I want a group regression, not separated regression...
>   >
>   
>   I'm not sure what you mean by a "group regression" 
> rather than "separated 
>   regressions." The multivarite linear regression that I 
> suggested will give
>   you 12 slopes and 12 intercepts. They are exactly what 
> you'd get from 12
>   individual least-squares regression of each Y on X, but 
> the multivariate 
>   regression can also give you, e.g., the covariances 
> among all of the
>   coefficients (if you want them).
>   
>   John
>   
>   > Thanks
>   >
>   >
>   >
>   >
>   > On 11/15/06, John Fox < [EMAIL PROTECTED] 
> <mailto:[EMAIL PROTECTED]> > wrote:
>   >
>   >   Dear Michael,
>   >
>   >   This looks like a multivariate simple regression --
>   > that is, 12 response
>   >   variables, one predictor. If the data are in the matrix 
>   > X, then lm(X[,1:12]
>   >   ~ X[,13]) should do the trick.
>   >
>   >   I hope this helps,
>   >   John
>   >
>   >   
>   >   John Fox
>   >   Department of Sociology 
>   >   McMaster University
>   >   Hamilton, Ontario
>   >   Canada L8S 4M4
>   >   905-525-9140x23604
>   >   http://socserv.mcmaster.ca/jfox
>   >   
>   >
>   >   > -Original Message-
>   >   > From: [EMAIL PROTECTED]
>   >   > [mailto: [EMAIL PROTECTED]
>   > <mailto:[EMAIL PROTECTED]> ] On Behalf 
> Of Michael
>   >   > Sent: Wednesday, November 15, 2006 12:23 AM
>   >   > To: R-help@stat.math.ethz.ch
>   >   > Subject: [R] how to create this design m

Re: [R] how to create this design matrix?

[John Fox]

Dear Michael,

This looks like a multivariate simple regression -- that is, 12 response
variables, one predictor. If the data are in the matrix X, then lm(X[,1:12]
~ X[,13]) should do the trick.

I hope this helps,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Michael
> Sent: Wednesday, November 15, 2006 12:23 AM
> To: R-help@stat.math.ethz.ch
> Subject: [R] how to create this design matrix?
> 
> Hi all,
> 
> I have a multiple-linear regression problem.
> 
> There are 13 columns of data, the whole data matrix is:  n x 
> 13, where n is the number of samples.
> 
> Now I want to regress  EACH of the first 12 columns onto the 
> 13th column, with 2-parameter linear model  y_i = b0 + b1 * 
> x_i, where i goes from 1 to n, and b0 is the intercept.
> 
> How do I create a design matrix to do the 12-column 
> regression collectively all at once using multiple linear regressions?
> 
> Thanks a lot
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Making a case for using R in Academia

[John Fox]

Dear Charilaos,

It's very difficult to give definitive answers to the questions that you
pose because we don't have any good data (at least as far as I know) about
how widely R is used. I recall a claim, I think on the r-help list, that R
is now second to SAS in use world-wide, but I'm not sure how one would
establish that. I do, however, have anecdotal evidence about R use that may
be of some help, and perhaps others can contribute their impressions,
especially if they differ from mine.

I think that it's fair to say that the S language is more or less the
standard among statisticians, and that R has now far surpassed S-PLUS in
use. For rough evidence one could look, for example, to the relative levels
of activity on the r-help and s-news email lists.

Among social scientists the picture is not as clear. My impression is that
SPSS is used very widely for low-levels methods courses taught to
undergraduates, and not very extensively in the best social-science graduate
programmes. I would expect that, at present, Stata use in social-science
graduate programmes exceeds R, and that SAS and R would also be used fairly
widely. In my opinion, these are the only reasonable choices -- I don't
think that SPSS is sufficiently capable to compete with R, Stata, or SAS.
There are, for example, several different packages used at the ICPSR Summer
Program in Quantitative Methods for Social Research, but several relatively
advanced courses now use R. Likewise, the Oxford Spring School, hosted by
the Department of Politics and International Relations at Oxford, has mostly
employed R and Stata.

Of course, my own preference is for R.

Regards,
 John

----
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of 
> Charilaos Skiadas
> Sent: Thursday, November 09, 2006 11:18 AM
> To: r-help@stat.math.ethz.ch
> Subject: Re: [R] Making a case for using R in Academia
> 
> As a addendum to all this, this is one of the responses I got 
> from one of my colleagues:
> 
> "The problem with R is that our students in many social science  
> fields, are expected to know SPSS when they go to graduate school.   
> Not having a background in SPSS would put these students at a 
> disadvantage."
> 
> Is this really the case? Does anyone have any such statistics?
> 
> Charilaos Skiadas
> Department of Mathematics
> Hanover College
> P.O.Box 108
> Hanover, IN 47243
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] correlation argument for lmer?

[John Fox]

Dear r-help members,

Can lmer() in the lme4 package fit models that have a specified within-group
correlation structure, as provided, for example, by the correlation argument
to lme() in the nlme package?

Thanks,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] poly() question

[John Fox]

Dear Jonathan,

An alternative is to use the raw argument to poly(), possibly after
centering the predictor. This probably will work OK as long as the degree of
the polynomial isn't high.

I hope this helps,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of 
> Jonathan Greenberg
> Sent: Wednesday, November 01, 2006 9:39 PM
> To: R-help
> Subject: [R] poly() question
> 
> Besides the primary citation, "Kennedy, W. J. Jr and Gentle, 
> J. E. (1980) Statistical Computing Marcel Dekker." (which is 
> $300 and my library doesn't have it), is there any other 
> documentation on how to take a poly() object and predict "by 
> hand" new data?  E.g. What do those coefficients actually 
> mean ("The orthogonal polynomial is summarized by the 
> coefficients, which can be used to evaluate it via the 
> three-term recursion...")?  We created a GLM model with a 
> poly() term and I'm trying to apply this model in another 
> program (grass gis via mapcalc) without requiring the direct 
> link with R if at all possible.  We'd like to avoid making a 
> lookup table if at all possible.  Thanks!
> 
> --j
> 
> --
> Jonathan A. Greenberg, PhD
> NRC Research Associate
> NASA Ames Research Center
> MS 242-4
> Moffett Field, CA 94035-1000
> Office: 650-604-5896
> Cell: 415-794-5043
> AIM: jgrn307
> MSN: [EMAIL PROTECTED]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] coefficient constraints

[John Fox]

Dear John,

You can subtract the variables in question from the left-hand-side of the
model [e.g., lm(y - x1 - x2 ~ 1)], or use the offset argument to lm [e.g.,
lm(y ~ 1, offset=x1 + x2)], or use offset() in the model formula [e.g., lm(y
~ offset(x1 + x2)]. All of these options are equivalent.

I hope this helps,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of John Gauss
> Sent: Tuesday, October 31, 2006 7:28 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] coefficient constraints
> 
> Hi,
> 
> Thank you for your help.  I'm trying to constrain the 
> coefficients on a linear regression model, <-lm(...), to all 
> equal one, except for the intercept.  I've searched 
> help(glm), help(model.fit), etc. but I can not find where to 
> add the constraint.  Your help would be greatly appreciated.
> 
> Thanks.
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Qualitative Data??(String command)

[John Fox]

Dear Davia,

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Davia Cox
> Sent: Friday, October 27, 2006 3:37 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Qualitative Data??(String command)
> 
> I am using the read.table function to load an Excel data set into R.
> It has a few variables with very long qualitative (free 
> response typically in sentences) response that I would like 
> to keep, but would like to limit the "length" of the response 
> that R shows. Is there some sort of string or column width 
> command I can include in the read.table function to limit the 
> length of words used in the variables that have long 
> responses?? I do know which variable names are the ones with 
> the qualitative responses. Will the Rcmdr program do this for 
> me? 

No. When you read an Excel file (or an ascii file) into R using the Rcmdr,
character data will be converted into factors.

> I know STATA has this function but I don't have it on my 
> computer to use.
>

If I understand you correctly, it's not hard to do what you want. The
following function takes a data frame and desired character length as
arguments and truncates both character variables and the levels of factors.
Note that you have to be a little careful with factors in that you could
truncate different levels to the same string, inadvertently combining them.

truncateValues <- function(Data, len=10){
for (i in 1:ncol(Data)){
var <- Data[[i]]
if (is.factor(var)) {
levels(var) <- sapply(levels(var), 
function(x) substr(x, 1, min(nchar(x), len)))
Data[[i]] <- var
}
else if (is.character(var)) {
Data[[i]] <- sapply(var, function(x) substr(x, 1, min(nchar(x),
len)))
}
}
Data
}

I hope this helps,
 John

> Davia
> 
> 
> --
> Davia S. Cox
> "Image creates desire. You will what you imagine."
> ~J.G. Gallimore~
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Error using "White-corrected" Anova (white.adjust="hc3")

[John Fox]

Dear Kilian,

This error was introduced when the linear.hypothesis function was modified,
and I've not yet fixed it because I'm not entirely sure what's the best way
to proceed. It's a bad idea simply to allow Anova() to throw a cryptic
error, so I'll try to resolve the issue some time in the next several days.

Sorry,
 John

----
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Kilian_Wasmer
> Sent: Friday, October 27, 2006 5:19 AM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Error using "White-corrected" Anova (white.adjust="hc3")
> 
> Hello.
> 
> I try to compute a "White-corrected" ANOVA. Therefore I use 
> the command "Anova" implemented in the car()-package):
> Anova(modelname, white.adjust="hc3")
> and receive the error message:
> "Fehler in SS[i] <- SS.term(names[i]) : nichts zu ersetzen" 
> (German) "Error in SS[i] <- SS.term(names[i]) : nothing to 
> replace" (probably being the message in the English version) 
> I also tried the example mentioned in Fox (2002, "An R and 
> S-Plus companion..."), with the "ornstein" data set but 
> received the same message.
> If I don't use the argument "white.adjust="hc3" the standard 
> anova is computed without any error messages..
> 
> Any hints?
> 
> Kilian
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Recoding categorical variables

[John Fox]

Dear Murray,

How about as.numeric(factor(y)) ?

I hope this helps,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Murray 
> Jorgensen
> Sent: Wednesday, October 25, 2006 7:13 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Recoding categorical variables
> 
> I want to recode an integer-valued variable y so that its 
> values become 1:length(y). I can do this using a loop but 
> maybe someone can suggest code without a loop. My code is this:
> 
> y <- round(20*runif(25))
> table(y)
> suy <- sort(unique(y))
> m <- length(suy)
> z <- y + max(suy)
> for(i in 1:m) z[y==suy[i]] <- i
> rbind(y,z)
> 
> (the recoded y is stored in z)
> 
> Murray Jorgensen
> -- 
> Dr Murray Jorgensen  http://www.stats.waikato.ac.nz/Staff/maj.html
> Department of Statistics, University of Waikato, Hamilton, New Zealand
> Email: [EMAIL PROTECTED]Fax 7 838 4155
> Phone  +64 7 838 4773 wkHome +64 7 825 0441Mobile 021 1395 862
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Help with "recode" and "factor" functions

[John Fox]

Dear Chris,

I hesitate to answer a question arising from an exam -- you should probably
address the question to the person who set the exam -- but it doesn't hurt,
I think, to point out the following: The recode function that you're using
is in the car package, which is associated with a book (see ?car) that has
several examples of the use of the function.

Regards,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Chris Linton
> Sent: Monday, October 23, 2006 4:16 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Help with "recode" and "factor" functions
> 
> I have a data set with seven inputs.  Four of which are 
> categorical.  For my midterm, my professor wants us to scale 
> all the inputs.  This means, I pressume, that I have to use 
> 'recode' or 'factor' to transform the categorical data in 
> numerical.  For example, one input variable is 
> 'race=(b,w,h,o)'.  I just want to assign a numerical value to 
> all 'b,w,h,o'.  I thought 'recode' should do this, but it 
> doesn't work.  Here's the code I'm using for recode:
> 
> recode(race, "b='1';w='2';h='3';o='4'")
> 
> this is the error I get:
> Error in eval(expr, envir, enclos) : object "o" not found
> 
> 
> It's not that there's no "o".  If I change the order or 
> combination of the variables, it always can't find one of them.
> 
> I could also use 'factor', from what I hear.  But, I looked 
> at the help section on this function and I ended up more confused.
> 
> 
> How do I code it so these variables take on numerial values?  
> I need to be able to use:
> 
> race.centered = race - mean(race)
> 
> 
> This scaling code doesn't really make sense if the values of 
> 'race' are non-numerical.  I might end up dividing by 2 SD's 
> as well.  But, I don't know if I need to.  I'll have to do 
> some more reading.
> 
> 
> Thank you for your help!
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] How to get multiple Correlation Coefficients

[John Fox]

Dear David and Kum-Hoe Hwang,

David has pointed out to me that polychor() will incorrectly return a
polychoric correlation when its arguments are length-one character vectors:

> polychor("a", "b")
[1] 0.1055909

Actually, the problem is more general, since polychor() will erroneously
compute a polychoric correlation in any event when the contingency table for
the data contains only one cell:

> polychor(1, 2)
[1] 0.1055909
> polychor(rep("a", 100), rep("b", 100))
[1] 0.1055909

I've fixed this bug so that polychor() and polyserial() report an error when
a categorical variable has just one level. I'll upload the new version of
the package to CRAN shortly.

Regards,
 John


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of David Barron
> Sent: Thursday, October 19, 2006 5:54 AM
> To: Kum-Hoe Hwang; r-help
> Subject: Re: [R] How to get multiple Correlation Coefficients
> 
> The problem is that in the expression  polychor(vars[i], 
> vars[j]), vars[i] and vars[j] refer to the names of the 
> variables, not the variables themselves.  So, use sdi[,i] and 
> sdi[,j] instead.
> 
> On 19/10/06, Kum-Hoe Hwang <[EMAIL PROTECTED]> wrote:
> > Hi
> >
> > I have used a polycor package for categorical correlation 
> coefficients.
> > I run the following script. But there were no results.
> >
> > Could you tell me how to correct the script?
> >
> > Thanks in advance,
> >
> > vars <- names(sdi)
> > for (i in 1:length(vars)) {
> > for (j in 1:length(vars)) {
> >   paste(vars[i]," and ", vars[j])
> >   polychor(vars[i], vars[j])
> >   # corr
> > }
> > }
> >
> >
> >
> > --
> > Kum-Hoe Hwang, Ph.D.Phone : 82-31-250-3516Email : [EMAIL PROTECTED]
> >
> > __
> > R-help@stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide 
> > http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
> 
> 
> --
> =
> David Barron
> Said Business School
> University of Oxford
> Park End Street
> Oxford OX1 1HP
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] regression analyses using a vector of means anda variance-covariance matrix

[John Fox]

Dear John,

I'm not aware of an existing function that will do what you want
conveniently, but it's not hard to write one:

reg <- function(V, means, y, x, n){
# V: covariance matrix of variables
# means: vector of means
# y: index of response in V
# x: indices of regressors in V
SSP <- (n - 1)*V + n*outer(means, means)
SSP <- cbind(n*means, SSP)
SSP <- rbind(c(n, n*means), SSP)
XtXinv <- solve(SSP[c(1,x+1), c(1,x+1)])
b <- as.vector(XtXinv %*% SSP[c(1,x+1), y+1])
R2 <- as.vector(V[y,x] %*% solve(V[x,x]) %*% V[x,y])/V[y,y]
s2 <- (1 - R2)*V[y,y]*(n - 1)/(n - length(x) - 1)
Vb <- s2*XtXinv
list(coef=b, vcov=Vb, R2=R2, s2=s2, SSP=SSP)
}

For example,

> library(car)
> data <- Duncan[,c("prestige", "income", "education")]
> reg(var(data), colMeans(data), 1, 2:3, nrow(data))
$coef
[1] -6.0646629  0.5987328  0.5458339

$vcov
  incomeeducation
  18.2494814 -0.15184501 -0.150706025
income-0.1518450  0.01432027 -0.008518551
education -0.1507060 -0.00851855  0.009653582

$R2
[1] 0.8281734

$s2
[1] 178.7309

$SSP
   prestige income education
45 2146   1884  2365
prestige  2146   146028 118229147936
income1884   118229 105148122197
education 2365   147936 122197163265

> summary(lm(prestige ~ income + education, data=Duncan)) # check

Call:
lm(formula = prestige ~ income + education, data = Duncan)

Residuals:
 Min   1Q   Median   3Q  Max 
-29.5380  -6.4174   0.6546   6.6051  34.6412 

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.064664.27194  -1.4200.163
income   0.598730.11967   5.003 1.05e-05 ***
education0.545830.09825   5.555 1.73e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 13.37 on 42 degrees of freedom
Multiple R-Squared: 0.8282, Adjusted R-squared:  0.82 
F-statistic: 101.2 on 2 and 42 DF,  p-value: < 2.2e-16

Some caveats: This function will work OK for well conditioned data, since
I'm inverting the matrix of sums of squares and products; one could take a
more sophisticated approach. As well, there's a bit of redundancy in the
computation of R2, but I didn't have time to figure out how to avoid it.
Finally, there's an obvious advantage to computing on the original data --
e.g., one can get residuals.

I hope this helps,
 John 


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of John Sorkin
> Sent: Saturday, October 14, 2006 3:27 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] regression analyses using a vector of means anda 
> variance-covariance matrix
> 
> R 2.2.0
> windows XP
> 
> How can I perform a regression analyses using a vector of 
> means, a variance-covariance matrix? I looked at the help 
> screen for lm and did not see any option for using the afore 
> mentioned structures as input to lm.
> Thanks,
> 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) [EMAIL PROTECTED]
> 
> Confidentiality Statement:
> This email message, including any attachments, is for the\...{{dropped}}

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Differences in estimates calculated in R and SPSS

[John Fox]

Dear Katja,

If the fitted values are the same and the coefficients differ, then the
parametrization of the models differ. I suspect that SPSS is using
sum-to-zero contrasts (contr.sum in R, coded -1, 0, 1), while R uses
dummy-coded contrasts (contr.treatment, coded 0,1) by default. If this is
the case, then you should get the same coefficients in R by setting
options(contrasts=c("contr.sum", "contr.poly")). There are other ways of
doing this as well. 

See ?options, ?contr.treatment, etc., and section 11.1.1 of the introductory
manual that comes with R (and perhaps read something more detailed about how
R handles factors in linear models).

I hope this helps,
 John

--------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Katja 
> Walentowitz
> Sent: Saturday, October 14, 2006 11:08 AM
> To: Dirk Eddelbuettel
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] Differences in estimates calculated in R and SPSS
> 
> I have compared the models, they are exactly the same. The 
> same for the data. The problem occurs most often when I use 
> the aov-function.
> When I calculate the fitted values, they are the same, but 
> the intercept is somewhat different. In R the intercept is 
> more or less the mean of the regressand values, as defined in 
> some statistical books - but still: it is not the exact mean. 
> In SPSS it is not.
> 
> Thanks again for helping me.
> Katja
> 
> On 10/14/06, Dirk Eddelbuettel <[EMAIL PROTECTED]> wrote:
> >
> >
> > On 14 October 2006 at 16:59, Katja Walentowitz wrote:
> > | I am a bit confused because I get different estimates for the 
> > | regression coefficients from R and SPSS.
> > | Does anyone have an explanation for this?  As I am not so 
> > | experienced in data analysis, I don't know why this happens.
> >
> > Maybe you are not estimating the same model -- one may have an 
> > intercept, and the other doesn't.  Or maybe something happens while 
> > you read the data -- start by comparing summary statistics, 
> > tabulations and plots of your data.
> >
> > But as you didn't tell what you estimated, and how, it is a 
> tad hard 
> > for use to guess.  There is a Posting Guide recommending 
> how to phrase 
> > questions and what to suplly ...
> >
> > Dirk
> >
> > --
> > Hell, there are no rules here - we're trying to accomplish 
> something.
> >   -- Thomas 
> A. Edison
> >
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.