similar to: Balanced incomplete block analysis

Bother! This cold has made me accident-prone. I meant to hit Reply-all. Clarification below. -------- Original Message -------- Subject: Re: [R] Parameterization puzzle Date: Fri, 21 Jul 2006 19:10:03 +1200 From: Murray Jorgensen <maj at waikato.ac.nz> To: Prof Brian Ripley <ripley at stats.ox.ac.uk> References: <44C063E5.3020703 at waikato.ac.nz>

If rb.lm is an lm-object, I can access the error mean square as s2 <- sum(rb.lm$residuals^2)/rb.lm$df.residual This seems a bit like hard work for such a commonly wanted quantity. Is there a better way to do this? Murray Jorgensen Dr Murray Jorgensen http://www.stats.waikato.ac.nz/Staff/maj.html Department of Statistics, University of Waikato, Hamilton, New Zealand Email: maj at

Does anyone understand this? I have a long vector parlist containing values of 138 parameters in 100 simulations. The factor parno is defined by trys <- 100 sim <- gl(trys,138) parno <- gl(138,1,trys*138) I want to extract and compare results for related groups of parameters but I get inconsistent results: > length(parlist[parno==37:42]) [1] 600 >

Has anyone got a cute way of compiling a table of the available plotting symbols? I've been repeatedly pasting in cno <- cno + 1 cno plot(x,y,pch=cno) which is tedious, and besides I have to note down the symbols by hand. Dr Murray Jorgensen http://www.stats.waikato.ac.nz/Staff/maj.html Department of Statistics, University of Waikato, Hamilton, New Zealand Email: maj at

G'day list! I'm thinking about programming likelihood and score intervals for generalized linear models in R based on the paper "On the computation of likelihood ratio and score test based confidence intervals in generalized linear models" by Juha Alho (1992) (Statistics in Medicine, 11, 923-930). Being lazy, I thought that I would ask if anyone else on the list has

I can ask my question using and example from Chapter 1 of Pinheiro & Bates. > # 1.4 An Analysis of Covariance Model > > OrthoFem <- Orthodont[ Orthodont$Sex == "Female", ] > fm1OrthF <- + lme( distance ~ age, data = OrthoFem, random = ~ 1 | Subject ) > summary( fm1OrthF ) Linear mixed-effects model fit by REML Data: OrthoFem AIC BIC

Having just update to R 2.6.2 on my old Windows laptop I notice that the number of packages is growing exponentially and my usual approach of get-em-all may not be viable much longer. Has any thought been given to dividing "contributed" binaries into a recommended set, perhaps a couple of hundred, and the remained. That way one could install the recommended ones routinely and add in

I forgot to mention that I am using Windows XP. -------- Original Message -------- Subject: Pasting data into scan() Date: Tue, 02 May 2006 11:55:03 +1200 From: Murray Jorgensen <maj at stats.waikato.ac.nz> To: r-help at stat.math.ethz.ch The file TENSILE.DAT from the Hand et al "Handbook of Small Data Sets" looks like this: [...] -- Dr Murray Jorgensen

How do you stop R from putting up a dialog box when you quit Rgui? (I use Windows and I never save workspaces that way) Murray -- Dr Murray Jorgensen http://www.stats.waikato.ac.nz/Staff/maj.html Department of Statistics, University of Waikato, Hamilton, New Zealand Email: maj at waikato.ac.nz Fax 7 838 4155 Phone +64 7 838 4773 wk +64 7 849 6486 home

It is not uncommon to want to coarsen a factor by grouping levels together. I have found one way to do this in R: > sites [1] F A A D A A B F C F A D E E D C F A E D F C E D E F F D B C Levels: A B C D E F > regions <- list(I = c("A","B","C"), II = "D", III = c("E","F")) > library(Epi) > region <-

In one saved workspace I have the results of a simulation experiment stored as an array "resarray". > dim(resarray) [1] 10 6 500 3 In another workspace I have a similar array from another run of the simulation. I want to combine the two arrays into a single array of dimensions 10, 6, 1000, 3 What's the best way to do this? Murray Jorgensen Dr Murray Jorgensen

I wanted a sort of matrix product of an array and a matrix. As there does not seem to be any array multiplication apart from outer() I proceeded as follows: lambda <- array(0, c(n,m,d)) # stuff omitted # zed is an n by m matrix # # \lamb.star_{ik} lamb.star <- matrix(0, nrow=n, ncol=d) for (i in 1:n) { for (k in 1:d) { for (j in 1:m) { lamb.star[i,k] = lamb.star[i,k] +

I assigned a class the first problem in Pinheiro & Bates, which uses the data set PBIB from the SASmixed package. I have recently downloaded 2.0.1 and its associated packages. On trying library(SASmixed) data(PBIB) library(nlme) plot(PBIB) I get a warning message Warning message: replacing previous import: coef in: namespaceImportFrom(self, asNamespace(ns)) after library(nlme) and a

The birth weight example from ?stepAIC in package MASS runs well as indeed it should. However when I change stepAIC() calls to step() calls I get warning messages that I don't understand, although the output is similar. Warning messages: 1: In model.response(m, "numeric") : using type="numeric" with a factor response will be ignored (and three more the same.) Checked

I got a shock a few days ago when I accidentally visited www.r-project.com . I thought that the r-project site had been hacked until I realised my mistake. There is also a site www.r-project.net. Both of these sites appear to be Japanese. Does anyone know anything about them? I suppose that it is not unusual for names close to those of popular sites to be used. It is good that they use a

I am exploring the result of clustering a large multivariate data set into a number of groups, represented, say, by a factor G. I wrote a function to see how categorical variables vary between groups: > ddisp <- function(dvar) { + csqt <- chisq.test(G,dvar) + print(csqt$statistic) + print(csqt$observed) + print(round(csqt$expected)) + round(csqt$residuals) + } > > x

I am wondering if stepAIC in the MASS library may be used for model selection in an overdispersed Poisson situation. What I thought of doing was to get an estimate of the overdispersion parameter phi from fitting a model with all or most of the available predictors (we have a large number of observations so this should not be problematical) and then use stepAIC with scale = phi. Should this

I had an error message while running a macro from Yudi Pawitan's web site: > source("ex2-13.r") Error in parse(file, n, text, prompt) : syntax error on line 2 Inspecting ex2-13.r I found that the error was generated by another source() command. Clearly R does not like nested source()s, which is fair enough when you think about it. Still it's something that you might want

We have a choice when calculating the Huber location estimate: > set.seed(221205) > y <- 7 + 3*rt(30,1) > library(MASS) > huber(y)$mu [1] 5.9117 > coefficients(rlm(y~1)) (Intercept) 5.9204 I was surprised to get two different results. The function huber() works directly with the definition whereas rlm() uses iteratively reweighted least squares. My surprise is

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.1 Comp.2 Comp.3 Comp.4 Standard deviation 1.2268300 0.9690865 0.7918504 0.55295970 Proportion of Variance 0.4456907 0.2780929