[R] glmmADMB error

2014-01-30 Thread Sebastián Daza
I ran an example using glmmADMB, and I got this error:

data(bacteria,package=MASS)
bacteria$present - as.numeric(bacteria$y)-1
(bfit -  glmmadmb(present ~ trt + I(week  2), random = ~ 1 | ID,
+ family = binomial, data = bacteria))

Error in paste0(symbol1, paste0(paste0(var, collapse = symbol2))) :
  argument symbol1 is missing, with no default

Any ideas?

R version 3.0.2 (2013-09-25)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] splines   stats graphics  grDevices utils datasets  methods
base

other attached packages:
[1] multcomp_1.3-1TH.data_1.0-3 survival_2.37-7   mvtnorm_0.9-9997
 data.table_1.8.10 glmmADMB_0.7.7R2admb_0.7.10
[8] MASS_7.3-29   foreign_0.8-59

loaded via a namespace (and not attached):
[1] grid_3.0.2  lattice_0.20-24 Matrix_1.1-1.1  nlme_3.1-113
 sandwich_2.3-0  tools_3.0.2 zoo_1.7-10

-- 
Sebastián Daza

[[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.


[R] CTM and survival analysis with heterogeneity

2012-12-03 Thread Sebastián Daza
Hello R experts,

I wonder if there is any package to estimate  this kind of models in R:

Multi-state Multi-spell Survival Models with Heterogeneity

One of the most powerful programs for survival models is CTM, the
Continuous Time Model, developed for NIH at NORC under the direction
of James J. Heckman. Generalizing the competing risks model, CTM
allows for transitions between any number of states, repeated spells
within a state, time varying covariates, person and state specific
heterogeneity, and arbitrary duration dependence specified in a highly
flexible hazard model. Although in all cases the baseline hazard is
fully parametric, it can be specified as a sufficiently rich function
of time to capture a very wide set of duration dependencies.

Any references?
Thanks!

-- 
Sebastián Daza

__
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.


[R] survey package question

2012-10-11 Thread Sebastián Daza
Hello,

I have got a cluster sample using an election dataset where I already
had the final results of a county-specific election. I am trying to
figure out what would be the best sampling design for my data.

The  structure of the dataset is:

1) polling station (in general schools where people vote, for a
county, for example, there are 15 polling stations)
2) inside each polling station, there are voting units, where people
actually vote (on average there are about 40 voting units for polling
station)
3) for each voting unit I have the total votes by candidate (e.g.,
candidate 1 =322, candidate 2=122, candidate 3= 89)

The initial sampling design is:
1) selection of 5 polling stations PPS (based on number of voters)
2) selection of 10 voting units (SRS)

I am interested in estimating the proportion of votes by candidate
(let's assume we have 3 candidates). My naive estimate would be:

votes for candidate 1 / all valid votes = proportion

e.g.

candidate 1= 2132 / 10874= .1906
candidate 2= 5323 / 10874= .4895
candidate 3= 3419 / 10874= .3144

In this case, the unit of analysis is voters (or votes).

 If I specify the sampling design using the survey package in this way...

design -svydesign(id=~station + unit  fpc=~probstation +probunit,
data=sample, pps=brewer)

svyciprop(~I(candidate1/totalVotes), design)

... I am assuming that the unit of analysis is the voting unit, right?
and I am estimating an average among voting units?

I should expand my database at individual level (voters), or I just
have to include a unit weight according to the number of voters for
voting unit? In other words, is there a way to estimate, for instance,
votes for candidate 1 / all valid votes = proportion, directly from
the survey package or I have to expand  the database at people level
(voters), and then estimate the proportion using svymean and the
respective design.

I would appreciate any advice or help.

Sebastian

__
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] survey package question

2012-10-11 Thread Sebastián Daza
Hello Thomas,

I use both svymean (with the expanded sample = people), and svyratio
(voting unit level), using the same design:

design -svydesign(id=~station + unit, fpc=~probstation+probunits,
data=sample, pps=brewer)

I got different results using the same sample:

svyratio (voting unit)

Ratio   2.5%  97.5% Result
Cand1  0.05252871 0.04537301 0.05968441 0.05181146
Cand20.47226973 0.45215097 0.49238849 0.49041590
Cand3   0.47520156 0.45460831 0.49579482 0.45777264

svymean (expanded sample, individuals or votes)

   Mean  SE  2.5 % 97.5 %Results
