Re: [R] lm: Displaying 95% confidence and prediction in tervals on scatterplots
jenny tan jennytimp at hotmail.com writes: May I know how does one superimpose the 95% confidence and prediction intervals on the linear regression line of a scatterplot? Check http://finzi.psych.upenn.edu/R/Rhelp02a/archive/8380.html Dieter __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Questions about extract-lme.cov
Li, Hua HUL at stowers-institute.org writes: I am using extract.lme.cov to extract the covariance matrix of lme. But the results are not expected. b - lme(travel~1,Rail,~1|Rail) The default correlation for lme is no correlation within groups. extract.lme.cov(b,Rail) The part of covariance matrix looks like: 123456 1 631.4778 615.3111 615.3111 0. 0. 0. 2 615.3111 631.4778 615.3111 0. 0. 0. 3 615.3111 615.3111 631.4778 0. 0. 0. 4 0. 0. 0. 631.4778 615.3111 615.3111 5 0. 0. 0. 615.3111 631.4778 615.3111 6 0. 0. 0. 615.3111 615.3111 631.4778 Where does the covariance come from? Maybe I'm missing sth here? Maybe VarCorr(b) gives all you need, because extract.lme.cov2(b,Rail,start.level=2) or even extract.lme.cov(b,Rail,start.level=2) are a bit redundant. Dieter __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] lm: Displaying 95% confidence and prediction intervals on scatterplots
May I know how does one superimpose the 95% confidence and prediction intervals on the linear regression line of a scatterplot? You could use ggplot: install.packages(ggplot) library(ggplot) qplot(wt, mpg, data=mtcars, type=c(point,smooth), method=lm) (which gives the 95% confidence interval) 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
[R] How to Interpret Results of Regression in R
- Howdy, Gurus I am appying R package for regression analysis as followings. A dependent variable is jhnet that means ratio of dividing internal trip with all trips in a traffic zone. There are many indepentent variables including factor or dummy varibles such as parkfee, ohouse, Devt2, corridor1. I have three questions. First, What are estimated for regression model? Second, Which results should I trust among results from lm, anova, Anova? Third, are there any differences among results from lm, anova, Anova? Thank you in advance, - Call: lm(formula = I(jhnet^1) ~ as.factor(parkfee) + fare + as.factor(ohouse) + tripMIN + as.factor(Devt2) + cityarea + coreDistance + density + as.factor(corridor1) + as.factor(ic) + OD_total + nfirm_dong + house1 + house2 + house3 + house4, data = sdi.data[w, ], na.action = na.omit, singular.ok = TRUE) Residuals: Min 1Q Median 3Q Max -0.27703 -0.02942 0.00552 0.02325 0.18894 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 5.68e-01 1.11e-02 51.15 2e-16 *** as.factor(parkfee)1 3.89e-03 6.40e-030.61 0.5430 as.factor(parkfee)2 3.84e-03 2.51e-031.53 0.1255 fare 7.40e-07 3.07e-060.24 0.8097 as.factor(ohouse)2 -3.49e-03 2.55e-03 -1.37 0.1704 as.factor(ohouse)34.49e-04 4.83e-030.09 0.9260 as.factor(ohouse)4 -2.16e-02 8.11e-03 -2.66 0.0078 ** tripMIN -6.14e-05 3.75e-05 -1.64 0.1020 as.factor(Devt2)1-2.57e-02 2.89e-03 -8.91 2e-16 *** as.factor(Devt2)2 2.29e-02 3.57e-036.42 1.5e-10 *** as.factor(Devt2)3 6.04e-02 3.34e-03 18.10 2e-16 *** cityarea -6.59e-05 1.17e-05 -5.61 2.2e-08 *** coreDistance 2.59e-07 2.39e-071.09 0.2773 density 7.55e-06 3.28e-07 23.02 2e-16 *** as.factor(corridor1)B2.54e-02 5.40e-034.70 2.7e-06 *** as.factor(corridor1)C7.86e-02 1.19e-026.58 5.5e-11 *** as.factor(corridor1)D 8.39e-02 7.76e-03 10.81 2e-16 *** as.factor(corridor1)E -2.72e-01 2.22e-02 -12.22 2e-16 *** as.factor(corridor1)F -4.39e-02 5.21e-03 -8.43 2e-16 *** as.factor(ic)11.51e-02 5.57e-032.70 0.0069 ** OD_total -1.36e-08 2.88e-08 -0.47 0.6375 nfirm_dong -9.77e-06 1.73e-06 -5.66 1.6e-08 *** house1 -4.41e-06 2.27e-07 -19.44 2e-16 *** house2 -8.13e-08 8.12e-08 -1.00 0.3165 house37.40e-06 3.22e-07 23.03 2e-16 *** house4 -3.83e-06 2.40e-07 -15.96 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0546 on 2974 degrees of freedom Multiple R-Squared: 0.584, Adjusted R-squared: 0.58 F-statistic: 167 on 25 and 2974 DF, p-value: 2e-16 Analysis of Variance Table Response: I(jhnet^1) Df Sum Sq Mean Sq F value Pr(F) as.factor(parkfee) 2 0.060.03 10.05 4.5e-05 *** fare1 0.030.03 10.77 0.0010 ** as.factor(ohouse) 3 0.080.038.90 7.2e-06 *** tripMIN 1 0.020.028.35 0.0039 ** as.factor(Devt2)3 1.770.59 197.95 2e-16 *** cityarea1 2.552.55 855.60 2e-16 *** coreDistance1 2.472.47 827.62 2e-16 *** density 1 0.790.79 265.23 2e-16 *** as.factor(corridor1)5 1.370.27 91.91 2e-16 *** as.factor(ic) 1 0.0035 0.00351.17 0.2787 OD_total1 0.030.03 10.05 0.0015 ** nfirm_dong 1 0.390.39 131.97 2e-16 *** house1 1 0.490.49 164.06 2e-16 *** house2 1 0.500.50 167.25 2e-16 *** house3 1 1.111.11 373.09 2e-16 *** house4 1 0.760.76 254.75 2e-16 *** Residuals2974 8.86 0.0030 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 [1] Anova Results Anova Table (Type II tests) Response: I(jhnet^1) Sum Sq Df F value Pr(F) as.factor(parkfee) 0.0121.19 0.3035 fare 0.0001710.06 0.8097 as.factor(ohouse) 0.0332.89 0.0343 * tripMIN 0.0112.68 0.1020 as.factor(Devt2)1.773 198.18 2e-16 *** cityarea0.091 31.50 2.2e-08 *** coreDistance 0.0035211.18 0.2773 density
[R] Find peaks in histograms / Analysis of cumulative frequency
Hello all, I have some histograms of amount of DNA in some cells (DU145 cells overexpressing Bax and Bcl-xL for those who wish to know). The histograms show not only two peaks as expected, but three, indicating that some cells have more than normal amounts of DNA. I am interested in knowing how much of the cell populations are in each peak as well as between. I am not really sure how to go about it; I have been considering fitting a gaussian distribution to each peak and integrate the part between the peaks as described by Watson et al (1987 Cytometry 8:1-8). A more straight forward and more visual approach appears to be plotting the cumulative frequencies. In either case, I should like to find the peaks in the histogram automatically, as well as getting proper information about the peaks. How would I go about finding peaks using R? Also I have really not been able to figure out how to fit a distribution. Is there a way to analyse the cumulative frequencies? the knots() function appears to return far too many knots. I am relatively new to R, but do have good programming experience, though I am mostly biologist. Thank you in advance for any inputs. PS. An example of the histogram can be found herehttp://photos1.blogger.com/blogger/7029/2724/1600/DU145-Bax3-Bcl-xL.png [[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
[R] termplot and ylim
Hi together, I always get an error message with using ylim in termplot(), like this: x-(1:10) y-(10:1) l-lm(y~x) termplot(l,ylim=c(1,2)) Is this a bug, or is there another possibility to do that? Especially, I would like to use term.plot() for gamlss objects. Thanks for your help! Andreas -- Echte DSL-Flatrate dauerhaft für 0,- Euro*! __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Error in Quantile Regression - Clear Message
[EMAIL PROTECTED] wrote: Dear Users, I loaded my dataset as following: presu - read.table(C:/_Ricardo/Paty/qtdata_f.txt, header=TRUE, sep=\t, na.strings=NA, dec=., strip.white=TRUE) dep-presu[,3]; exo-presu[,4:92]; When I try: rq(dep ~ exo, ...) or mle.stepwise(dep ~ exo, ...) I got the same error: rq(dep ~ exo) Error in model.frame(formula, rownames, variables, varnames, extras, extranames, : invalid variable type for 'exo' Any hint in how to fix it? I think this is due my data format. Maybe, but you will have to tell us the class() and mode() of exo, some output of the str() function would be fine. Uwe Ligges Thanks, Ricardo Gonçalves Silva, M. Sc. Apoio aos Processos de Modelagem Matemática Econometria Inadimplência Serasa S.A. (11) - 6847-8889 [EMAIL PROTECTED] __ Hi, I load my data set and separate it as folowing: presu - read.table(C:/_Ricardo/Paty/qtdata_f.txt, header=TRUE, sep=\t, na.strings=NA, dec=., strip.white=TRUE) dep-presu[,3]; exo-presu[,4:92]; Now, I want to use it using the wls and quantreg packages. How I change the data classes for mle and rq objects? Thanks a lot, Ricardo __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] termplot and ylim
It looks like a bug or at least an omission.
Try the following. It creates a new termplot function which
in turn defines a plot function which looks for a ylims= arg
and, if present, replaces the ylim= arg with ylims= and
then calls the real plot from the graphics package with the
new set of arguments.
It then defines a new proto object, i.e. an environment,
whose parent is the environment within the new termplot
we are setting up. It places a copy of termplot from the
stats package in that proto object naming the copy f.
This has the side effect of resetting f''s parent to
the proto object so that when f calls plot it finds
the plot we just defined instead of the usual plot.
Note that one uses ylims= instead of ylim in the call.
termplot - function(...) {
plot - function(...) { # replace ylim= with ylims=
args - list(...)
if (ylims %in% names(args)) {
args$ylim - args$ylims
args$ylims - NULL
}
do.call(graphics::plot, args)
}
proto(f = stats::termplot)[[f]](...)
}
# test
library(proto)
L - lm(y ~ x, data.frame(x = 1:10, y = 10:1))
termplot(L, ylims = 1:2)
On 7/15/06, Andreas Beyerlein [EMAIL PROTECTED] wrote:
Hi together,
I always get an error message with using ylim in termplot(), like this:
x-(1:10)
y-(10:1)
l-lm(y~x)
termplot(l,ylim=c(1,2))
Is this a bug, or is there another possibility to do that? Especially, I
would like to use term.plot() for gamlss objects.
Thanks for your help!
Andreas
--
Echte DSL-Flatrate dauerhaft für 0,- Euro*!
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Excel to R
Hi peolple! I have a many excel tables with mode than 100 variables. And I want use R to analize that. But I have a problem, a group of this variables (more than 50) in any table is a factor and other part is a number. Tha factors variables have tha values enconde this form (1=Yes,2=No and 9 = NA) Well I use this scripts to import the database require(RODBC) channel - odbcConnectExcel(f:/teste.xls) data - sqlFetch(channel, Sheet1) summary(data) qw ee Min. :1.000 Min. :1.000 1st Qu.:1.000 1st Qu.:1.500 Median :1.000 Median :2.000 Mean :1.333 Mean :2.429 3rd Qu.:1.750 3rd Qu.:3.500 Max. :2.000 Max. :4.000 NA's :1.000 But qw is a factor (and is colnum type isvtext) Is possible modify my script for this utcome summary(data) qw ee 1 :4 Min. :1.000 2 :2 1st Qu.:1.500 NA's:1 Median :2.000 Mean :2.429 3rd Qu.:3.500 Max. :4.000 Thanks in advance Bernardo Rangel Tura, MD, MSc National Institute of Cardiology Laranjeiras Rio de Janeiro Brazil __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] names() function and lmer()
Hello All, I would like to retrieve some of the results from the lmer(...) function in library lme4. If I run a model, say fm.1 - lmer(y ~ 1 + (1 | x), data = dog) and try names(fm.1), I get NULL. Is there anyway to retrieve the information? Thanks __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Multiple tests on 2 way-ANOVA
comments in line Grathwohl, Dominik, LAUSANNE, NRC-BAS wrote: Dear r-helpers, I have a question about multiple testing. Here an example that puzzles me: All matrixes and contrast vectors are presented in treatment contrasts. 1. example: library(multcomp) n-60; sigma-20 # n = sample size per group # sigma standard deviation of the residuals cov1 - matrix(c(3/4,-1/2,-1/2,-1/2,1,0,-1/2,0,1), nrow = 3, ncol=3, byrow=TRUE, dimnames = list(c(A, B, C), c(C.1, C.2, C.3))) # cov1 = variance covariance matrix of the beta coefficients of a # 2x2 factorial design (see Piantadosi 2005, p. 509) cm1 - matrix(c(0, 1, 0, 0, 0, 1), nrow = 2, ncol=3, byrow=TRUE, dimnames = list(c(A, B), c(C.1, C.2, C.3))) # cm1 = contrast matrix for main effects v1 - csimint(estpar=c(100, 6, 5), df=4*n-3, covm=cov1*sigma^2/n, cmatrix=cm1, conf.level=0.95) summary(v1) The adjusted p-values are almost the Bonferroni p-values. If I understood right: You need not to adjust for multiple testing on main effects in a 2x2 factorial design assuming the absence of interaction. SG: Where did you get this idea? A p value of 0.05 says that if the null hypothesis of no effect is true, a result at least as extreme as that observed work actually occur with probability 0.05. Thus, with 2 independent tests, the probability of getting a result at least that extreme in one or both of the tests is 1-(1-0.05)^2 = 0.0975, which is almost 2*0.05. Thus, if I were to consider only main effects in a 2x2 factorial design, this is what I would get from Bonferroni. I do not think that there is a bug, I want to understand, why multcomp does adjust for multiple tests having all information about the design of the trial (variance covariance matrix)? Or do I have to introduce somehow more information? 2. example: And I have second question: How do I proper correct for multiple testing if I want to estimate in the presence of interaction the two average main effects. Can some one point me to some literature where I can learn these things? Here the example, 2x2 factorial with interaction, estimation of average main effects: cov2 - matrix( c(1,-1,-1, 1, -1, 2, 1,-2, -1, 1, 2,-2, 1,-2,-2, 4) , nrow=4, ncol=4, byrow=TRUE) cm2 - matrix(c(0, 1, 0, 1/2, 0, 0, 1, 1/2), nrow = 2, ncol=4, byrow=TRUE, dimnames = list(c(A, B), c(C.1, C.2, C.3, C.4))) v2 - csimint(estpar=c(100, 6, 5, 2), df=4*n-4, covm=cov2*sigma^2/n, cmatrix=cm2, conf.level=0.95) summary(v2) SG: The Bonferroni p-value is the observed times the number of rows in the contrast matrix. The number of columns is irrelevant to Bonferroni. SG: I'm not sure, but I believe that the the adjusted p value would likely be close to (if not exactly) the rank of the contrast matrix; given the rank, it is (I think) independent of the number of rows and columns. SG: These two assertions are consistent with the following example, where I increase the number of dimensions by a factor of 4 without changing the rank. The Bonferroni p value increased by a factor of 4 while the adjusted p value did not change, as predicted. cm2.4 - rbind(cm2, cm2, cm2, cm2) v2.4 - csimint(estpar=c(100, 6, 5, 2), df=4*n-4, covm=cov2*sigma^2/n, cmatrix=cm2.4, conf.level=0.95) summary(v2.4) I do not believe that this is the most efficient way for doing this, since I made already bad experience with the first example. SG: I hope this reply converts the bad experience to good. As for efficiency, you did very well by including such simple but elegant examples. Your post might have more efficiently elicited more and more elegant responses sooner with a more carefully chosen Subject, perhaps like Multiple Comparisons Questions. However, the selection of a possible better subject might rely on information you didn't have. Hope this helps. Spencer Graves My R.version: platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor2.1 year 2005 month12 day 20 svn rev 36812 language R __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] names() function and lmer()
On 7/15/06, A.R. Criswell [EMAIL PROTECTED] wrote: Hello All, I would like to retrieve some of the results from the lmer(...) function in library lme4. If I run a model, say fm.1 - lmer(y ~ 1 + (1 | x), data = dog) and try names(fm.1), I get NULL. Is there anyway to retrieve the information? Yes. The recommended way of retrieving information form a fitted model object like an lmer object is with extractor functions. See ?lmer-class for a list of such methods. That help page also documents the slots in the object. The lmer class is an S4 class and uses typed slots instead of named components. The str function displays the structure of pretty well any type of R object, including S3 classed objects or S4 classed objects. That is my favorite way of checking the structure of an object. Please remember that it is risky to count on being able to reach in to an object and pull out slots or components and operate on them. The names and contents of slots are not guaranteed to stay constant. The lme4 and Matrix packages have been under development for a long time and should continue to be regarded as under development. When we change the internal representation we do change the extractor functions accordingly. It is a bug if an internal change causes an extractor function in the package to fail to return correct results. It is not a bug if an internal change causes your code that assumes a particular, obsolete representation to break. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] names() function and lmer()
lme4 is S4 package.So you should use slotNames to see what slots an object has,and use @ instead of $ to extract the slot you want. library(lme4) Loading required package: Matrix Loading required package: lattice Loading required package: lattice example(lmer) slotNames(fm2) [1] assign frametermsflistZt X [7] ywts wrkres method useScale family [13] call cnames nc Gp XtX ZtZ [19] ZtX Zty Xty OmegaLRZX [25] RXX rZy rXy devComp deviance fixef [31] ranefRZXinv bVar gradComp status [EMAIL PROTECTED] [1] 1.5127787 -40.3739988 -39.1811143 24.5188772 22.9144206 9.2219771 [7] 17.1561300 -7.4517287 0.5786303 34.7680264 -25.7543295 -13.8649853 [13] 4.9159565 20.9290870 3.2586571 -26.4758005 0.9056393 12.4217767 [19] 9.3234692 -8.5991469 -5.3877753 -4.9686374 -3.1939301 -0.3084938 [25] -0.2872083 1.1159883 -10.9059430 8.6275943 1.2806874 6.7563873 [31] -3.0751268 3.5122002 0.8730485 4.9837782 -1.0052909 1.2583989 2006/7/15, A.R. Criswell [EMAIL PROTECTED]: Hello All, I would like to retrieve some of the results from the lmer(...) function in library lme4. If I run a model, say fm.1 - lmer(y ~ 1 + (1 | x), data = dog) and try names(fm.1), I get NULL. Is there anyway to retrieve the information? Thanks __ 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 -- 黄荣贵 Department of Sociology Fudan 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
Re: [R] storing the estimates from lmer
The structure of 'lmer' objects and helper functions is outlined in the 'lmer' and 'lmer-class' help pages. The latter mentions 'vcov 'signature(object = mer)': Calculate variance-covariance matrix of the _fixed_ effect terms, see also 'vcov'. Thus, sqrt(diag(vcov(lmer.object))) should give you standard errors of the fixed effects. The parameter estimates can be obtained using 'VarCorr'. However, characterizing their random variability is harder, because their distribution is not as simply summarized. The 'lmer-class' help page says that an 'lmer' object includes a slot, 'Omega': A list of positive-definite matrices stored as 'dpoMatrix' objects that are the relative precision matrices of the random effects associated with each of the grouping factors. However, I don't know how to use this. For problems like this, the 'lme4' and 'Matrix' packages include a function 'simulate', which is what I might use, at least until I figured out how to get essentially the same answer from the Omega slot. Hope this helps, Spencer Graves prabhu bhaga wrote: Dear all, I'm trying to store/extract the mean standard error of the fixed effects parameter and the variance of the random effects parameter from lmer procedure from mlmre4 package developed by bates n pinheiro. while storing fixed effects parameter is straight forward, the same is not true for storing the variance parameter of the random effects. kindly help me ~prabhu [[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 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] storing the estimates from lmer
p.s. I intended to include the following extension to an example from the 'lmer' help page: fm1 - lmer(Reaction ~ Days + (Days|Subject), sleepstudy) (fm1.fix - fixef(fm1)) (fm1.fix.se - sqrt(diag(vcov(fm1 fm1.fix/fm1.fix.se fm1.ran - VarCorr(fm1) diag(fm1.ran$Subject) The structure of 'lmer' objects and helper functions is outlined in the 'lmer' and 'lmer-class' help pages. The latter mentions 'vcov 'signature(object = mer)': Calculate variance-covariance matrix of the _fixed_ effect terms, see also 'vcov'. Thus, sqrt(diag(vcov(lmer.object))) should give you standard errors of the fixed effects. The parameter estimates can be obtained using 'VarCorr'. However, characterizing their random variability is harder, because their distribution is not as simply summarized. The 'lmer-class' help page says that an 'lmer' object includes a slot, 'Omega': A list of positive-definite matrices stored as 'dpoMatrix' objects that are the relative precision matrices of the random effects associated with each of the grouping factors. However, I don't know how to use this. For problems like this, the 'lme4' and 'Matrix' packages include a function 'simulate', which is what I might use, at least until I figured out how to get essentially the same answer from the Omega slot. Hope this helps, Spencer Graves prabhu bhaga wrote: Dear all, I'm trying to store/extract the mean standard error of the fixed effects parameter and the variance of the random effects parameter from lmer procedure from mlmre4 package developed by bates n pinheiro. while storing fixed effects parameter is straight forward, the same is not true for storing the variance parameter of the random effects. kindly help me ~prabhu [[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 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] storing the estimates from lmer
Hi Spencer, On 7/15/06, Spencer Graves [EMAIL PROTECTED] wrote: p.s. I intended to include the following extension to an example from the 'lmer' help page: fm1 - lmer(Reaction ~ Days + (Days|Subject), sleepstudy) (fm1.fix - fixef(fm1)) (fm1.fix.se - sqrt(diag(vcov(fm1 fm1.fix/fm1.fix.se fm1.ran - VarCorr(fm1) diag(fm1.ran$Subject) I'm confident that you are aware that sqrt(diag(vcov(fm1))) and sqrt(diag(fm1.ran$Subject)) refer to different parameters the model. However, some readers of your reply may not see the distinction. Allow me to try to clarify. The summary for this fitted model is lmer (fm1 - lmer(Reaction ~ Days + (Days | Subject), sleepstudy)) Linear mixed-effects model fit by REML Formula: Reaction ~ Days + (Days | Subject) Data: sleepstudy AIC BIClogLik MLdeviance REMLdeviance 1753.628 1769.593 -871.8141 1751.986 1743.628 Random effects: Groups NameVariance Std.Dev. Corr Subject (Intercept) 612.090 24.7405 Days 35.072 5.9221 0.066 Residual 654.941 25.5918 number of obs: 180, groups: Subject, 18 Fixed effects: Estimate Std. Error t value (Intercept) 251.4051 6.8246 36.838 Days 10.4673 1.5458 6.771 Correlation of Fixed Effects: (Intr) Days -0.138 The fixed effects described in the lower part of this summary are a familiar type of parameter in a statistical model. They are coefficients in a linear predictor. We produce estimates of these parameters and also provide a measure of the precision of these estimates - their standard errors. The vcov generic function returns an estimate of the precision of the estimated parameters (typically parameters that are coefficients in a linear predictor). Thus sqrt(diag(vcov(fm1))) [1] 6.824558 1.545789 provides the standard errors of the fixed effects estimates. The random effects, although they also appear in the linear predictor, are not formally parameters in the model. They are unobserved random variables whose variance covariance matrix has a known form but unknown value. It is the values that determine the variance-covariance matrix that are parameters in the model. These parameters are returned by VarCorr VarCorr(fm1) $Subject 2 x 2 Matrix of class dpoMatrix (Intercept) Days (Intercept) 612.09032 9.60428 Days9.60428 35.07165 attr(,sc) [1] 25.59182 In other words, the 2 x 2 matrix shown above is the estimate of the variance-covariance matrix of the random effects associated with the grouping factor Subject. Thus sqrt(diag(VarCorr(fm1)$Subject)) [1] 24.740459 5.922132 gives the estimated standard deviations of the random effects. These are estimates of parameters in the model. They are *not* standard errors of parameter estimates. The lmer function and related software does not return standard errors of the estimates of variance components. This is intentional. Below I give my familiar rant on why I think returning such standard errors is a bad practice. I encourage users of lmer who wish to determine the precision of the estimates of the variance components to create a Markov chain Monte Carlo sample of the parameters and evaluate the HPDintervals. sm1 - mcmcsamp(fm1, 5) library(coda) Warning message: use of NULL environment is deprecated HPDinterval(sm1) lowerupper (Intercept) 236.6518363 266.5465536 Days 7.0136243 13.8947993 log(sigma^2) 6.25500826.7295329 log(Sbjc.(In))5.49282057.5751372 log(Sbjc.Days)2.81975234.6337518 atanh(Sbj.(I).Dys) -0.69886320.9836688 deviance 1752.5158501 1766.6461469 attr(,Probability) [1] 0.95 rant Some software, notably SAS PROC MIXED, does produce standard errors for the estimates of variances and covariances of random effects. In my opinion this is more harmful than helpful. The only use I can imagine for such standard errors is to form confidence intervals or to evaluate a z-statistic or something like that to be used in a hypothesis test. However, those uses require that the distribution of the parameter estimate be symmetric, or at least approximately symmetric, and we know that the distribution of the estimate of a variance component is more like a scaled chi-squared distribution which is anything but symmetric. It is misleading to attempt to summarize our information about a variance component by giving only the estimate and a standard error of the estimate. /rant The structure of 'lmer' objects and helper functions is outlined in the 'lmer' and 'lmer-class' help pages. The latter mentions 'vcov 'signature(object = mer)': Calculate variance-covariance matrix of the _fixed_ effect terms, see also 'vcov'. Thus, sqrt(diag(vcov(lmer.object))) should give you standard errors of the fixed effects. The parameter estimates can be obtained using
[R] Some problems with latex(ftable)
Dear Users, I've got some problems with LaTeX from the Hmisc package. I write a table to LaTeX from ftable. The table is correctly converted but the headers vanish. What should I do to copy the headers of the table into LaTeX? Thanks, Wilfred __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] intersect() question
Hi, I have a matrix containing ID numbers in each column. I would like to program function which calculate common number of ID numbers between each pair of columns. Suppose: 5 6 7 1 5 3 6 7 2 Then the result should be: 0 2 0 2 0 1 0 1 0 The main problem is how to implement intersect() function to walk through each pair of columns and write result to result matrix. Thanks in advance for any suggestion, Andrej __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] intersect() question
Define a generalized crossproduct and then apply it with the indicated function. Multiply the diagonal elements by zero as the sample output seems to be forcing them that way. mm - matrix(c(5, 1, 6, 6, 5, 7, 7, 3, 2), 3) # test matrix # generalized crossproduct inner - function(a,b=a,f=crossprod) apply(b,2,function(x)apply(a,2,function(y)f(x,y))) inner(mm, f = function(x,y) length(intersect(x,y))) * !diag(ncol(mm)) On 7/15/06, Andrej Kastrin [EMAIL PROTECTED] wrote: Hi, I have a matrix containing ID numbers in each column. I would like to program function which calculate common number of ID numbers between each pair of columns. Suppose: 5 6 7 1 5 3 6 7 2 Then the result should be: 0 2 0 2 0 1 0 1 0 The main problem is how to implement intersect() function to walk through each pair of columns and write result to result matrix. Thanks in advance for any suggestion, Andrej __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] intersect() question
Gabor Grothendieck wrote: Define a generalized crossproduct and then apply it with the indicated function. Multiply the diagonal elements by zero as the sample output seems to be forcing them that way. mm - matrix(c(5, 1, 6, 6, 5, 7, 7, 3, 2), 3) # test matrix # generalized crossproduct inner - function(a,b=a,f=crossprod) apply(b,2,function(x)apply(a,2,function(y)f(x,y))) inner(mm, f = function(x,y) length(intersect(x,y))) * !diag(ncol(mm)) On 7/15/06, Andrej Kastrin [EMAIL PROTECTED] wrote: Hi, I have a matrix containing ID numbers in each column. I would like to program function which calculate common number of ID numbers between each pair of columns. Suppose: 5 6 7 1 5 3 6 7 2 Then the result should be: 0 2 0 2 0 1 0 1 0 The main problem is how to implement intersect() function to walk through each pair of columns and write result to result matrix. Thanks in advance for any suggestion, Andrej __ 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 Thanks for fine solution. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Some problems with latex(ftable)
Bi-Info (http://members.home.nl/bi-info) wrote: Dear Users, I've got some problems with LaTeX from the Hmisc package. I write a table to LaTeX from ftable. The table is correctly converted but the headers vanish. What should I do to copy the headers of the table into LaTeX? Thanks, Wilfred I'm not sure what you mean by ftable, and you did not include the code you are trying to run. I suspect you are trying to latex( ) something other than a list, matrix, or data frame and that you will have to use arguments to latex.default to get the desired result. You may have to write your own latex method. As an example print latex.summary.formula.cross. -- 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
Re: [R] AICc vs AIC for model selection
Regarding AIC.c, have you tried RSiteSearch(AICc) and RSiteSearch(AIC.c)? This produced several comments that looked to me like they might help answer your question. Beyond that, I've never heard of the forecast package, and I got zero hits for RSiteSearch(best.arima), so I can't comment directly on your question. Do you have only one series or multiple? If you have only one, I think it would be hard to justify more than a simple AR(1) model. Almost anything else would likely be overfitting. If you have multiple series, have you considered using 'lme' in the 'nlme' package? Are you familiar with Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer)? If not, I encourage you to spend some quality time with this book. My study of it has been amply rewarded, and I believe yours will likely also. Best Wishes, Spencer Graves Sachin J wrote: Hi, I am using 'best.arima' function from forecast package to obtain point forecast for a time series data set. The documentation says it utilizes AIC value to select best ARIMA model. But in my case the sample size very small - 26 observations (demand data). Is it the right to use AIC value for model selection in this case. Should I use AICc instead of AIC. If so how can I modify best.arima function to change the selection creteria? Any pointers would be of great help. Thanx in advance. Sachin - [[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 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] DTW - dynamic time warping - and time series in R
I found references on Google but not RSiteSearch for dynamic time warping. What do you want to do? Unless you've received a conflicting reply that I haven't seen, it looks like you would have to code it yourself. If you can find something written in another language, you could either translate that into R or like from R to it. Why do you want dynamic time warping? If you want to do research in it, then I suspect you would be best creating your own code. If you want to try it as a solution to some other problem, I suggest you reconsider the other problem: What are your real objectives? If you'd like further help from this listserve, please submit another post. First, however, I suggest you read the posting guide! www.R-project.org/posting-guide.html. There is substantial logic and anecdotal evidence to suggest that posts more consistent with that guide are more likely to get more useful replies quicker. Hope this helps. Spencer Graves Ondrej Vyborny wrote: Hello, can anybody tell me if there exists functions for DTW in R? I didn't find anything at CRAN's search page... Also any information about packages for time series preprocessing (for pattern matching) would be useful... Thanks a lot, ondra __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Some problems with latex(ftable)
The ftable structure is not an ordinary matrix. Instead, it has the
body of the table with several cbind- and rbind-ed rows and columns of
label information. The example in ?ftable has two row factors and two
column factors.
Continuing with the example in ?ftable, enter
tmp - ftable(mtcars$cyl, mtcars$vs, mtcars$am, mtcars$gear, row.vars = c(2,
4),
dnn = c(Cylinders, V/S, Transmission, Gears))
print.default(tmp)
To get what you are looking for, you will need to intercept write.ftable
with, for example,
trace(write.ftable, exit=recover)
then do
3
tmp.latex - latex(t(x))
print.default(tmp.latex)
Now open up t.latex and prepend
\documentstyle{article}
\begin{document}
and append
\end{document}
then latex it.
This gets you close to what you want and you can work with the generated
t.tex file to get the rest of the detail. Or you can work with the
numerous arguments we built into latex (see ?latex) to get some of them
automatically generated.
tmp2.latex - latex(t(x), col.just=rep(c(l,r), c(3,6)),
n.rgroup=c(3,6), file=t2.tex)
Now open up t2.latex and pre- and append the latex controls to it.
This works well for one or two examples. To do many, then you will
need to follow Frank's suggestion and build all of this into a method.
Once the method works well, send it to Frank and he will consider
including it in the next release (Frank, I hope that's a fair offer I
make for you to do).
This raises a question from me to the R developers. I would have
written write.ftable to return t(x), not ox. print is constrained to
return its argument. I thought write had more freedom. Then I would
respecify
print.ftable - function (x, digits = getOption(digits), ...) {
write.ftable(x, quote = FALSE, digits = digits)
invisible(x)
}
Had this been done then the current task would simplify to
latex(write(tmp))
Te question is: why was write.ftable designed to follow the
print.ftable constraint on returned value.
ps. I just designed the method. Take write.ftable, drop the last line
and give it a new name. then you can latex the output of that new
function.
Rich
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] last OSX 3.9 release
I am looking for ready to install binaries of the last, and I assume most stable, version of R for OSX 3.9. I believe this is 2.2.1. I have downloaded the sources, but to not have to time or resources to compile them. -.-. --.- Roger Coppock [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
Re: [R] AICc vs AIC for model selection
Spencer Graves wrote: Regarding AIC.c, have you tried RSiteSearch(AICc) and RSiteSearch(AIC.c)? This produced several comments that looked to me like they might help answer your question. Beyond that, I've never heard of the forecast package, and I got zero hits for RSiteSearch(best.arima), so I can't comment directly on your question. The forecast package is at http://www-personal.buseco.monash.edu.au/~hyndman/Rlibrary/forecast/ Kjetil Do you have only one series or multiple? If you have only one, I think it would be hard to justify more than a simple AR(1) model. Almost anything else would likely be overfitting. If you have multiple series, have you considered using 'lme' in the 'nlme' package? Are you familiar with Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer)? If not, I encourage you to spend some quality time with this book. My study of it has been amply rewarded, and I believe yours will likely also. Best Wishes, Spencer Graves Sachin J wrote: Hi, I am using 'best.arima' function from forecast package to obtain point forecast for a time series data set. The documentation says it utilizes AIC value to select best ARIMA model. But in my case the sample size very small - 26 observations (demand data). Is it the right to use AIC value for model selection in this case. Should I use AICc instead of AIC. If so how can I modify best.arima function to change the selection creteria? Any pointers would be of great help. Thanx in advance. Sachin - [[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 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] AICc vs AIC for model selection
Hi, Kjetil: Thanks. Spencer Graves Kjetil Brinchmann Halvorsen wrote: Spencer Graves wrote: Regarding AIC.c, have you tried RSiteSearch(AICc) and RSiteSearch(AIC.c)? This produced several comments that looked to me like they might help answer your question. Beyond that, I've never heard of the forecast package, and I got zero hits for RSiteSearch(best.arima), so I can't comment directly on your question. The forecast package is at http://www-personal.buseco.monash.edu.au/~hyndman/Rlibrary/forecast/ Kjetil Do you have only one series or multiple? If you have only one, I think it would be hard to justify more than a simple AR(1) model. Almost anything else would likely be overfitting. If you have multiple series, have you considered using 'lme' in the 'nlme' package? Are you familiar with Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer)? If not, I encourage you to spend some quality time with this book. My study of it has been amply rewarded, and I believe yours will likely also. Best Wishes, Spencer Graves Sachin J wrote: Hi, I am using 'best.arima' function from forecast package to obtain point forecast for a time series data set. The documentation says it utilizes AIC value to select best ARIMA model. But in my case the sample size very small - 26 observations (demand data). Is it the right to use AIC value for model selection in this case. Should I use AICc instead of AIC. If so how can I modify best.arima function to change the selection creteria? Any pointers would be of great help. Thanx in advance. Sachin - [[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 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] put R on a web server
Dear R People: Has anyone put R on a web server any time, recently, please? (Red Hat Linux) The University of Montana put a version up in 2003, but I was wondering if anyone had done so, please? Also, where would I find information on such an installation, please? thanks, Sincerely, Erin Hodgess Associate Professor Department of Computer and Mathematical Sciences University of Houston - Downtown mailto: [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] princomp and eigen
Consider the following output [R2.2.0; Windows XP] set.seed(160706) X - matrix(rnorm(40),nrow=10,ncol=4) Xpc - princomp(X,cor=FALSE) summary(Xpc,loadings=TRUE, cutoff=0) Importance of components: Comp.1Comp.2Comp.3 Comp.4 Standard deviation 1.2268300 0.9690865 0.7918504 0.55295970 Proportion of Variance 0.4456907 0.2780929 0.1856740 0.09054235 Cumulative Proportion 0.4456907 0.7237836 0.9094576 1. Loadings: Comp.1 Comp.2 Comp.3 Comp.4 [1,] -0.405 -0.624 0.466 0.479 [2,] -0.199 -0.636 -0.346 -0.660 [3,] 0.884 -0.443 0.023 0.148 [4,] 0.122 0.099 0.814 -0.559 eigen(var(X)) $values [1] 1.6723465 1.0434763 0.6966967 0.3397382 $vectors [,1] [,2][,3] [,4] [1,] -0.4048158 0.6240510 0.46563382 0.4794473 [2,] -0.1994853 0.6361009 -0.34634256 -0.6600213 [3,] 0.8839775 0.4429553 0.02261302 0.1478618 [4,] 0.1221215 -0.0986234 0.81407655 -0.5591414 I would have expected the princomp component standard deviations to be the square roots of the eigen() $values and they clearly are not. 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
Re: [R] Prediction interval of Y using BMA
From what I've heard, the primary reason for using Bayesian Model Averaging (BMA) is to get better predictions, both better point estimates and better estimates of the uncertainties. I could not find a function to compute that. However, if you've got point estimates of predictions, I can give you a formula for the variance of the prediction error: var(Y|x) = var(E(Y|x,i)over i)+E(var(Y|x,i)over i), where Y is the response variable, and 'i' indicates which model. This is a companion to the formula for the predictions: E(Y|x) = E(E(Y|x,i)over i). If you've got a function to compute E(Y|x,i) and compute their expectation over i, it should not be too difficult to modify that function to also compute var(Y|x,i) and then to compute both E(var(Y|x,i)over i) and var(E(Y|x,i)over i). Then if distributions are roughly normal, you've got what you need. Does this make sense? Hope this helps. Spencer Graves p.s. I've copied the maintainer for the BMA package on this email. I hope he will correct any misstatement in these comments and update us on any relevant capabilities we may have missed. [EMAIL PROTECTED] wrote: Hello everybody, In order to predict income for different time points, I fitted a linear model with polynomial effects using BMA (bicreg(...)). It works fine, the results are consistent with what we are looking for. Now, we would like to predict income for a future time point t_next and of course draw the prediction interval around the estimated value for this point t_next. I've found the formulae for the ponctual estimation of t_next but I cannot succeed in finding the formula to compute the prediction interval based on BMA sets of models, which requires the adequate variance. The paper of Hoeting et al. on BMA gives the posterior variance of a parameter but what I want is the posterior variance of a predicted Y and not the posterior variance of a beta, since the goal here is not to build an explanatory model but an accurately predictive model. Is there a way to get around it with the function of the BMA package or does anybody has indication on where to find the formulae ? I searched the web for hours but nothing like this is to be found. Any help will be appreciated. Thank's in advance, Nicolas Meyer Département de Santé Publique Unité de Biostatistique et Méthodologie 03 88 11 63 58 [EMAIL PROTECTED] mailto:[EMAIL PROTECTED] La conjecture de Feynman : To report a significant result and reject the null in favor of an alternative hypothesis is meaningless unless the alternative hypothesis has been stated before the data was obtained The meaning of it all, 1998. [[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 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] princomp and eigen
Murray Jorgensen wrote: set.seed(160706) X - matrix(rnorm(40),nrow=10,ncol=4) Xpc - princomp(X,cor=FALSE) summary(Xpc,loadings=TRUE, cutoff=0) Importance of components: Comp.1Comp.2Comp.3 Comp.4 Standard deviation 1.2268300 0.9690865 0.7918504 0.55295970 [...] I would have expected the princomp component standard deviations to be the square roots of the eigen() $values and they clearly are not. It's an 1/n vs. 1/(n-1) thing: eX - eigen(var(X)) sqrt(eX$values) [1] 1.2931924 1.0215069 0.8346836 0.5828707 sqrt(9/10 * eX$values) [1] 1.2268300 0.9690865 0.7918504 0.5529597 -- Bjørn-Helge Mevik __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Some problems with latex(ftable)
Richard M. Heiberger wrote:
The ftable structure is not an ordinary matrix. Instead, it has the
body of the table with several cbind- and rbind-ed rows and columns of
label information. The example in ?ftable has two row factors and two
column factors.
Continuing with the example in ?ftable, enter
tmp - ftable(mtcars$cyl, mtcars$vs, mtcars$am, mtcars$gear, row.vars = c(2,
4),
dnn = c(Cylinders, V/S, Transmission, Gears))
print.default(tmp)
To get what you are looking for, you will need to intercept write.ftable
with, for example,
trace(write.ftable, exit=recover)
then do
3
tmp.latex - latex(t(x))
print.default(tmp.latex)
Now open up t.latex and prepend
\documentstyle{article}
\begin{document}
and append
\end{document}
then latex it.
This gets you close to what you want and you can work with the generated
t.tex file to get the rest of the detail. Or you can work with the
numerous arguments we built into latex (see ?latex) to get some of them
automatically generated.
tmp2.latex - latex(t(x), col.just=rep(c(l,r), c(3,6)),
n.rgroup=c(3,6), file=t2.tex)
Now open up t2.latex and pre- and append the latex controls to it.
This works well for one or two examples. To do many, then you will
need to follow Frank's suggestion and build all of this into a method.
Once the method works well, send it to Frank and he will consider
including it in the next release (Frank, I hope that's a fair offer I
make for you to do).
Yes, working with Charles Thomas Dupont we'll be glad to add such a
contribution to latex. Thanks for your excellent ideas above Rich.
-Frank
This raises a question from me to the R developers. I would have
written write.ftable to return t(x), not ox. print is constrained to
return its argument. I thought write had more freedom. Then I would
respecify
print.ftable - function (x, digits = getOption(digits), ...) {
write.ftable(x, quote = FALSE, digits = digits)
invisible(x)
}
Had this been done then the current task would simplify to
latex(write(tmp))
Te question is: why was write.ftable designed to follow the
print.ftable constraint on returned value.
ps. I just designed the method. Take write.ftable, drop the last line
and give it a new name. then you can latex the output of that new
function.
Rich
__
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
