Re: [R] QQ plotting of various distributions...

2009-09-25 Thread Sunil Suchindran
#same shape

some_data - rgamma(500,shape=6,scale=2)
test_data - rgamma(500,shape=6,scale=2)
plot(sort(some_data),sort(test_data))
# You can also use qqplot(some_data,test_data)
abline(0,1)

# different shape

some_data - rgamma(500,shape=6,scale=2)
test_data - rgamma(500,shape=4,scale=2)
plot(sort(some_data),sort(test_data))
abline(0,1)

It is helpful to assess the sampling variability, by
creating repeated sets of test_data, and plotting
all of these along with your observations to create
a confidence envelope.

The SuppDists provides Inverse Gauss.


On Thu, Sep 17, 2009 at 11:46 AM, Petar Milin pmi...@ff.uns.ac.rs wrote:

 Hello!
 I am trying with this question again:
 I would like to test few distributional assumptions for some behavioral
 response data. There are few theories about true distribution of those data,
 like: normal, lognormal, gamma, ex-Gaussian (exponential-Gaussian), Wald
 (inverse Gaussian) etc. The best way would be via qq-plot, to show to
 students differences. First two are trivial:
 qqnorm(dat$X)
 qqnorm(log(dat$X))
 Then, things are getting more hairy. I am not sure how to make plots for
 the rest. I tried gamma with:
 qqmath(~ X, data=dat, distribution=function(X)
   qgamma(X, shape, scale))
 Which should be the same as:
 plot(qgamma(ppoints(dat$X), shape, scale), sort(dat$X))
 Shape and scale parameters I got via mhsmm package that has gammafit() for
 shape and scale parameters estimation.
 Am I on right track? Does anyone know how to plot the rest: ex-Gaussian
 (exponential-Gaussian), Wald (inverse Gaussian)?

 Thanks,
 PM

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


[[alternative HTML version deleted]]

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


Re: [R] plotting least-squares regression against x-axis

2009-09-19 Thread Sunil Suchindran
x - seq(50)

y - 10 + x * 2 + rnorm(50,0,10)

plot(y~x)

mylm = lm(y~x)

# Use str(mylm) to see how to get the residuals

plot(x,mylm$residuals)


On Sat, Sep 19, 2009 at 8:35 PM, Jason Priem pr...@email.unc.edu wrote:

 Hi,
 I want to plot the residuals of a least-squares regression.

 plot(lm(y~x), which=1)

 does this, but it plots the y-axis of my data on the x-axis of the
 residuals plot.  That is, it plots the residual for each y-value in the
 data.  Can I instead use the x-axis of my data as the x-axis of the
 residuals plot, showing the residual for a given x?

 Thanks!

 Jason Priem
 University of North Carolina at Chapel Hill
 School of Information and Library Science

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


[[alternative HTML version deleted]]

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


Re: [R] What does model.matrix() return?

2009-09-17 Thread Sunil Suchindran
This describes the way in which categorical variables are coded in
the model. In this case, we see treatment coding,
although this can go by different names in textbooks (for example, reference
cell coding).

attr(,assign)
[1] 0 1 2 2

means that the first column of the model matrix corresponds to the
intercept,
the second column corresponds to the first variable A,
and the third and fourth columns correspond to the second variable B.

A has 2 levels so it can be represented by 1 column.
Similarly, B has 3 levels so it is represented by 2 columns.
For variable A, the regression coefficient represents an estimated contrast
(difference)
between the two levels when B is held constant.

On Thu, Sep 17, 2009 at 11:13 AM, Peng Yu pengyu...@gmail.com wrote:

 Hi,

 I don't understand what the meaning of the following lines returned by
 model.matrix(). Can somebody help me understand it? What can they be
 used for?

 attr(,assign)
 [1] 0 1 2 2
 attr(,contrasts)
 attr(,contrasts)$A
 [1] contr.treatment

 attr(,contrasts)$B
 [1] contr.treatment

 Regards,
 Peng

  a=2
  b=3
  n=4
  A = rep(sapply(1:a,function(x){rep(x,n)}),b)
  B = as.vector(sapply(sapply(1:b, function(x){rep(x,n)}),
 function(x){rep(x,a)}))
  Y = A + B + rnorm(a*b*n)
  fr = data.frame(Y=Y,A=as.factor(A),B=as.factor(B))
  afit=aov(Y ~ A + B,fr)
  model.matrix(afit)
   (Intercept) A2 B2 B3
 11  0  0  0
 21  0  0  0
 31  0  0  0
 41  0  0  0
 51  1  0  0
 61  1  0  0
 71  1  0  0
 81  1  0  0
 91  0  1  0
 10   1  0  1  0
 11   1  0  1  0
 12   1  0  1  0
 13   1  1  1  0
 14   1  1  1  0
 15   1  1  1  0
 16   1  1  1  0
 17   1  0  0  1
 18   1  0  0  1
 19   1  0  0  1
 20   1  0  0  1
 21   1  1  0  1
 22   1  1  0  1
 23   1  1  0  1
 24   1  1  0  1
 attr(,assign)
 [1] 0 1 2 2
 attr(,contrasts)
 attr(,contrasts)$A
 [1] contr.treatment

 attr(,contrasts)$B
 [1] contr.treatment

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