Cand1   0.0528433 0.004562755 0.04390047 0.06178614 0.05181146
Cand2 0.4717504 0.010201398 0.45175605 0.49174480 0.49041590
Cand30.4754063 0.010429222 0.45496538 0.49584718 0.45777264

Point estimators are different, and confidence intervals are more
narrow using svyratio.
Could you give me any clue about what is going on?

Thank you in advance.
Sebastian

On Thu, Oct 11, 2012 at 7:50 PM, Sebastián Daza
sebastian.d...@gmail.com wrote:
 Hello Thomas,

 I use both svymean (with the expanded sample = people), and svyratio
 (voting unit level), using the same design:

 design -svydesign(id=~station + unit, fpc=~probstation+probunits,
 data=sample, pps=brewer)

 I got different results using the same sample:

 svyratio (voting unit)

 Ratio   2.5%  97.5% Result
 Cand1  0.05252871 0.04537301 0.05968441 0.05181146
 Cand20.47226973 0.45215097 0.49238849 0.49041590
 Cand3   0.47520156 0.45460831 0.49579482 0.45777264

 svymean (expanded sample, individuals or votes)

Mean  SE  2.5 % 97.5 %Results
 Cand1   0.0528433 0.004562755 0.04390047 0.06178614 0.05181146
 Cand2 0.4717504 0.010201398 0.45175605 0.49174480 0.49041590
 Cand30.4754063 0.010429222 0.45496538 0.49584718 0.45777264

 Point estimators are different, and confidence intervals are more
 narrow using svyratio.
 Could you give me any clue about what is going on?

 Thank you in advance.
 Sebastian

 On Thu, Oct 11, 2012 at 3:56 PM, Sebastián Daza
 sebastian.d...@gmail.com wrote:
 Thank you Thomas!

 On Thu, Oct 11, 2012 at 2:33 PM, Thomas Lumley tlum...@uw.edu wrote:
 On Fri, Oct 12, 2012 at 6:56 AM, Sebastián Daza
 sebastian.d...@gmail.com wrote:
 Hello,

 I have got a cluster sample using an election dataset where I already
 had the final results of a county-specific election. I am trying to
 figure out what would be the best sampling design for my data.

 The  structure of the dataset is:

 1) polling station (in general schools where people vote, for a
 county, for example, there are 15 polling stations)
 2) inside each polling station, there are voting units, where people
 actually vote (on average there are about 40 voting units for polling
 station)
 3) for each voting unit I have the total votes by candidate (e.g.,
 candidate 1 =322, candidate 2=122, candidate 3= 89)

 The initial sampling design is:
 1) selection of 5 polling stations PPS (based on number of voters)
 2) selection of 10 voting units (SRS)

 I am interested in estimating the proportion of votes by candidate
 (let's assume we have 3 candidates). My naive estimate would be:

 votes for candidate 1 / all valid votes = proportion

 e.g.

 candidate 1= 2132 / 10874= .1906
 candidate 2= 5323 / 10874= .4895
 candidate 3= 3419 / 10874= .3144

 In this case, the unit of analysis is voters (or votes).

  If I specify the sampling design using the survey package in this way...

 design -svydesign(id=~station + unit  fpc=~probstation +probunit,
 data=sample, pps=brewer)

 svyciprop(~I(candidate1/totalVotes), design)

 ... I am assuming that the unit of analysis is the voting unit, right?
 and I am estimating an average among voting units?


 You want a ratio estimator

 svyratio(~candidate1, ~totalVotes, design)


-thomas

 --
 Thomas Lumley
 Professor of Biostatistics
 University of Auckland



 --
 Sebastián Daza



 --
 Sebastián Daza



-- 
Sebastián Daza

__
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.


[R] Loop question

2012-05-30 Thread Sebastián Daza
Hello,
I have a dataframe (Lx) with 5 Lb, and 5 Lw variables. I want to
create several variables according to the loop presented below (data
attached).

Lx - read.csv(Lx.csv, header=T, sep=,)

for (i in 1:20) {
  Lx$sb1[i]  - Lx$Lb1[i+1]/Lx$Lb1[i]
  Lx$sb2[i]  - Lx$Lb2[i+1]/Lx$Lb2[i]
  Lx$sb3[i]  - Lx$Lb3[i+1]/Lx$Lb3[i]
  Lx$sb4[i]  - Lx$Lb4[i+1]/Lx$Lb4[i]
  Lx$sb5[i]  - Lx$Lb5[i+1]/Lx$Lb5[i]
  Lx$sw1[i]  - Lx$Lw1[i+1]/Lx$Lw1[i]
  Lx$sw2[i]  - Lx$Lw2[i+1]/Lx$Lw2[i]
  Lx$sw3[i]  - Lx$Lw3[i+1]/Lx$Lw3[i]
  Lx$sw4[i]  - Lx$Lw4[i+1]/Lx$Lw4[i]
  Lx$sw5[i]  - Lx$Lw5[i+1]/Lx$Lw5[i]
}

How I could also create the variable names (letters b and w, and
numbers from 1 to 5 in s and L variables) using a loop in R?
In Stata I can use this command:

 foreach r in b w {;
   foreach s of numlist 0(5)40 {;
   foreach ed of numlist 1/5 {;
   local f = `s'+5;
   scalar define rL`r'`ed'_`f'over`s' = L`r'`ed'_`f' / L`r'`ed'_`s';
   };
 };
 };

but I do not how to do it in R.
Thank you in advance!

-- 
Sebastián Daza
__
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] trellis plot

