1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
|
cmvglm <- function(formula, model, ndim,data=NULL, fact=NULL){
toBuildFormula<-function(Xnames,sepp="+"){
lng<-length(Xnames)
rhs<-NULL
if (lng!=0){
if(lng==1){
rhs=Xnames
}else{
for (j in 1:(lng-1)){
rhs<-paste(rhs,as.name(Xnames[[j]]))
rhs<-paste(rhs,sepp)
}
rhs<-paste(rhs,Xnames[[lng]])
}
}
return (rhs)
}
tt<-terms(formula)
attr(tt,"systEqns")<-names(formula)
p<-make.parameters(tt,shape="matrix")
vars<-rownames(p)
cm<-vector("list", length(vars))
names(cm)<-vars
for(i in 1:length(cm))
cm[[i]]<-diag(1, ndim)
constrain<-attr(tt,"constraints")
if(!is.logical(constrain)){
tmp <- sort(colnames(constrain))
for (i in 1:length(tmp)) {
ci<-constrain[,i]
if (is.null(na.omit(ci)) || length(unique(na.omit(ci)))!=1)
stop("invalid input for constrain")
minj <- match(FALSE, is.na(ci))
whatvar <- pmatch(unique(na.omit(ci)), names(cm))
for (j in 1:3)
if (!is.na(ci[j])) {
cm[[whatvar]][j,j]<-0
cm[[whatvar]][j,minj]<-1
}
}
}
for(i in rownames(p)){
for(j in 1:ncol(p)){
if(is.na(p[i,j]))
cm[[i]][j,j]<-0
}
}
# if(!is.null(constant))
# for(i in 1:length(constant))
# for(j in 1:length(cm))
# if(names(cm)[j]!="(Intercept)")
# cm[[j]][constant[i],]<-matrix(0, ncol=ncol(cm[[j]]))
for(i in 1:length(cm))
cm[[i]]<-as.matrix(cm[[i]][,apply(cm[[i]], 2, sum)!=0])
rhs<-toBuildFormula(attr(tt,"indVars"))
if(!(is.null(rhs)))
rhs<-(paste("~",rhs))
else
rhs<-"~1"
Ynames<-unlist(attr(tt,"depVars"))
if(!is.null(fact))
lhs<-fact
else{
if(length(Ynames)>1){
lhs<-toBuildFormula(Ynames,",")
if (!(is.null(lhs))){
lhs<-paste("cbind(",lhs)
lhs<-paste(lhs,")")
}
}else{
lhs=Ynames
}
}
formula<-as.formula(paste(lhs,rhs))
list("formula"=formula, "constraints"=cm)
}
|
