Nothing
mywhiskers<-function(x, y, nclass=10, limits=NA,
add=FALSE, se=TRUE, main="", xlab="x", ylab="y",
ylim=NA, lwd=1, highlight="red") {
away<-is.na(x+y)
x<-x[!away]
y<-y[!away]
means<-ses<-NA
nclass<-nclass+1
if (is.na(limits[1])) {
while (sum(is.na(means+ses))&nclass>1) { # decrease the number of classes if there are too few observations
nclass<-nclass-1
# limits<-seq(min(x),max(x)+1e-10,length=nclass+1)
limits<-quantile(x,probs=seq(0,1,length=nclass+1))
limits[nclass+1]<-limits[nclass+1]+1e-8
means<-sapply(1:nclass,function(i) mean(y[x>=limits[i]&x<limits[i+1]]))
if (se) {
ses<-sapply(1:nclass,function(i) sd(y[x>=limits[i]&x<limits[i+1]])/
sqrt(sum(x>=limits[i]&x<limits[i+1])))
} else {
ses<-sapply(1:nclass,function(i) sd(y[x>=limits[i]&x<limits[i+1]]))
}
}
} else {
nclass=length(limits)-1
means<-sapply(1:nclass,function(i) mean(y[x>=limits[i]&x<limits[i+1]]))
if (se) {
ses<-sapply(1:nclass,function(i) sd(y[x>=limits[i]&x<limits[i+1]])/
sqrt(sum(x>=limits[i]&x<limits[i+1])))
} else {
ses<-sapply(1:nclass,function(i) sd(y[x>=limits[i]&x<limits[i+1]]))
}
}
lb<-means-1.96*ses
ub<-means+1.96*ses
xclass<-1/2*(limits[-1]+limits[-nclass-1])
if (add) {
points(xclass,means)
} else {
if(is.na(ylim[1])) ylim<-c(min(lb),max(ub))
plot(xclass,means,ylim=ylim,main=main,xlab=xlab,ylab=ylab,xlim=range(x))
}
color<-rep("black",nclass)
if (sum(ub*lb>0)) color[ub*lb>0]<-highlight
sapply(1:nclass,function(i) lines(xclass[c(i,i)],c(lb[i],ub[i]),lwd=lwd,col=color[i],lend="butt"))
list(x=xclass,m=means,s=ses,lb,ub)
}
# This is for plotting the response against age and connecting the observations of the same subject
# in a longitudinal study.
# in a longitudinal study.
#linesplot<-function(x, y, group,
# xlab="x", ylab="y", main="", cex=0.5, pch=19,
# col=1, col.lin=1, lw=FALSE, ylim=NULL, xlim=NULL,
# add=FALSE, lty="solid",lwd=1) {
# group<-group[order(x)]
# y<-y[order(x)]
#
# col<-rep(col,length(x))
# if (length(col.lin)==1) col.lin=rep(col.lin,length(x))
# col.lin<-col.lin[order(x)]
#
# if (length(lty)==1) lty=rep(lty,length(x))
# lty<-lty[order(x)]
#
# if (length(lwd)==1) lwd=rep(lwd,length(x))
# lwd<-lwd[order(x)]
#
# x<-x[order(x)]
# if (!add) {
# plot(x,y,type="n",xlab=xlab,ylab=ylab,main=main,cex=cex,pch=pch,col=col,ylim=ylim,xlim=xlim)
# }
# apu<-unique(group)
# for (i in 1:length(apu)) {
# print(col.lin[group==apu[i]])
# lines(x[group==apu[i]],y[group==apu[i]],col=col.lin[group==apu[i]],lty=lty[group==apu[i]],
# lwd=lwd[group==apu[i]])
# }
# points(x,y,col=col,cex=cex)
# if (lw) {
# apu<-unique(col)
# for (i in 1:length(apu))
# lines(lowess(x[col==apu[i]],y[col==apu[i]],f=0.5,iter=5),col=apu[i],lwd=3)
# }
#}
# Mikan paivitys 3.4.2018.
# Laurin paivitys 18.1.2018 (split-funktion kaytto nopeutti toimintaa)
linesplot<-function (x, y, group, xlab = "x", ylab = "y", main = "", cex = 0.5,
pch = 19, col = 1, col.lin = 1, lw = FALSE, ylim = NULL,
xlim = NULL, add = FALSE, lty = "solid", lwd = 1)
{
xorig<-x
yorig<-y
group <- group[order(x)]
y <- y[order(x)]
# if (length(col) == 1)
# col <- rep(col, length(x))
# col <- col[order(x)]
if (length(col.lin) == 1)
col.lin = rep(col.lin, length(x))
col.lin <- col.lin[order(x)]
if (length(lty) == 1)
lty = rep(lty, length(x))
lty <- lty[order(x)]
if (length(lwd) == 1)
lwd = rep(lwd, length(x))
lwd <- lwd[order(x)]
x <- x[order(x)]
if (!add) {
plot(xorig, yorig, type = "n", xlab = xlab, ylab = ylab, main = main,ylim = ylim, xlim = xlim)
}
sx<-split(x,group)
sy<-split(y,group)
# scol<-split(col,group)
scol.lin<-split(col.lin,group)
slty<-split(lty,group)
slwd<-split(lwd,group)
for (i in 1:length(sx)) {
lines(sx[[i]], sy[[i]], col = scol.lin[[i]], lty = slty[[i]], lwd = slwd[[i]])
}
points(xorig, yorig, col = col, cex = cex, pch = pch)
if (lw) {
apu <- unique(col)
for (i in 1:length(apu)) lines(lowess(x[col == apu[i]],
y[col == apu[i]], f = 0.5, iter = 5), col = apu[i], lwd = 3)
}
}
circle<-function(x,y,r,border="black",lty="solid",lwd=1,fill=NULL) {
xapu<-sin(seq(0,pi,length=50)-pi/2)
for (i in 1:length(x)) {
xv1<-x[i]+xapu*r[i]
xv2<-x[i]-xapu*r[i]
yv1<-sqrt(pmax(0,r[i]^2-(xv1-x[i])^2))+y[i]
yv2<--sqrt(pmax(0,r[i]^2-(xv2-x[i])^2))+y[i]
yv<-c(yv1,yv2)
xv<-c(xv1,xv2)
polygon(xv,yv,border=border,col=fill,lty=lty,lwd=lwd)
}
}
# Newton-Raphson algorithm for n multidimensional functions using numeric differentiation
NRnum<-function(init, fnlist, crit=6,...) {
par<-init
# sapply(par,function(x) cat(x," "))
# cat("\n")
value<-sapply(fnlist,function(f) f(par,...))
j<-1
while (round(sum(abs(value)),crit)!=0&j<100) {
# cat(j,"\n")
grad<-t(sapply(fnlist,function(f) attributes(numericDeriv(quote(f(par)),c("par")))$gradient))
deltax<-solve(grad,-value)
par<-par+deltax
# sapply(par,function(x) cat(x," "))
# cat("\n")
value<-sapply(fnlist,function(f) f(par,...))
j<-j+1
}
if (j<100) {
list(par=par,value=value)
} else {
list(par=NA,value=NA)
}
}
# An up-and-down algorithm
# Solves fn(x)=0 for x between l (lower bound) and u (upper bound).
updown<-function(l, u, fn, crit=6) {
fnu<-fn(u)
fnl<-fn(l)
if (fnl*fnu>0) return(NA)
value<-fn((u+l)/2)
while (round(value,crit)!=0&(u-l)>10^(-crit)) {
if (fnu*value>0) {
u<-(u+l)/2
fnu<-value
} else {
l<-(u+l)/2
fnl<-value
}
value<-fn((u+l)/2)
# cat(l,fnl,u,fnu,"\n")
}
nollakohta<-(l+u)/2
nollakohta
}
# Solve unidimensional equation using Newton-Rapson algorithm
# init=Initial guess,
# fn=function,
# gr=gradient
# crit= convergence criteria
# range= range of the possible values
NR<-function(init, fn, gr, crit=6, range=c(-Inf,Inf)) {
par<-init
value<-fn(par)
j<-1
while (round(value,crit)!=0&j<100000) {
grad<-gr(par)
a<-value-grad*par
par<--a/grad
par<-min(max(par,range[1]),range[2])
value<-fn(par)
cat(par,value,"\n")
j<-j+1
}
if (j<100000) {
list(par=par,value=value)
} else {
list(par=NA,value=NA)
}
}
HDnaslund<-function(d,a,b,bh=1.3) {
d^2/(a+b*d)^2+bh
}
startHDnaslund<-function(d,h,bh=1.3) {
start<-coef(lm(I(d/sqrt(h-bh))~d))
names(start)<-c("a","b")
start
}
# Michailoff 1943. H=a exp{-bD^{-1}}
# Lundqvist-Korf with c=1
HDmichailoff<-function(d,a,b,bh=1.3) {
a*exp(-b*d^(-1))+bh
}
startHDmichailoff<-function(d,h,bh=1.3) {
a<-1.3*max(h-bh)
mod<-lm(log((h-bh))~I(1/d))
start<-c(exp(coef(mod)[1]),-coef(mod)[2])
names(start)<-c("a","b")
start
}
# Curtis says it can be found in Prodan(book, 1965) or Strand (article, 1951).
HDcurtis<-function(d,a,b,bh=1.3) {
a*(d/(1+d))^b+bh
}
startHDcurtis<-function(d,h,bh=1.3) {
start<-coef(lm(log(h)~log(d/(1+d))))
start[1]<-exp(start[1])
names(start)<-c("a","b")
start
}
# Meyer: H=a(1-exp{-bD}) Meyer 1940
HDmeyer<-function(d,a,b,bh=1.3) {
a*(1-exp(-b*d))+bh
}
startHDmeyer<-function(d,h,bh=1.3) {
a<-1.3*max(h-bh)
b<--coef(lm(log(1-(h-bh)/a)~d-1))
start<-c(a,b)
names(start)<-c("a","b")
start
}
# Power function
# Stoffels 1950
HDpower<-function(d,a,b,bh=1.3) {
a*d^b+bh
}
startHDpower<-function(d,h,bh=1.3) {
start<-coef(lm(log(h-bh)~log(d))) # includes starting values estimates of a and b
start[1]<-exp(start[1])
names(start)<-c("a","b")
start
}
HDnaslund2<-function(d,a,b,bh=1.3) {
d^2/(a+exp(b)*d)^2+bh
}
startHDnaslund2<-function(d,h,bh=1.3) {
start<-coef(lm(I(d/sqrt(h-bh))~d))
start[2]<-log(start[2])
names(start)<-c("a","b")
start
}
HDnaslund3<-function (d, a, b, bh = 1.3) {
d^2/(exp(a) + b * d)^2 + bh
}
startHDnaslund3<-function (d, h, bh = 1.3) {
start <- coef(lm(I(d/sqrt(h - bh)) ~ d))
start[1]<-log(max(start[1],0.1))
names(start) <- c("a", "b")
start
}
HDnaslund4<-function (d, a, b, bh = 1.3) {
d^2/(exp(a) +exp( b) * d)^2 + bh
}
startHDnaslund4<-function (d, h, bh = 1.3) {
start <- coef(lm(I(d/sqrt(h - bh)) ~ d))
start[1]<-log(max(start[1],0.1))
start[2]<-log(max(start[2],0.1))
names(start) <- c("a", "b")
start
}
# Bates and Watts aD/(b+d) (Calama and Montero 2004)
HDmicment<-function(d,a,b,bh=1.3) {
a*d/(b+d)+bh
}
startHDmicment<-function(d,h,bh=1.3) {
tmp<-coef(lm(I(d/(h-bh))~d))
start<-c(1/tmp[2],tmp[1]/tmp[2])
names(start)<-c("a","b")
start
}
# Bates and Watts aD/(b+d) (Calama and Montero 2004)
HDmicment2<-function(d,a,b,bh=1.3) {
d/(a+b*d)+bh
}
startHDmicment2<-function(d,h,bh=1.3) {
start<-coef(lm(I(d/(h-bh))~d))
names(start)<-c("a","b")
start
}
HDwykoff<-function(d,a,b,bh=1.3) {
exp(a+b/(d+1))
}
startHDwykoff<-function(d,h,bh=1.3) {
mod<-lm(log((h-bh))~I(1/d))
start<-c(coef(mod)[1],coef(mod)[2])
names(start)<-c("a","b")
start
}
# 3- parameter models
HDprodan<-function(d,a,b,c,bh=1.3) {
d^2/(a+b*d+c*d^2)+bh
}
startHDprodan<-function(d,h,bh=1.3) {
start<-coef(lm(I(d^2/(h-bh))~d+I(d^2)))
names(start)<-c("a","b","c")
start
}
HDlogistic<-function(d,a,b,c,bh=1.3) {
a/(1+b*exp(-c*d))+bh
}
startHDlogistic<-function(d,h,bh=1.3) {
a<-1.3*max(h-bh)
start<-c(a,coef(lm(log((a-h+bh)/(h-bh))~d)))
start[2]<-exp(start[2])
start[3]<--start[3]
names(start)<-c("a","b","c")
start
}
HDrichards<-function(d,a,b,c,bh=1.3) {
a*(1-exp(-b*d))^c+bh
}
startHDrichards<-function(d,h,bh=1.3,b=0.04) {
start<-coef(lm(log(h-bh)~log(1-exp(-b*d))))
start<-c(exp(start[1]),b,start[2])
names(start)<-c("a","b","c")
start
}
HDweibull<-function(d,a,b,c,bh=1.3) {
a*(1-exp(-b*d^c))+bh
}
startHDweibull<-function(d,h,bh=1.3) {
a<-1.3*max(h-bh)
start<-c(a,coef(lm(log(-log(1-(h-bh)/a))~log(d))))
start[2]<-exp(start[2])
names(start)<-c("a","b","c")
start
}
HDgomperz<-function(d,a,b,c,bh=1.3) {
a*exp(-b*exp(-c*d))+bh
}
startHDgomperz<-function(d,h,bh=1.3) {
a<-1.3*max(h-bh)
start<-c(a,coef(lm(log(-log((h-bh)/a))~d)))
start[2]<-exp(start[2])
start[3]<--start[3]
names(start)<-c("a","b","c")
start
}
HDsibbesen<-function(d,a,b,c,bh=1.3) {
a*d^(b*d^(-c))+bh
}
startHDsibbesen<-function(d,h,bh=1.3,a=0.5) {
start<-c(a,coef(lm(log(log((h-bh)/a)/log(d))~log(d))))
start[2]<-exp(start[2])
start[3]<--start[3]
names(start)<-c("a","b","c")
start
}
# Korf: H=a exp{-bD^{-c}}
HDkorf<-function(d,a,b,c,bh=1.3) {
val<-a*exp(-b*d^(-c))+bh
# print(a,b,c)
# points(d,val,col="green")
val
}
startHDkorf<-function(d,h,bh=1.3) {
a<-1.3*max(h-bh)
y<-log(-log((h-bh)/a))
startmod<-lm(y[is.finite(y)]~log(d)[is.finite(y)],na.action=na.omit)
start<-c(a,exp(coef(startmod)[1]),-coef(startmod)[2])
names(start)<-c("a","b","c")
start
}
# Korf: H=a exp{-bD^{-c}}
HDratkowsky<-function(d,a,b,c,bh=1.3) {
val<-a*exp(-b/(d+c))+bh
# print(a,b,c)
# points(d,val,col="green")
val
}
startHDratkowsky<-function(d,h,bh=1.3,c=5) {
start<-c(coef(lm(log(h-bh)~I(1/(d+c)))),c)
start[1]<-exp(start[1])
start[2]<--start[2]
names(start)<-c("a","b","c")
start
}
## Korf: H=a exp{-bD^{-c}}
#HDhossfeldIV<-function(d,a,b,c,bh=1.3) {
# val<-a/(1+1/b*d^(-c))+bh
# val
# }
#
#startHDhossfeldIV<-function(d,h,bh=1.3,c=5) {
# a<-1.3*max(h-bh)
# start<-c(a,coef(lm(log(a/(h-bh)-1)~log(d))))
# start[2]<-1/exp(start[2])
# start[3]<--start[3]
# names(start)<-c("a","b","c")
# start
# }
HDhossfeldIV<-function(d,a,b,c,bh=1.3) {
val<-a/(1+1/(b*d^c))+bh
val
}
startHDhossfeldIV<-function(d,h,bh=1.3,c=5) {
a<-1.3*max(h-bh)
start<-c(a,coef(lm(log(1/(a/h-1))~log(d))))
start[2]<-exp(start[2])
names(start)<-c("a","b","c")
start
}
HDmehtatalo<-function(d,a,b,c,bh=1.3) {
val<-d^c/(a+exp(b)*d)^c+bh
points(d,val,col="red")
val
}
startHDmehtatalo<-function(d,h,bh=1.3) {
start<-c(coef(lm(I(d/sqrt(h-bh))~d)),2)
start[2]<-log(start[2])
names(start)<-c("a","b","c")
start
}
#startHDkorf<-function(d,h,bh=1.3) {
# a<-1.3*max(h-bh)
# mod<-lm(log((h-bh)/a)~I(1/d)-1)
# start<-c(a,-coef(mod))
# names(start)<-c("a","b")
# start
# }
#startHDkorf3<-function(d,h,bh=1.3) {
# a<-max(h-bh)
# b<--coef(lm(log((h-bh)/a)~I(1/d)-1))
# start<-c(a,b,1)
# names(start)<-c("a","b","c")
# start
# }
# Lineaariset
# h=a+b*log(D)
# a*D/(D+1)+b*D
# 2. asteen polynomi
# spliniregressio
# Fits one of the commonly used models to H-D data and returns an
# object of class hdmod.
fithd<-function(d, h, plot=c(), modelName="naslund",
nranp=2, random=NA, varf=0, na.omit=TRUE,
start=NA, bh=1.3, control = list(), SubModels=NA, vfstart=0) {
if (!is.na(list(random))) nranp<-1
varf<-as.numeric(varf)
if (na.omit) {
sel<-is.na(h*d)
d<-d[!sel]
plot<-plot[!sel]
h<-h[!sel]
}
# if (!is.na(SubModels)[1]) {
dmean<-tapply(d,plot,mean)
dmean<-dmean[match(plot,names(dmean))]
dsd<-tapply(d,plot,sd)
dsd<-dsd[match(plot,names(dsd))]
dstd<-(d-dmean)/dsd
# }
if (varf==1) w<-d else if (varf==2) w<-pmax(1,dstd+3)
# Linear model fitting
if (inherits(modelName,"formula")) {
if (!is.na(SubModels)[1]) warning("Argument 'SubModels' not implemented for linear model fit")
modEq<-modelName
if (!is.na(random)) {
# If a formula is provided in ranp, use it.
ranEq=as.formula(paste(deparse(random),"|plot"))
} else if (nranp>0) {
# Otherwise take the given number of trems of model from the beginning.
# Constant is always assumed to be the first if model just has it.
ranTerms<-paste("+",attributes(terms(modelName))$term.labels,sep="")
if (attributes(terms(modelName))$intercept) {
ranTerms<-c("1",ranTerms)
ranp<-1:nranp
} else {
ranTerms<-c("-1",ranTerms)
ranp<-1:(nranp+1)
}
ranEq<-"~"
for (j in ranp) {
ranEq<-paste(ranEq,ranTerms[j],sep="")
}
ranEq<-as.formula(paste(ranEq,"|plot",sep=""))
}
if (!nranp) {
if (varf) {
mod<-gls(modEq,
weights=varPower(vfstart,~w),
control=control)
} else {
mod<-gls(modEq,
control=control)
}
} else {
if (varf) {
mod<-lme(fixed=modEq,
random=ranEq,
weights=varPower(vfstart,~w),
control=control)
} else {
mod<-lme(fixed=modEq,
random=ranEq,
control=control)
}
}
class(mod)<-c("hdmod",class(mod))
attributes(mod)$model<-modelName
attributes(mod)$d<-d
attributes(mod)$h<-h
attributes(mod)$plot<-plot
mod$call[[2]]<-modEq
if (nranp) mod$call[[3]]<-ranEq
} else {
# Nonlinear model fitting
if (is.na(start[1])) {
start<-eval(call(paste("startHD",modelName,sep=""),d=d,h=h,bh=bh))
if (!is.na(SubModels)[1]) {
start2<-c()
for (i in 1:length(SubModels)) {
start2<-c(start2,start[i])
start2<-c(start2,0[as.character(SubModels[i])!="1"])
}
start<-start2
}
}
# # The # of random parameters should not exceed the # of parameters
# ranp<-min(ranp,length(formals(paste("HD",model,sep="")))-2)
if (!is.na(list(random))) {
ranEq<-random
} else if (nranp) {
ranEq<-lapply(paste(names(formals(paste("HD",modelName,sep="")))[1+1:nranp],"~1",sep=""),as.formula)
}
# if (nranp==1) {
# ranEq<-list(a~1)
# } else if (nranp==2) {
# ranEq<-list(a~1,b~1)
# } else if (nranp==3) {
# ranEq<-list(a~1,b~1,c~1)
# }
if (length(formals(paste("HD",modelName,sep="")))==5) {
# if (modelName=="korf3") {
modEq<-as.formula(paste("h~HD",modelName,"(d,a,b,c,bh=",bh,")",sep=""))
if (is.na(SubModels)[1]) {
fixEq<-list(as.formula("a~1"),as.formula("b~1"),as.formula("c~1"))
} else if (length(SubModels==3)) {
fixEq<-lapply(paste(c("a~","b~","c~"),SubModels),as.formula)
# fixEq<-list(as.formula("a~1"),as.formula("b~1"),as.formula("c~1"))
} else {
stop("length of SubModels is not 3")
}
} else {
modEq<-as.formula(paste("h~HD",modelName,"(d,a,b,bh=",bh,")",sep=""))
if (is.na(SubModels[1])) {
fixEq<-list(as.formula("a~1"),as.formula("b~1"))
} else if (length(SubModels==2)) {
fixEq<-lapply(paste(c("a~","b~"),SubModels),as.formula)
# fixEq<-list(as.formula("a~1"),as.formula("b~1"),as.formula("c~1"))
} else {
stop("length of SubModels is not 2")
}
}
if (!nranp) {
if (!varf) {
mod<-gnls(modEq,
params=fixEq,
start=start,
control=control
)
} else {
mod<-gnls(modEq,
params=fixEq,
start=start,
weights=varPower(vfstart,~w),
control=control
)
}
} else {
if (!varf) {
mod<-nlme(modEq,
fixed=fixEq,
random=ranEq,
groups=~plot,
start=start,
control=control)
} else {
mod<-nlme(modEq,
fixed=fixEq,
random=ranEq,
groups=~plot,
start=start,
weights=varPower(vfstart,~w),
control=control)
}
}
class(mod)<-c("hdmod",class(mod))
attributes(mod)$model<-modelName
attributes(mod)$d<-d
attributes(mod)$h<-h
attributes(mod)$plot<-plot
mod$call[[2]]<-modEq
mod$call[[3]]<-fixEq
if (nranp) mod$call[[4]]<-ranEq
}
mod
}
# Plot a height-diameter model
plot.hdmod<-function(x,col.point="blue",highlight="red",standd=TRUE,cex=1,corD=FALSE,ask=TRUE,...) {
model<-x
res<-resid(model,type="p")
d<-attributes(model)$d
plot<-attributes(model)$plot
plots<-unique(plot)
dstd<-d
# compute diameter that is standardized for each plot
if (standd) {
for (i in 1:length(plots)) {
# thisplot will include observations only from i:th plot
dthis<-d[plot==plots[i]]
dstd[plot==plots[i]]<-(dthis-mean(dthis))/sd(dthis)
}
plot(dstd,
res,
col=col.point,cex=cex,
main=paste(attributes(model)$model,"s.e.=",round(sd(resid(model)),3)),
xlab="Relative diameter",
ylab="Standardized residual")
} else {
plot(dstd,
res,
col=col.point,cex=cex,
main=paste(attributes(model)$model,"s.e.=",round(sd(resid(model)),3)),
xlab="Diameter, cm",
ylab="Standardized residual")
}
mywhiskers(dstd,res,add=TRUE,highlight="black",se=FALSE)
mywhiskers(dstd,res,add=TRUE,highlight=highlight,lwd=3)
abline(0,0)
if (inherits(model,"nlme")&corD) {
re<-ranef(model)
re$plot<-as.numeric(rownames(ranef(model)))
dmean<-tapply(d,plot,mean)
re$dmean<-dmean[match(re$plot,names(dmean))]
nranp<-dim(ranef(model))[2]
for (i in 1:nranp) {
devAskNewPage(ask)
plot(log(re$d),re[,i],xlab="log(dmean)",ylab=paste(i,"th random effect",sep=""))
abline(lm(re[,i]~log(re$d)))
}
}
}
# Plot a height-diameter model
qqplotHD<-function(model,startnew=TRUE) {
res<-resid(model,type="p")
if (sum(class(model)=="lme")) {
rem<-ranef(model)
nplots<-dim(rem)[2]+1
if (startnew) {
par(mfcol=c(2,ceiling(nplots/2)))
}
qqnorm(res,main="N-QQ-plot of standardized residuals")
qqline(res)
sapply(1:dim(rem)[2],function(x) {qqnorm(rem[,x],main=paste("QQ-plot of",colnames(rem)[x])); qqline(rem[,x])})
} else {
qqnorm(res,main="N-QQ-plot of standardized residuals")
qqline(res)
}
}
ImputeHeights<-function(d, h, plot=c(),modelName="naslund",nranp=2, varf=TRUE,
addResidual=FALSE, makeplot=TRUE, level=1,
start=NA,bh=1.3,control=list(),random=NA) {
data<-data.frame(d=d,h=h,plot=plot)
data1<-data[!is.na(data$h*data$d),]
mod<-fithd(data1$d,data1$h,data1$plot,
modelName=modelName,nranp=nranp,random=random,
varf=varf,start=start,bh=bh,control=control)
hplots<-unique(data$plot[!is.na(data$h*data$d)])
type<-rep(2,dim(data)[1])
for (i in 1:length(hplots)) {
type[data$plot==hplots[i]]<-1
}
nohplots<-unique(data$plot[type==2])
if (makeplot) plot(mod)
# To initialize, compute the predictions using the fixed part only.
# This will remain as the final predioction
# a) if nranp=0 (no randiom effects) and addResidual==0
# b) if level=0 (fixed part prediction is requested) and addResidual==0.
predall<-predict(mod,newdata=data, level=0)
# If a higher level prediction is requested, then the realized random effects are used.
# If Some plots did not have measured heights for the preiction of random effects,
# Then random effects are sampled from among the predictions. (Not using MVN(0,D)!)
if (level>0 && nranp>0) {
predall[type==1]<-predict(mod,newdata=data[type==1,])
D<-as.matrix(mod$modelStruct$reStruct[[1]])*mod$sigma^2
if (length(nohplots)) {
for (i in 1:length(nohplots)) {
# thiscoef<-fixef(mod)+mvrnorm(1,mu=rep(0,dim(D)[1]),Sigma=D)
thiscoef<-coef(mod)[sample(1:mod$dims$ngrps[1],1),]
if (length(formals(paste("HD",modelName,sep="")))==5) {
predall[data$plot==nohplots[i]]<-
eval(call(paste("HD",modelName,sep=""),
d=data$d[data$plot==nohplots[i]],a=thiscoef[[1]],b=thiscoef[[2]],c=thiscoef[[3]],bh=bh))
} else {
predall[data$plot==nohplots[i]]<-
eval(call(paste("HD",modelName,sep=""),
d=data$d[data$plot==nohplots[i]],a=thiscoef[[1]],b=thiscoef[[2]],bh=bh))
}
# predall[data$plot==nohplots[i]] <- HDnaslund(d=data$d[data$plot==nohplots[i]],a=thiscoef[[1]],b=thiscoef[[2]])
}
}
}
# If level = 0 and AddResidual=True, then random effects are randomly assigned for each plot
# before adding the tree level residual
if (level==0 & nranp>0 & addResidual) {
plots<-unique(data$plot)
for (i in 1:length(plots)) {
# thiscoef<-fixef(mod)+mvrnorm(1,mu=rep(0,dim(D)[1]),Sigma=D)
thiscoef<-coef(mod)[sample(1:mod$dims$ngrps[1],1),]
if (length(formals(paste("HD",modelName,sep="")))==5) {
predall[data$plot==plots[i]]<-
eval(call(paste("HD",modelName,sep=""),
d=data$d[data$plot==plots[i]],a=thiscoef[[1]],b=thiscoef[[2]],c=thiscoef[[3]],bh=bh))
} else {
predall[data$plot==plots[i]]<-
eval(call(paste("HD",modelName,sep=""),
d=data$d[data$plot==plots[i]],a=thiscoef[[1]],b=thiscoef[[2]],bh=bh))
}
# predall[data$plot==plots[i]] <- HDnaslund(d=data$d[data$plot==plots[i]],a=thiscoef[[1]],b=thiscoef[[2]])
}
}
# If addResidual==True, a random residual is added from a normal distribution
# using the estimated variance function.
if (addResidual) {
if (varf) {
ressd<-mod$sigma*data$d^mod$modelStruct$varStruct
} else if (!varf) {
ressd<-mod$sigma
}
predall<-predall+rnorm(length(data$d),0,ressd)
}
type[!is.na(h)]<-0
h[is.na(h)]<-predall[is.na(h)]
list(h=h,imputed=is.na(data$h),model=mod,predType=type,hpred=predall)
}
# Predicts tree volumes (liters) using the functions of Laasasenaho.
# pl=species (1= pine)
# malli mean model. Value of 1 uses only d; balue of 2 uses d and h.
predvol<-function(species,d,h=0,model=1) {
pl<-species
malli<-model
pl2<-pl
pl2[pl2>4|pl2<1]<-5
if (malli==1) {
pars<-matrix(c(-5.39417,3.48060,0.039884,-5.39934,3.46468,0.0273199,
-5.41948,3.57630,0.0395855,-5.41948,3.57630,0.0395855,
rep(NA,3)),ncol=3,byrow=TRUE)
exp(pars[pl2,1]+pars[pl2,2]*log(2+1.25*d)-pars[pl2,3]*d)
} else if (malli==2) {
pars<-matrix(c(0.036089,2.01395,0.99676,2.07025,-1.07209,
0.022927,1.91505,0.99146,2.82541,-1.53547,
0.011197,2.10253,0.98600,3.98519,-2.65900,
0.011197,2.10253,0.98600,3.98519,-2.65900,
rep(NA,5)),ncol=5,byrow=TRUE)
pars[pl2,1]*d^pars[pl2,2]*pars[pl2,3]^d*h^pars[pl2,4]*(h-1.3)^pars[pl2,5]
}
}
# The stand table example
# The multinomial distribution
# Linear interpolation.
# Assuming given values of X and Y, the function
# returns the interpolated values at points x and the first derivative.
# In the case of percentile-based distribution, these are the distribution function and density.
interpolate<-function(x,X,Y) {
n<-length(X) # the number of points, the number of intervals is n-1
iv<-1:n
# define the interval at which each of x values is located
# i is a vector including values from 1 to n-1
i<-sapply(x,function(x) max(c(1,iv[x>=X])))
i[x>=X[n]]<-n-1 # extrapolation at the upper end
i[x<X[1]]<-1 # extrapolation at the lower end
# compute the n-1 long vector of slope coefficients
slope<-(Y[iv+1]-Y[iv])/(X[iv+1]-X[iv])
slope<-slope[-length(slope)]
# compute the value for each value in x by picking the correct constant and slope
y<-Y[i]+slope[i]*(x-X[i])
list(y=y,dy=slope[i],const=Y[-length(slope)]-slope*X[-length(slope)],slope=slope)
}
ddcomp<-function(d,density="dweibull",power=0,
limits=seq(0,100),limitsd=limits,plot=FALSE,...) {
n<-length(limits)-1
if (is.character(density)) {
f <- match.fun(density)
pdf <- function(x) f(x, ...)
cdf<-function(x) sapply(x,function(xi) integrate(pdf,0,xi)$value)
if (cdf(limits[1])>1e-6|cdf(limits[n+1])<(1-1e-6)) {
warning(paste("The provided limits seems not cover
the whole range of estimated pdf:", cdf(limits[1]),cdf(limits[n+1]))) }
EX<-integrate(function(x) x*pdf(x),0,Inf)$value
EX2<-integrate(function(x) x^2*pdf(x),0,Inf)$value
} else if (is.numeric(density)) {
nd<-length(limitsd)-1
meansd<-(limitsd[-1]+limitsd[-(nd+1)])/2
if (length(density)!=nd) {
stop("Incompatible number of diameter classes in the estimated density") }
if (abs(sum(density)-1)>1e-8) {
stop("Estimated densities should sum up to unity") }
f<-density
EX<-sum(f*meansd)
EX2<-sum(f*meansd^2)
cdf<-function(x) {
F<-c(0,sapply(1:length(f),function(i) sum(f[1:i])))
interpolate(x,limitsd,F)$y
}
} else {
stop("Argument density should be either character or numeric.")
}
varX<-EX2-EX^2
meanD<-mean(d)
sdD<-sd(d)
mudif<-meanD-EX
vardif<-sdD^2-varX
sddif<-sdD-sqrt(varX)
cdf2<-function(y) cdf(sqrt(varX)/sdD*(y-meanD)+EX)
if (plot) {
plot(ecdf(d))
yapu<-seq(min(d),max(d),length=100)
lines(yapu,cdf2(yapu))
lines(yapu,cdf(yapu),lty="dashed")
legend("topleft",legend=c("scaled","original"),lty=c("solid","dashed"))
}
means<-(limits[-1]+limits[-(n+1)])/2
if (min(d)<limits[1]|max(d)>limits[n+1]) {
stop("The provided limits do not cover all observed diameters.") }
dp<-d^power
sumdp<-sum(dp)
ftrue<-sapply(1:n,function(i) sum(dp[d>=limits[i]&d<limits[i+1]]))/sumdp
Fpred<-cdf(limits)
fpred<-Fpred[-1]-Fpred[-n]
fpred<-fpred*means^power/sum(fpred*means^power)
Fpred2<-cdf2(limits)
fpred2<-Fpred2[-1]-Fpred2[-n]
fpred2<-fpred2*means^power/sum(fpred2*means^power)
list(mudif=mudif,vardif=vardif,
sddif=sddif,ei1=sum(abs(fpred-ftrue)),
ei2=sum(abs(fpred2-ftrue)))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.