[[alternative HTML version deleted]]

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


Re: [R] lean text label below barplot table

2009-09-10 Thread Sunil Suchindran
#Example data

df - data.frame(trt = factor(c(A long label, Another long \n label)),
outcome = c(1,4))

#Install ggplot2 if needed

library(ggplot2)
p - ggplot(df, aes(y=outcome, x=trt))
p - p + geom_bar(position=dodge, stat=identity)
p - p + opts(axis.text.x = theme_text(angle = 45, hjust=1))
p



On Mon, Sep 7, 2009 at 5:03 PM, Xiaogang Yang gavinxy...@gmail.com wrote:

 Hi, everyone:
 I am plotting an graph with bar plot, but the label after every bar is too
 long, I wanna if I can draw the label lean to an angle
 thanks

 --
 Xiaogang Yang
 Sensorweb Research Laboratory
 http://sensorweb.vancouver.wsu.edu/
 Washington State University Vancouver

[[alternative HTML version deleted]]

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


[[alternative HTML version deleted]]

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


Re: [R] standard error associated with correlation coefficient

2009-08-30 Thread Sunil Suchindran
Here are some options for confidence intervals.
#by hand

sample_r - .5
n_sample - 100
df - n_sample-1
fisher_z - 0.5*(log(1+sample_r)-log(1-sample_r))
se_z - 1/sqrt(n_sample-3)
t_crit - qt(.975,df ,lower.tail=TRUE)
z_lci - fisher_z - (t_crit * se_z)
z_uci - fisher_z + (t_crit * se_z)

(r_lci - tanh(z_lci))
(r_uci - tanh(z_uci))
# Compare to the CIr function in psychometric library
#that uses standard normal instead of t

library(psychometric)
CIr(sample_r,n_sample)


And see

http://www.stat.umn.edu/geyer/5601/examp/bse.html

for bootstrap options.

On Thu, Aug 27, 2009 at 3:03 PM, Donald Braman dbra...@law.gwu.edu wrote:

 I want the standard error associated with a correlation.  I can calculate
 using cor  var, but am wondering if there are libraries that already
 provide this function.

[[alternative HTML version deleted]]

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


[[alternative HTML version deleted]]

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


Re: [R] Select values at random by id value

2009-07-01 Thread Sunil Suchindran
#Highlight the text below (without the header)
# read the data in from clipboard

df - do.call(data.frame, scan(clipboard, what=list(id=0,
date=,loctype=0 ,haptype=0)))

# split the data by date, sample 1 observation from each split, and rbind

sampled_df - do.call(rbind, lapply(split(df,
df$date),function(x)x[sample(1:nrow(x), 1),]))

On Mon, Jun 29, 2009 at 9:11 AM, James Martin just.strut...@gmail.comwrote:

 All,

 I have data that looks like below.  For each id there may be more than one
 value per day.  I want to select a random value for that day for that id.
 The end result would hopefully be a matrix with the id as rows, date as
 columns and populated by the random hab value.  Thanks to someone on here
 (Jim) I know how to do the matrix, but now realize I need to randomly
 select
 some of my values.  All help is appreciated. jm
 id, date, loctype, habtype
  50022 1/25/2006 0 6  50022 1/31/2006 0 6  50022 2/8/2006 0 6  50022
 2/13/2006 0 6  50022 2/15/2006 0 6  50022 2/24/2006 0 6  50022 3/2/2006 0 6
 50022 3/9/2006 0 6  50022 3/16/2006 0 6  50022 3/24/2006 0 6  50022
 4/9/2006
 0 3  50022 4/18/2006 0 6  50022 4/27/2006 0 3  50022 5/23/2006 1 3  50022
 5/23/2006 1 6  50022 5/24/2006 1 3  50022 5/24/2006 1 3  50022 5/24/2006 1
 3
 50022 5/24/2006 1 3  50022 5/25/2006 1 3  50022 5/25/2006 1 3  50022
 5/25/2006 1 3  50022 5/26/2006 1 5  50022 5/26/2006 1 3  50022 5/27/2006 1
 5
 50022 5/27/2006 1 3  50022 5/28/2006 1 3  50022 5/29/2006 1 3  50022
 5/30/2006 1 5  50022 5/30/2006 1 3  50022 5/31/2006 1 3  50022 5/31/2006 1
 3
 50022 6/1/2006 1 3  50022 6/2/2006 1 3  50022 6/3/2006 1 3  50022 6/4/2006
 1
 3  50022 6/5/2006 1 3  50022 6/6/2006 1 5  50022 6/6/2006 1 5  50022
 6/6/2006 1 5  50022 6/6/2006 1 3  50022 6/6/2006 1 3  50022 6/7/2006 1 5
 50022 6/7/2006 1 3  50022 6/7/2006 1 3  50022 6/7/2006 1 3  50022 6/7/2006
 1
 3  50022 6/8/2006 1 3  50022 6/8/2006 1 3  50022 6/8/2006 1 3  50022
 6/9/2006 1 5  50022 6/10/2006 1 3  50022 6/11/2006 1 5