2012-03-26 Thread Sebastián Daza
Thank you Peter. The problem with you solution is that it doesn't
represent the actual values of disruption. Look this example:

person  - rep(1:2, each=4)
income  - c(100, 120, 150, 200, 90, 100,120, 150)
disruption  - c(0,0,0,1,0,1,1,0)
time  - rep(c(1:4),2)
dat  - as.data.frame(cbind(person,time, income, disruption))

mycols - c(2, 5)
library(lattice)
xyplot(income~time|as.factor(person),data=dat,
       type=c(p,g,o), col.line=black, col.symbol =
mycols[dat$disruption+1],
       xlab=Time,
       ylab=Familiar Income)

I am looking for a way to get the colors of dots according to actual
disruption values, and at the same time to draw the lines between all
dots for each person. Thank you!



 Hi everyone,
 I am just trying to figure out how to do a xyplot where in addition to
  dots and lines I can change dots' colors according to an individual
 variable (e.g., marital disruption across time, a dummy 0/1). When I
 use groups specification (see below), I get two different lines for
 each individual based on groups, and what I want is to get one line
 connecting dots, and different dots' colors according to marital
 disruption.

 Any ideas about how to do that?

 person- rep(1:2, each=4)
 income- c(100, 120, 150, 200, 90, 100,120, 150)
 disruption- rep(c(0,1), 4)
 time- rep(c(1:4),2)
 dat- as.data.frame(cbind(person,time, income, disruption))


 library(lattice)
 xyplot(income~time|as.factor(person),data=dat,
        type=c(p,g,o), col.line=black,
        xlab=Time,
        ylab=Familiar Income)

 # I just want to change dots' colors according to disruption, not to
 get two different lines:

 xyplot(income~time|as.factor(person),data=dat,
        type=c(p,g,o), col.line=black, groups=disruption,
        xlab=Time,
        ylab=Familiar Income)


 If I understand correctly what you're after, try adding the
 'col.symbol' argument instead of the 'groups' argument.

     col.symbol = disruption + 1

 or, better, for arbitrary colours:

    mycols - c(2, 5)
    xyplot(,
           col.symbol = mycols[disruption + 1],
           )

 Peter Ehlers



 Thank you in advance!





 --
 Sebastián Daza



-- 
Sebastián Daza

__
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.


[R] trellis plot

2012-03-25 Thread Sebastián Daza
Hi everyone,
I am just trying to figure out how to do a xyplot where in addition to
 dots and lines I can change dots' colors according to an individual
variable (e.g., marital disruption across time, a dummy 0/1). When I
use groups specification (see below), I get two different lines for
each individual based on groups, and what I want is to get one line
connecting dots, and different dots' colors according to marital
disruption.

Any ideas about how to do that?

person  - rep(1:2, each=4)
income  - c(100, 120, 150, 200, 90, 100,120, 150)
disruption  - rep(c(0,1), 4)
time  - rep(c(1:4),2)
dat  - as.data.frame(cbind(person,time, income, disruption))


library(lattice)
xyplot(income~time|as.factor(person),data=dat,
   type=c(p,g,o), col.line=black,
   xlab=Time,
   ylab=Familiar Income)

# I just want to change dots' colors according to disruption, not to
get two different lines:

xyplot(income~time|as.factor(person),data=dat,
   type=c(p,g,o), col.line=black, groups=disruption,
   xlab=Time,
   ylab=Familiar Income)



