[R] R2WinBUGS with Multivariate Logistic Regression
Dear R-User,
I have written a simple code to analyze some data using Bayesian logistic
regression via the R2WinBUGS package. The code when run in WinBUGS stops
WinBUGS from running it and using the package returns no results also.
I attach herewith, the code and a sample of the dataset.
Any suggestion will be greatly appreciated.
Chris Guure
Biostatistics Department
University of Ghana
library(R2WinBUGS)
library(coda)model1<-function(){
for (i in 1:N) {
# likelihood function
ms[i] ~ dbin( p [ i ], N )
logit(p [ i ] ) <- alpha + bage*age[ i ] + bpam*pam[i ] + bpah*pah[i]
}
### prior for intercept
alpha ~ dnorm(0,0.0001)
# prior for slopes
bage ~ dnorm(0,0.0001)
bpam ~ dnorm(0,0.0001)
bpah ~ dnorm(0,0.0001)
# OR for alpha
or.age<-exp(bage)
# OR for hbp
or.pam <- exp(bpam)
# OR for fdm
or.pah <- exp(bpah)
}
data=cbind(ms=c(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0),
age=c(77, 83, 75, 78, 75, 83, 85, 80, 80, 85, 76, 77, 80, 76, 88, 77, 81, 78,
85, 81),
pam=c(0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1),
pah=c(1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0))
N=20
initial values
###
lineinits <- function(){
list(alpha=1,bage= 0.05, bpam=0.05, bpah=0.05)
}
## CODA output
lineout1 <- bugs(data, lineinits, c("bage", "alpha","bpam", "bpah", "or.age",
"or.pam", "or.pah"), model1,n.iter = 11000, n.burnin = 1000, n.chains = 2,
codaPkg = T, DIC = TRUE)
### Posterior summaries
line.coda <- read.bugs(lineout1)
summary(line.coda)
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Clarification on Simulation and Iteration
Thanks Dave.
What I actually want is to obtain say 10, different sets of (n=50) data for
every 10,000 iterations I run. You will realise that the current code produces
one set of data (n=50). I want 10 different sets of 50 observations at one run.
I hope this makes sense.
Chris Guure
On Saturday, August 1, 2015 3:32 AM, David Winsemius dwinsem...@comcast.net
wrote:
On Jul 31, 2015, at 6:36 PM, Christopher Kelvin via R-help wrote:
Dear All,
I am performing some simulations for a new model. I run about 10,000
iterations with a sample of 50 datasets and this returns one set of 50
simulated data.
Now, what I need to obtain is 10 sets of the 50 simulated data out of the
10,000 iterations and not just only 1 set. The model is the Copas selection
model for publication bias in Mete-analysis. Any one who knows this model has
any suggestion for the improvement of my code is most welcome.
Below is my code.
Kind regards
Chris Guure
University of Ghana
install.packages(msm)
library(msm)
rho1=-0.3; tua=0.020; n=50; d=-0.2; rr=1; a=-1.3; b=0.06
si-rtnorm(n, mean=0, sd=1, lower=0, upper=0.2)# I used this to generate
standard errors for each study
set.seed(2) ## I have stored the data and the output in this seed
for( i in 1:rr){
mu-rnorm(n,d,tua^2) # prob. of each effect estimate
rho-si*rho1/sqrt(tua^2 + si^2) # estimate of the correlation coefficient
mu0- a + b/si # mean of the truncated normal model (Copas selection
model)
y1-rnorm(mu,si^2)# observed effects zise
z-rnorm(mu0,1) # selection model
rho2-cor(y1, z)
select-pnorm((mu0 + rho*(y1-mu)/sqrt(tua^2 + si^2))/sqrt(1-rho^2))
probselect-ifelse(selectz, y1, NA)# the prob that the study is selected
probselect
data-data.frame(probselect,si)# this contains both include and exclude
data
data
data1-data[complete.cases(data),] # Contains only the included data for
analysis
data1
}
OK. The code runs without error. So what exactly is the problem? (I have
no experience with the Copas selection model if in fact that is what is being
exemplified.)
--
David Winsemius
Alameda, CA, USA
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Clarification on Simulation and Iteration
Dear All,
I am performing some simulations for a new model. I run about 10,000 iterations
with a sample of 50 datasets and this returns one set of 50 simulated data.
Now, what I need to obtain is 10 sets of the 50 simulated data out of the
10,000 iterations and not just only 1 set. The model is the Copas selection
model for publication bias in Mete-analysis. Any one who knows this model has
any suggestion for the improvement of my code is most welcome.
Below is my code.
Kind regards
Chris Guure
University of Ghana
install.packages(msm)
library(msm)
rho1=-0.3; tua=0.020; n=50; d=-0.2; rr=1; a=-1.3; b=0.06
si-rtnorm(n, mean=0, sd=1, lower=0, upper=0.2)# I used this to generate
standard errors for each study
set.seed(2) ## I have stored the data and the output in this seed
for( i in 1:rr){
mu-rnorm(n,d,tua^2) # prob. of each effect estimate
rho-si*rho1/sqrt(tua^2 + si^2) # estimate of the correlation coefficient
mu0- a + b/si # mean of the truncated normal model (Copas selection
model)
y1-rnorm(mu,si^2)# observed effects zise
z-rnorm(mu0,1) # selection model
rho2-cor(y1, z)
select-pnorm((mu0 + rho*(y1-mu)/sqrt(tua^2 + si^2))/sqrt(1-rho^2))
probselect-ifelse(selectz, y1, NA)# the prob that the study is selected
probselect
data-data.frame(probselect,si)# this contains both include and exclude
data
data
data1-data[complete.cases(data),] # Contains only the included data for
analysis
data1
}
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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 with Newton_Raphson
Hello,
I have being trying to estimate the parameters of the generalized exponential
distribution. The random number generation for the GE distribution
is x-(-log(1-U^(1/p1))/b), where U stands for uniform dist. The data i have
generated to estimate the parameters is right censored and the code is given
below; The problem is that, the newton-Raphson approach isnt working and i do
not know what is wrong. Can somebody suggest something or help in identifying
what the problem might be.
p1-0.6;b-2
n=20;rr=5000
U-runif(n,0,1)
for (i in 1:rr){
x-(-log(1-U^(1/p1))/b)
meantrue-gamma(1+(1/p1))*b
meantrue
d-meantrue/0.30
cen- runif(n,min=0,max=d)
s-ifelse(x=cen,1,0)
q-c(x,cen)
z-function(data, p){
shape-p[1]
scale-p[2]
log1-n*sum(s)*log(p[1])+
n*sum(s)*log(p[2])+(p[1]-1)*sum(s)*log(1-((exp(-(p[2])*sum(x)
-(p[2])*sum(t) + (p[1])*log((exp(-(p[2])*sum(x-
(p[1])*sum(s)*log((exp(-(p[2])*sum(x
return(-log1)
}
}
start - c(1,1)
zz-optim(start,fn=z,data=q,hessian=T)
zz
m1-zz$par[2]
p-zz$par[1]
Thank you
Chris Kelvin
INSPEM. UPM
__
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] Problem with Newton_Raphson
By your suggestion if i understand you, i have changed (p[1]) and (p[2]) and
have also corrected the error sum(t), but if i run it, the parameters estimates
are very large.
Can you please run it and help me out? The code is given below.
p1-0.6;b-2
n=20;rr=5000
U-runif(n,0,1)
for (i in 1:rr){
x-(-log(1-U^(1/p1))/b)
meantrue-gamma(1+(1/p1))*b
meantrue
d-meantrue/0.30
cen- runif(n,min=0,max=d)
s-ifelse(x=cen,1,0)
q-c(x,cen)
}
z-function(data, p){
shape-p[1]
scale-p[2]
log1-n*sum(s)*log(shape)+
n*sum(s)*log(scale)+(shape-1)*sum(s)*log(1-((exp(-(scale)*sum(x)
-(scale)*sum(x) + (shape)*log((exp(-(scale)*sum(x-
(shape)*sum(s)*log((exp(-(scale)*sum(x
return(-log1)
}
start - c(1,1)
zz-optim(start,fn=z,data=q,hessian=T)
zz
Thank you
Chris Kelvin
- Original Message -
From: Berend Hasselman b...@xs4all.nl
To: Christopher Kelvin chris_kelvin2...@yahoo.com
Cc: r-help@r-project.org r-help@r-project.org
Sent: Thursday, September 20, 2012 8:52 PM
Subject: Re: [R] Problem with Newton_Raphson
On 20-09-2012, at 13:46, Christopher Kelvin wrote:
Hello,
I have being trying to estimate the parameters of the generalized exponential
distribution. The random number generation for the GE distribution is
x-(-log(1-U^(1/p1))/b), where U stands for uniform dist. The data i have
generated to estimate the parameters is right censored and the code is given
below; The problem is that, the newton-Raphson approach isnt working and i do
not know what is wrong. Can somebody suggest something or help in identifying
what the problem might be.
Newton-Raphson? is not a method for optim.
p1-0.6;b-2
n=20;rr=5000
U-runif(n,0,1)
for (i in 1:rr){
x-(-log(1-U^(1/p1))/b)
meantrue-gamma(1+(1/p1))*b
meantrue
d-meantrue/0.30
cen- runif(n,min=0,max=d)
s-ifelse(x=cen,1,0)
q-c(x,cen)
z-function(data, p){
shape-p[1]
scale-p[2]
log1-n*sum(s)*log(p[1])+
n*sum(s)*log(p[2])+(p[1]-1)*sum(s)*log(1-((exp(-(p[2])*sum(x)
-(p[2])*sum(t) + (p[1])*log((exp(-(p[2])*sum(x-
(p[1])*sum(s)*log((exp(-(p[2])*sum(x
return(-log1)
}
}
start - c(1,1)
zz-optim(start,fn=z,data=q,hessian=T)
zz
m1-zz$par[2]
p-zz$par[1]
Running your code given above gives an error message:
Error in sum(t) : invalid 'type' (closure) of argument
Where is object 't'?
Why are you defining function z within the rr loop? Only the last definition is
given to optim.
Why use p[1] and p[2] explicitly in the calculation of log1 in the body of
function z when you can use shape and scale defined in the lines before log1 -.
Berend
__
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] Problem with Newton_Raphson
Thank you very much for everything. Your suggestions were very helpful.
Chris
- Original Message -
From: Berend Hasselman b...@xs4all.nl
To: Christopher Kelvin chris_kelvin2...@yahoo.com
Cc: r-help@r-project.org r-help@r-project.org
Sent: Thursday, September 20, 2012 10:06 PM
Subject: Re: [R] Problem with Newton_Raphson
On 20-09-2012, at 15:17, Christopher Kelvin wrote:
By your suggestion if i understand you, i have changed (p[1]) and (p[2]) and
have also corrected the error sum(t), but if i run it, the parameters
estimates are very large.
Can you please run it and help me out? The code is given below.
p1-0.6;b-2
n=20;rr=5000
U-runif(n,0,1)
for (i in 1:rr){
x-(-log(1-U^(1/p1))/b)
meantrue-gamma(1+(1/p1))*b
meantrue
d-meantrue/0.30
cen- runif(n,min=0,max=d)
s-ifelse(x=cen,1,0)
q-c(x,cen)
}
z-function(data, p){
shape-p[1]
scale-p[2]
log1-n*sum(s)*log(shape)+
n*sum(s)*log(scale)+(shape-1)*sum(s)*log(1-((exp(-(scale)*sum(x)
-(scale)*sum(x) + (shape)*log((exp(-(scale)*sum(x-
(shape)*sum(s)*log((exp(-(scale)*sum(x
return(-log1)
}
start - c(1,1)
zz-optim(start,fn=z,data=q,hessian=T)
zz
1. You think you are using a (quasi) Newton method. But the default method is
Nelder-Mead. You should specify method explicitly for example method=BFGS.
When you do that you will see that the results are completely different. You
should also carefully inspect the return value of optim. In your case
zz$convergence is 1 which means indicates that the iteration limit maxit had
been reached..
When you use method=BFGS you will get zz$ convergence is 0.
This may do what you want
zz-optim(start, fn=z, data=q, method=BFGS, hessian=T)
zz
2. when providing examples such as yours it is a good idea to issue
set.seed(number) at the start of your script to ensure reproducibility.
3. The function z does not calculate what you want: since fully formed
expressions are terminated by a newline and since the remainder of the line
after log1- is a complete expression, log1 does include the value of the two
following lines.
See the difference between
a - 1
b - 11
a
-b
and
a-
b
So you could write this
log1 - n*sum(s)*log(shape)+ n*sum(s)*log(scale)+
(shape-1)*sum(s)*log(1-((exp(-(scale)*sum(x) -(scale)*sum(x) +
(shape)*log((exp(-(scale)*sum(x -
(shape)*sum(s)*log((exp(-(scale)*sum(x
and you will see rather different results.
You will have to determine if they are what you expect/desire.
A final remark on function z:
- do not calculate things like n*sum(s) repeatedly: doing something like
A-n*sum(s) and reusing A is more efficient.
- same thing for log((exp(-(scale)*sum(x (recalculated four times)
- same thing for sum(x)
See below for a partly rewritten function z and results with method=BFGS.
I have put a set.seed(1) at the start of your code.
good luck,
Berend
z-function(p,data){
shape-p[1]
scale-p[2]
log1 - n*sum(s)*log(shape)+ n*sum(s)*log(scale)+
(shape-1)*sum(s)*log(1-((exp(-(scale)*sum(x) -(scale)*sum(x) +
(shape)*log((exp(-(scale)*sum(x -
(shape)*sum(s)*log((exp(-(scale)*sum(x
return(-log1)
}
start - c(1,1)
z(start, data=q)
[1] -138.7918
zz-optim(start, fn=z, data=q, method=BFGS, hessian=T)
There were 34 warnings (use warnings() to see them)
zz
$par
[1] 1.009614e+11 8.568498e+01
$value
[1] -1.27583e+15
$counts
function gradient
270 87
$convergence
[1] 0
$message
NULL
$hessian
[,1] [,2]
[1,] -62500 0
[2,] 0 0
__
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] Optim Problem
Hello,
I want to estimate the exponential parameter by using optim with the following
input, where t contains 40% of the data and q contains 60% of the data within
an interval. In implementing the code command for optim i want it to contain
both the t and q data so i can obtain the correct estimate. Is there any
suggestion as to how this can be done. I have tried h-c(t,q) but it is not
working because q lies within an interval.
rate-15;n-100;a-40;b-60;rr-1000
ms11=ms22=0
bia11=bia22=0
t-rexp(a,rate)
for(i in 1:rr){
C1-runif(b,0,rate)
C2-rexp(b,rate)
f2 - function(C1, C2) {
r - pmax(C1 , C2 + C1)
cbind(C1, r)
}
m-f2(C1,C2)
x[1:b]-(m[,1])
u-x[1:b]
x[1:b]-(m[,2])
v-x[1:b]
q-cbind(u,v)
h-c(t,q)
z-function(data ){
rate-p[2]
log1--(n/log(p[2]))-sum(t/(p[2]))+sum(log(exp(-(u/(p[2])))-exp(-(v/(p[2])
return(-log1)
}
}
start - c(1,1)
zz-optim(start,fn=z,data=h,hessian=T)
m1-zz$par[2]
thank you
chris b guure
researcher
institute for mathematical research
upm
__
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 with ifelse
Dear all,
The code below is used to generate interval censored data but unfortunately
there is an error with the ifelse which i am not able to rectify.
Can somebody help correct it for me.
Thank you
t-rexp(20,0.2)
v-c(0,m,999)
y-function(t,v){
z-numeric(length(t ((
s-numeric(length(t ((
for(i in 1:length(t)){
for(j in 1:length(v-1))
{ ifelse ((t[i]v[j] t v[j+1] ),{z[i]-v[j];s[i]-v[j+1]},NA)}}
return(cbind(z,s))}
y(t,v)
Chris Kelvin
Institute for Mathematical Research
UPM
__
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] standard error
Dear all, I want to determine the standard error or the mean squared error for the parameter estimate for beta and eta base on the real data. Any help on how to obtain these estimated errors. library(survival) d - data.frame(ob=c(149971, 70808, 133518, 145658, 175701, 50960, 126606, 82329), state=1) s - Surv(d$ob,d$state) sr - survreg(s~1,dist=weibull) beta-1/sr$scale p1=(beta) p1 eta-exp(sr$coefficients[1]) b=(eta) b Thank you Chris Guure Researcher Institute for Mathematical Research UPM __ 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 replace NA with zero (0)
Hello,
When i generate data with the code below there appear NA as part of the
generated data, i prefer to have zero (0) instead of NA on my data.
Is there a command i can issue to replace the NA with zero (0) even if it is
after generating the data?
Thank you
library(survival)
p1-0.8;b-1.5;rr-1000
for(i in 1:rr){
r-runif(45,min=0,max=1)
t-rweibull(45,p1,b)
w=Surv(r,t,type=interval2)
x[1:45]-(w[,1])
u-x[1:45]
y[1:45]-(w[,2])
v-y[1:45]
}
w
u
v
Chris G
Researcher
Institute for Mathematical Research
UPM
__
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] interval censoring
Hello,
May i know whether it is possible to generate data twice from Weibull
distribution and use one as the start time and the
other as the end time, below is my code.
Any suggestion on how to estimate the parameters of Weibull distribution with
interval data will be highly appreciated.
Thank you
shape=2;scale=4;rr=5000
for(i in 1:rr){
x-rweibull(50,shape,scale)
y=rweibull(50,shape,scale)
w=Surv(y,x,type=interval2)
}
w
Chris Guure
Researcher
Institute for Mathematical Research
UPM
__
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] standard error
Hello,
I have tried obtaining the value of standard error from the code below but i
get different values when i compare it with the
standard error obtained from the hessian matrix. Can somebody help me out?
Thank you
n=100;rr=1000
p1=1.2;b=1.5
sq11=sq21=0
for (i in 1:rr){
t-rweibull(n,shape=p1,scale=b)
meantrue-gamma(1+(1/p1))*b
meantrue
d-meantrue/0.40
cen- runif(n,min=0,max=d)
s-ifelse(t=cen,1,0)
q-c(t,cen)
z-function(data, p){
beta-p[1]
eta-p[2]
log1-(n*sum(s)*log(p[1])-n*sum(s)*(p[1])*log(p[2])+sum(s)*(p[1]-1)*sum(log(t))-n*sum((t/(p[2]))^(p[1])))
return(-log1)
}
start - c(1,1)
zz-optim(start,fn=z,data=q,hessian=T)
zz
m1-zz$par[2]
p-zz$par[1]
sq11-sq11+(1/rr*(sum((q-m1)^2)))
sq21-sq21+(1/rr*(sum((q-Lm1)^2)))
}
se11-sqrt(sq11)/(n-1)
se11
se21-sqrt(sq21)/(n-1)
se21
f-solve(zz$hessian)
se-sqrt(diag(f))
se
Chris Guure
Researcher
Institute for Mathematical Research
UPM
__
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] Standard error
Hello,
I have tried obtaining the value of standard error from the code below but i
get different values when i compare it with the
standard error obtained from the hessian matrix. Can somebody help me out?
Thank you
n=100;rr=1000
p1=1.2;b=1.5
sq11=sq21=0
for (i in 1:rr){
t-rweibull(n,shape=p1,scale=b)
meantrue-gamma(1+(1/p1))*b
meantrue
d-meantrue/0.40
cen- runif(n,min=0,max=d)
s-ifelse(t=cen,1,0)
q-c(t,cen)
z-function(data, p){
beta-p[1]
eta-p[2]
log1-(n*sum(s)*log(p[1])-n*sum(s)*(p[1])*log(p[2])+sum(s)*(p[1]-1)*sum(log(t))-n*sum((t/(p[2]))^(p[1])))
return(-log1)
}
start - c(1,1)
zz-optim(start,fn=z,data=q,hessian=T)
zz
m1-zz$par[2]
p-zz$par[1]
sq11-sq11+(1/rr*(sum((q-m1)^2)))
sq21-sq21+(1/rr*(sum((q-Lm1)^2)))
}
se11-sqrt(sq11)/(rr-1)
se11
se21-sqrt(sq21)/(rr-1)
se21
f-solve(zz$hessian)
se-sqrt(diag(f))
se
Chris Guure
Researcher
Institute for Mathematical Research
UPM
__
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] Interval censorin
Hello,
May i know whether it is possible to generate data twice from Weibull
distribution and use one as the start time and the
other as the end time, below is my code.
Any suggestion on how to estimate the parameters of Weibull distribution with
interval data will be highly appreciated.
Thank you
shape=2;scale=4;rr=5000
for(i in 1:rr){
x-rweibull(50,shape,scale)
y=rweibull(50,shape,scale)
w=Surv(y,x,type=interval2)
}
w
Chris Guure
Researcher
Institute for Mathematical Research
UPM
__
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] R: Help; error in optim
Hello,
When i run the code below from Weibull distribution with 30% censoring by using
optim i get an error form R, which states that
Error in optim(start, fn = z, data = q, hessian = T) :
objective function in optim evaluates to length 25 not 1
can somebody help me remove this error. Is my censoring approach correct.
n=25;rr=1000
p=1.5;b=1.2
for (i in 1:rr){
q-c(t,cen)
t-rweibull(25,shape=p,scale=b)
meantrue-gamma(1+(1/p))*b
meantrue
d-meantrue/0.30
cen- runif(25,min=0,max=d)
cen
s-ifelse(t=cen,1,0)
z-function(data,p){
beta-p[1]
eta-p[2]
log1-(n*cen*log(p[1])-n*cen*(p[1])*log(p[2])+cen*(p[1]-1)*sum(log(t))-n*sum((t/(p[2]))^(p[1])))
return(-log1)
}
start -c(0.5,0.5)
zz-optim(start,fn=z,data=q,hessian=T)
m1-zz$par[2]
p-zz$par[1]
}
m1
p
Thank you
Chris Guure
Researcher
Institute for Mathematical Research
UPM
__
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] R-Help: Censoring data
Hello,
I want to estimate weibull parameters with 30% censored data. I have below the
code for the censoring
but how it must be put into the likelihood equation to obtain the desire
estimate is where i have a problem with,
can some body help?
My likelihood equation is for a random type-I censoring where time for the
censored units is different for each unit.
n=50;r=35
p=0.8;b=1.5
t-rweibull(50,shape=p,scale=b)
meantrue-gamma(1+(1/p))*b
meantrue
d-meantrue/0.30
cen- runif(50,min=0,max=d)
cen
s-ifelse(t=cen,1,0)
s
z-function(p){
shape-p[1]
scale-p[2]
log1-(r*log(p[1])-r*(p[1])*log(p[2])+(p[1]-1)*sum(log(t))-sum((t/(p[2]))^(p[1])
)-((n-r)*(sum(cen)/(p[2]))^(p[1])))
return(-log1)
}
start - c(1,1)
zz-optim(start,fn=z,hessian=T)
zz
Thanks in anticipation
Chris Guure
Researcher
Institute for Mathematical Research
UPM
__
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] R-help; generating censored data
Hello, can i implement this as 10% censored data where t gives me failure and x censored. Thank you p=2;b=120 n=50 set.seed(132); r-sample(1:50,45) t-rweibull(r,shape=p,scale=b) t set.seed(123); cens - sample(1:50, 5) x-runif(cens,shape=p,scale=b) x Chris Guure Researcher, Institute for Mathematical Research UPM __ 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] R-help; Censoring
Hello, I wish to censor 10% of my sample units of 50 from a Weibull distribution. Below is the code for it. I will need to know whether what i have done is correct and if not, can i have any suggestion to improve it? Thank you p=2;b=120 n=50 r=45 t-rweibull(r,shape=p,scale=b) meantrue-gamma(1+(1/p))*b meantrue cen- runif(n-r,min=0,max=meantrue) cen Chris Guure Researcher, Institute for Mathematical Research UPM __ 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] simulation
Hello,
i need to simulate 100 times, n=40 ,
the distribution has 90% from X~N(0,1) + 10% from X~N(20,10)
Is my loop below correct?
Thank you
n=40
for(i in 1:100){
x-rnorm(40,0,1) # 90% of n
z-rnorm(40,20,10) # 10% of n
}
x+z
__
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] matrix with Loop
Hello!
I got something to ask..whether you can help me with the R program...i got this
for example 5x4 matrix..and i want to find:
i) mean for each row of the matrix
ii) median for each column of the matrix
and i need to do this using a loop function...below is my program..u try to
check it for me as the output that i got is not what i desired...thanks..
data-rnorm(20,0,1)
dim(data)-c(5,4)
is.matrix(data)
data
a-function(x)
{
for(i in 1:nrow(x))
{
for(j in 1:ncol(x))
{
med-median(x[,j])
mean-mean(x[i,])
print(c(med=med,mean=mean))
}
}
}
a(data)
Chris Guure
Researcher
Institute for Mathematical Research
UPM
__
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] R- Fisher Information
Dear All,
Can you help me, with the code below how do I obtain the fisher information
from it.
Is my q-replicate(1000,x) the right way to do simulation.
thank you.
x-rweibull(100,0.8,1.5)
q-replicate(1000,x)
z-function(p){
beta-p[1]
eta-p[2]
log1-(n*log(beta)-n*beta*log(eta)+(beta-1)*sum(log(x))-sum((x/eta)^beta))
return(-log1)
}
zz-optim(c(0.5,0.5),z)
zz
Chris Guure
postgraduate researcher/tutor
Institute for Mathematical Research
Universiti Putra Malaysia
[[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] Fisher Imformation
Hello,
i have used the code below to estimate the parameters of weibull distribution
and i want to obtain the fisher information
by providing the the next code but i receive errors anytime i try to, what do i
do?
by the way is my replication correct and is it placed at the right position for
replicating x to obtain the estimates
thank you
n=100
library(survival)
x-rweibull(n,0.8,1.5)
q-replicate(1000,x)
z-function(p){
beta-p[1]
eta-p[2]
log1-(n*log(beta)-n*beta*log(eta)+(beta-1)*sum(log(x))-sum((x/eta)^beta))
return(-log1)
}
zz-optim(c(0.5,0.5),z)
zz
library(MASS)
out - nlm(z,zz,x=x hessian = TRUE)
fish - out$hessian
fish
solve(fish)
Chris Guure
postgraduate researcher/tutor
Institute for Mathematical Research
Universiti Putra Malaysia
__
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] r-help; fisher information
Hello,
i have used the code below to estimate the parameters of weibull distribution
and i want to obtain the fisher information
by providing the the next code but i receive errors anytime i try to, what do i
do?
by the way is my replication correct and is it placed at the right position for
replicating x to obtain the estimates
thank you
n=100
library(survival)
x-rweibull(n,0.8,1.5)
q-replicate(1000,x)
z-function(p){
beta-p[1]
eta-p[2]
log1-(n*log(beta)-n*beta*log(eta)+(beta-1)*sum(log(x))-sum((x/eta)^beta))
return(-log1)
}
zz-optim(c(0.5,0.5),z)
zz
library(MASS)
out - nlm(z,zz,x=x hessian = TRUE)
fish - out$hessian
fish
solve(fish)
[[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] parameter estimate
I need help, the codes below estimates the weibull parameters with complete failure, my question is how do i change the state to include some censoring (may be right, type-I or type-II) to generate and estimate the parameters. thank you x=rweibull(10,2,2) library(survival) d-data.frame(ob=c(x),state=1) s - Surv(d$ob,d$state) sr - survreg(s~1,dist=weibull) print(paste(beta =,1/sr$scale)) print(paste(eta =,exp(sr$coefficients[1]))) or library(MASS) set.seed(123) m - replicate(1000, coef(fitdistr(rweibull(50, 0.8, 2), weibull))) summary(t(m)) [[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] r-help; parameter estimate
I need help, the codes below estimates the weibull parameters with complete failure, my question is how do i change the state to include some censoring (may be right, type-I or type-II) to generate and estimate the parameters. thank you x=rweibull(10,2,2) library(survival) d-data.frame(ob=c(x),state=1) s - Surv(d$ob,d$state) sr - survreg(s~1,dist=weibull) print(paste(beta =,1/sr$scale)) print(paste(eta =,exp(sr$coefficients[1]))) or library(MASS) set.seed(123) m - replicate(1000, coef(fitdistr(rweibull(50, 0.8, 2), weibull))) summary(t(m)) [[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] r-help; weibull parameter estimate
Hello,
If i write a function as below using log of weibull distribution i do not get
the required
results in estimating the parameters what do i do, please
a/b * (t/b)^a-1 * exp(-t/b)^a
n=500
x-rweibull(n,2,2)
z-function(p) {(-n*log(p[1])+n*log(p[2])-
(p[1]-1)*sum(log(x))+(p[1]-1)*log(p[2])+(sum(x/p[2])^(p[1])) )}
zz-optim(c(0.5,0.5),z)
zz
[[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] r-help; weibull distribution
Please, Help me, How do I generate data from the weibull distribution if the data contain both failure and interval censored, For example, I want to generate n=100, shape=2 and scale =4 with 30% interval censored. What about right censoring Thank you [[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] R-HELP
please help; I want to know how to generate an interval-censored data of about 20% and a right censored data of about 30% using the weibull distribution of say, x=rweibull(100,shape=1.2,scale=1.5) [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] R-Help
Can somebody help me, How do I generate data from the weibull distribution if the data contain both failure and interval censored, For example, I want to generate n=100, shape=2 and scale =4 with 30% interval censored. Thank you [[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] R-help mailing list submissions
R-help mailing list submissions [[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.