[[alternative HTML version deleted]]

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


[[alternative HTML version deleted]]

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


Re: [R] text() to label points in ggplot

2009-05-19 Thread Sunil Suchindran
# Here are two options:

p - ggplot(mtcars, aes(wt, mpg)) + geom_point() + geom_text(aes(x = 5, y =
30, label = A Label))

#or

response - c(2,4)
xvar - c(1,2)
label - response;
myData - data.frame(response,xvar,label)
p - ggplot(myData, aes(y=response, x=xvar))
p + geom_bar(position=dodge, stat=identity) +
geom_text(aes(y=response+1,label=label))

On Thu, May 14, 2009 at 1:22 PM, stephenb sten...@go.com wrote:


 is there a way to label points in a graph using text(locator(1),text)
 after ggplot() or qplot() ?

  qplot(date, psavert, data = economics, geom = line,main=jhdjd)-p
  p+opts(text(locator(1),),new=T)

 does not work.
 --
 View this message in context:
 http://www.nabble.com/text%28%29-to-label-points-in-ggplot-tp23545135p23545135.html
 Sent from the R help mailing list archive at Nabble.com.

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


[[alternative HTML version deleted]]

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


Re: [R] how to control lattice plot parameters

2009-04-18 Thread Sunil Suchindran
We can see the plotting options using trellis.par.get(). For example, one
listed parameter is $superpose.line for which we can set col, lty, and lwd.
Load the mlmRev package to obtain the Gcsemv data (used in Lattice book).
One way to set the parameters is

data(Gcsemv, package = mlmRev)

my.fonts  - list(font = 2,cex = 2)

my.theme - list(
superpose.symbol = list(col=c(1,3)),
superpose.line=list(col=c(1,3),lwd=c(1,2)),
par.xlab.text = my.fonts,
par.ylab.text = my.fonts
)

xyplot(written ~ course , data = Gcsemv, groups = gender, type = c(p,r),

   par.settings = my.theme,
   auto.key = list(lines = TRUE, points = FALSE))
On Wed, Apr 15, 2009 at 5:40 AM, Paulo E. Cardoso pecard...@netcabo.ptwrote:

 I'm not being able to control all parameters of a xyplot.

 dataset example:
  tabp[1:10,]
   time  n X id   name
 1 1 95  0.00  1 Coral reef
 2 1 93  0.00  2 Coral reef
 3 1 92  0.00  3 Coral reef
 4 1 90  0.00  4 Coral reef
 5 1 87  8.994321  5 Coral reef
 6 1 86 12.580143  6 Coral reef
 7 1 84 17.004030  7 Coral reef
 8 1 83 18.469500  8 Coral reef
 9 1 82 37.919033  9 Coral reef
 101 81 39.059352 10 Coral reef
 ...

 plot code:
 xyplot(
  X~n,groups=name,data=tabp[tabp$time %in% c(1,4,8),],
  ylab=% of target,
  xlab=PU selection Frequency,
  lty=c(1,3,5), #! not responding
  xlim=c(-10,110),
  scales=list(cex=0.7,at=seq(0,360,by=20),labels=seq(0,360,by=20)), #!
 not responding
  panel=function(x,y,groups,subscripts)
  {
panel.xyplot(x,y,
subscripts=subscripts,
groups=groups,
type=l)
panel.abline(h=30,lty=2,col=grey50)
  },
  key=list(
space=top,
columns=1,
text=list(levels(tabp$name)[c(2,5,11)]),
lines=list(col=c(grey,blue,green))
)
  )

 A few controls are not responding and I don't know how to:
 1) control cex for xlab and ylab
 2) control position of legend
 3) control lty for each group level
 4) control plot groups colors
 5) match legend colors with graph colors
 Any help will be very appreciated.

 
 Paulo E. Cardoso

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


