[R] principal component analysis
Hello, I have used "prcomp" and the variances for the first 3 PC's are 2.65, 1.97 and 0.38. When I plot the principal component values for each data point I can see that the points lie in a plane as one might expect from the variances. But this plane is diagonal through the 3D space of the first 3 PC's. My question is: is it usual to follow a PC analysis with an eigen analysis to align this plane with with 2 principal directions, or is there a way of incorporating this into "prcomp". Thanks, David __ 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] thickness of tick marks
Hi, I am unable to get the tick marks to appear thicker in plot. I have tried things like par(lw=2) but this only seems to affect other line thicknesses. The use of axes directly fixes the problem because lw = 2 applies to both the axis and the ticks. Is there is way of feeding a single parameter to plot or setting a par parameter to do this? I am using R 1.9.0 on a windows 2000 platform. Thanks, David __ [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
Re: [R] Getting Pr from Summary(lm)
Thanks for all the relies. I recently discovered "names" and applied it to "lm" objects but did not think to apply it to the "summary" object. Cheers, David __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Specifying xlevels in effects library
Andy and John,
I looked at typical when xlevels did not work but when I saw that it was
a function I went no further. Setting the function to a constant was a
good idea.
John's method seems to require that I change the model:
> eff <-effect("sex*age",mod,xlevel=(Age=c(120,120)))
Error in all(!(factors[, term1] & (!factors[, term2]))) :
subscript out of bounds
In addition: Warning message:
sex:age does not appear in the model in: effect("sex*age", mod, xlevel =
(age = c(120, 120)))
Andy's method works as suggested for this simple case.
Thanks for your time.
Cheers, David
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Getting Pr from Summary(lm)
Hello, I am trying to get the P values from the output of a summary for lm. lm <- lm(y ~ age + sex) s <- summary(lm) I thought that I might be able to get them using a combination of scan, grep and sub. But I got stuck on the first step - being able to process "s" as a text string. I could perhaps write it to file than scan it back but there is probably an easier way to do the whole thing. Help would be welcome, David __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Specifying xlevels in effects library
library(effects)
mod <- lm(Measurement ~ Age + Sex, data=d)
e <-effect("Sex",mod)
The effect is evaluated at the mean age.
> e
Sex effect
Sex
FM
43.33083 44.48531
>
> e$model.matrix
(Intercept) Age SexM
11 130.58590
23 1 130.58591
To evaluate the effect at Age=120 I tried:
e <-effect("Sex",mod,xlevels=list(Age=c(120)))
but the effect was still evaluated at 130.5859.
Is this an incorrect usage of xlevels?
Thanks, David
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] loess prediction limits
Hello I am plotting a loess curve with confidence limits as below. How do I create the prediction limits? Is multiplying the standard errors by sqrt(n) appropriate? data <- mndata lo <- loess(data[[variableName]] ~ Age, data, span=1.0, control = loess.control(surface = "direct")) xPoints <- data.frame(age = seq(1,240,1)) lo1 <- predict(lo, xPoints, se = TRUE) age <- xPoints$age lines(age,lo1$fit, col=4) # now do +/- 2 std errors lo1p <- lo1$fit + 2*lo1$se.fit lo1m <- lo1$fit - 2*lo1$se.fit lines(age,lo1p, col=4) lines(age,lo1m, col=4) Thanks, David __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Getting the same values of adjusted mean and standard errors as SAS
Thanks for the help. Both the "Design" package and the "effects" package look as though they are what I need although it will probably take me a while to get on top of both. I have had a brief go at the Design package and the contrast function is particularly useful. A question on the Design package: There are 5 types for factor "group", one is the reference - call it "a". f <- ols(y ~ age + sex + group, data=dd) contrast(f, list(group='a'), list(group='b')) I can use this to contrast pairs but can I use this to contrast b against c,d, and e as a group. Also "a" against the rest? Thanks, David Frank E Harrell Jr wrote: On Thu, 27 May 2004 16:34:58 +0930 "David J. Netherway" <[EMAIL PROTECTED]> wrote: Hello, I am trying to get the same values for the adjusted means and standard errors using R that are given in SAS for the following data. The model is Measurement ~ Age + Gender + Group. I can get the adusted means at the mean age by using predict. I do not know how to get the appropriate standard errors at the adjusted means for Gender using values from predict. So I attempted to get them directly from the residuals as follows. The data is at the end of the email. While there is a match for the males there is a large difference for the females indicating that what I am doing is wrong. # meanAge <- mean(dd$Age) meanAgeM <- mean(dd$Age[d$Gender=="M"]) meanAgeF <- mean(dd$Age[d$Gender=="F"]) . . . . By using sex-specific means of age you are not getting adjusted estimates in the usual sense. I prefer to think of effects as differences in predicted values rather than as complex SAS-like contrasts. The Design package's contrast function makes this easy (including SEs and confidence limits): library(Design) # also requires Hmisc d <- datadist(dd); options(datadist='d') f <- ols(y ~ age + sex + group, data=dd) contrast(f, list(sex='M'), list(sex='F')) # usual adjusted difference M vs F contrast(f, list(sex='M',age=mean(dd$age[dd$sex=='M']), list(sex='F',age=mean(dd$age[dd$sex=='F')) # M vs F not holding age constant You can also experiment with specifying age=tapply(age, sex, mean, na.rm=TRUE) using some of the contrast.Design options. --- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt University __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Getting the same values of adjusted mean and standard errors as SAS
Hello,
I am trying to get the same values for the adjusted means and standard
errors using R that are given in SAS for the
following data. The model is Measurement ~ Age + Gender + Group. I can
get the adusted means at the mean age
by using predict. I do not know how to get the appropriate standard
errors at the adjusted means for Gender
using values from predict. So I attempted to get them directly from the
residuals as follows. The data is at the end
of the email. While there is a match for the males there is a large
difference for the females indicating that what I am doing is wrong.
#
meanAge <- mean(dd$Age)
meanAgeM <- mean(dd$Age[d$Gender=="M"])
meanAgeF <- mean(dd$Age[d$Gender=="F"])
# determine adjusted means for the males and females at meanAge using
predict
# set up data frame to get predicted values at meanAge
evalDF <- data.frame(Age = meanAge, Gender = c("F","M"),
Group = c("1NC","1NC", "2UCLP", "2UCLP", "3BCLP", "3BCLP",
"4ICP", "4ICP", "5CLPP", "5CLPP"))
mod <-lm(Measurement ~ Age + Gender + Group, data=dd)
pred <- predict(mod, evalDF, se.fit = TRUE)
adjDF <- data.frame(evalDF,fit = pred$fit, se = pred$se.fit)
adjMeanMale <- mean(adjDF$fit[adjDF$Gender=="M"]); # match: 3.889965 cf
3.88996483
adjMeanFemale <- mean(adjDF$fit[adjDF$Gender=="F"]); # match: 3.9 cf
3.9036
# Try to get standard errors at the adjusted means for Gender as follows:
ddr <- data.frame(dd, res = residuals(mod))
nM <- summary(ddr$Gender)[["M"]]
seRegM <- sqrt(mean( ddr$res[ddr$Gender=="M"]**2 ))
sxxM <- sum((dd$Age[d$Gender=="M"]-meanAgeM)**2)
syM <- seRegM * sqrt(1/nM + (meanAge-meanAgeM)**2/sxxM); #0.1103335 cf
0.11032602 - matches to 5 decimal places
nF <- summary(ddr$Gender)[["F"]]
seRegF <- sqrt(mean( ddr$res[ddr$Gender=="F"]**2 ))
sxxF <- sum((dd$Age[d$Gender=="F"]-meanAgeF)**2)
syF <- seRegF * sqrt(1/nF + (meanAge-meanAgeF)**2/sxxF); # wrong:
0.07279221 cf 0.14256466
> dd
Measurement Age Gender Group
1 3.8 94 M 3BCLP
2 2.7 88 F 3BCLP
3 3.0 155 M 3BCLP
4 2.7 33 M 3BCLP
5 4.6 109 M 5CLPP
6 5.1 325 M 5CLPP
7 3.9 79 M 5CLPP
8 4.2 126 M 5CLPP
9 3.9 77 F 5CLPP
10 4.0 61 F 5CLPP
11 3.6 49 F 5CLPP
12 3.7 14 F 4ICP
13 4.2 160 F 4ICP
14 3.9 60 M 4ICP
15 5.0 61 M 4ICP
16 3.9 222 F 4ICP
17 3.8 82 F 4ICP
18 4.8 340 F 4ICP
19 3.2 206 M 4ICP
20 3.8 19 M 1NC
21 4.9 166 M 1NC
22 3.8 93 M 1NC
23 3.6 142 M 1NC
24 4.8 241 M 1NC
25 3.9 81 M 1NC
26 4.5 41 M 1NC
27 5.1 244 F 1NC
28 4.6 100 M 1NC
29 5.1 122 F 1NC
30 4.7 194 F 1NC
31 5.1 297 M 1NC
32 3.9 69 M 2UCLP
33 2.5 141 M 2UCLP
34 3.2 104 M 2UCLP
35 3.8 90 M 2UCLP
36 3.8 92 M 2UCLP
37 3.6 149 F 2UCLP
38 3.8 53 F 2UCLP
39 4.7 111 M 2UCLP
40 3.8 116 F 2UCLP
41 3.3 81 M 2UCLP
>
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Setting up contrasts
I am using the following model: lm <- lm(mydata[[variableName]] ~ Age + Gender + Group, data=mydata) There are 5 groups in "Group": nonc (the control), c1,c2,c3 and c4. How do I contrast nonc vs the others? and How do I contrast c1 vs other c's (ie c2,c3,c4 as a subgroup)? I have looked at the contrasts option in lm and model.matrix and am really none the wiser. Though it looks like I need to read up on model matrices again. Any help would be appreciated. Thanks, David __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