Thank you in advance!

-- 
Sebastián Daza

__
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.


[R] calculations combining values from different rows

2012-02-09 Thread Sebastián Daza
Hi everyone,
I looking for functions or systematic ways to do calculations between
different columns and rows in R, as one can do easily in Excel. For
example, I have two variables, a and b, where a1 represents an a value
in row 1, and b2 represents a b value in row 2, etc.

a  - c(4,3,5,5,6,7,3,2,1,4)
b  - c(2,4,1,2,5,3,1,8,7,5)
data  - cbind(a,b)

I have to calculate something like this:

x1 = NA
x2 = -b1 /24 * a1 + b2 /2 * a2 + b3 /24 * a3
x3 = -b2 /24 *a2 + b3 /2 * a3 + b4 /24 * a4
x4 = -b3 /24 *a3 + b4 /2 * a4 + b5 /24 * a5
...
x9 = -b8 /24* a8 + b9 /2 * a9 + b10 /24 * a10
x10= NA

For example, x2 would be equal to: -2/24*4 +4/2*3 + 1/24 *5

Any ideas?
Thank you in advance.

-- 
Sebastián Daza

__
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.


[R] data management question

2012-02-07 Thread Sebastián Daza
Hi everyone,
I would like to have a function to compute some values in a dataset.
First, I have to define a value for the lx variable in row 1 (e.g.,
100,000), npx is a given proportion. lx of row 2 is equal to lx of row
1 times npx of row 1. I can do this row by row...

data[1,lx]  - 10
data[2,lx]  - data[1,lx]*data[1,npx]
data[3,lx]  - data[2,lx]*data[2,npx]
data[4,lx]  - data[3,lx]*data[3,npx]
...
data[19,lx]  - data[18,lx]*data[18,npx]

Any ideas about how to define this in a function or in a more systematic way?
Thank you in advance.

-- 
Sebastián Daza

__
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.


[R] raking weighting

2011-11-09 Thread Sebastián Daza
Hi everyone,
Does anyone know if there is a package to compute raking weights using
R? What I need is to create a variable with weights base in some
demographic variables (e.g. sexo, age group, area) using the raking
procedure.

Thank you in advance!

-- 
Sebastián Daza

__
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] raking weighting

2011-11-09 Thread Sebastián Daza
I did. However, I am not sure that I can create a weight variable in
my database with the rake function you mention. That is why I am
asking.

El día 9 de noviembre de 2011 13:44, David Winsemius
dwinsem...@comcast.net escribió:

 On Nov 9, 2011, at 2:27 PM, Sebastián Daza wrote:

 Hi everyone,
 Does anyone know if there is a package to compute raking weights using
 R? What I need is to create a variable with weights base in some
 demographic variables (e.g. sexo, age group, area) using the raking
 procedure.

 There is a `rake` function in package survey.

 In the future you ought to do this at the R console _before_ sending a
 question to R-help:

 RSiteSearch(raking)


 Thank you in advance!

 --
 Sebastián Daza

 __
 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.

 David Winsemius, MD
 West Hartford, CT





-- 
Sebastián Daza

__
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.


[R] identifying groups

2011-05-20 Thread Sebastián Daza

Hi everyone,
Does anyone know if there is a package to perform a  Moody's Crowds 
routine to identify groups using R, or other algorithms designed to 
search groups by maximizing modularity scores? Otherwise, other packages 
to identify groups, or other programs to do it.


I have about 80 networks clustered in schools. So, I have to identify 
groups for each school.


Thank  you in advance!
--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] data network format and grouping analysis

2011-05-18 Thread Sebastián Daza

Hi everyone,
I have a dataset of friendship with this format:

 ego alter
47461 2
97421 3
14738   1NA
47472NA
974323
14739   21
4748313
97443   5
14740   314
47494NA
97454NA
14741   4NA
47505NA
9746513
14742   510
47516   12
97476 7
...


NA means that individuals don't select any friend. Does anyone know how 
to format this dataset to use sna or igraph packages? I don't know how 
to convert it into a matrix or a edgelist in R without losing isolated 
individuals .


Next question, anyone knows if there is a package to perform a  Moody's 
Crowds routine to identify groups using R, or other algorithms designed 
to search groups by maximizing modularity scores?


Thank  you in advance!

--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] setCoefTemplate

2011-04-03 Thread Sebastián Daza

Hi everyone,
I am trying to build a table putting standard errors horizontally. I 
haven't been able to do it.