[[alternative HTML version deleted]]

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


Re: [R] loglm fitting

2009-01-14 Thread Sunil Suchindran
If we check the data, we see a mistake in the entry.

 tbl.8.3
  mara cig alc count
1  Yes Yes Yes   911
2   No  No  No   538
3  Yes Yes Yes44
4   No  No  No   456
5  Yes Yes Yes 3
6   No  No  No43
7  Yes Yes Yes 2
8   No  No  No   279

Try:

table8.3 - read.table(textConnection(alc cig mar count
Yes Yes Yes 911
Yes Yes No  538
Yes No  Yes 44
Yes No  No  456
No  Yes Yes 3
No  Yes No  43
No  No  Yes 2
No  No  No  279),header=TRUE)
closeAllConnections()
myglm - glm(count ~ mar+cig+alc, family=poisson, data=table8.3)
myglm$fitted.values

 myglm$fitted.values
1 2 3 4 5 6
7 8
539.98258 740.22612 282.09123 386.70007  90.59739 124.19392  47.32880
64.87990

This reproduces the fitted values for your example (Agresti, pg 323)

On Wed, Jan 14, 2009 at 7:23 AM, Gerard M. Keogh gmke...@justice.ie wrote:


 Dear all,

 sorry to bother you all with this but I've been trying to use the loglm in
 MASS package (v2.8.0) and cannot get any sensible output.
 I'm wondering am I doing something very foolish or missing something
 obvious.

 For example, I tried the documentation help(loglm) example - here's the
 code

 # Case 1: frequencies specified as an array.
 sapply(minn38, function(x) length(levels(x)))
 ## hs phs fol sex f
 ##  3   4   7   2 0
 minn38a - array(0, c(3,4,7,2), lapply(minn38[, -5], levels))
 minn38a[data.matrix(minn38[,-5])] - minn38$fol
 fm - loglm(~1 + 2 + 3 + 4, minn38a)  # numerals as names.
 deviance(fm)
  [1] 0

 The deviance is zero.
 I tried other examples as well (Laura Thompson) and get pretty much
 nonsense.

 E.g.

  library(MASS)
  options(contrasts=c(contr.treatment, contr.poly))
  tbl.8.3 = data.frame(mara
  =factor(c(Yes,No),levels=c(Yes,No)),
cig
  =factor(c(Yes,No),levels=c(No,Yes)),
alc
  =factor(c(Yes,No),levels=c(Yes,No)),
count=c(911,538,44,456,3,43,2,279) )
  fit.acm = loglm( count ~ alc*cig*mara,data=tbl.8.3, param=T, fit=T)


 Any advice greatly appreciated.

 Gerard



 **
 The information transmitted is intended only for the person or entity to
 which it is addressed and may contain confidential and/or privileged
 material. Any review, retransmission, dissemination or other use of, or
 taking of any action in reliance upon, this information by persons or
 entities other than the intended recipient is prohibited. If you received
 this in error, please contact the sender and delete the material from any
 computer.  It is the policy of the Department of Justice, Equality and Law
 Reform and the Agencies and Offices using its IT services to disallow the
 sending of offensive material.
 Should you consider that the material contained in this message is
 offensive you should contact the sender immediately and also mailminder[at]
 justice.ie.

 Is le haghaidh an duine nó an eintitis ar a bhfuil sí dírithe, agus le
 haghaidh an duine nó an eintitis sin amháin, a bheartaítear an fhaisnéis a
 tarchuireadh agus féadfaidh sé go bhfuil ábhar faoi rún agus/nó faoi
 phribhléid inti. Toirmisctear aon athbhreithniú, atarchur nó leathadh a
 dhéanamh ar an bhfaisnéis seo, aon úsáid eile a bhaint aisti nó aon ghníomh
 a dhéanamh ar a hiontaoibh, ag daoine nó ag eintitis seachas an faighteoir
 beartaithe. Má fuair tú é seo trí dhearmad, téigh i dteagmháil leis an
 seoltóir, le do thoil, agus scrios an t-ábhar as aon ríomhaire. Is é beartas
 na Roinne Dlí agus Cirt, Comhionannais agus Athchóirithe Dlí, agus na
 nOifígí agus na nGníomhaireachtaí a úsáideann seirbhísí TF na Roinne,
 seoladh ábhair cholúil a dhícheadú.
 Más rud é go measann tú gur ábhar colúil atá san ábhar atá sa
 teachtaireacht seo is ceart duit dul i dteagmháil leis an seoltóir
 láithreach agus le mailminder[ag]justice.ie chomh maith.

 ***



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


[[alternative HTML version deleted]]

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