[R] Plackett-Dale Model in R
Dear R users, Can someone inform me about a library/function in R that fits a Plackett-Dale model ? Thanks in advance Pryseley - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] lme convergence
Dear R-Users, Is it possible to get the covariance matrix from an lme model that did not converge ? I am doing a simulation which entails fitting linear mixed models, using a "for loop". Within each loop, i generate a new data set and analyze it using a mixed model. The loop stops When the "lme function" does not converge for a simulated dataset. I want to inquire if there is a method to suppress the error message from the lme function, or better still, a way of going about this issue of the loop ending once the lme function does not converge. Thanks in advance, Pryseley - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Help: lme
Good day, Yeah i tried most of the suggestions and still had the same problem. Some people also suggested that i supplied a data set so they could verify what i am saying. That is why i am still asking the same question with more details as required. My main focus now is to calculate the values in the "unaccessible"-printed output from the accessible ones as mentioned towards the end of my latest mail. Kind regards Pryseley "Doran, Harold" <[EMAIL PROTECTED]> wrote: Pryseley: You asked this question last week and received answers to that question along with other suggestions from Doug Bates and myself. Is there some reason the post by Andrew Robinson is not working? Harold > -Original Message- > From: [EMAIL PROTECTED] > [mailto:[EMAIL PROTECTED] On Behalf Of Pryseley Assam > Sent: Thursday, June 01, 2006 6:25 AM > To: r-help@stat.math.ethz.ch > Subject: [R] Help: lme > > Good day R-Users, > > I have a problem accessing some values in the output from > the summary of an lme fit. > > > The structure of my data is as shown below (I have attached a > copy of the full data). > > > id trials endp Z.sas ST > 1 1 -1 -1 42.42884 > 1 1 1 -1 48.12007 > 2 1 -1 -1 43.42878 > 2 1 1 -1 46.82817 > 3 1 -1 -1 45.56736 > ... > 598 10 -1 1 49.59715 > 598 10 1 1 55.46324 > 599 10 -1 1 46.80873 > 599 10 1 1 53.52223 > 600 10 -1 1 48.58257 > 600 10 1 1 56.78674 > > I used the following code to fit my model of interest: > > ggg <- lme (ST~ -1 + as.factor(endp):Z.sas + > as.factor(endp), data=dat4a, > random=~-1 + as.factor(endp) + > as.factor(endp):Z.sas|as.factor(trials), > correlation = > corSymm(form=~1|as.factor(trials)/as.factor(id)), > weights=varIdent(form=~1|endp)) > > hh <- summary(ggg) > hh > > Below is the following part of the output of interest I > wish to access: > > Correlation Structure: General > Formula: ~1 | as.factor(trials)/as.factor(id) Parameter estimate(s): > Correlation: > 1 > 2 0.785 > Variance function: > Structure: Different standard deviations per stratum > Formula: ~1 | endp > Parameter estimates: > -1 1 > 1.000 0.9692405 > > I wish to access the value of the correlation (0.785) and > the vector of the variance function estimates (1, 0.969). I > know these can be done through the intervals function, but > sometimes when the estimated Hessian matrix is not positive > definite or something like that (i am not quite sure), the > intervals function delivers an error message. > > Thus, i will like to ask if there is another way to access > these values. I tried using the following code: > > hh$modelStruct$corStruct[1] > hh$modelStruct$varStruct[1] > > Rather the output was: > > > hh$modelStruct$corStruct[1] > [1] -1.308580 > > hh$modelStruct$varStruct[1] > [1] -0.03124255 > > I presume there is a way to calculate the correlation and > variance function coefficients using these values. > > Could you tell me how to access those values (without using > the intervals function) or better still how to calculate the > values from the last output values. > > Kind regards > Pryseley > > > > __ > > > __ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Help: lme
Good day R-Users, I have a problem accessing some values in the output from the summary of an lme fit. The structure of my data is as shown below (I have attached a copy of the full data). id trials endp Z.sas ST 1 1 -1 -142.42884 1 11 -1 48.12007 2 1 -1 -143.42878 2 11 -1 46.82817 3 1 -1 -145.56736 . 598 10 -11 49.59715 598 10 11 55.46324 599 10 -1 1 46.80873 599 101 1 53.52223 600 10 -1 1 48.58257 600 101 1 56.78674 I used the following code to fit my model of interest: ggg <- lme (ST~ -1 + as.factor(endp):Z.sas + as.factor(endp), data=dat4a, random=~-1 + as.factor(endp) + as.factor(endp):Z.sas|as.factor(trials), correlation = corSymm(form=~1|as.factor(trials)/as.factor(id)), weights=varIdent(form=~1|endp)) hh <- summary(ggg) hh Below is the following part of the output of interest I wish to access: Correlation Structure: General Formula: ~1 | as.factor(trials)/as.factor(id) Parameter estimate(s): Correlation: 1 2 0.785 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | endp Parameter estimates: -1 1 1.000 0.9692405 I wish to access the value of the correlation (0.785) and the vector of the variance function estimates (1, 0.969). I know these can be done through the intervals function, but sometimes when the estimated Hessian matrix is not positive definite or something like that (i am not quite sure), the intervals function delivers an error message. Thus, i will like to ask if there is another way to access these values. I tried using the following code: hh$modelStruct$corStruct[1] hh$modelStruct$varStruct[1] Rather the output was: > hh$modelStruct$corStruct[1] [1] -1.308580 > hh$modelStruct$varStruct[1] [1] -0.03124255 I presume there is a way to calculate the correlation and variance function coefficients using these values. Could you tell me how to access those values (without using the intervals function) or better still how to calculate the values from the last output values. Kind regards Pryseley __ "id","trials","endp","Z.sas","ST" "1",1,1,-1,-1,42.4288378507851 "2",1,1,1,-1,48.1200688213539 "3",2,1,-1,-1,43.4287823712575 "4",2,1,1,-1,46.8281725658322 "5",3,1,-1,-1,45.5673614208191 "6",3,1,1,-1,47.8126803355643 "7",4,1,-1,-1,44.5853106980775 "8",4,1,1,-1,48.3336576286004 "9",5,1,-1,-1,42.8799097392805 "10",5,1,1,-1,45.9461944592287 "11",6,1,-1,-1,41.9963773500656 "12",6,1,1,-1,46.7681426907703 "13",7,1,-1,-1,45.0950139016658 "14",7,1,1,-1,47.1927110682639 "15",8,1,-1,-1,44.2809181126174 "16",8,1,1,-1,49.3343251826364 "17",9,1,-1,-1,43.6697149526949 "18",9,1,1,-1,47.0506641010342 "19",10,1,-1,-1,45.3395919514918 "20",10,1,1,-1,48.6456098943332 "21",11,1,-1,-1,44.7191120116416 "22",11,1,1,-1,47.6787837956721 "23",12,1,-1,-1,42.574882443111 "24",12,1,1,-1,46.5033799702558 "25",13,1,-1,-1,45.5942908784639 "26",13,1,1,-1,46.6163292019490 "27",14,1,-1,-1,41.3926590716193 "28",14,1,1,-1,45.9854614499554 "29",15,1,-1,-1,46.5254842348644 "30",15,1,1,-1,49.2734985879016 "31",16,1,-1,-1,45.6755380933392 "32",16,1,1,-1,51.9223064720151 "33",17,1,-1,-1,43.4295061026270 "34",17,1,1,-1,46.4532673614053 "35",18,1,-1,-1,42.5950564348038 "36",18,1,1,-1,45.0631534193601 "37",19,1,-1,-1,44.3736086293811 "38",19,1,1,-1,48.5882675006395 "39",20,1,-1,-1,44.1957249299453 "40",20,1,1,-1,46.4500267334107 "41",21,1,-1,-1,48.1591332470503 "42",21,1,1,-1,50.8229744661494 "43",22,1,-1,-1,44.5018680448579 "44",22,1,1,-1,46.4587625288197 "45",23,1,-1,-1,44.7030305150597 "46",23,1,1,-1,48.6532707122201 "47",24,1,-1,-1,43.5459931026717 "48",24,1,1,-1,47.6356177630651 "49",25,1,-1,-1,41.8399591786956 "50",25,1,1,-1,46.8069821855383 "51",26,1,-1,-1,44.1651824914932 "52",26,1,1,-1,47.5448386968692 "53",27,1,-1,-1,40.3324483325434 "54",27,1,1,-1,44.825438710612 "55",28,1,-1,-1,46.8166603188741 "56",28,1,1,-1,49.0892308043607 "57",29,1,-1,-1,44.3224523205635 "58",29,1,1,-1,47.270775923145 "59",30,1,-1,-1,48.8009313925536 "60",30,1,1,-1,49.4285855123831 "61",31,1,-1,1,49.5652653624918 "62",31,1,1,1,56.644638206898 "63",32,1,-1,1,48.8987022278928 "64",32,1,1,1,53.4153955207579 "65",33,1,-1,1,50.1327002968025 "66",33,1,1,1,56.5215711024782 "67",34,1,-1,1,47.0853189927834 "68",34,1,1,1,54.4921427423021 "69",35,1,-1,1,47.6719233208383 "70",35,1,1,1,52.8554367239123 "71",36,1,-1,1,49.111827773968 "72",36,1,1,1,56.493182088657 "73",37,1,-1,1,48.1116745312135 "74",37,1,1,1,55.0787395153315 "75",38,1,-1,1,49.0183004093922 "76",38,1,1,1,55.6325519444112 "77",39,1,-1,1,49.6556999705135 "78",39,1,1,1,55.2358303483103 "79",40,1,
[R] Query: lme output
Dear R-Users I have a problem accessing some values in the output from the summary of an lme fit. I fit the model below: ggg <- lme (ST~ -1 + as.factor(endp):Z.sas + as.factor(endp), data=dat4a, random=~-1 + as.factor(endp) + as.factor(endp):Z.sas|as.factor(trials), correlation = corSymm(form=~1|as.factor(trials)/as.factor(id)), weights=varIdent(form=~1|endp)) hh <- summary(ggg) hh Below is the following part of the output of interest: Correlation Structure: General Formula: ~1 | as.factor(trials)/as.factor(id) Parameter estimate(s): Correlation: 1 2 0.785 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | endp Parameter estimates: -1 1 1.000 0.9692405 I wish to access the value of the correlation (0.785) and the vector of the variance function estimates (1,0.969). I know these can be done throught the intervals function, but sometimes when the estimated Hessian matrix is not positive definite or something like that (i am not quite sure), the intervals function delivers an error message. Thus, i will like to ask if there is another way to access these values. I tried using the following code: hh$modelStruct$corStruct[1] hh$modelStruct$varStruct[1] Rather the output was: > hh$modelStruct$corStruct[1] [1] -1.308580 > hh$modelStruct$varStruct[1] [1] -0.03124255 I presume there is a way to calculate the correlation and variance function coefficients using these values. Could someone tell me how to access those values (without using the intervals function) or better still how to calculate the values from the last output values. Kind regards Pryseley - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] query: lme
Good R-users, I have difficulties accessing the variance components for an lme fit when the variance covariance matrix of the random effects is not positive definite. For example, i fit the following model: ggg <- lme (ST~ -1 + as.factor(endp):Z.sas + as.factor(endp), data=dat2a, random=~-1 + as.factor(endp) + as.factor(endp):Z.sas|as.factor(trials), correlation = corSymm(form=~1|as.factor(trials)/as.factor(id)), weights=varIdent(form=~1|endp)) intervals(ggg, which="var-cov") when i try to access the variance components using the 'intervals' function i get the following error message: "Error in intervals.lme(ggg, which = "var-cov") : Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance" Is there a way out of this? or better still Is there another function through which i can access these variance components other than the intervals function? Kind regards Pryseley - Sneak preview the all-new Yahoo.com. It's not radically different. Just radically better. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] query: lme
Dear R Users I have difficulties accessing the variance components for an lme fit when the variance covariance matrix of the random effects is not positive definite. Can anyone inform me on how to get by this ? Thanks in advance Pryseley - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Pairewise Likelihood
Dear R-users Can anyone inform me of a library or more specifically functions that can maximise (or calculate) a Pairwsie likelihood from a data. Better still, i would like to know if there is a function (library) that fits regression models based on pairwise likelihoods. Thanks Pryseley - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Log Cholesky parametrization in lme
Dear R-Users I used the nlme library to fit a linear mixed model (lme). The random effect standard errors and correlation reported are based on a Log-Cholesky parametrization. Can anyone tell me how to get the Covariance matrix of the random effects, given the above mentioned parameters based on the Log-Cholesky parametrization?? Thanks in advance Pryseley - Find great deals to the top 10 hottest destinations! [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] lme and gls : accessing values from correlation structure and variance functions
Dear R-users I am relatively new to R, i hope my many novice questions are welcome. I have problems accessing some objects (specifically the random effects, correlation structure and variance function) from an object of class gls and lme. I used the following models: yah <- gls (outcome~ -1 + as.factor(Trial):as.factor(endpoint)+ as.factor(Trial):as.factor(endpoint):trt, data=datt[datt$Trial<4,], correlation = corSymm(form=~1|as.factor(Trial)/as.factor(subject)), weights=varIdent(form=~1|endpoint)) bm <- lme (outcome~ -1 + as.factor(endpoint)+ as.factor(endpoint):trt, data=datt[datt$Trial<4,], random=~-1 + as.factor(endpoint) + as.factor(endpoint):trt|as.factor(Trial), correlation = corSymm(form=~1|as.factor(Trial)/as.factor(subject)), weights=varIdent(form=~1|endpoint)) When i print the object "bm" i get the following output: -- > bm Linear mixed-effects model fit by REML Data: datt[datt$Trial < 4, ] Log-restricted-likelihood: -52.23147 Fixed: outcome ~ -1 + as.factor(endpoint) + as.factor(endpoint):trt as.factor(endpoint)-1 as.factor(endpoint)1 as.factor(endpoint)-1:trt -3.663087 -1.772427 -3.661823 as.factor(endpoint)1:trt -3.209671 Random effects: Formula: ~-1 + as.factor(endpoint) + as.factor(endpoint):trt | as.factor(Trial) Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr as.factor(endpoint)-1 2.05744327 as.()-1 as.()1 a.()-1: as.factor(endpoint)1 0.08400874 -0.976 as.factor(endpoint)-1:trt 1.90318009 0.975 -0.967 as.factor(endpoint)1:trt 3.25432832 -0.992 0.982 -0.972 Residual 6.48819860 Correlation Structure: General Formula: ~1 | as.factor(Trial)/as.factor(subject) Parameter estimate(s): Correlation: 1 2 0.812 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | endpoint Parameter estimates: 1 -1 1.00 1.764878 Number of Observations: 18 Number of Groups: 3 Can somebody tell me how to access the values in blue above? Also, when i tried accessing these values i obtained the following bm$modelStruct corStruct parameters: [1] -1.394879 varStruct parameters: [1] 0.5680815 What do these values represent. Thanks in advance.. Pryseley sample data: RowNames Trial subject VISUAL0 TRT VISUAL24 VISUAL52 TREAT outcome endpoint trt 4 11003 65 4 65 55 2 01 1 8 11007 67 1 64 68 2 -31 -1 12 21110 59 4 53 42 2 -61 1 14 2 64 1 72 65 2 81 -1 16 21112 39 1 37 37 2 -21 -1 18 21115 59 4 54 58 2 -51 1 24 31806 46 4 27 24 2 -191 1 26 31813 31 4 33 48 2 21 1 28 31815 64 1 67 64 2 31 -1 4 11003 65 4 65 55 2 -10 -1 1 8 11007 67 1 64 68 2 1 -1 -1 12 21110 59 4 53 42 2 -17 -1 1 14 2 64 1 72 65 2 1 -1 -1 16 21112 39 1 37 37 2 -2 -1 -1 18 21115 59 4 54 58 2 -1 -1 1 24 31806 46 4 27 24 2 -22 -1 1 26 31813 31 4 33 48 2 17 -1 1 28 31815 64 1 67 64 2 0 -1 -1 - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Three level linear mixed models
Hello R-users Is it possible to fit a three level linear mixed effect model in R? If anyone has an idea or sample code, i will appreciate it very much if i can receive it. I am reading the book by Pinheiro and Bates but have not come across that yet! Kind regards Pryseley - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Linear mixed model with correlated residual
Dear R- users Does any one know how to fit a linear mixed model such that the residuals ( grouped by a variable say "gender") are correlated and have a covariance matrix (in this case a 2 by 2 covariance matrix). Thanks in advance Pryseley - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Help with lme and correlated residuals
Dear R - Users I have some problems fitting a linear mixed effects model using the lme function (nlme library). A sample data is as shown at the bottom of this mail. I fit my linear mixed model using the following R code: bmr <-lme (outcome~ -1 + as.factor(endpoint)+ as.factor(endpoint):trt, data=datt, random=~-1 + as.factor(endpoint) + as.factor(endpoint):trt|as.factor(Trial), correlation = corSymm(form=~subject|as.factor(endpoint)), weights=varIdent (form=~subject|endpoint)) With this code, i want to obtain random effects for each Trial. ALso my residuals are correlated, and the correlated residuals are heteroscedastic over endpoint. That is, each endpoint has a different constant variance. However, I recieve the following error message using the code above: Error in lme.formula(outcome ~ -1 + as.factor(endpoint) + as.factor(endpoint):trt, : Incompatible formulas for groups in "random" and "correlation" May someone kindly inform me of how to correct my code. Does this imply that the grouping factor in the correlation formula must be the same grouping factor in the random formula ? If so, is there a way of getting pass this restriction ? In essence, I wish to know how to fit a linear mixed model with correlated residuals (based on a variable in my model) and to obtain the correlation matrix. Best regards Pryseley Sample data: RowNames Trial subject VISUAL0 TRT VISUAL24 VISUAL52 TREAT outcome endpoint trt 4 11003 65 4 65 55 2 01 1 8 11007 67 1 64 68 2 -31 -1 12 21110 59 4 53 42 2 -61 1 14 2 64 1 72 65 2 81 -1 16 21112 39 1 37 37 2 -21 -1 18 21115 59 4 54 58 2 -51 1 24 31806 46 4 27 24 2 -191 1 26 31813 31 4 33 48 2 21 1 28 31815 64 1 67 64 2 31 -1 4 11003 65 4 65 55 2 -10 -1 1 8 11007 67 1 64 68 2 1 -1 -1 12 21110 59 4 53 42 2 -17 -1 1 14 2 64 1 72 65 2 1 -1 -1 16 21112 39 1 37 37 2 -2 -1 -1 18 21115 59 4 54 58 2 -1 -1 1 24 31806 46 4 27 24 2 -22 -1 1 26 31813 31 4 33 48 2 17 -1 1 28 31815 64 1 67 64 2 0 -1 -1 Thanks - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Help with lme code
Dear R - Users I have some problems fitting a linear mixed effects model using the lme function (from the library nlme). A sample data is as shown at the bottom of this mail. I fit my linear mixed model using the following R code: bmr <-lme (outcome~ -1 + as.factor(endpoint)+ as.factor(endpoint):trt, data=datt, random=~-1 + as.factor(endpoint) + as.factor(endpoint):trt|as.factor(Trial), correlation = corSymm(form=~subject|as.factor(endpoint)), weights=varIdent(form=~subject|endpoint)) With this code, i want to obtain random effects for each Trial. ALso my residuals are correlated, and the correlated residuals are heteroscedastic over endpoint. That is, each endpoint has a different constant variance. However, I recieve the following error message using the code above: Error in lme.formula(outcome ~ -1 + as.factor(endpoint) + as.factor(endpoint):trt, : Incompatible formulas for groups in "random" and "correlation" May someone kindly inform me how to correct this code. Does this imply that the grouping factor in the correlation formula must be the same grouping factor in the random formula ? If so, is there a way of getting pass this restriction ? Best regards Pryseley Sample data: RowNames Trial subject VISUAL0 TRT VISUAL24 VISUAL52 TREAT outcome endpoint trt 4 11003 65 4 65 55 2 01 1 8 11007 67 1 64 68 2 -31 -1 12 21110 59 4 53 42 2 -61 1 14 2 64 1 72 65 2 81 -1 16 21112 39 1 37 37 2 -21 -1 18 21115 59 4 54 58 2 -51 1 24 31806 46 4 27 24 2 -191 1 26 31813 31 4 33 48 2 21 1 28 31815 64 1 67 64 2 31 -1 4 11003 65 4 65 55 2 -10 -1 1 8 11007 67 1 64 68 2 1 -1 -1 12 21110 59 4 53 42 2 -17 -1 1 14 2 64 1 72 65 2 1 -1 -1 16 21112 39 1 37 37 2 -2 -1 -1 18 21115 59 4 54 58 2 -1 -1 1 24 31806 46 4 27 24 2 -22 -1 1 26 31813 31 4 33 48 2 17 -1 1 28 31815 64 1 67 64 2 0 -1 -1 - Bring photos to life! New PhotoMail makes sharing a breeze. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Help with classes and generic function
Hello R-Users
I am pretty new to R and forgive me if my questions are "childish" ,
nevertheless i need help.
I have some problems while trying to grasp the idea of a generic function. I
understand that i have to build methods with naming syntax "function.class",
where class is the CLASS of an (first) object passed to the generic function,
and is used to decide which method shall be invoked.
Suppose i have a function say,
unis <- function(x,y,...){ ..function body...}
and i call the function, assigning the result to an object say,
results <-unis(outcome, trial, ..)
now i want the object, results, to have a particular CLASS say "surr" so that
i can define methods like
summary.surr
plot.surr
My main trouble is:
what do i need to do, so that any object assigned the value from a call of
the function unis has CLASS as "surr"
Thanks
Pryseley
-
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Help with functions
Dear R-users
I intend to create a function which calls some smaller other functions in
return. Some of these smaller functions all call some functions. I do not know
a good way to do this. I tried using the source() function to include the
smaller functions within the main functions before they are called. This does
not work, or maybe i am not doing the right thing.
For example:
the following quantities are computed in the main function (only a part of
the main function is shown below)
###
# quantities for bootstrap confidence intervals
###
r.value<- cor(x,y)
npb.B <- mean(r)-r.value
npb.mean <- mean(r)
npb.var <- var(r)
r.star <- sort(r)
z.star <- (r.star - mean(r.star))/sqrt(var(r.star))
after which the main function calls the function (i.s.conf.int) shown below,
which is saved in a different script file
###
# Confidence Interval
###
i.s.conf.int <- function(type) {
switch(type,
a = nci(),
b = inci(),
c = bbi(),
d = pi(),
e = si())
}
Also, this 'sub' function (i.s.conf.int) may call one of 5 other smaller
functions saved in yet another file, or the same file as the function
"i.s.conf.int".
However, these 5 functions need the quantities calculated in the main
function. One of these functions is shown below:
###
# nonparametric bootstrap normal confidence interval
###
nci <-function () {
npb.n.L <- r.value - 1.96*sqrt(npb.var)
npb.n.U <- r.value + 1.96*sqrt(npb.var)
npb.n <- array(data = c( npb.n.L, npb.n.U),
dim = c(1,2),
dimnames = list("Estimates" ,c("LowerBound","UpperBound")))
npb.n.CI <- list( Type = "Normal Confidence Interval",
ConfidenceInterval= npb.n)
npb.n.CI
}
I use the source() function to include these "secondary/sub" functions, just
before they are called the main function. I get an Error message when any of
the 5 small functions is executing : that it can not find the quantities
calculated in the main function.
So my question boils down to how can i call functions from different files in
the main function, so that the recognise quantities already calculated in the
main function or passed as arguements to the main function ?
Thanks
Pryseley.
-
Bring words and photos together (easily) with
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Help with R: functions
Dear R-users Thanks for all the help/suggestions, they have been most helpful. Best regards Pryseley - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Help with R: functions
Hello R-users
I am new to R and trying to write some functions. I have problems writing
functions that takes a data set as an arguement and uses variables in the data.
I illustrate my problem with a small example below:
sample data #--
visual24<-rnorm(30,3,5)
visual52<-rt(30,7)
dats<- data.frame(cbind(visual24,visual52))
remove(visual24, visual52)
# first code
#--
st <-function(data,x,y){
rcc<-coef(lm(y~x))
plot(x,y)
abline(rcc[1],rcc[2])
}
st(data=dats,x=dats$visual24,y=dats$visual52)
This code works fine, but with such a code the data as an arguement to the
funtion is not necessary.
However, i wish to write a function that reads the variables from the data
directly.
I tried using the function below but it does not work.
# second code
#--
st <-function(data,x,y){
rcc<-coef(lm(data$y~data$x))
plot(data$x,data$y)
abline(rcc[1],rcc[2])
}
st(dats,visual24,visual52)
I wish to inquire if any one has an idea of what i need to adjust in the
function so that it works.
I believe that the referencing $x or $y in the function is not doing the
correct thing.
Better still, will it be a problem if i code the functions as in the first
code above?
I mean given that they will be used to create a library
Best regards
Pryseley
-
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Help with mixed effects models
Dear R-users I have problems using lme The model i want to fit can be viewed as a two-level bivariate model Two-level bivariate: bivariate (S coded as -1,T coded as 1) endpoint within trial OR It can equivalently be considered as a three-level model.Three-level: endpoint within patient, patient within trial. My code tries to model the levels through a RANDOM statement and a within group correlation structure. -- Now then, i used the following R code: bm <- lme (outcome~ -1 + as.factor(endpoint)+ as.factor(endpoint):trt, data=datt, random=~-1 + as.factor(endpoint) + as.factor(endpoint):trt |as.factor(Trial), corr = corSymm(form~-1+as.factor(endpoint)|trial/subject)) I beleive the fixed effects part of the code is okay. My intention for the random effects part is to estimate an intercept and treatment effect for each endpoint at the trial level. The correlation structure should produce a within correlation matrix for the enpoints at the subject level. Thus the random effects matrix is 4 by 4 and the within correlation matrix is 2 by 2 When i run the code in R, i get the following error message "Error in Initialize.corSymm(X[[2]], ...) : Initial value for corSymm parameters of wrong dimension" I hope someone will correct my codes . Kind regards Pryseley - Ring in the New Year with Photo Calendars. Add photos, events, holidays, whatever. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