library(memisc)
berkeley - aggregate(Table(Admit,Freq)~.,data=UCBAdmissions)

berk0 - glm(cbind(Admitted,Rejected)~1,data=berkeley,family=binomial)
berk1 - 
glm(cbind(Admitted,Rejected)~Gender,data=berkeley,family=binomial)
berk2 - 
glm(cbind(Admitted,Rejected)~Gender+Dept,data=berkeley,family=binomial)


setCoefTemplate(est.se=c(est = ($est:#)($se:#)))

mtable(berk0,berk1,berk2,
+   coef.style=est.se,
+   summary.stats=c(Deviance,AIC,N))
Error in dim(ans) - newdims :
  dims [product 1] do not match the length of object [2]

Thank you in advance.
--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] setCoefTemplate

2011-04-03 Thread Sebastián Daza

Hi everyone,
I am trying to build a table putting standard errors horizontally. I 
haven't been able to do it.


library(memisc)
berkeley - aggregate(Table(Admit,Freq)~.,data=UCBAdmissions)

berk0 - glm(cbind(Admitted,Rejected)~1,data=berkeley,family=binomial)
berk1 - 
glm(cbind(Admitted,Rejected)~Gender,data=berkeley,family=binomial)
berk2 - 
glm(cbind(Admitted,Rejected)~Gender+Dept,data=berkeley,family=binomial)


setCoefTemplate(est.se=c(est = ($est:#)($se:#)))

mtable(berk0,berk1,berk2,
+   coef.style=est.se,
+   summary.stats=c(Deviance,AIC,N))
Error in dim(ans) - newdims :
  dims [product 1] do not match the length of object [2]

Thank you in advance.
--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] coefficients style

2011-04-03 Thread Sebastián Daza

Hi everyone,
I am trying to build a table putting standard errors horizontally. I 
haven't been able to do it.


library(memisc)
berkeley - aggregate(Table(Admit,Freq)~.,data=UCBAdmissions)

berk0 - glm(cbind(Admitted,Rejected)~1,data=berkeley,family=binomial)
berk1 - 
glm(cbind(Admitted,Rejected)~Gender,data=berkeley,family=binomial)
berk2 - 
glm(cbind(Admitted,Rejected)~Gender+Dept,data=berkeley,family=binomial)


setCoefTemplate(est.se=c(est = ($est:#)($se:#)))

mtable(berk0,berk1,berk2,
+ coef.style=est.se,
+ summary.stats=c(Deviance,AIC,N))
Error in dim(ans) - newdims :
  dims [product 1] do not match the length of object [2]

Thank you in advance.

--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] style question

2011-04-03 Thread Sebastián Daza

Hi everyone,
I am trying to build a table putting standard errors horizontally. I 
haven't been able to do it.


library(memisc)
berkeley - aggregate(Table(Admit,Freq)~.,data=UCBAdmissions)

berk0 - glm(cbind(Admitted,Rejected)~1,data=berkeley,family=binomial)
berk1 - 
glm(cbind(Admitted,Rejected)~Gender,data=berkeley,family=binomial)
berk2 - 
glm(cbind(Admitted,Rejected)~Gender+Dept,data=berkeley,family=binomial)


setCoefTemplate(est.se=c(est = ($est:#)($se:#)))

mtable(berk0,berk1,berk2,
+ coef.style=est.se,
+ summary.stats=c(Deviance,AIC,N))
Error in dim(ans) - newdims :
  dims [product 1] do not match the length of object [2]

Thank you in advance.

--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] simple if question

2011-03-26 Thread Sebastián Daza

Hi everyone,
I have just got different samples from a dataframe (independent and 
exclusive, there aren't common elements among them). I want to create a 
variable that indicate the sampling selection of the elements in the 
original dataframe (for example, 0 = no selected, 1= sample 1, 2=sample 
2, etc.).


I have tried to do it with ifelse command, but the problem is that the 
second line replaces the values of the first line, and I haven't been 
able to do it with the if command (I got this error: In if (data$school 
%in% sample1) { :

  the condition has length  1 and only the first element will be used)

data$selection - ifelse(data$school %in% sample1, 1, 0)
data$selection - ifelse(data$school %in% sample2, 2, 0)

Any ideas?
Thank you in advance.

--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] pattern of panel waves using R

2011-03-18 Thread Sebastián Daza
Finally, I solved my problem using the following procedure and a 
database called ej


ej$a - 1

head(ej)

   ano nunico a
1 2008  1 1
2 2009  1 1
3 2008  2 1
4 2009  2 1
5 2008  3 1
6 2009  3 1

library(reshape)
dej - cast(ej, nunico ~ ano, sum, margins = FALSE)

head(dej)

  nunico 2008 2009 2010
1  1110
2  2110
3  3110
4  4111
5  5101
6  6111

dej[dej==0] - NA

library(Hmisc)
na.pattern(dej[,c(2:4)])


pattern
 000  001  010  011  110
3385 1073  203  338  573

That way I can review the pattern of my panel data.

Sebastian

__
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.


[R] pattern of panel waves using R

2011-03-17 Thread Sebastián Daza

Hi everyone,

Is there any command to identify the pattern of responses of a database 
with this format:


year  id
20081
20091
20082
20092
20083
20093
20084
20094
20104


I just need the frequency of the patterns grouped by id:

2008 2009 2010 = 80
2009 2010 = 30
2008 2009 = 10
and so on

Thank you in advance!

--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] lavaan diagram

2011-03-08 Thread Sebastián Daza

Hi everyone,
I got the following error when I tried to diagram a confirmatory factor 
analysis using lavaan:


lavaan.diagram(conf)
Error in strwidth(xvars) : plot.new has not been called yet

Does anyone have any idea about what is happening?
Thank you in advance.

--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] lavaan diagram

2011-03-06 Thread Sebastián Daza

Hi everyone,
I get the following message when I try to get a diagram of a CFA:

Error in lavaan.diagram(fit) :
  no slot of name GLIST for this object of class Fit

Does anyone know what the problem is?

Regards.
--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] aggregate in R

2011-02-22 Thread Sebastián Daza

plyr is very useful to aggregate data... I strongly recommend it.

On 2/22/2011 5:59 PM, Jeff Newmiller wrote:

Use ?plyr::ddply

---
Jeff Newmiller The . . Go Live...
DCN:jdnew...@dcn.davis.ca.us  Basics: ##.#. ##.#. Live Go...
Live: OO#.. Dead: OO#.. Playing
Research Engineer (Solar/Batteries O.O#. #.O#. with
/Software/Embedded Controllers) .OO#. .OO#. rocks...1k
---
Sent from my phone. Please excuse my brevity.

Hongwei Dongpdxd...@gmail.com  wrote:

Hi, R users, I'm wondering how I can aggregate data in R with different functions for different columns. For example: x-rep(1:5,3) 
y-cbind(x,a=1:15,b=21:35) y-data.frame(y) I want to aggregate a and b in y by x. With a, I want to use 
function mean; with b, I want to use function sum. I tried:  aggregate(y,x,mean(y$a),sum(y$b)) But I got the error: Error 
in match.fun(FUN) : 'mean(y$a)' is not a function, character or symbol Anyone can tell me how to fix this problem? Thanks. Gary [[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.


[[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.



--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] id number by group and correlative

2011-02-16 Thread Sebastián Daza

Hello everyone,
I am new in R and I am trying to create a id number (a correlative 
sequence of numbers) by group, and a correlative sequence of numbers 
inside each group (my idea is to get statistics by group  without having 
to aggregate the database). Here an example:


group   id_groupcorrelative_group
A   1   1
A   1   2
A   1   3
A   1   4
B   2   1
B   2   2
B   2   3
C   3   1
C   3   2
C   3   3
C   3   4
C   3   5

Unfortunately, I have been able to find an explicit lag function to get 
id_group (I know I can get it using aggregate and merge but I'm just 
wondering if there is another way to do it). With regard to the 
correlative_group, I don't have any clue about how to do it.


PD: Last question. Is there any way to save in a variable the internal 
that R uses (those numbers that appear on left side of a dataframe)?


Thank you,
Sebastian.
--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] very basic HLM question

2011-02-06 Thread Sebastián Daza

Thank you for your reply and sorry for my ambiguity.

I computed:

summary(anova1 - aov(math ~ as.factor(schoolid), data=nels88))
ICC1(anova1)

ICC1 comes from the multilevel package. I have found an article where it 
is pointed out that with this formula:


  34.011/(34.011+72.256)
[1] 0.3200523

I should get practically an identical result than using ICC1. I have 
used the database of that article and I have got an identical result. 
But with the dataset I am using (schools) it doesn't work.


This is the example of the article mentioned above:

library(multilevel)
base(bh1996)

summary(lmer(WBEING ~ 1 + (1|GRP), data=bh1996))

0.035801/(0.035801+0.789497)

summary(test - aov(WBEING~as.factor(GRP),data=bh1996))
ICC1(test)


I have attached the small database where this equality doesn't exist.
Thank you in advance.


On 2/5/2011 11:07 PM, Paul Johnson wrote:

2011/2/5 Sebastián Dazasebastian.d...@gmail.com:

Hi everyone,

I need to get a between-component variance (e.g. random effects Anova), but
using lmer I don't get the same results (variance component) than using
random effects Anova. I am using a database of students, clustered on
schools (there is not the same number of students by school).

According to the ICC1 command, the interclass correlation is .44


ICC1(anova1)

[1] 0.4414491


If you don't tell us exactly what model you are calculating in
anova1, how would we guess if there is something wrong?

Similarly, I get this

ICC1

Error: object 'ICC1' not found

so it must mean you've loaded a package or written a function, which
you've not shown us.

I googled my way to a package called multilevel that has ICC1, and
its code for ICC1 shows a formula that does not match the one you used
to calculate ICC from lmer.

function (object)
{
 MOD- summary(object)
 MSB- MOD[[1]][1, 3]
 MSW- MOD[[1]][2, 3]
 GSIZE- (MOD[[1]][2, 1] + (MOD[[1]][1, 1] + 1))/(MOD[[1]][1,
 1] + 1)
 OUT- (MSB - MSW)/(MSB + ((GSIZE - 1) * MSW))
 return(OUT)
}

I'm not saying that's right or wrong, just not obviously identical to
the formula you proposed.



However, I cannot get the same ICC from the lmer output:


anova2- lmer(math ~ 1 + (1|schoolid), data=nels88)
summary(anova2- lmer(math ~ 1 + (1|schoolid), data=nels88))




Instead, do this (same thing, fits model only once):


anova2- lmer(math ~ 1 + (1|schoolid), data=nels88)
summary(anova2)


Note that lmer is going to estimate a normally distributed random
effect for each school, as well as an individual observation random
effect (usual error term) that is assumed independent of the
school-level effect.  What is anova1 estimating?




--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] very basic HLM question

2011-02-05 Thread Sebastián Daza

Hi everyone,

I need to get a between-component variance (e.g. random effects Anova), 
but using lmer I don't get the same results (variance component) than 
using random effects Anova. I am using a database of students, clustered 
on schools (there is not the same number of students by school).


According to the ICC1 command, the interclass correlation is .44

 ICC1(anova1)
[1] 0.4414491

However, I cannot get the same ICC from the lmer output:

 anova2 - lmer(math ~ 1 + (1|schoolid), data=nels88)
 summary(anova2 - lmer(math ~ 1 + (1|schoolid), data=nels88))

Linear mixed model fit by REML
Formula: math ~ 1 + (1 | schoolid)
   Data: nels88
  AIC  BIC logLik deviance REMLdev
 1878 1888 -935.8 18751872
Random effects:
 Groups   NameVariance Std.Dev.
 schoolid (Intercept) 34.011   5.8319
 Residual 72.256   8.5003
Number of obs: 260, groups: schoolid, 10

Fixed effects:
Estimate Std. Error t value
(Intercept)   48.861  1.927   25.36

The intercept random effect is 34.011. If I compute the ICC manually I get:

 34.011/(34.011+72.256)
[1] 0.3200523

According to my Anova analysis, the between-component variance is 59.004.
Does anyone know what the problem is? How can I get the 59.004 figure 
using R?



--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] grey scale graphs

2011-02-02 Thread Sebastián Daza

Hi everyone,
Does anyone know how to get black and white theme (grey scale,, I 
would say) graphs using lattice or ggplot2, as it is shown in this 
webpage: http://lmdvr.r-forge.r-project.org/figures/figures.html?


I am using Sweave, and I cannot get that color configuration.
I have added the following option: trellis.device(color=FALSE) but I got 
a pdf file with color graphs.


Thank in advance.
--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Latent Class Logit Models in discrete choice experiments

2011-01-31 Thread Sebastián Daza

See: https://www.msu.edu/~chunghw/downloads.html
Maybe you can find something useful there!

Regards

On 1/31/2011 12:35 PM, Daniel Vecchiato wrote:

Dear R users,

I would like to perform Latent Class Logit Models for the analysis of choice 
experiments in environmental valuation.

This kind of analysis is usually performed with NLogit Software 
(http://www.limdep.com).

I attach the results I usually obtain using NLogit and NLogit model 
specifications.

For Random parameter models and Logit Models I usually perform my analysis with 
the package mlogit ( http://cran.r-project.org/web/packages/mlogit/index.html ).

The models I would like to run are presented in this ppt:

http://pages.stern.nyu.edu/~wgreene/.../Lectures/Part13-LatentClassModels.ppt

and an overview is given in this paper by Hoyos:

http://dx.doi.org/10.1016/j.ecolecon.2010.04.011

Anybody has any package to suggest for this kind of analysis? (poLCA does not 
provide me the same estimates)


Thanks in advance

Daniel


__
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.



--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] rj packages and eclipse

2011-01-30 Thread Sebastián Daza

Hi everyone!

Does anyone know how to install the rj package in Windows (7)? I am 
trying to set up eclipse, but I got [INFO] The R package 'rj' is not 
available, R-StatET tools cannot be initialized.


Thank you in advance!

__
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.


[R] LTA

2011-01-24 Thread Sebastián Daza

Hi everyone,
Does anyone know if there is a package to run Latent Transitional 
Analysis using R?


Regards!

--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Using summaryBy with weighted data

2011-01-17 Thread Sebastián Daza

Hi everyone,
I am trying to run Sweave.bat (batchfiles_0.6-1) from the command line 
on Windows, but I get this error:


C:\batchfiles_0.6-1Sweave.bat Sweave-test-1
Error: rterm.exe not found

I don't know how to set up the path if this one were the problem... I 
ran rcmd.bat and I got this... so I don't know if it is a path problem.


C:\batchfiles_0.6-1Rcmd,bat
R_ARCH=/x64
R_ARCH0=x64
R_ARCH0=x64
cmdpath=C:\R\R-2.12.1\bin\x64\Rcmd.exe
args=,bat
'bat' is not recognized as an internal or external command,
operable program or batch file.

the path of rterm.exe in my computer is: C:\R\R-2.12.1\bin\x64
thank you in advance!

--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] sweave.bat

2011-01-17 Thread Sebastián Daza

Hi everyone,
I am trying to run Sweave.bat (batchfiles_0.6-1) from the command line 
on Windows, but I get this error:


C:\batchfiles_0.6-1Sweave.bat Sweave-test-1
Error: rterm.exe not found

I don't know how to set up the path if this one were the problem... I 
ran rcmd.bat and I got this... so I don't know if it is a path problem.


C:\batchfiles_0.6-1Rcmd,bat
R_ARCH=/x64
R_ARCH0=x64
R_ARCH0=x64
cmdpath=C:\R\R-2.12.1\bin\x64\Rcmd.exe
args=,bat
'bat' is not recognized as an internal or external command,
operable program or batch file.

the path of rterm.exe in my computer is: C:\R\R-2.12.1\bin\x64
thank you in advance!

--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] easy loop question

2011-01-12 Thread Sebastián Daza

Hi everyone,

I am new in R and programming. I have tried to remove the values out of 
range in some variables using a loop:


1)

var  - names(est8vo[, 77:83])   # I got the variable names

 var
[1] p16.1 p16.2 p16.3 p16.4 p16.5 p16.6 p16.7

for (i in 1:7)  {
var.i  - var[i]
est8vo$var.i[ est8vo$var.i==3] - 99
  }

I got this error:

Error in `$-.data.frame`(`*tmp*`, var.i, value = numeric(0)) :
  replacement has 0 rows, data has 215700


2) The second step would be to define the factors, but I got the same error:

for (i in 1:7)  {
var.i  - var[i]
est8vo$var.i- factor(est8vo$var.i,
  levels=c(0, 1, 2, 99),
  labels=c(vacío, sí, no, doble marca)
  )
  }


I don't know how to do it.
Thank you in advance!
Sebastian

__
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.


[R] integration Sweave and TexMakerX

2011-01-05 Thread Sebastián Daza

Hi,

Does anyone know how to integrate texmakerx and sweave on Windows? I 
mean, to run .rnw files directly from texmakerx  and get a pdf or dvi file.


Thank you in advance,

--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.


[R] two-part growth analysis

2010-12-21 Thread Sebastián Daza

Hi everyone!

Does anyone know if there is a package to do two-part growth analysis 
with R?


Regards,
Sebastian


--
Sebastián Daza
sebastian.d...@gmail.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.html
and provide commented, minimal, self-contained, reproducible code.