[R] Three-dimensional contingency table
Hi,
I am trying to assemble a three-way contingency table examining the
presence/absence of mussels, water depth (Depth1 and Depth 2) and water
velocity (Flow vs. No Flow). I have written the following code listed below;
however, when run the glm I get the following message, "Error in
model.frame.default(formula = Count ~ MP + wd + wv, drop.unused.levels = TRUE)
: variable lengths differ (found for 'MP')". This may be something simple, if
so I apologize. Any help would be greatly appreciated.
Best,
C.R.
numbers <- c(1134,956,328,529,435,599,27,99)
dim(numbers) <- c(2,2,2)
numbers
dimnames(numbers)[[3]] <-list("Mussels", "No Mussels")
dimnames(numbers)[[2]] <- list("Flow", "No Flow")
dimnames(numbers)[[1]] <- list("Depth1", "Depth2")
ftable(numbers)
as.data.frame.table(numbers)
frame <- as.data.frame.table(numbers)
names(frame) <- c("wd", "wv", "MP", "Count")
frame
attach(frame)
model1 <- glm(Count~MP+wd+wv,poisson)
__
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] Question regarding significance of a covariate in a coxme survival model
On Sat, Aug 28, 2010 at 8:38 PM, David Winsemius wrote: > > On Aug 27, 2010, at 4:32 PM, Teresa Iglesias wrote: > > Christopher David Desjardins umn.edu> writes: >> >> >>> Hi, >>> I am running a Cox Mixed Effects Hazard model using the library coxme. I >>> am trying to model time to onset (age_sym1) of thought problems (e.g. >>> hearing voices) (sym1). As I have siblings in my dataset, I have >>> decided to account for this by including a random effect for family >>> (famid). My covariate of interest is Mother's diagnosis where a 0 is >>> bipolar, 1 is control, and 2 is major depression. I am trying to fit >>> the following model. >>> >>> thorp1 <- coxme(Surv(age_sym1, sym1) ~ lifedxm + (1 | famid), data = >>> bip.surv) >>> >>> Which provides the following output: >>> >>> - >>> thorp1 >>> Cox mixed-effects model fit by maximum likelihood >>> Data: bip.surv >>> events, n = 99, 189 (3 observations deleted due to missingness) >>> Iterations= 10 54 >>>NULL Integrated Penalized >>> Log-likelihood -479.0372 -467.3549 -435.2096 >>> Chisqdf p AIC BIC >>> Integrated loglik 23.36 3.00 3.3897e-05 17.36 9.58 >>> Penalized loglik 87.66 47.27 3.2374e-04 -6.88 -129.54 >>> Model: Surv(age_sym1, sym1) ~ lifedxm + (1 | famid) >>> Fixed coefficients >>> coef exp(coef) se(coef) z p >>> lifedxmCONTROL -1.0576257 0.3472794 0.3872527 -2.73 0.0063 >>> lifedxmMAJOR -0.6329957 0.5309987 0.3460847 -1.83 0.0670 >>> Random effects >>> Group Variable Std DevVariance >>> famid Intercept 0.980 0.9757033 >>> >>> -- >>> >>> So it appears that there is a difference in hazard rates associated with >>> Mother's diagnosis between BIPOLAR and CONTROL. Just to be safe, I fit a >>> model without this covariate present and decided to compare the models >>> with AIC and see if fit decreased with this covariate not in the model. >>> >>> -- >>> thorp1.n <- coxme(Surv(age_sym1, sym1) ~ (1 | famid), data = bip.surv) >>> thorp1.n >>> Cox mixed-effects model fit by maximum likelihood >>> Data: bip.surv >>> events, n = 99, 189 (3 observations deleted due to missingness) >>> Iterations= 10 54 >>>NULL Integrated Penalized >>> Log-likelihood -479.0372 -471. -436.0478 >>> Chisqdf pAIC BIC >>> Integrated loglik 15.41 1.00 0.8663 13.41 10.81 >>> Penalized loglik 85.98 49.48 0.00099915 -12.97 -141.37 >>> Model: Surv(age_sym1, sym1) ~ (1 | famid) >>> Random effects >>> Group Variable Std Dev Variance >>> famid Intercept 1.050949 1.104495 >>> >>> >>> The problem is that the AIC for the model w/o lifedxm is 13.41 and the >>> model w/ lifedxm is 17.36. So fit actually improved without lifedxm in >>> the model but really the fit is indistinguishable if I use Burnham & >>> Anderson, 2002's criteria. >>> >>> So my conundrum is this - The AIC and the p-values appear to disagree. >>> The p-value would suggest that lifedxm should be included in the model >>> and the AIC says it doesn't improve fit. In general, I tend to favor the >>> AIC (I usually work with large sample sizes and multilevel models) but I >>> am just beginning with survival models and I am curious if a survival >>> model guru might suggest whether lifedxm needs to be in the model or not >>> based on my results here? Would people generally use the AIC in this >>> situation? Also, I can't run anova() on this models because of the >>> random effect. >>> >>> I am happy to provide the data if necessary. >>> >>> Please cc me on a reply. >>> >>> Thanks, >>> Chris >>> >> >> Hi Chris, >> Did you get an answer to why the p-value and AIC score disagreed? >> I am getting the same results with my own data. They're not small >> disagreements either. The AIC score is telling me the opposite of >> what the p-value and the parameter coef are saying. >> The p-value and the coef for the predictor variable are in agreement. >> >> I've also noticed in some published papers with tables containing >> p-values and AIC scores that the chisq p-value and AIC score >> are in opposition too. This makes me think I'm missing something >> and that this all actually makes sense. >> >> However given that AIC = â 2 (log L) + 2K (where K is the number of >> parameters) >> and seeing as how the log-likelihood scores below are in the hundreds, >> shouldn't >> the AIC score be in the hundreds as well?? >> > > That is different than my understanding of AIC. I thought that the AIC and > BIC both took as input the difference in -2LL and then adjusted those > differences for the differences in number of degrees of freedom. That > perspective would inform one that the absolute value of the deviance ( = > -2LL) is immaterial and o
Re: [R] odd subset behavior
Dear Erin, -0.50 is greater than -1.00, this means that your logical test returns all FALSE answers. Now consider the results of: x[FALSE] you are selecting none of 'x', hence numeric(0). Perhaps you mean: x[x <= -0.50 & x >= -1.0] HTH, Josh On Sat, Aug 28, 2010 at 7:47 PM, Erin Hodgess wrote: > Dear R People: > > I just produced the following: > > x <- seq(-3,3,.01) > x[x>= -0.50 & x<= -1.0] > numeric(0) > x[x>= -0.50 && x<= -1.0] > numeric(0) > > What am I doing wrong, please? This seems strange. > > Thanks, > Erin > > > -- > Erin Hodgess > Associate Professor > Department of Computer and Mathematical Sciences > University of Houston - Downtown > mailto: erinm.hodg...@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. > -- Joshua Wiley Ph.D. Student, Health Psychology University of California, Los Angeles http://www.joshuawiley.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] Please ignore previous stupid question.
Sorry. -- Erin Hodgess Associate Professor Department of Computer and Mathematical Sciences University of Houston - Downtown mailto: erinm.hodg...@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] odd subset behavior
Dear R People: I just produced the following: x <- seq(-3,3,.01) x[x>= -0.50 & x<= -1.0] numeric(0) x[x>= -0.50 && x<= -1.0] numeric(0) What am I doing wrong, please? This seems strange. Thanks, Erin -- Erin Hodgess Associate Professor Department of Computer and Mathematical Sciences University of Houston - Downtown mailto: erinm.hodg...@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] Question regarding significance of a covariate in a coxme survival model
On Aug 27, 2010, at 4:32 PM, Teresa Iglesias wrote: Christopher David Desjardins umn.edu> writes: Hi, I am running a Cox Mixed Effects Hazard model using the library coxme. I am trying to model time to onset (age_sym1) of thought problems (e.g. hearing voices) (sym1). As I have siblings in my dataset, I have decided to account for this by including a random effect for family (famid). My covariate of interest is Mother's diagnosis where a 0 is bipolar, 1 is control, and 2 is major depression. I am trying to fit the following model. thorp1 <- coxme(Surv(age_sym1, sym1) ~ lifedxm + (1 | famid), data = bip.surv) Which provides the following output: - thorp1 Cox mixed-effects model fit by maximum likelihood Data: bip.surv events, n = 99, 189 (3 observations deleted due to missingness) Iterations= 10 54 NULL Integrated Penalized Log-likelihood -479.0372 -467.3549 -435.2096 Chisqdf p AIC BIC Integrated loglik 23.36 3.00 3.3897e-05 17.36 9.58 Penalized loglik 87.66 47.27 3.2374e-04 -6.88 -129.54 Model: Surv(age_sym1, sym1) ~ lifedxm + (1 | famid) Fixed coefficients coef exp(coef) se(coef) z p lifedxmCONTROL -1.0576257 0.3472794 0.3872527 -2.73 0.0063 lifedxmMAJOR -0.6329957 0.5309987 0.3460847 -1.83 0.0670 Random effects Group Variable Std DevVariance famid Intercept 0.980 0.9757033 -- So it appears that there is a difference in hazard rates associated with Mother's diagnosis between BIPOLAR and CONTROL. Just to be safe, I fit a model without this covariate present and decided to compare the models with AIC and see if fit decreased with this covariate not in the model. -- thorp1.n <- coxme(Surv(age_sym1, sym1) ~ (1 | famid), data = bip.surv) thorp1.n Cox mixed-effects model fit by maximum likelihood Data: bip.surv events, n = 99, 189 (3 observations deleted due to missingness) Iterations= 10 54 NULL Integrated Penalized Log-likelihood -479.0372 -471. -436.0478 Chisqdf pAIC BIC Integrated loglik 15.41 1.00 0.8663 13.41 10.81 Penalized loglik 85.98 49.48 0.00099915 -12.97 -141.37 Model: Surv(age_sym1, sym1) ~ (1 | famid) Random effects Group Variable Std Dev Variance famid Intercept 1.050949 1.104495 The problem is that the AIC for the model w/o lifedxm is 13.41 and the model w/ lifedxm is 17.36. So fit actually improved without lifedxm in the model but really the fit is indistinguishable if I use Burnham & Anderson, 2002's criteria. So my conundrum is this - The AIC and the p-values appear to disagree. The p-value would suggest that lifedxm should be included in the model and the AIC says it doesn't improve fit. In general, I tend to favor the AIC (I usually work with large sample sizes and multilevel models) but I am just beginning with survival models and I am curious if a survival model guru might suggest whether lifedxm needs to be in the model or not based on my results here? Would people generally use the AIC in this situation? Also, I can't run anova() on this models because of the random effect. I am happy to provide the data if necessary. Please cc me on a reply. Thanks, Chris Hi Chris, Did you get an answer to why the p-value and AIC score disagreed? I am getting the same results with my own data. They're not small disagreements either. The AIC score is telling me the opposite of what the p-value and the parameter coef are saying. The p-value and the coef for the predictor variable are in agreement. I've also noticed in some published papers with tables containing p-values and AIC scores that the chisq p-value and AIC score are in opposition too. This makes me think I'm missing something and that this all actually makes sense. However given that AIC = − 2 (log L) + 2K (where K is the number of parameters) and seeing as how the log-likelihood scores below are in the hundreds, shouldn't the AIC score be in the hundreds as well?? That is different than my understanding of AIC. I thought that the AIC and BIC both took as input the difference in -2LL and then adjusted those differences for the differences in number of degrees of freedom. That perspective would inform one that the absolute value of the deviance ( = -2LL) is immaterial and only differences are subject to interpretation. The p-values are calculated as Wald tests, which are basically first order approximations and have lower credence than model level comparisons. The OP also seemed to have ignored the fact that the penalized estimates of AIC and BIC went in the opposite direction from what he might have hoped. -- model0 (i
Re: [R] R.matlab package help
Hi.
On Sat, Aug 28, 2010 at 3:50 PM, michael wrote:
> Henrik,
> OK, finally I got the problem: I have an apostrophe in my
> windowns 7 user name. That mess up the file name. So I logged in using a
> guest account and it works:
>
> Received cmd: 1
> "eval" string: "B"
> B =
> -0.1347
> Sent byte: 0
> Received cmd: 1
> "eval" string: "variables = {'B'};"
> Sent byte: 0
> Received cmd: 2
> save tempname-V6 B
> answer=0
>
> Thanks a lot for your patience and help.
Good that you found a workaround.
> One final question, the variable B I got from Matlab is not just a numeric
> value, it is:
>
>> B
> $B
> [,1]
> [1,] -0.1346952
> attr(,"header")
> attr(,"header")$description
> [1] "MATLAB 5.0 MAT-file, Platform: PCWIN, Created on: Sat Aug 28 18:40:26
> 2010 "
> attr(,"header")$version
> [1] "5"
> attr(,"header")$endian
> [1] "little"
> How can I get only the numeric value, that ' -0.1346952' part?
The getVariable() function returns a named list where each element
contains the value of a particular Matlab object. You can request more
than one variable in each call. You need to use standard R list
operators to extract the element you'd like. For example:
data <- getVariable(matlab, c("B", "A"));
B <- data$B;
A <- data$A;
Thus, instead of writing:
B <- getVariable(matlab, "B");
B <- B$B;
it is less confusing if you write:
data <- getVariable(matlab, "B");
B <- data$B;
/Henrik
>
>
> Thanks again,
>
> Michael
> On Fri, Aug 27, 2010 at 7:59 AM, michael wrote:
>
>> Hi,all
>> I have a problem running R.matlab package
>> (under 2.10.1 version). I can set up the matlab server under local
>> machine(run the MatlabServer.m), "
>>
>>
>> And I can use setVariable and evaluate matlab functions in R. But when I
>> ask
>> Matlab to send the value back to R using getVariable function it
>> always returns an error:
>> "
>> ??? Error: A MATLAB string constant is not terminated properly.
>>
>> Error in ==> MatlabServer at 197
>> eval(expr);
>> "
>>
>> it seems matlab have put the data into a temporary file, so my remote
>> option is actually FALSE? (how to set it to be true?), or otherwise
>> what could be the possible problem since I can send data to matlab
>> from R but not vice versa.
>>
>>
>>
>> > sessionInfo()
>> R version 2.10.1 (2009-12-14)
>> i386-pc-mingw32
>>
>> locale:
>> [1] LC_COLLATE=English_United States.1252
>> [2] LC_CTYPE=English_United States.1252
>> [3] LC_MONETARY=English_United States.1252
>> [4] LC_NUMERIC=C
>> [5] LC_TIME=English_United States.1252
>>
>> attached base packages:
>> [1] stats graphics grDevices utils datasets methods base
>>
>> other attached packages:
>> [1] R.utils_1.5.0 R.matlab_1.3.1 R.oo_1.7.3 R.methodsS3_1.2.0
>>
>>
>> > traceback()
>> 5: file(con, open = "rb")
>> 4: readMat.default(filename)
>> 3: readMat(filename)
>> 2: getVariable.Matlab(matlab, "B")
>> 1: getVariable(matlab, "B")
>>
>> Thanks,
>>
>> Michael
>>
>
> [[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-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] Question regarding significance of a covariate in a coxme survival model
Christopher David Desjardins umn.edu> writes: > > Hi, > I am running a Cox Mixed Effects Hazard model using the library coxme. I > am trying to model time to onset (age_sym1) of thought problems (e.g. > hearing voices) (sym1). As I have siblings in my dataset, I have > decided to account for this by including a random effect for family > (famid). My covariate of interest is Mother's diagnosis where a 0 is > bipolar, 1 is control, and 2 is major depression. I am trying to fit > the following model. > > thorp1 <- coxme(Surv(age_sym1, sym1) ~ lifedxm + (1 | famid), data = > bip.surv) > > Which provides the following output: > > - > > thorp1 > Cox mixed-effects model fit by maximum likelihood >Data: bip.surv >events, n = 99, 189 (3 observations deleted due to missingness) >Iterations= 10 54 > NULL Integrated Penalized > Log-likelihood -479.0372 -467.3549 -435.2096 >Chisqdf p AIC BIC > Integrated loglik 23.36 3.00 3.3897e-05 17.36 9.58 > Penalized loglik 87.66 47.27 3.2374e-04 -6.88 -129.54 > Model: Surv(age_sym1, sym1) ~ lifedxm + (1 | famid) > Fixed coefficients > coef exp(coef) se(coef) z p > lifedxmCONTROL -1.0576257 0.3472794 0.3872527 -2.73 0.0063 > lifedxmMAJOR -0.6329957 0.5309987 0.3460847 -1.83 0.0670 > Random effects > Group Variable Std DevVariance > famid Intercept 0.980 0.9757033 > > -- > > So it appears that there is a difference in hazard rates associated with > Mother's diagnosis between BIPOLAR and CONTROL. Just to be safe, I fit a > model without this covariate present and decided to compare the models > with AIC and see if fit decreased with this covariate not in the model. > > -- >thorp1.n <- coxme(Surv(age_sym1, sym1) ~ (1 | famid), data = bip.surv) > > thorp1.n > Cox mixed-effects model fit by maximum likelihood >Data: bip.surv >events, n = 99, 189 (3 observations deleted due to missingness) >Iterations= 10 54 > NULL Integrated Penalized > Log-likelihood -479.0372 -471. -436.0478 >Chisqdf pAIC BIC > Integrated loglik 15.41 1.00 0.8663 13.41 10.81 > Penalized loglik 85.98 49.48 0.00099915 -12.97 -141.37 > Model: Surv(age_sym1, sym1) ~ (1 | famid) > Random effects > Group Variable Std Dev Variance > famid Intercept 1.050949 1.104495 > > > The problem is that the AIC for the model w/o lifedxm is 13.41 and the > model w/ lifedxm is 17.36. So fit actually improved without lifedxm in > the model but really the fit is indistinguishable if I use Burnham & > Anderson, 2002's criteria. > > So my conundrum is this - The AIC and the p-values appear to disagree. > The p-value would suggest that lifedxm should be included in the model > and the AIC says it doesn't improve fit. In general, I tend to favor the > AIC (I usually work with large sample sizes and multilevel models) but I > am just beginning with survival models and I am curious if a survival > model guru might suggest whether lifedxm needs to be in the model or not > based on my results here? Would people generally use the AIC in this > situation? Also, I can't run anova() on this models because of the > random effect. > > I am happy to provide the data if necessary. > > Please cc me on a reply. > > Thanks, > Chris > > Hi Chris, Did you get an answer to why the p-value and AIC score disagreed? I am getting the same results with my own data. They're not small disagreements either. The AIC score is telling me the opposite of what the p-value and the parameter coef are saying. The p-value and the coef for the predictor variable are in agreement. I've also noticed in some published papers with tables containing p-values and AIC scores that the chisq p-value and AIC score are in opposition too. This makes me think I'm missing something and that this all actually makes sense. However given that AIC = − 2 (log L) + 2K (where K is the number of parameters) and seeing as how the log-likelihood scores below are in the hundreds, shouldn't the AIC score be in the hundreds as well?? -- model0 (intercept only model) NULL Integrated Penalized Log-likelihood -119.8470 -117.9749 -113.1215 Chisq dfp AICBIC Integrated loglik 3.74 1.00 0.052989 1.74 0.08 Penalized loglik 13.45 7.06 0.063794 -0.68 -12.43 Random effects Group Variable Std Dev Variance SiteIntercept0.6399200 0.4094977 _ model1 (before vs after) NULL Integrated Penalized Log-likelihood -119.8470 -112.1598 -108.1663
[R] maxNR in maxLik package never stops
Greetings, I use maxNR function under maxLik package to find the REML estimates of the parameters of variance components in heteroskedastic linear regression models. I assume that in the model there is additive/multiplicative/mixed heteroskedasticity and I need estimate the respective parameters of additive/multiplicative/mixed variance components. For my research purposes I make a batch of simulations and estimate for the given design matrix but different (simulated) response vectors the variance component parameters. However, it happens very frequently, usually after the first simulation that the log-likelihood function maximization step (for which maxNR is used) R "freezes". It appears nothing wrong with R, it is working but the maximization never stops (frankly, I cannot say that it is really "never", I did not allow it more time than one day). Therefore, my simulations crash - I cannot even perform one simulation. I see two possible solutions: 1. Using some additional parameters of maxNR (such as number of iterations or something). 2. Using some function that could help me control the length of the maximization step - if it were to take to long, it would be evaluated as an unsuccessful simulation. Since I do not have much experience with R, it is clear to me that I cannot resolve the thing myself. Maybe there are other viable solutions. I shall be much obliged for your help and advice. Martin BoÄa (Matej Bel University / Faculty of Natural Sciences / Banská Bystrica / Slovakia). [[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] R.matlab package help
OK, I think it is just B[1].
Thanks!
Michael
On Sat, Aug 28, 2010 at 6:50 PM, michael wrote:
> Henrik,
>OK, finally I got the problem: I have an apostrophe in my
> windowns 7 user name. That mess up the file name. So I logged in using a
> guest account and it works:
>
>
> Received cmd: 1
> "eval" string: "B"
> B =
>-0.1347
> Sent byte: 0
> Received cmd: 1
> "eval" string: "variables = {'B'};"
> Sent byte: 0
> Received cmd: 2
> save tempname-V6 B
> answer=0
>
> Thanks a lot for your patience and help.
> One final question, the variable B I got from Matlab is not just a numeric
> value, it is:
>
> > B
> $B
>[,1]
> [1,] -0.1346952
> attr(,"header")
> attr(,"header")$description
> [1] "MATLAB 5.0 MAT-file, Platform: PCWIN, Created on: Sat Aug 28 18:40:26
> 2010 "
> attr(,"header")$version
> [1] "5"
> attr(,"header")$endian
> [1] "little"
> How can I get only the numeric value, that ' -0.1346952' part?
>
>
> Thanks again,
>
> Michael
>
>
>
[[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] R.matlab package help
Henrik,
OK, finally I got the problem: I have an apostrophe in my
windowns 7 user name. That mess up the file name. So I logged in using a
guest account and it works:
Received cmd: 1
"eval" string: "B"
B =
-0.1347
Sent byte: 0
Received cmd: 1
"eval" string: "variables = {'B'};"
Sent byte: 0
Received cmd: 2
save tempname-V6 B
answer=0
Thanks a lot for your patience and help.
One final question, the variable B I got from Matlab is not just a numeric
value, it is:
> B
$B
[,1]
[1,] -0.1346952
attr(,"header")
attr(,"header")$description
[1] "MATLAB 5.0 MAT-file, Platform: PCWIN, Created on: Sat Aug 28 18:40:26
2010 "
attr(,"header")$version
[1] "5"
attr(,"header")$endian
[1] "little"
How can I get only the numeric value, that ' -0.1346952' part?
Thanks again,
Michael
On Fri, Aug 27, 2010 at 7:59 AM, michael wrote:
> Hi,all
> I have a problem running R.matlab package
> (under 2.10.1 version). I can set up the matlab server under local
> machine(run the MatlabServer.m), "
>
>
> And I can use setVariable and evaluate matlab functions in R. But when I
> ask
> Matlab to send the value back to R using getVariable function it
> always returns an error:
> "
> ??? Error: A MATLAB string constant is not terminated properly.
>
> Error in ==> MatlabServer at 197
> eval(expr);
> "
>
> it seems matlab have put the data into a temporary file, so my remote
> option is actually FALSE? (how to set it to be true?), or otherwise
> what could be the possible problem since I can send data to matlab
> from R but not vice versa.
>
>
>
> > sessionInfo()
> R version 2.10.1 (2009-12-14)
> i386-pc-mingw32
>
> locale:
> [1] LC_COLLATE=English_United States.1252
> [2] LC_CTYPE=English_United States.1252
> [3] LC_MONETARY=English_United States.1252
> [4] LC_NUMERIC=C
> [5] LC_TIME=English_United States.1252
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] R.utils_1.5.0 R.matlab_1.3.1R.oo_1.7.3R.methodsS3_1.2.0
>
>
> > traceback()
> 5: file(con, open = "rb")
> 4: readMat.default(filename)
> 3: readMat(filename)
> 2: getVariable.Matlab(matlab, "B")
> 1: getVariable(matlab, "B")
>
> Thanks,
>
> Michael
>
[[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] extracting columns
On Aug 28, 2010, at 1:24 PM, Laetitia Schmid wrote: Hi, Can anybody show me how to extract all columns in my dataset that are polymorphic? Or phrased in another way I would like to delete all columns that have no more than one letter in it (that are monomorphic). Assuming you read this in with read.table: > table(sapply(txt.in, function(x) length(levels(x 0 1 2 3 4 18 112 31 10 1 I am further assuming you do not want the 18 columns with no levels (presumably numeric or logical, ooops only ones with all "T"'s, and got converted to logical, so no matter ) then: > names(txt.in)[ sapply(txt.in, function(x) length(levels(x))) >1] [1] "V1" "V2" "V5" "V23" "V27" "V29" "V30" "V34" "V45" "V47" "V48" "V49" [13] "V50" "V53" "V56" "V62" "V63" "V64" "V66" "V67" "V71" "V73" "V81" "V90" [25] "V95" "V104" "V107" "V109" "V114" "V116" "V121" "V124" "V126" "V131" "V133" "V139" [37] "V156" "V158" "V163" "V166" "V168" "V172" so: newdata <- txt.in[ , ] -- David Winsemius, MD West Hartford, CT __ 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] extracting columns
Wow. That was fast. And it works. Thank you! Laetitia Am 29.08.2010 um 00:35 schrieb Jorge Ivan Velez: > index <- apply(d, 2, function(x) length(table(x)) > 1) > d[, index] [[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] extracting columns
Hi, Can anybody show me how to extract all columns in my dataset that are polymorphic? Or phrased in another way I would like to delete all columns that have no more than one letter in it (that are monomorphic). Thank you. Laetitia "V1" "V2" "V3" "V4" "V5" "V6" "V7" "V8" "V9" "V10" "V11" "V12" "V13" "V14" "V15" "V16" "V17" "V18" "V19" "V20" "V21" "V22" "V23" "V24" "V25" "V26" "V27" "V28" "V29" "V30" "V31" "V32" "V33" "V34" "V35" "V36" "V37" "V38" "V39" "V40" "V41" "V42" "V43" "V44" "V45" "V46" "V47" "V48" "V49" "V50" "V51" "V52" "V53" "V54" "V55" "V56" "V57" "V58" "V59" "V60" "V61" "V62" "V63" "V64" "V65" "V66" "V67" "V68" "V69" "V70" "V71" "V72" "V73" "V74" "V75" "V76" "V77" "V78" "V79" "V80" "V81" "V82" "V83" "V84" "V85" "V86" "V87" "V88" "V89" "V90" "V91" "V92" "V93" "V94" "V95" "V96" "V97" "V98" "V99" "V100" "V101" "V102" "V103" "V104" "V105" "V106" "V107" "V108" "V109" "V110" "V111" "V112" "V113" "V114" "V115" "V116" "V117" "V118" "V119" "V120" "V121" "V122" "V123" "V124" "V125" "V126" "V127" "V128" "V129" "V130" "V131" "V132" "V133" "V134" "V135" "V136" "V137" "V138" "V139" "V140" "V141" "V142" "V143" "V144" "V145" "V146" "V147" "V148" "V149" "V150" "V151" "V152" "V153" "V154" "V155" "V156" "V157" "V158" "V159" "V160" "V161" "V162" "V163" "V164" "V165" "V166" "V167" "V168" "V169" "V170" "V171" "V172" "V1" "G" "C" "C" "G" "G" "C" "A" "G" "C" "C" "A" "G" "A" "C" "A" "G" "G" "C" "G" "T" "C" "C" "A" "G" "A" "C" "G" "G" "A" "A" "G" "T" "G" "C" "G" "C" "G" "C" "C" "C" "C" "G" "C" "T" "C" "A" "G" "G" "G" "G" "T" "A" "G" "T" "C" "G" "G" "C" "T" "C" "G" "C" "C" "G" "T" "C" "G" "A" "A" "A" "G" "A" "C" "G" "C" "A" "C" "C" "A" "C" "T" "G" "G" "G" "C" "A" "A" "T" "C" "C" "G" "T" "A" "G" "C" "A" "G" "G" "C" "T" "T" "G" "C" "A" "T" "G" "G" "G" "C" "G" "A" "C" "G" "T" "T" "C" "A" "T" "G" "G" "G" "A" "G" "C" "G" "G" "C" "C" "A" "T" "T" "T" "G" "G" "C" "G" "A" "C" "T" "C" "G" "T" "A" "G" "A" "A" "C" "A" "G" "G" "A" "G" "A" "G" "G" "G" "G" "C" "A" "G" "C" "A" "G" "G" "C" "G" "A" "C" "C" "A" "T" "A" "V2" "C" "C" "C" "G" "G" "C" "A" "G" "C" "C" "A" "G" "A" "C" "A" "G" "G" "C" "G" "T" "C" "C" "A" "G" "A" "C" "C" "G" "A" "A" "G" "T" "G" "C" "G" "C" "G" "C" "C" "C" "C" "G" "C" "T" "C" "A" "G" "C" "G" "G" "T" "A" "G" "T" "C" "G" "G" "C" "T" "C" "G" "A" "C" "G" "T" "C" "G" "A" "A" "A" "C" "A" "C" "G" "C" "A" "C" "C" "A" "C" "T" "G" "G" "G" "C" "A" "A" "T" "C" "C" "G" "T" "A" "G" "C" "A" "G" "G" "C" "T" "T" "G" "C" "A" "T" "G" "G" "G" "C" "G" "A" "C" "G" "T" "T" "C" "A" "T" "G" "G" "T" "A" "G" "C" "G" "G" "C" "C" "A" "T" "T" "T" "G" "G" "C" "G" "A" "C" "T" "C" "G" "T" "A" "G" "A" "A" "C" "A" "G" "G" "A" "G" "A" "G" "G" "G" "G" "C" "A" "G" "C" "A" "G" "G" "C" "G" "A" "C" "C" "A" "T" "A" "V3" "C" "C" "C" "G" "G" "C" "A" "G" "C" "C" "A" "G" "A" "C" "A" "G" "G" "C" "G" "T" "C" "C" "A" "G" "A" "C" "C" "G" "A" "A" "G" "T" "G" "C" "G" "C" "G" "C" "C" "C" "C" "G" "C" "T" "C" "A" "G" "C" "G" "G" "T" "A" "G" "T" "C" "G" "G" "C" "T" "C" "G" "G" "C" "G" "T" "C" "G" "A" "A" "A" "C" "A" "C" "G" "C" "A" "C" "C" "A" "C" "C" "G" "G" "G" "C" "A" "A" "T" "C" "C" "G" "T" "A" "G" "C" "A" "G" "G" "C" "T" "T" "G" "C" "A" "T" "G" "G" "G" "C" "G" "A" "C" "G" "T" "T" "C" "A" "T" "G" "G" "T" "A" "G" "C" "G" "G" "C" "C" "A" "T" "T" "T" "C" "G" "C" "G" "A" "C" "T" "C" "G" "T" "A" "G" "A" "A" "C" "A" "G" "G" "A" "G" "A" "G" "G" "G" "G" "C" "A" "G" "C" "A" "G" "G" "C" "G" "A" "C" "C" "A" "T" "A" "V4" "C" "C" "C" "G" "G" "C" "A" "G" "C" "C" "A" "G" "A" "C" "A" "G" "G" "C" "G" "T" "C" "C" "A" "G" "A" "C" "C" "G" "A" "A" "G" "T" "G" "C" "G" "C" "G" "C" "C" "C" "C" "G" "C" "T" "C" "A" "G" "C" "G" "G" "T" "A" "G" "T" "C" "G" "G" "C" "T" "C" "G" "G" "C" "G" "T" "C" "G" "A" "A" "A" "C" "A" "C" "G" "C" "A" "C" "C" "A" "C" "G" "G" "G" "G" "C" "A" "A" "T" "C" "C" "G" "T" "A" "G" "C" "A" "G" "G" "C" "T" "T" "G" "C" "A" "T" "G" "G" "G" "C" "G" "A" "C" "G" "T" "T" "C" "A" "T" "G" "G" "A" "A" "G" "C" "G" "G" "C" "C" "A" "T" "T" "T" "C" "G" "C" "G" "A" "C" "T" "C" "G" "T" "A" "G" "A" "A" "C" "A" "G" "G" "A" "G" "A" "G" "G" "G" "G" "C" "A" "G" "C" "A" "G" "G" "C" "G" "A" "A" "C" "A" "T" "C" "V5" "A" "C" "C" "G" "G" "C" "A" "G" "C" "C" "A" "G" "A" "C" "A" "G" "G" "C" "G" "T" "C" "C" "A" "G" "A" "C" "C" "G" "A" "A" "G" "T" "G" "C" "G" "C" "G" "C" "C" "C" "C" "G" "C" "T" "C" "A" "G" "C" "G" "G" "T" "A" "G" "T" "C" "G" "G" "C" "T" "C" "G" "G" "C" "G" "T" "C" "G" "A" "A" "A" "C" "A" "C" "G" "C" "A" "C" "C" "A" "C" "G" "G" "G" "G" "C" "A" "A" "T" "C" "C" "G" "T" "A" "G" "C" "A" "G" "G" "C" "T" "T" "G" "C" "A" "T" "G" "G" "G" "C" "G" "A" "C" "G" "T" "T" "C" "A" "T" "G" "G" "A" "A" "G" "C" "G" "G" "C" "C" "A" "T" "T" "T" "C" "G" "C" "G" "A" "C" "T" "C" "G" "T" "A" "G" "A" "A" "C" "A" "G" "G" "A" "G" "A" "G" "G" "G" "G" "C" "A" "G" "C" "A" "G" "G" "C" "G" "A" "A" "C" "A" "T" "C" "V6" "A" "C" "C" "G" "C" "C" "A" "G" "C" "C" "A" "G" "A" "C" "A" "G" "G" "C" "G" "T" "C" "C" "A" "G" "A" "C" "G" "G" "A" "A" "G" "T" "G" "G" "
Re: [R] conditionally sum
See ?table for details. > tmp <- ' yx + 1 1 2008 + 2 1 2008 + 3 0 2008 + 4 0 2009 + 5 1 2009' > data <- read.table(textConnection(tmp), header=TRUE) > data yx 1 1 2008 2 1 2008 3 0 2008 4 0 2009 5 1 2009 > table(data$x, data$y) 0 1 2008 1 2 2009 1 1 > __ 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] binomial distribution
Perfect! Thank You! 2010/8/28 Peter Dalgaard > On 08/28/2010 10:23 PM, tamas barjak wrote: > > Hello! > > > > I need some help. > > How I know it to draw the formula of the binomial distribution? > > > > expr<-expression(P(xi == k) == choose(n, k)* p^k*(1-p)^(n-k)) ---> not > good > > > > on the screen the "choose(n, k)" not the Binomial Formula, but "choose(n, > > k)" > > Check demo(plotmath), 4th-to-last example. > > -- > Peter Dalgaard > Center for Statistics, Copenhagen Business School > Phone: (+45)38153501 > Email: pd@cbs.dk Priv: pda...@gmail.com > [[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] nlme help
I would like to fit a nonlinear model with "nlme."
The model is a nonlinear growth model with five time points: y = X*b +
e, where design matrix X is defined as
X= | 1 0 |
| 1 1 |
| 1 k |
| 1 2k |
| 1 3k |
and parameter vector b = (b0, b1). And where "k" is a parameter to be
estimated. Of course I also want to estimate the intercept (b0), slope
(b1). Error variances (e) with 5 free parameters and random effects
covariance matrix (2x2: for b0 and b1).
I am not certain how to get the design matrix in the code so that the
"nlme" procedure will be able to
estimate parameter "k" as well as the intercept and slope.
I have tried the following:
meanfunc <- function(x,b0,b1,k){
meangrad <- array(0,c(length(x),2),list(NULL,c("b0","b1")))
for(i in 1:length(x)){
if(x[i]==0){
err[i] <- b0 + b1*x[i]
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- x[i]}
if(x[i]==1){
err[i] <- b0 + b1*x[i]
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- x[i]}
if(x[i]==2){
err[i] <- b0 + b1*(x[i]-1)*k
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- (x[i]-1)*k}
if(x[i]==3){
err[i] <- b0 + b1*2*(x[i]-1)*k
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- (x[i]-1)*k}
if(x[i]==4){
err[i] <- b0 + b1*3*(x[i]-1)*k
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- (x[i]-1)*k}
}
# compute analytical derivatives
attr(err,"gradient") <- meangrad
err
}
data.mlfit <- nlme(y ~ meanfunc(x,b0,b1,k),
fixed=list(b0 ~ 1, b1 ~ 1, k ~ 1),
random=list(b0 ~ 1, b1 ~ 1),
groups = ~subj,
data=nd,
start=list(fixed=c(300,50,1)),
method="ML",verbose=T)
I get the following error message:
Error in meangrad[i, "b1"] <- (x[i] - 1) * k :
number of items to replace is not a multiple of replacement length
In addition: Warning messages:
1: In err[i] <- b0 + b1 * x[i] :
number of items to replace is not a multiple of replacement length
2: In err[i] <- b0 + b1 * x[i] :
number of items to replace is not a multiple of replacement length
3: In err[i] <- b0 + b1 * (x[i] - 1) * k :
number of items to replace is not a multiple of replacement length
Any suggestions or insights would be greatly appreciated.
Thank you,
Jeff
--
**
Jeffrey R. Harring, Assistant Professor
Department of Measurement, Statistics & Evaluation (EDMS)
1230 Benjamin Building
University of Maryland
College Park, MD 20742-1115
Phone: 301.405.3630
Fax:301.314.9245
Email: harr...@umd.edu
Web:http://www.education.umd.edu/EDMS/fac/Harring/webpage.html
__
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] conditionally sum
Hi, Suppose I have a dataframe like: > data yx 1 1 2008 2 1 2008 3 0 2008 4 0 2009 5 1 2009 Now I want to know how much "1" corresponds with 2008 and how much with 2009, the answers are 3 and 1. How can I achieve this? Thanks for the anwer. André [[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] R.matlab package help
Michael, could you please give me *verbatim* details on what messages you are seeing. In your previous reply you did *not* report seeing "save(tmpname, '-V6', 'B');" and now you say you get it. Please do not abbreviate what you are getting (e.g. "save(temp,-v6,B)"), because that will not be useful and only confuse me more. Can you please show me exactly what you do and what you get in both Matlab and R, otherwise there is no way I can help you solve this. I also do not know whether or not you installed that new version that I suggest. You are the first one having this problem and it has to do with you having an apostrophe in your Windows user name. /Henrik On Sat, Aug 28, 2010 at 2:22 PM, michael wrote: > Henrik, > > Yes I replaced the MatlabServer.m and restarted Matlab and R, it seems the > problem is still there. I see save(temp,-v6,B) there too. I use > MatlabServer.m mannually in Matlab, actually I couldn't use it in R, Matlab > will pop up a command window and never respond. > > > Thanks, > > Michael > > [[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-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] expression() and plot title
On Aug 28, 2010, at 4:48 PM, Sancar Adali wrote:
My function is like this
sim.res<-gaussian_simulation(p=3, r=4, q=3, c=0.1,d=2,
Wchoice = "avg",
pre.scaling = TRUE,
oos = TRUE,
alpha = NULL,
n = 100, m = 100, nmc = 100)
which is defined as
gaussian_simulation <- function(p, r, q, c,
d = p-1,
pprime1 = p+q, # cca arguments
pprime2 = p+q, # cca arguments
Wchoice = "avg",
pre.scaling = TRUE,
oos = TRUE,
alpha = NULL,
n = 100, m = 100, nmc = 100)
Not much of a definition...no body???
and I want to title the plot after I invoke the gaussian_simulation
function
sim.res<-gaussian_simulation(p=3, r=4, q=3, c=0.1,d=2,
Wchoice = "avg",
pre.scaling = TRUE,
oos = TRUE,
alpha = NULL,
n = 100, m = 100, nmc = 100)
plot(sim.res)
title("p=3, r=4, q=3, c=0.1,d=2")
The following assumes you want the values of those parameters for the
title to be picked up from the function's arguments, because the code
you offer for the title is syntactically correct. You first need to re-
read Bill Venables response regarding bquote which you show no signs
of adapting, and then you need to study:
?plotmath
... where you will learn the proper syntax for "=" ("==") and for
comma-separated lists. You should not need the quotes inside the
bquote call. The "~" (space) and "*" (no space) are syntactic
separators inside expressions (in the limited scope of the plotmath
perspective). The plotmath interpretation of these is not well
described in the documentation (IMO). The help page could easily lead
you to think that commas or spaces might be valid separators, or that
list() was the same as list() in regular R, and it's only by repeated
errors and reading worked examples on r-help that I have learned
otherwise.
Consider this and apply lessons learned:
plot(1,1); p=3; r=4; q=3; c=0.1; d=2
title(main= bquote( list(p==.(p), r==.(r), q==.(q), c==.(c), d==.
(d)) ) )
# you might get away with not using the "main" argument name, but I
think it's bad practice.
--
David.
On Sat, Aug 28, 2010 at 12:53 AM, Sancar Adali
wrote:
What I want to do is put the arguments I supply to a function into
the
title of a plot
Say I'm calling func.1
func.1(a=4,b=4)
plot(,..., title("a=4, b=4"))
If I'm calling func.1 with different arguments, I want the plot
title to
reflect that.
A small detail is that func.1 might have an argument with a default
like
c=a+b. I tried using expression but couldn't get it to work.
Is there a way to do this using expression() ?
--
Sancar Adali
David Winsemius, MD
West Hartford, CT
__
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] R.matlab package help
Henrik, Yes I replaced the MatlabServer.m and restarted Matlab and R, it seems the problem is still there. I see save(temp,-v6,B) there too. I use MatlabServer.m mannually in Matlab, actually I couldn't use it in R, Matlab will pop up a command window and never respond. Thanks, Michael [[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] expression() and plot title
My function is like this
sim.res<-gaussian_simulation(p=3, r=4, q=3, c=0.1,d=2,
Wchoice = "avg",
pre.scaling = TRUE,
oos = TRUE,
alpha = NULL,
n = 100, m = 100, nmc = 100)
which is defined as
gaussian_simulation <- function(p, r, q, c,
d = p-1,
pprime1 = p+q, # cca arguments
pprime2 = p+q, # cca arguments
Wchoice = "avg",
pre.scaling = TRUE,
oos = TRUE,
alpha = NULL,
n = 100, m = 100, nmc = 100)
and I want to title the plot after I invoke the gaussian_simulation function
sim.res<-gaussian_simulation(p=3, r=4, q=3, c=0.1,d=2,
Wchoice = "avg",
pre.scaling = TRUE,
oos = TRUE,
alpha = NULL,
n = 100, m = 100, nmc = 100)
plot(sim.res)
title("p=3, r=4, q=3, c=0.1,d=2")
On Sat, Aug 28, 2010 at 12:53 AM, Sancar Adali wrote:
> What I want to do is put the arguments I supply to a function into the
> title of a plot
> Say I'm calling func.1
> func.1(a=4,b=4)
> plot(,..., title("a=4, b=4"))
> If I'm calling func.1 with different arguments, I want the plot title to
> reflect that.
> A small detail is that func.1 might have an argument with a default like
> c=a+b. I tried using expression but couldn't get it to work.
>
> Is there a way to do this using expression() ?
>
> --
> Sancar Adali
>
>
--
Sancar Adali
Johns Hopkins University
Graduate Student
[[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] binomial distribution
On 08/28/2010 10:23 PM, tamas barjak wrote: > Hello! > > I need some help. > How I know it to draw the formula of the binomial distribution? > > expr<-expression(P(xi == k) == choose(n, k)* p^k*(1-p)^(n-k)) ---> not good > > on the screen the "choose(n, k)" not the Binomial Formula, but "choose(n, > k)" Check demo(plotmath), 4th-to-last example. -- Peter Dalgaard Center for Statistics, Copenhagen Business School Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@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] binomial distribution
try: ?pbinom -- View this message in context: http://r.789695.n4.nabble.com/binomial-distribution-tp2398568p2398705.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.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R.matlab package help
Hi.
On Fri, Aug 27, 2010 at 7:08 PM, michael wrote:
> Henrik,
> The line before that is:
>
> Received cmd: 2
> save
> C:\Users\FAN'S~1\AppData\Local\Temp\tpe2b4012b_f9ed_402d_af0f_f21ebd8116a6.mat
> -V6 B
You should see have seen something like:
save(tmpname, '-V6', 'B');
For some reason Matlab is not picking up the patched version of MatlabServer.m.
Q1. Did you download that file and *replaced* the existing
MatlabServer.m with the downloaded one?
Q2. Did you shut down Matlab/the MatlabServer after replacing MatlabServer.m?
Q3. How do you start the MatlabServer? Manually or via R?
If you don't manage to figure out how you replace MatlabServer.m, then
try to install the preliminary R.matlab v1.3.2:
source("http://www.braju.com/R/hbLite.R";);
installPackages("http://www.braju.com/R/patches/R.matlab/20100827/R.matlab_1.3.2.zip";);
and then restart R and the MatlabServer.
/Henrik
> [[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-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] binomial distribution
Hello! I need some help. How I know it to draw the formula of the binomial distribution? expr<-expression(P(xi == k) == choose(n, k)* p^k*(1-p)^(n-k)) ---> not good on the screen the "choose(n, k)" not the Binomial Formula, but "choose(n, k)" Thanx! [[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] AIC using nls function
Bert Gunter gene.com> writes: > > John: > > 1. As always, and as requested (see posting guide), a small > reproducible example might help. Bert is right that things aren't well defined. However, AIC is still *widely* used for nonlinear models. For the sloppy folks among us, here are some useful (?) pieces of information: 1. nls counts the variance of the residual error as a parameter 2. As long as you compute AIC *within the same framework* (e.g. comparing nls fits to each other, or glm fits, or nlme fits ...) these decisions are generally made consistently. Comparing across model types requires attention to detail to make sure that parameters are counted using similar rules, and that additive constants are consistently included or not. __ 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] R.matlab package help
michael gmail.com> writes:
>
> C:\Users\FAN'S~1\AppData\Local\Temp\tpe2b4012b_f9ed_402d_af0f_f21ebd8116a6.mat
> -V6 B
I bet the problem is with the single quote (') in the path.
__
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] Non-standard sorts on vectors
The following code shows that the unlisted data frame assigns an index to
each member. When one sorts the data frame based on ULST, he in fact uses
the (implicit) indices of ULST but not the actual values! Therefore, your
guess is correct.
> test.vec=data.frame(c("A","F","C"))
> ULST=unlist(test.vec)
> ULST
c..AFC..1 c..AFC..2 c..AFC..3
A F C
Levels: A C F
> sort(as.vector(ULST),index.return=T)
$x
[1] "A" "C" "F"
$ix
[1] 1 3 2
# also the rank and he index are different
> rank(ULST)
c..AFC..1 c..AFC..2 c..AFC..3
1 3 2
--
View this message in context:
http://r.789695.n4.nabble.com/Non-standard-sorts-on-vectors-tp2340431p2384591.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.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] About plot graphs
Gavin gave some problems with relying attaching data, here is another example,
somewhat artificial, but not unrealistic (I had similar to this happen to me
before I learned better):
attach(women)
# do some stuff
library(lattice)
attach(singer)
# do some more stuff
# now we want to go back and look at the women data
plot(weight,height)
#or even worse
detach()
attach(singer[1:15,])
plot(weight,height)
# what conclusions do we draw from the plot?
detach()
detach()
--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111
> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of Stephen Liu
> Sent: Friday, August 27, 2010 11:14 PM
> To: gavin.simp...@ucl.ac.uk
> Cc: r-help@r-project.org
> Subject: Re: [R] About plot graphs
>
> Hi Gavin,
>
>
> Thanks for your advice and the examples explaining plotting settings.
>
> The steps on your examples work on my test.
>
>
> > 2) Don't attach() things, you are asking for trouble
>
> > If a function has a formula method (which plot does) then use it like
> > this: plot(Draft_No. ~ Day_of_year, data = Test01)
>
> > If the function doesn't have a formula method, then wrap it in a
> > with()
> > call:
>
> > with(Test01, plot(Day_of_year, Draft_No.))
>
> > No need for attach.
>
>
> Noted and thanks. What will be the problem caused by "attach()"?
>
>
> > dev.new(height = 6, width = 12)
> > layout(matrix(1:2, ncol = 2))
> > op <- par(pty = "s") ## this is the important bit
> > plot(runif(100), rnorm(100))
> > plot(runif(100), rnorm(100), col = "red")
> > par(op) ## now reset the pars
> > layout(1)
>
> What is the function of layout(1) ? Tks
>
>
> B.R.
> satimis
>
>
>
>
> - Original Message
> From: Gavin Simpson
> To: Stephen Liu
> Cc: "r-help@r-project.org"
> Sent: Fri, August 27, 2010 5:38:40 PM
> Subject: Re: [R] About plot graphs
>
> On Fri, 2010-08-27 at 02:05 -0700, Stephen Liu wrote:
> > Hi Gavin,
> >
> > Thanks for your advice which works for me.
> >
> >
> > (rectangular window)
> > dev.new(height = 6, width = 12)
> > layout(matrix(1:2, nrow=1))
> > plot(Test01$Day_of_year, Test01$Draft_No.)
> > attach(Test01)
> > plot(Day_of_year,Draft_No.)
>
> 1) I can't reproduce this; where/what is Test01? But don;t bother
> sending, see my example below
> 2) Don't attach() things, you are asking for trouble
>
> If a function has a formula method (which plot does) then use it like
> this: plot(Draft_No. ~ Day_of_year, data = Test01)
>
> If the function doesn't have a formula method, then wrap it in a with()
> call:
>
> with(Test01, plot(Day_of_year, Draft_No.))
>
> No need for attach.
>
> >
> > (rectangular window in vertical position)
> > dev.new(height = 12, width = 4)
> > layout(matrix(1:2, nrow=2))
> > plot(Test01$Day_of_year, Test01$Draft_No.)
> > plot(Day_of_year,Draft_No.)
> >
> > (height = 12, width = 6) can't work. The graphs plotted are
> distorted off
> > square shape. I must reduce "width = 4"
> >
> > Why? TIA
>
> Because you don't appreciate that the dimensions of the device are not
> the same as the dimensions of the plotting region *on* the device. Most
> notably, the margins on the device are given by par("mar") for example
> and are not square:
>
> > par("mar")
> [1] 5.1 4.1 4.1 2.1
>
> So more space is set aside on the bottom then anywhere else, and the
> margin on the right is quite small.
>
> You have already been provided with an answer that you dismissed
> because
> you didn't appear to appreciate what you were being told.
>
> Compare this:
>
> dev.new(height = 6, width = 12)
> layout(matrix(1:2, ncol = 2))
> plot(runif(100), rnorm(100))
> plot(runif(100), rnorm(100), col = "red")
> layout(1)
>
> with this:
>
> dev.new(height = 6, width = 12)
> layout(matrix(1:2, ncol = 2))
> op <- par(pty = "s") ## this is the important bit
> plot(runif(100), rnorm(100))
> plot(runif(100), rnorm(100), col = "red")
> par(op) ## now reset the pars
> layout(1)
>
> So now the regions are square, we have the asymmetric margins like all
> R
> plots and we have drawn this on a device that has ~ appropriate
> dimensions.
>
> If you want to fiddle more with this, then you'll need to make the
> margins equal, but you don't want to do that really as you need more
> space in certain areas to accommodate axis labels and tick labels etc.
>
> For the vertical one, this works for me, though note that because of
> the
> margins, pty = "s" is making the individual plotting regions smaller to
> respect the square plotting region you asked for.
>
> dev.new(height = 12, width = 6)
> layout(matrix(1:2, ncol = 1))
> op <- par(pty = "s") ## this is the important bit
> plot(runif(100), rnorm(100))
> plot(runif(100), rnorm(100), col = "red")
> par(op) ## now reset the pars
> layout(1)
>
> This is because you have 5.1 plus 4.1 lines of margin in the vertical
> direction per plot (so 18.4 lines in total) versus 4.1 + 2.1 = 6.2
> lines
> of margin in the horizo
Re: [R] how to un-crosstabulate data without changing numeric values to text?
On Aug 28, 2010, at 11:31 AM, David Winsemius wrote:
On Aug 28, 2010, at 6:07 AM, Nevil Amos wrote:
I have a large amount of data read in from over 140 excel files in
the format of x. r1 to r5 are repeat measures for a given
Wavelength and ANWC_NO.
I need to rearrange x into 3 columns, ANWC_NO,Wavelegth, value ie
ANWC_NOWavelength r1
ANWC_NOWavelength,r2
ANWC_NOWavelength r3
etc...
I can rearrange the data using the code below, however all the
columns end up as strings, not numeric values. I cannot then
summaries the data, ( whcih I need to do in bins of wavelanght for
each ANWC_NO)
> x
Wavelength r1 r2 r3 r4 r5 ANWC_NO
1300 0.003126 0.005382 0.001094 0.012529 0.005632 39239
2302 0.004924 0.006280 0.002366 0.015234 0.006204 39239
3304 0.004769 0.005960 0.002759 0.015856 0.006804 39239
4306 0.005181 0.006717 0.004033 0.017380 0.007675 39239
5308 0.005872 0.008083 0.004429 0.018334 0.008504 39239
6310 0.007164 0.010775 0.005949 0.019952 0.009594 39239
Here are two other ways (the first perhaps less explicit and relying
on argument recycling in the data.frame function for repeating the
values for the first and last columns):
> data.frame(Wave=x$Wavelength,
ANWC =x$ANWC_NO,
values = unlist( x[ ,grep("^r",names(x))] ) )
And of course:
> require(reshape)
> melt(x, measure.vars=2:6)
The last one arguable the cleanest.
--
David.
Try:
reshape(x, idvar=c("Wavelength", "ANWC_NO"),
varying=2:6, v.names="r",
direction="long")
I think that melt in package reshape would also work.
--
David
> y =NULL
> rows<-nrow(x)
> for(r in 1:rows){
+ for(c in 2:6){
+ row<-c(c(x[r,7]),as.numeric(c(x[r,1])),as.numeric(c(x[r,c])))
+ y<-rbind(y,row)
+ }}
> colnames(y)<-c("ANWC_NO","WAVELENGTH","VALUE")
> head (y)
ANWC_NO WAVELENGTH VALUE
row "39239" "300" "0.003126"
row "39239" "300" "0.005382"
row "39239" "300" "0.001094"
row "39239" "300" "0.012529"
row "39239" "300" "0.005632"
row "39239" "302" "0.004924"
> mean(y$VALUE)
Error in y$VALUE : $ operator is invalid for atomic vectors
how do I get the data arranged in three columns but maintaining
WavelENGTH and the values as numeric in a data.frame?
Many thanks
Nevil Amos
Monash University
__
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] About plot graphs
On Aug 28, 2010, at 9:15 AM, Stephen Liu wrote:
Hi Gavin,
Lot of thank for your detail explanation.
Just looked at;
?with
?within
Both "Evaluate an Expression in a Data Environment"
Usage: Both;
with(data, expr, ...)
within(data, expr, ...)
Details:
.
‘within’ is similar, except that it examines the environment after
the evaluation of ‘expr’ and makes the corresponding modifications
to ‘data’ (this may fail in the data frame case if objects are
created which cannot be stored in a data frame), and returns it.
‘within’ can be used as an alternative to ‘transform’.
What does it mean ". examines the environment after the
evaluation of 'expr'
I take it to mean that the results of expr (applied to "data" and
anything else used as arguments) are then used to update "data", but
only if possible, i.e., if those results are congruent with the
structure of "data". I read that phrase as covering the "if possible"
assessment process.
."
Tks
B.R.
Stephen
B.R.
Stephen L
- Original Message
From: Gavin Simpson
To: Stephen Liu
Cc: "r-help@r-project.org"
Sent: Sat, August 28, 2010 3:28:34 PM
Subject: Re: [R] About plot graphs
On Fri, 2010-08-27 at 22:14 -0700, Stephen Liu wrote:
Hi Gavin,
Thanks for your advice and the examples explaining plotting settings.
The steps on your examples work on my test.
2) Don't attach() things, you are asking for trouble
If a function has a formula method (which plot does) then use it
like
this: plot(Draft_No. ~ Day_of_year, data = Test01)
If the function doesn't have a formula method, then wrap it in a
with()
call:
with(Test01, plot(Day_of_year, Draft_No.))
No need for attach.
Noted and thanks. What will be the problem caused by "attach()"?
If you change the underlying data, this is not reflected in the
attached
copy, because it is just that, a "copy"[1] created at the point at
which
you attached the object. E.g.
## Some data, which we attach
dat <- data.frame(A = 1:10, B = letters[1:10])
attach(dat)
## Look at A
A
## Change or dat object by altering the A component
dat <- within(dat, A <- LETTERS[1:10])
## Look at A
A
## Look at what A really is
with(dat, A)
Using with() and within() etc have two key advantages over attach: i)
only one version of the data/object exists, ii) the intention of code
using:
with(dat, "do something with A")
is much more clear than
"do something with A"
A could be anything, anywhere. More info is on the ?attach help page.
[1] ?attach contains the details of what I mean by "copy"
dev.new(height = 6, width = 12)
layout(matrix(1:2, ncol = 2))
op <- par(pty = "s") ## this is the important bit
plot(runif(100), rnorm(100))
plot(runif(100), rnorm(100), col = "red")
par(op) ## now reset the pars
layout(1)
What is the function of layout(1) ? Tks
The opposite of
layout(matrix(1:4, ncol = 2))
for example. It ( layout(1) ) says create a layout with a single
plotting region. So we use it to reset the current device back to
normal. I find it is good working practice to tidy up after doing
plots
like this. In the code above, we change both the layout() and the
plotting parameters (via par() ), and the last two lines of code in
that
example reset these changes.
G
B.R.
satimis
- Original Message
From: Gavin Simpson
To: Stephen Liu
Cc: "r-help@r-project.org"
Sent: Fri, August 27, 2010 5:38:40 PM
Subject: Re: [R] About plot graphs
On Fri, 2010-08-27 at 02:05 -0700, Stephen Liu wrote:
Hi Gavin,
Thanks for your advice which works for me.
(rectangular window)
dev.new(height = 6, width = 12)
layout(matrix(1:2, nrow=1))
plot(Test01$Day_of_year, Test01$Draft_No.)
attach(Test01)
plot(Day_of_year,Draft_No.)
1) I can't reproduce this; where/what is Test01? But don;t bother
sending, see my example below
2) Don't attach() things, you are asking for trouble
If a function has a formula method (which plot does) then use it like
this: plot(Draft_No. ~ Day_of_year, data = Test01)
If the function doesn't have a formula method, then wrap it in a
with()
call:
with(Test01, plot(Day_of_year, Draft_No.))
No need for attach.
(rectangular window in vertical position)
dev.new(height = 12, width = 4)
layout(matrix(1:2, nrow=2))
plot(Test01$Day_of_year, Test01$Draft_No.)
plot(Day_of_year,Draft_No.)
(height = 12, width = 6) can't work. The graphs plotted are
distorted off
square shape. I must reduce "width = 4"
Why? TIA
Because you don't appreciate that the dimensions of the device are
not
the same as the dimensions of the plotting region *on* the device.
Most
notably, the margins on the device are given by par("mar") for
example
and are not square:
par("mar")
[1] 5.1 4.1 4.1 2.1
So more space is set aside on the bottom then anywhere else, and the
margin on the right is quite small.
You have already been provided with an answer that you dismissed
because
you didn't appear to appreciate what you were being told.
Re: [R] How to define new matrix based on an elementary row oper
On Aug 28, 2010, at 11:32 AM, Cuckovic Paik wrote: Thank you very much, David; for row swapping: R2<==>R3 A=diag(1:4) A [,1] [,2] [,3] [,4] [1,]1000 [2,]0200 [3,]0030 [4,]0004 A1=A[c(1,3,2,4),] A1 [,1] [,2] [,3] [,4] [1,]1000 [2,]0030 [3,]0200 [4,]0004 Yes, obviously. My point was illustration of construction of equivalent matrix operators (which was what I assumed was the educational goal motivating your insistence on a one step process for a translation operation, but maybe your goal was merely compact code.) The manipulation of indices does not generalize to rotations very well, either. -- David. __ 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 plot an expression-label with variable text (lattice variant)
For the record: lattice behaves slightly differently, it requires an
"as.expression" with bquote.
Dieter
plotExp <- function(what) {
plot.new()
lab = bquote(Estimated~t[50]~ from ~.(what) ) ;
text(0.5,0.2,lab)
}
plotExp("tgv")
library(lattice)
plotLattice <- function(what) {
lab = bquote(Estimated~t[50]~ from ~.(what) ) ;
# Note that ylab is wrong, xlab is correct
xyplot(1~1,xlab=as.expression(lab), ylab=lab)
}
plotLattice("tgv")
--
View this message in context:
http://r.789695.n4.nabble.com/How-to-plot-an-expression-label-with-variable-text-tp2341465p2367176.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.html
and provide commented, minimal, self-contained, reproducible code.
[R] NLME question...
I would like to fit a nonlinear model with "nlme."
The model is a nonlinear growth model with five time points: y = X*b +
e, where design matrix X is defined as
X= | 1 0 |
| 1 1 |
| 1 k |
| 1 2k |
| 1 3k |
and parameter vector b = (b0, b1). And where "k" is a parameter to be
estimated. Of course I also want to estimate the intercept (b0), slope
(b1). Error variances (e) with 5 free parameters and random effects
covariance matrix (2x2: for b0 and b1).
I am not certain how to get the design matrix in the code so that the
"nlme" procedure will be able to
estimate parameter "k" as well as the intercept and slope.
I have tried the following:
meanfunc <- function(x,b0,b1,k){
meangrad <- array(0,c(length(x),2),list(NULL,c("b0","b1")))
for(i in 1:length(x)){
if(x[i]==0){
err[i] <- b0 + b1*x[i]
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- x[i]}
if(x[i]==1){
err[i] <- b0 + b1*x[i]
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- x[i]}
if(x[i]==2){
err[i] <- b0 + b1*(x[i]-1)*k
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- (x[i]-1)*k}
if(x[i]==3){
err[i] <- b0 + b1*2*(x[i]-1)*k
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- (x[i]-1)*k}
if(x[i]==4){
err[i] <- b0 + b1*3*(x[i]-1)*k
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- (x[i]-1)*k}
}
# compute analytical derivatives
attr(err,"gradient") <- meangrad
err
}
data.mlfit <- nlme(y ~ meanfunc(x,b0,b1,k),
fixed=list(b0 ~ 1, b1 ~ 1, k ~ 1),
random=list(b0 ~ 1, b1 ~ 1),
groups = ~subj,
data=nd,
start=list(fixed=c(300,50,1)),
method="ML",verbose=T)
I get the following error message:
Error in meangrad[i, "b1"] <- (x[i] - 1) * k :
number of items to replace is not a multiple of replacement length
In addition: Warning messages:
1: In err[i] <- b0 + b1 * x[i] :
number of items to replace is not a multiple of replacement length
2: In err[i] <- b0 + b1 * x[i] :
number of items to replace is not a multiple of replacement length
3: In err[i] <- b0 + b1 * (x[i] - 1) * k :
number of items to replace is not a multiple of replacement length
Any suggestions or insights would be greatly appreciated.
Thank you,
Jeff
--
**
Jeffrey R. Harring, Assistant Professor
Department of Measurement, Statistics & Evaluation (EDMS)
1230 Benjamin Building
University of Maryland
College Park, MD 20742-1115
Phone: 301.405.3630
Fax:301.314.9245
Email: harr...@umd.edu
Web:http://www.education.umd.edu/EDMS/fac/Harring/webpage.html
__
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] star models
Hi, I am traying to implement an STAR model, but I have some problems. I am following the instruction of the model, that they are in: http://bm2.genes.nig.ac.jp/RGM2/R_current/library/tsDyn/man/star.html that they are from: http://bm2.genes.nig.ac.jp/RGM2/pkg.php?p=tsDyn The model is: star(x, m=2, noRegimes, d = 1, steps = d, series, rob = FALSE, mTh, thDelay, thVar, sig=0.05, trace=TRUE, control=list(), ...) I have implemented the example: mod.star <- star(log10(lynx), m=2, mTh=c(0,1), control=list(maxit=3000)) mod.star However, when I put my variables from my data.frame, R does not find my variables, so, for this reason I have try to create a matrix with my data variables, and to implement the model with the matrix, but I can not create the matrix, either. Please, could someone help me to implement this model? Thanks. Mercedes. __ 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 define new matrix based on an elementary row oper
Thank you very much, David; for row swapping: R2<==>R3 > A=diag(1:4) > A [,1] [,2] [,3] [,4] [1,]1000 [2,]0200 [3,]0030 [4,]0004 > A1=A[c(1,3,2,4),] > A1 [,1] [,2] [,3] [,4] [1,]1000 [2,]0030 [3,]0200 [4,]0004 > -- View this message in context: http://r.789695.n4.nabble.com/How-to-define-new-matrix-based-on-an-elementary-row-operation-in-a-single-step-tp2341768p2364398.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.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to un-crosstabulate data without changing numeric values to text?
On Aug 28, 2010, at 6:07 AM, Nevil Amos wrote:
I have a large amount of data read in from over 140 excel files in
the format of x. r1 to r5 are repeat measures for a given
Wavelength and ANWC_NO.
I need to rearrange x into 3 columns, ANWC_NO,Wavelegth, value ie
ANWC_NOWavelength r1
ANWC_NOWavelength,r2
ANWC_NOWavelength r3
etc...
I can rearrange the data using the code below, however all the
columns end up as strings, not numeric values. I cannot then
summaries the data, ( whcih I need to do in bins of wavelanght for
each ANWC_NO)
> x
Wavelength r1 r2 r3 r4 r5 ANWC_NO
1300 0.003126 0.005382 0.001094 0.012529 0.005632 39239
2302 0.004924 0.006280 0.002366 0.015234 0.006204 39239
3304 0.004769 0.005960 0.002759 0.015856 0.006804 39239
4306 0.005181 0.006717 0.004033 0.017380 0.007675 39239
5308 0.005872 0.008083 0.004429 0.018334 0.008504 39239
6310 0.007164 0.010775 0.005949 0.019952 0.009594 39239
Try:
reshape(x, idvar=c("Wavelength", "ANWC_NO"),
varying=2:6, v.names="r",
direction="long")
I think that melt in package reshape would also work.
--
David
> y =NULL
> rows<-nrow(x)
> for(r in 1:rows){
+ for(c in 2:6){
+ row<-c(c(x[r,7]),as.numeric(c(x[r,1])),as.numeric(c(x[r,c])))
+ y<-rbind(y,row)
+ }}
> colnames(y)<-c("ANWC_NO","WAVELENGTH","VALUE")
> head (y)
ANWC_NO WAVELENGTH VALUE
row "39239" "300" "0.003126"
row "39239" "300" "0.005382"
row "39239" "300" "0.001094"
row "39239" "300" "0.012529"
row "39239" "300" "0.005632"
row "39239" "302" "0.004924"
> mean(y$VALUE)
Error in y$VALUE : $ operator is invalid for atomic vectors
how do I get the data arranged in three columns but maintaining
WavelENGTH and the values as numeric in a data.frame?
Many thanks
Nevil Amos
Monash University
__
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-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 define new matrix based on an elementary row oper
On Aug 28, 2010, at 10:12 AM, (Ted Harding) wrote: On 28-Aug-10 13:15:47, Cuckovic Paik wrote: Thank all help help. Ted's intuitive single step definition is what I want. I try to teach elementary Linear Algebra using R to manupilate matrices. Since my students have no programming experience at all, any fancy and muliple step definition in matrix row operation will confuse them. Again, I appreciate the suggestions from all of you! -- Just to avoid any misunderstanding: My contribution was an explanatory comment to Gabor's solution. It was Gabor who gave the solution! But, while I am at it, if you are teaching elementary Linear Algebra and matrix manipulation, note that the desired result (rows 1, 2 & 4 of A, with row 3 = 2*(row 1 of A) + (row 3 of A)) can be obtained using a matrix product: 1 0 0 0 %*% 1 5 9 13 0 1 0 02 6 10 14 2 0 1 03 7 11 15 0 0 0 14 8 12 16 Which was exactly what my solution was doing, since it becomes in expanded form : D2 <- I %*% A2 + matrix (c(0,0,2,rep(0,13)), 4) %*% A2 = 1 0 0 0 %*% 1 5 9 13 + 0 0 0 0 %*% 1 5 9 13 0 1 0 02 6 10 14 0 0 0 02 6 10 14 0 0 1 03 7 11 15 2 0 0 03 7 11 15 0 0 0 14 8 12 16 0 0 0 04 8 12 16 And A %*% C + B %*% C = (A+B) %*% C and 1 0 0 0 + 0 0 0 0 = 1 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 2 0 0 0 2 0 1 0 0 0 0 1 0 0 0 0 0 0 0 1 There are other elementary operations that can be expressed as matrix operations: Swapping rows 2 and 3 for instance by pre-multiplication by: 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 == matrix (c(1,0,0,0,0,0,1,0,0,1,0,0,0,0,0,1), 4) Or swapping cols 2 and 3 by post multiplication with the same matrix. -- David. i.e. (in R): (diag(nrow(A)) + c(0,0,2,0, 0,0,0,0, 0,0,0,0, 0,0,0,0))%*%A # [,1] [,2] [,3] [,4] # [1,]159 13 # [2,]26 10 14 # [3,]5 17 29 41 # [4,]48 12 16 (as before). Ted. E-Mail: (Ted Harding) Fax-to-email: +44 (0)870 094 0861 Date: 28-Aug-10 Time: 15:12:03 -- XFMail -- __ 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-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] nlme question
I would like to fit a nonlinear model with "nlme."
The model is a nonlinear growth model with five time points: y = X*b +
e, where design matrix X is defined as
X= | 1 0 |
| 1 1 |
| 1 k |
| 1 2k |
| 1 3k |
and parameter vector b = (b0, b1). And where "k" is a parameter to be
estimated. Of course I also want to estimate the intercept (b0), slope
(b1). Error variances (e) with 5 free parameters and random effects
covariance matrix (2x2: for b0 and b1).
I am not certain how to get the design matrix in the code so that the
"nlme" procedure will be able to
estimate parameter "k" as well as the intercept and slope.
I have tried the following:
meanfunc <- function(x,b0,b1,k){
meangrad <- array(0,c(length(x),2),list(NULL,c("b0","b1")))
for(i in 1:length(x)){
if(x[i]==0){
err[i] <- b0 + b1*x[i]
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- x[i]}
if(x[i]==1){
err[i] <- b0 + b1*x[i]
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- x[i]}
if(x[i]==2){
err[i] <- b0 + b1*(x[i]-1)*k
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- (x[i]-1)*k}
if(x[i]==3){
err[i] <- b0 + b1*2*(x[i]-1)*k
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- (x[i]-1)*k}
if(x[i]==4){
err[i] <- b0 + b1*3*(x[i]-1)*k
meangrad[i,"b0"] <- 1
meangrad[i,"b1"] <- (x[i]-1)*k}
}
# compute analytical dervivatives
attr(err,"gradient") <- meangrad
err
}
data.mlfit <- nlme(y ~ meanfunc(x,b0,b1,k),
fixed=list(b0 ~ 1, b1 ~ 1, k ~ 1),
random=list(b0 ~ 1, b1 ~ 1),
groups = ~subj,
data=nd,
start=list(fixed=c(300,50,1)),
method="ML",verbose=T)
I get the following error message:
Error in meangrad[i, "b1"] <- (x[i] - 1) * k :
number of items to replace is not a multiple of replacement length
In addition: Warning messages:
1: In err[i] <- b0 + b1 * x[i] :
number of items to replace is not a multiple of replacement length
2: In err[i] <- b0 + b1 * x[i] :
number of items to replace is not a multiple of replacement length
3: In err[i] <- b0 + b1 * (x[i] - 1) * k :
number of items to replace is not a multiple of replacement length
Any suggestions or insights would be greatly appreciated.
Thank you,
Jeff
--
**
Jeffrey R. Harring, Assistant Professor
Department of Measurement, Statistics & Evaluation (EDMS)
1230 Benjamin Building
University of Maryland
College Park, MD 20742-1115
Phone: 301.405.3630
Fax:301.314.9245
Email: harr...@umd.edu
Web:http://www.education.umd.edu/EDMS/fac/Harring/webpage.html
__
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 define new matrix based on an elementary row oper
On 28-Aug-10 13:15:47, Cuckovic Paik wrote: > Thank all help help. Ted's intuitive single step definition > is what I want. > I try to teach elementary Linear Algebra using R to manupilate > matrices. > Since my students have no programming experience at all, any fancy and > muliple step definition in matrix row operation will confuse them. > Again, I appreciate the suggestions from all of you! > -- Just to avoid any misunderstanding: My contribution was an explanatory comment to Gabor's solution. It was Gabor who gave the solution! But, while I am at it, if you are teaching elementary Linear Algebra and matrix manipulation, note that the desired result (rows 1, 2 & 4 of A, with row 3 = 2*(row 1 of A) + (row 3 of A)) can be obtained using a matrix product: 1 0 0 0 %*% 1 5 9 13 0 1 0 02 6 10 14 2 0 1 03 7 11 15 0 0 0 14 8 12 16 i.e. (in R): (diag(nrow(A)) + c(0,0,2,0, 0,0,0,0, 0,0,0,0, 0,0,0,0))%*%A # [,1] [,2] [,3] [,4] # [1,]159 13 # [2,]26 10 14 # [3,]5 17 29 41 # [4,]48 12 16 (as before). Ted. E-Mail: (Ted Harding) Fax-to-email: +44 (0)870 094 0861 Date: 28-Aug-10 Time: 15:12:03 -- XFMail -- __ 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] Band-wise Sum
Dear David Sir,
Thanks a lot for your guidance. You reply besides helping, also taught me the
importance of sharing your knowledge. It also helped me understand where do I
stand. I am a starter in R and I have started going through at least some mails
everyday whenever possible so that I can learn something from THE WISE like you.
Thanks once again Sir. Your help was great and it means a lot to me and for
other freshers like me.
Regards
Vincy Pyne
--- On Fri, 8/27/10, David Winsemius wrote:
From: David Winsemius
Subject: Re: [R] Band-wise Sum
To: "Vincy Pyne"
Cc: r-help@r-project.org
Received: Friday, August 27, 2010, 2:36 PM
On Aug 27, 2010, at 9:49 AM, Vincy Pyne wrote:
> Hi
>
> I have a large credit portfolio (exceeding 5 borrowers). For particular
> process I need to add up the exposures based on the bands. I am giving a
> small test data below.
I would think that cut() would be the accepted method for defining a factor
variable based on specified cutpoints. If you then wanted to see what the
cumsum() was across the range of possible levels, that to would be a fairly
simple task.
df$ead.cat <- cut(df$ead, breaks=c(0, 10, 50, 100, 200, 500
, 1000, 1) )
df
with(df, tapply(ead.cat, rating, length))
# A AA AAA B BB BBB
# 10 8 2 1 4 7
with(df, tapply(ead.cat, rating, table))
# returns a list of table objects by bond rating
lapply( with(df, tapply(ead.cat, rating, table)) , cumsum)
#returns the cumsum of those tables
# sapply gives a more compact output of that result:
sapply( with(df, tapply(ead.cat, rating, table)) , cumsum)
A AA AAA B BB BBB
(0,1e+05] 4 2 1 0 3 1
(1e+05,5e+05] 8 2 1 1 3 1
(5e+05,1e+06] 9 2 1 1 3 1
(1e+06,2e+06] 9 4 2 1 4 3
(2e+06,5e+06] 9 5 2 1 4 4
(5e+06,1e+07] 10 5 2 1 4 7
(1e+07,1e+08] 10 8 2 1 4 7
Loops, you say we need loops? We don't need no stinkin' loops.
--David.
>
> rating <- c("A", "AAA", "A", "BBB","AA","A","BB", "BBB", "AA", "AA", "AA",
> "A", "A", "AA","BB","BBB","AA", "A", "AAA","BBB","BBB", "BB", "A", "BB", "A",
> "AA", "B","A", "AA", "BBB", "A", "BBB")
>
> ead <- c(169229.93,100, 5877794.25, 9530148.63, 75040962.06, 21000, 1028360,
> 600, 17715000, 14430325.24, 1180946.57, 15, 167490, 81255.16,
> 54812.5, 3000, 1275702.94, 9100, 1763142.3, 3283048.61, 120, 11800,
> 3000, 96894.02, 453671.72, 7590, 106065.24, 940711.67, 2443000, 950,
> 39000, 1501939.67)
>
> ## First I have sorted the data rating-wise as
>
> df <- data.frame(rating, ead)
>
> df_sorted <-
> df[order(df$rating),]
>
> df_sorted_AAA <- subset(df_sorted, rating=="AAA")
> df_sorted_AA <- subset(df_sorted, rating=="AA")
> df_sorted_A <- subset(df_sorted, rating=="A")
> df_sorted_BBB <- subset(df_sorted, rating=="BBB")
> df_sorted_BB <- subset(df_sorted, rating=="BB")
> df_sorted_B <- subset(df_sorted, rating=="B")
> df_sorted_CCC <- subset(df_sorted, rating=="CCC")
>
> ## we begin with BBB rating. The R output for df_sorted_BBB is as follows
>
>> df_sorted_BBB
> rating ead
> 4 BBB 9530149
> 8 BBB 600
> 16 BBB 3000
> 20 BBB 3283049
> 21 BBB 120
> 30 BBB 950
> 32 BBB 1501940
>
> My problem is I need to totals of eads falling in the respective bands
>
> I
> am defining bands in millions as
>
> seq_BBB <- seq(100, max(df_sorted_BBB$ead), by = 100)
>
> # The output is
> [1] 1e+06 2e+06 3e+06 4e+06 5e+06 6e+06 7e+06 8e+06 9e+06
>
> So for the sub data pertaining to Rating "BBB", I want corresponding ead
> totals i.e. I want ead totals where ead < 1e+06, then I want ead totals where
> 1+e06 < ead < 2e+06, 2e+06 < ead < 3e+06 ...and so on.
>
> I have tried the following code
>
> s_BBB <- NULL
>
> for (i in 1:length(s_BBB))
> {
> s_BBB[i] = sum(subset(df_sorted_BBB$ead, df_sorted_BBB$ead < s_BBB[i]))
> }
>
> I was trying to find totals ofads < 1e+06, ead < 2e+06, ead<3e+06and so on.
>
> but the result is
>
>> s_BBB
> [1] 0
>
>
> I apologize if I am not able to express my problem properly. My only
> objective is first to sort the whole portfolio rating-wise and then within
> each of these rating-wise sorted data, I wish to find out total of eads based
> on various bands starting <100, 100 - 20, 200 - 300,
> 300 - 400 and so on. Since the database contains more than 5
> records, various ead amounts ranging from few 000's to billion are available.
>
> Please guide
>
> Thanking you all in advance
>
> Vincy
>
>
>
>
>
>
>
>
>
>
>
>
> [[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.
David Winsemius, MD
West Hartford, CT
Re: [R] How to define new matrix based on an elementary row operation in a single step?
Thank all help help. Ted's intuitive single step definition is what I want. I try to teach elementary Linear Algebra using R to manupilate matrices. Since my students have no programming experience at all, any fancy and muliple step definition in matrix row operation will confuse them. Again, I appreciate the suggestions from all of you! -- View this message in context: http://r.789695.n4.nabble.com/How-to-define-new-matrix-based-on-an-elementary-row-operation-in-a-single-step-tp2341768p2362791.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.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] About plot graphs
Hi Gavin,
Lot of thank for your detail explanation.
Just looked at;
?with
?within
Both "Evaluate an Expression in a Data Environment"
Usage: Both;
with(data, expr, ...)
within(data, expr, ...)
Details:
.
‘within’ is similar, except that it examines the environment after
the evaluation of ‘expr’ and makes the corresponding modifications
to ‘data’ (this may fail in the data frame case if objects are
created which cannot be stored in a data frame), and returns it.
‘within’ can be used as an alternative to ‘transform’.
What does it mean ". examines the environment after the evaluation of
'expr'
."
Tks
B.R.
Stephen
B.R.
Stephen L
- Original Message
From: Gavin Simpson
To: Stephen Liu
Cc: "r-help@r-project.org"
Sent: Sat, August 28, 2010 3:28:34 PM
Subject: Re: [R] About plot graphs
On Fri, 2010-08-27 at 22:14 -0700, Stephen Liu wrote:
> Hi Gavin,
>
>
> Thanks for your advice and the examples explaining plotting settings.
>
> The steps on your examples work on my test.
>
>
> > 2) Don't attach() things, you are asking for trouble
>
> > If a function has a formula method (which plot does) then use it like
> > this: plot(Draft_No. ~ Day_of_year, data = Test01)
>
> > If the function doesn't have a formula method, then wrap it in a
> > with()
> > call:
>
> > with(Test01, plot(Day_of_year, Draft_No.))
>
> > No need for attach.
>
>
> Noted and thanks. What will be the problem caused by "attach()"?
If you change the underlying data, this is not reflected in the attached
copy, because it is just that, a "copy"[1] created at the point at which
you attached the object. E.g.
## Some data, which we attach
dat <- data.frame(A = 1:10, B = letters[1:10])
attach(dat)
## Look at A
A
## Change or dat object by altering the A component
dat <- within(dat, A <- LETTERS[1:10])
## Look at A
A
## Look at what A really is
with(dat, A)
Using with() and within() etc have two key advantages over attach: i)
only one version of the data/object exists, ii) the intention of code
using:
with(dat, "do something with A")
is much more clear than
"do something with A"
A could be anything, anywhere. More info is on the ?attach help page.
[1] ?attach contains the details of what I mean by "copy"
>
> > dev.new(height = 6, width = 12)
> > layout(matrix(1:2, ncol = 2))
> > op <- par(pty = "s") ## this is the important bit
> > plot(runif(100), rnorm(100))
> > plot(runif(100), rnorm(100), col = "red")
> > par(op) ## now reset the pars
> > layout(1)
>
> What is the function of layout(1) ? Tks
The opposite of
layout(matrix(1:4, ncol = 2))
for example. It ( layout(1) ) says create a layout with a single
plotting region. So we use it to reset the current device back to
normal. I find it is good working practice to tidy up after doing plots
like this. In the code above, we change both the layout() and the
plotting parameters (via par() ), and the last two lines of code in that
example reset these changes.
G
>
> B.R.
> satimis
>
>
>
>
> - Original Message
> From: Gavin Simpson
> To: Stephen Liu
> Cc: "r-help@r-project.org"
> Sent: Fri, August 27, 2010 5:38:40 PM
> Subject: Re: [R] About plot graphs
>
> On Fri, 2010-08-27 at 02:05 -0700, Stephen Liu wrote:
> > Hi Gavin,
> >
> > Thanks for your advice which works for me.
> >
> >
> > (rectangular window)
> > dev.new(height = 6, width = 12)
> > layout(matrix(1:2, nrow=1))
> > plot(Test01$Day_of_year, Test01$Draft_No.)
> > attach(Test01)
> > plot(Day_of_year,Draft_No.)
>
> 1) I can't reproduce this; where/what is Test01? But don;t bother
> sending, see my example below
> 2) Don't attach() things, you are asking for trouble
>
> If a function has a formula method (which plot does) then use it like
> this: plot(Draft_No. ~ Day_of_year, data = Test01)
>
> If the function doesn't have a formula method, then wrap it in a with()
> call:
>
> with(Test01, plot(Day_of_year, Draft_No.))
>
> No need for attach.
>
> >
> > (rectangular window in vertical position)
> > dev.new(height = 12, width = 4)
> > layout(matrix(1:2, nrow=2))
> > plot(Test01$Day_of_year, Test01$Draft_No.)
> > plot(Day_of_year,Draft_No.)
> >
> > (height = 12, width = 6) can't work. The graphs plotted are distorted off
> > square shape. I must reduce "width = 4"
> >
> > Why? TIA
>
> Because you don't appreciate that the dimensions of the device are not
> the same as the dimensions of the plotting region *on* the device. Most
> notably, the margins on the device are given by par("mar") for example
> and are not square:
>
> > par("mar")
> [1] 5.1 4.1 4.1 2.1
>
> So more space is set aside on the bottom then anywhere else, and the
> margin on the right is quite small.
>
> You have already been provided with an answer that you dismissed because
> you didn't appear to appreciate what you were being told.
>
> Compare this:
>
> dev.new(height = 6, width = 12)
> layout(matrix(1:2, ncol = 2))
> plot(runif(100), rnorm(100))
> pl
Re: [R] How to define new matrix based on an elementary row oper
On 28-Aug-10 11:27:30, Gabor Grothendieck wrote: > On Sat, Aug 28, 2010 at 1:32 AM, Cheng Peng > wrote: >> >> Sorry for possible misunderstanding: >> >> I want to define a matrix (B) based on an existing matrix (A) in a >> single >> step and keep A unchanged: >> >>> #Existing matrix >>> A=matrix(1:16,ncol=4) >>> A >> [,1] [,2] [,3] [,4] >> [1,] 159 13 >> [2,] 26 10 14 >> [3,] 37 11 15 >> [4,] 48 12 16 >>> # New matrix B is defined to be the submatrix after row1 and column1 >>> # are deleted. >>> B=A[-1,-1]# this single step deletes row1 nad column 1 and >>> assigns the name to the resulting submatrix. >>> B # check the new matrix B >> [,1] [,2] [,3] >> [1,] 6 10 14 >> [2,] 7 11 15 >> [3,] 8 12 16 >>> A# check the original matrix A >> [,1] [,2] [,3] [,4] >> [1,] 159 13 >> [2,] 26 10 14 >> [3,] 37 11 15 >> [4,] 48 12 16 >> >> Question: How can I do define a new matrix (D) by adding 2*row1 >> to row3 in A in a single step as what was done in the above example? >> >> If you do: _A[3,]=2*A[1,]+A[3,], _the new A is not the original A; >> if you D=A first, then D[3,]=2*D[1,]+D[3,], you used two step! >> >> Hope this clarifies my original question. Thanks again. > > Try using replace: > >D <- replace(A, row(A) == 3, 2 * A[1,] + A[3,]) That's neat -- and subtle! It's the sort of thing one wouldn't think of doing until one has seen it done! ?replace: Replace Values in a Vector Description: 'replace' replaces the values in 'x' with indices given in 'list' by those given in 'values'. If necessary, the values in 'values' are recycled. Usage: replace(x, list, values) Arguments: x: vector list: an index vector values: replacement values Value: A vector with the values replaced. Note: 'x' is unchanged: remember to assign the result. One needs to realise that a matrix A is a vector with a dimension attribute, hence is accessed "sequentially" like any vector (though displayed as an array). So, with the above matrix A: row(A)==3 # [,1] [,2] [,3] [,4] # [1,] FALSE FALSE FALSE FALSE # [2,] FALSE FALSE FALSE FALSE # [3,] TRUE TRUE TRUE TRUE # [4,] FALSE FALSE FALSE FALSE A[row(A)==3] # [1] 3 7 11 15 Then 'replace(A, row(A) == 3, 2 * A[1,] + A[3,])' replaces that part of the vector A as indexed by TRUE with the given expression 2 * A[1,] + A[3,] = 5 17 29 41 Hence replace(A, row(A) == 3, 2 * A[1,] + A[3,]) # [,1] [,2] [,3] [,4] # [1,]159 13 # [2,]26 10 14 # [3,]5 17 29 41 # [4,]48 12 16 is then assigned to D. Hoping this helps to clarify! Ted. E-Mail: (Ted Harding) Fax-to-email: +44 (0)870 094 0861 Date: 28-Aug-10 Time: 13:11:08 -- XFMail -- __ 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 define new matrix based on an elementary row operation in a single step?
On Sat, Aug 28, 2010 at 1:32 AM, Cheng Peng wrote: > > Sorry for possible misunderstanding: > > I want to define a matrix (B) based on an existing matrix (A) in a single > step and keep A unchanged: > >> #Existing matrix >> A=matrix(1:16,ncol=4) >> A > [,1] [,2] [,3] [,4] > [1,] 1 5 9 13 > [2,] 2 6 10 14 > [3,] 3 7 11 15 > [4,] 4 8 12 16 >> # New matrix B is defined to be the submatrix after row1 and column1 are >> deleted. >> B=A[-1,-1] # this single step deletes row1 nad column 1 and assigns the >> name to the resulting submatrix. >> B # check the new matrix B > [,1] [,2] [,3] > [1,] 6 10 14 > [2,] 7 11 15 > [3,] 8 12 16 >> A # check the original matrix A > [,1] [,2] [,3] [,4] > [1,] 1 5 9 13 > [2,] 2 6 10 14 > [3,] 3 7 11 15 > [4,] 4 8 12 16 > > > Question: How can I do define a new matrix (D) by adding 2*row1 to row3 in A > in a single step as what was done in the above example? > > If you do: A[3,]=2*A[1,]+A[3,], the new A is not the original A; if you > D=A first, then D[3,]=2*D[1,]+D[3,], you used two step! > > Hope this clarifies my original question. Thanks again. > Try using replace: D <- replace(A, row(A) == 3, 2 * A[1,] + A[3,]) -- Statistics & Software Consulting GKX Group, GKX Associates Inc. tel: 1-877-GKX-GROUP email: ggrothendieck at 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] how to un-crosstabulate data without changing numeric values to text?
I have a large amount of data read in from over 140 excel files in the
format of x. r1 to r5 are repeat measures for a given Wavelength and
ANWC_NO.
I need to rearrange x into 3 columns, ANWC_NO,Wavelegth, value ie
ANWC_NOWavelength r1
ANWC_NOWavelength,r2
ANWC_NOWavelength r3
etc...
I can rearrange the data using the code below, however all the columns
end up as strings, not numeric values. I cannot then summaries the
data, ( whcih I need to do in bins of wavelanght for each ANWC_NO)
> x
Wavelength r1 r2 r3 r4 r5 ANWC_NO
1300 0.003126 0.005382 0.001094 0.012529 0.005632 39239
2302 0.004924 0.006280 0.002366 0.015234 0.006204 39239
3304 0.004769 0.005960 0.002759 0.015856 0.006804 39239
4306 0.005181 0.006717 0.004033 0.017380 0.007675 39239
5308 0.005872 0.008083 0.004429 0.018334 0.008504 39239
6310 0.007164 0.010775 0.005949 0.019952 0.009594 39239
> y =NULL
> rows<-nrow(x)
> for(r in 1:rows){
+ for(c in 2:6){
+ row<-c(c(x[r,7]),as.numeric(c(x[r,1])),as.numeric(c(x[r,c])))
+ y<-rbind(y,row)
+ }}
> colnames(y)<-c("ANWC_NO","WAVELENGTH","VALUE")
> head (y)
ANWC_NO WAVELENGTH VALUE
row "39239" "300" "0.003126"
row "39239" "300" "0.005382"
row "39239" "300" "0.001094"
row "39239" "300" "0.012529"
row "39239" "300" "0.005632"
row "39239" "302" "0.004924"
> mean(y$VALUE)
Error in y$VALUE : $ operator is invalid for atomic vectors
how do I get the data arranged in three columns but maintaining
WavelENGTH and the values as numeric in a data.frame?
Many thanks
Nevil Amos
Monash University
__
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 define new matrix based on an elementary row operation in a single step?
Hi:
Here are some functions for computing elementary matrices so that you can do
Gauss elimination the hard way (the easy way is the LU decomposition, but I
digress) I wouldn't use these for serious work, since there are no
checks for matrix instability, but the essential ideas are there. It is
assumed that you know the relationship between elementary matrices, the
original matrix and its Gauss reduced form.
Eij <- function(p, i, j) {
# Permutation matrix that obtains by exchanging rows i and j
E <- diag(rep(1, p))
E[j, i] <- E0[i,j] <- 1
E[j, j] <- E0[i, i] <- 0
E
}
EiC <- function(p, i, c) {
# Elementary matrix obtained by multiplying i-th row of I by c
if(c == 0L) stop ('cannot multiply row by zero')
E <- diag(rep(1, p))
E[i, ] <- c * E0[i, ]
E
}
EijC <- function(p, i, j, c) {
# Elementary matrix obtained by multiplying c * row j and
# adding it to the i-th row of E
E <- diag(rep(1, p))
E[i, ] <- E[i, ] + c * E[j, ]
E
}
To answer your question,
> A=matrix(1:16,ncol=4)
> A
[,1] [,2] [,3] [,4]
[1,]159 13
[2,]26 10 14
[3,]37 11 15
[4,]48 12 16
> EijC(4, 3, 1, 2) %*% A # Add 2 * R1 of A to R3
[,1] [,2] [,3] [,4]
[1,]159 13
[2,]26 10 14
[3,]5 17 29 41
[4,]48 12 16
To reduce the first column of A:
> EijC(4, 2, 1, -2) %*% EijC(4, 3, 1, -3) %*% EijC(4, 4, 1, -4) %*% A
[,1] [,2] [,3] [,4]
[1,]159 13
[2,]0 -4 -8 -12
[3,]0 -8 -16 -24
[4,]0 -12 -24 -36
HTH,
Dennis
On Fri, Aug 27, 2010 at 10:32 PM, Cheng Peng wrote:
>
> Sorry for possible misunderstanding:
>
> I want to define a matrix (B) based on an existing matrix (A) in a single
> step and keep A unchanged:
>
> > #Existing matrix
> > A=matrix(1:16,ncol=4)
> > A
> [,1] [,2] [,3] [,4]
> [1,]159 13
> [2,]26 10 14
> [3,]37 11 15
> [4,]48 12 16
> > # New matrix B is defined to be the submatrix after row1 and column1 are
> > deleted.
> > B=A[-1,-1]# this single step deletes row1 nad column 1 and assigns
> the
> > name to the resulting submatrix.
> > B # check the new matrix B
> [,1] [,2] [,3]
> [1,]6 10 14
> [2,]7 11 15
> [3,]8 12 16
> > A # check the original matrix A
> [,1] [,2] [,3] [,4]
> [1,]159 13
> [2,]26 10 14
> [3,]37 11 15
> [4,]48 12 16
>
>
> Question: How can I do define a new matrix (D) by adding 2*row1 to row3 in
> A
> in a single step as what was done in the above example?
>
> If you do: A[3,]=2*A[1,]+A[3,], the new A is not the original A; if you
> D=A first, then D[3,]=2*D[1,]+D[3,], you used two step!
>
> Hope this clarifies my original question. Thanks again.
>
>
>
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/How-to-define-new-matrix-based-on-an-elementary-row-operation-in-a-single-step-tp2341768p2352538.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.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.
[R] problem after repackaging
Hy, I had a mistake on a function of a package i have created! I have solved it and then i repackaged and installed the modified package. I use to launch R from Excel! And so when i launch R, and next call my function from the workspace, i still find the problem on my function. And when i read on my workspace, the source code of my function, i find the old version of my function (the one from the precedent version of my package)! If one of you have ever met this kind of problem, could you help me. Thanks Samy Khezami M.S. Quantitative Finance EM Lyon Business School 0033625430829 samy-khezami@em-lyon.com [[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 define new matrix based on an elementary row operation in a single step?
On Aug 28, 2010, at 2:54 AM, Joshua Wiley wrote: Is this sufficiently single steppish for you? D <- A <- matrix(1:16, 4) D[3, ] <- 2 * D[1, ] + D[3, ] # Alternately, you could do this # but it is much messier, and I do not see how # two steps is really an issue # you want to end up with two matrices anyways # so it's not like you save memory by only making one assignment A2 <- matrix(1:16, 4) D2 <- rbind(A2[1:2, ], 2 * A2[1, ] + A2[3, ], A2[4, ]) You can do it the "matrix-ey" way as well by creating a row extractor/ multiplier operator. This matrix multiplication copies 2 times the 1st row to the third row: The operation: > matrix( c(0,0,2,rep(0,13)), 4) %*% A2 [,1] [,2] [,3] [,4] [1,]0000 [2,]0000 [3,]2 10 18 26 [4,]0000 So it can be used thusly: D2 <- A2 + matrix (c(0,0,2,rep(0,13)), 4) %*% A2 D2 [,1] [,2] [,3] [,4] [1,]159 13 [2,]26 10 14 [3,]5 17 29 41 [4,]48 12 16 This method does not pull the "donor" matrix apart. -- David all.equal(A, A2) all.equal(D, D2) Cheers, Josh On Fri, Aug 27, 2010 at 10:32 PM, Cheng Peng wrote: Sorry for possible misunderstanding: I want to define a matrix (B) based on an existing matrix (A) in a single step and keep A unchanged: #Existing matrix A=matrix(1:16,ncol=4) A [,1] [,2] [,3] [,4] [1,]159 13 [2,]26 10 14 [3,]37 11 15 [4,]48 12 16 # New matrix B is defined to be the submatrix after row1 and column1 are deleted. B=A[-1,-1]# this single step deletes row1 nad column 1 and assigns the name to the resulting submatrix. B # check the new matrix B [,1] [,2] [,3] [1,]6 10 14 [2,]7 11 15 [3,]8 12 16 A # check the original matrix A [,1] [,2] [,3] [,4] [1,]159 13 [2,]26 10 14 [3,]37 11 15 [4,]48 12 16 Question: How can I do define a new matrix (D) by adding 2*row1 to row3 in A in a single step as what was done in the above example? If you do: A[3,]=2*A[1,]+A[3,], the new A is not the original A; if you D=A first, then D[3,]=2*D[1,]+D[3,], you used two step! Hope this clarifies my original question. Thanks again. -- View this message in context: http://r.789695.n4.nabble.com/How-to-define-new-matrix-based-on-an-elementary-row-operation-in-a-single-step-tp2341768p2352538.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.html and provide commented, minimal, self-contained, reproducible code. -- Joshua Wiley Ph.D. Student, Health Psychology University of California, Los Angeles http://www.joshuawiley.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-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] Shifting of Principal amount while calculating Present Value
Matrix operations work "the other way" to what you expect on rows and
columns. Try
## Your data:
time <- c(1, 2, 3) # Please don't use `t' as a variable name
cash_flow <- c(7, 7, 107)
zero_rate <- data.frame(rating=factor(c("AAA","AA","A","B")),
year1=c(3.60, 3.65, 3.72, 4.10),
year2=c(4.17,4.22,4.32,4.67),
year3=c(4.73,4.78,4.93,5.25))
zero_rate1 <- zero_rate[, -1]
## Make clear we are dealing with matrices (not really needed)
zero_rate2 <- as.matrix(zero_rate1)
## The calculation
colSums(cash_flow/(t(1 + zero_rate2 / 100)^time))
# [1] 106.3549 106.2122 105.7969 104.8871
Here t() is the transpose function. Try to step through each part of
the calculation to figure it out. You may also want to look at c(1, 2,
3) / matrix(1, 3, 3) and matrix(sqrt(2), 3, 3)^c(1, 2, 4) for smaller
examples.
Hope this helps a little.
Allan.
On 20/08/10 12:24, Sarah Sanchez wrote:
> Dear R Helpers
>
> I have following data -
>
> cash_flow = c(7, 7, 107) # 107 = Principle 100 + interest of 7%
> t = c(1,2,3)
>
> and zero rate table as
>
> rating year1 year2 year3
> AAA3.604.17 4.73
> AA 3.654.22
> 4.78
> A 3.72 4.32 4.93
> BBB4.104.67 5.25
>
>
> For each of these ratings I need to calculate the Present Value as (say e.g.
> for AAA)
>
> PV(AAA) = 7/(1+3.60/100) + 7/(1+4.17/100)^2 + 107/(1+4.73/100)^3 which is
> equal to 106.3549
>
> Similarly when used the respective rates, PV(A) = 106.2122, PV(A) = 105.7969
> and PV(BBB) = 104.8871
>
> #
>
> ## My problem
>
> I have tried the following R code.
>
> zero_rate =
> read.csv('zero_rate_table.csv')
> zero_rate1 = zero_rate[, -1]
>
> cash_flow = c(7, 7, 107)
>
> t = c(1,2,3)
>
>
> PV_table = cash_flow / (1+zero_rate1/100)^t
>
> ## Then using rowSums, I should get the required result. However, I am
> getting following output as
>
>
>> PV_table
>>
> year_1year_2 year_3
> [1,] 6.756757 6.45078 93.147342
> [2,] 6.515675 94.521493 6.680664
> [3,] 95.8950646.710123 6.357680
> [4,] 6.724304 6.389305 91.773536
>
>
> rowSums(PV_table) gives
>
> [1] 106.35489 107.71783 108.96287 104.88714.
>
> Thus, the result is correct only in first case. In second case onwards, there
> is a shift of final cashflow e.g. in 2nd case, pertaining to year 2, value is
> 94.521493 which means it includes principal of 100.
>
> Likewise
> in third case, the principal is included in the first cash-flow itself.
> Again the fourth result is correct. So there is pre-shift of capital.
>
> I am not sure whether I am making any mistake in reading the zero_rate csv
> file or my formula is incorrect. Please guide.
>
> Thanking you all in advance
>
> Sarah
>
>
>
>
>
>
>
>
>
> [[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.
Re: [R] Calculating p and q values with R
On 08/28/2010 03:10 AM, ndar wrote: > > Hi, > I have a huge dataset (53 million records). I have to calculate the p and q > values of my data. How can I do it in R or perl? I have downloaded R (I'm > completely new to R). and the package qvalue but I don't understand how can > I call/use qvalue package with R. When I type library(qvalue), it gives me > an error that this package doesn't exist. What should I do? > Thanks! Looks like you didn't _install_ the package. Ch 13 in "An Introduction to R", which comes with R, and which you probably should at least skim through anyway. Apart from that, there is a tutorial to qvalue on CRAN which you should probably read; likely also references therein. It is not clear from your message how good your theoretical understanding is. For matters of this character, I'd say that if it is not very good, you need professional assistance beyond what can be expected from the mailing lists. -- Peter Dalgaard Center for Statistics, Copenhagen Business School Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@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] About plot graphs
On Fri, 2010-08-27 at 22:14 -0700, Stephen Liu wrote:
> Hi Gavin,
>
>
> Thanks for your advice and the examples explaining plotting settings.
>
> The steps on your examples work on my test.
>
>
> > 2) Don't attach() things, you are asking for trouble
>
> > If a function has a formula method (which plot does) then use it like
> > this: plot(Draft_No. ~ Day_of_year, data = Test01)
>
> > If the function doesn't have a formula method, then wrap it in a
> > with()
> > call:
>
> > with(Test01, plot(Day_of_year, Draft_No.))
>
> > No need for attach.
>
>
> Noted and thanks. What will be the problem caused by "attach()"?
If you change the underlying data, this is not reflected in the attached
copy, because it is just that, a "copy"[1] created at the point at which
you attached the object. E.g.
## Some data, which we attach
dat <- data.frame(A = 1:10, B = letters[1:10])
attach(dat)
## Look at A
A
## Change or dat object by altering the A component
dat <- within(dat, A <- LETTERS[1:10])
## Look at A
A
## Look at what A really is
with(dat, A)
Using with() and within() etc have two key advantages over attach: i)
only one version of the data/object exists, ii) the intention of code
using:
with(dat, "do something with A")
is much more clear than
"do something with A"
A could be anything, anywhere. More info is on the ?attach help page.
[1] ?attach contains the details of what I mean by "copy"
>
> > dev.new(height = 6, width = 12)
> > layout(matrix(1:2, ncol = 2))
> > op <- par(pty = "s") ## this is the important bit
> > plot(runif(100), rnorm(100))
> > plot(runif(100), rnorm(100), col = "red")
> > par(op) ## now reset the pars
> > layout(1)
>
> What is the function of layout(1) ? Tks
The opposite of
layout(matrix(1:4, ncol = 2))
for example. It ( layout(1) ) says create a layout with a single
plotting region. So we use it to reset the current device back to
normal. I find it is good working practice to tidy up after doing plots
like this. In the code above, we change both the layout() and the
plotting parameters (via par() ), and the last two lines of code in that
example reset these changes.
G
>
> B.R.
> satimis
>
>
>
>
> - Original Message
> From: Gavin Simpson
> To: Stephen Liu
> Cc: "r-help@r-project.org"
> Sent: Fri, August 27, 2010 5:38:40 PM
> Subject: Re: [R] About plot graphs
>
> On Fri, 2010-08-27 at 02:05 -0700, Stephen Liu wrote:
> > Hi Gavin,
> >
> > Thanks for your advice which works for me.
> >
> >
> > (rectangular window)
> > dev.new(height = 6, width = 12)
> > layout(matrix(1:2, nrow=1))
> > plot(Test01$Day_of_year, Test01$Draft_No.)
> > attach(Test01)
> > plot(Day_of_year,Draft_No.)
>
> 1) I can't reproduce this; where/what is Test01? But don;t bother
> sending, see my example below
> 2) Don't attach() things, you are asking for trouble
>
> If a function has a formula method (which plot does) then use it like
> this: plot(Draft_No. ~ Day_of_year, data = Test01)
>
> If the function doesn't have a formula method, then wrap it in a with()
> call:
>
> with(Test01, plot(Day_of_year, Draft_No.))
>
> No need for attach.
>
> >
> > (rectangular window in vertical position)
> > dev.new(height = 12, width = 4)
> > layout(matrix(1:2, nrow=2))
> > plot(Test01$Day_of_year, Test01$Draft_No.)
> > plot(Day_of_year,Draft_No.)
> >
> > (height = 12, width = 6) can't work. The graphs plotted are distorted off
> > square shape. I must reduce "width = 4"
> >
> > Why? TIA
>
> Because you don't appreciate that the dimensions of the device are not
> the same as the dimensions of the plotting region *on* the device. Most
> notably, the margins on the device are given by par("mar") for example
> and are not square:
>
> > par("mar")
> [1] 5.1 4.1 4.1 2.1
>
> So more space is set aside on the bottom then anywhere else, and the
> margin on the right is quite small.
>
> You have already been provided with an answer that you dismissed because
> you didn't appear to appreciate what you were being told.
>
> Compare this:
>
> dev.new(height = 6, width = 12)
> layout(matrix(1:2, ncol = 2))
> plot(runif(100), rnorm(100))
> plot(runif(100), rnorm(100), col = "red")
> layout(1)
>
> with this:
>
> dev.new(height = 6, width = 12)
> layout(matrix(1:2, ncol = 2))
> op <- par(pty = "s") ## this is the important bit
> plot(runif(100), rnorm(100))
> plot(runif(100), rnorm(100), col = "red")
> par(op) ## now reset the pars
> layout(1)
>
> So now the regions are square, we have the asymmetric margins like all R
> plots and we have drawn this on a device that has ~ appropriate
> dimensions.
>
> If you want to fiddle more with this, then you'll need to make the
> margins equal, but you don't want to do that really as you need more
> space in certain areas to accommodate axis labels and tick labels etc.
>
> For the vertical one, this works for me, though note that because of the
> margins, pty = "s" is making the individual plotting regions smaller to
> respect the square plotting
