title: When using glmnet how to report p-value significance to claim significance of predictors?

text: I have a large set of predictors (more than 43,000) for predicting a dependent variable which can take 2 values (0 or 1). The number of observations is more than 45,000. Most of the predictors are unigrams, bigrams and trigrams of words, so there is high degree of collinearity among them. There is a lot of sparsity in my dataset as well. I am using the logistic regression from the glmnet package, which works for the kind of dataset I have. My problem is how can I report p-value significance of the predictors. I do get the beta coefficient, but is there a way to claim that the beta coefficients are statistically significant? Here is my code:               library('glmnet')     data <- read.csv('datafile.csv', header=T)     mat = as.matrix(data)     X = mat[,1:ncol(mat)-1]      y = mat[,ncol(mat)]     fit <- cv.glmnet(X,y, family="binomial")      Another question is: I am using the default alpha=1, lasso penalty which causes the additional problem that if two predictors are collinear the lasso will pick one of them at random and assign zero beta weight to the other. I also tried with ridge penalty (alpha=0) which assigns similar coefficients to highly correlated variables rather than selecting one of them. However, the model with lasso penalty gives me a much lower deviance than the one with ridge penalty. Is there any other way that I can report both predictors which are highly collinear?


title: Intuition behind tensor product interactions in GAMs (MGCV package in R)

text: Generalized additive models are those where $$ y = \alpha + f_1(x_1) + f_2(x_2) + e_i$$ for example. the functions are smooth, and to be estimated. Usually by penalized splines. MGCV is a package in R that does so, and the author (Simon Wood) writes a book about his package with R examples. Ruppert et al (2003) write a far more accessible book about simpler versions of the same thing. My question is about interactions within these sorts of models. What if I want to do something like the following: $$ y = \alpha + f_1(x_1) + f_2(x_2) + f_3(x_1\times x_2) + e_i$$ if we were in OLS land (where the f is just a beta), I'd have no problem with interpreting $\hat{f}_3$. If we estimate via penalized splines, I also have no problem with interpretation in the additive context. But the MGCV package in GAM has these things called "tensor product smooths". I google "tensor product" and my eyes immediately glaze over trying to read the explanations that I find. Either I'm not smart enough or the math isn't explained very well, or both. Instead of coding               normal = gam(y~s(x1)+s(x2)+s(x1*x2))      a tensor product would do the same (?) thing by               what = gam(y~te(x1,x2))      when I do               plot(what)      or               vis.gam(what)      I get some really cool output. But I have no idea what is going on inside the black box that is te(), nor how to interpret the aforementioned cool output. Just the other night I had a nightmare that I was giving a seminar. I showed everyone a cool graph, they asked me what it meant, and I didn't know. Then I discovered that I had no clothes on. Could anyone help both me, and posterity, by giving a bit of mechanics and intuition on what is going on underneath the hood here? Ideally by saying a bit about the difference between the normal additive interaction case and the tensor case? Bonus points for saying everything in simple English before moving on to the math.


title: Issues when using neural network

text: I am having an issue with using neural networks. I started with something simple. I just used nntool with one hidden layer (with one neuron) with linear activation function. For the output also, I used the linear activation function. I just fed my Xs and Ys to the neural network tool and got the results on the testing set. I compared that with normal ridge function in MATLAB. I could see that neural network one performed much worse than ridge function. The first reason is that there are lots of negative values in the predictions, when my target is only positive. Ridge regression gave about 800 -ve values while nn gave around 5000 -ve values which totally ruined the accuracy of nntool. Since nntool can be used to perform linear regression, why nntool is not performing as well as ridge regression? What about the negative values what's the reason behind it and how to enforce positive values?


title: How do you visualize binary outcomes versus a continuous predictor?

text: I have some data I need to visualize and am not sure how best to do so. I have some set of base items $Q = \\{ q_1, \cdots, q_n \\}$ with respective frequencies $F = \\{f_1, \cdots, f_n \\}$ and outcomes $O \in \\{0,1\\}^n$. Now I need to plot how well my method "finds" (i.e., a 1-outcome) the low frequency items. I initially just had an x-axis of frequency and a y axis of 0-1 with point-plots, but it looked horrible (especially when comparing data from two methods). That is, each item $q \in Q$ is has an outcome (0/1) and is ordered by its frequency. > Here is an example with a single method's results: ![enter image description here](http://i.stack.imgur.com/ELfGN.png) My next idea was to divide the data into intervals and compute a local sensitivity over the intervals, but the problem with that idea is the frequency distribution is not necessarily uniform. So how should I best pick the intervals? > Does anyone know of a better/more useful way to visualize these sort of data > to portray the effectiveness of finding rare (i.e., very low-frequency) > items? EDIT: To be more concrete, I am showcasing the ability of some method to reconstruct biological sequences of a certain population. For validation using simulated data, I need to show the ability to reconstruct variants regardless of its abundance (frequency). So in this case I am visualizing the missed and found items, ordered by their frequency. This plot will not include reconstructed variants that are not in $Q$.


title: Mean regression factor scores and attributes cross tabulation yields. Why are all expected signs reversed?

text: I ran a factor analysis on 20 reasons for purchasing 4 different goods. These are ranked on a Likert scale from 1 to 5 with 5 being "extremely important", 4 "important", etc. I extracted 4 factors which I am able to label based on the group of reasons that loaded well into each of these factors. I saved these factors as regression factor scores in SPSS. Then I did a cross-tabulation with each of the goods as rows and regression factor scores as columns. I used the mean for the regression factor scores. I expected the sign of the factor score and each of the goods to be positive if there was a relationship with that factor score, or negative otherwise. However, the signs are exactly reversed. So for those goods that I thought would have the highest positive mean factor score, I am actually getting the highest negative, etc. Am I completely wrong in the approach I am taking to cross-tabulate this, or am I interpreting these wrong?


title: Re-parametrize posterior distribution

text: How can I write the unnormalized posterior $ f(p_1, p_2 | Y) = (z_1-1)*log(p_1) + (n_1-z_1-1)*log(1-p_1) + (z_2-1)*log(p_2) + (n_2-z_2-1)*log(1-p_2) $ in terms of the log odds-ratio $\alpha$ and the log odds-product $\eta$? $\alpha = log \left(\frac{p_2/(1-p_2)}{p_1/(1-p_1}\right) $ and $\eta = log \left(\frac{p_2}{1-p_2} \times \frac{p_1}{1-p_1}\right)$? where z_1, n_1, z_2, and n_2 are constants.


title: Ridge Regression

text: I have pretty much the same question as in this post: I have a situation where I want to regress 50 daily observations of interest rate swap yields on four factors. Actually it's the Nelson Siegel Svensson model http://en.wikipedia.org/wiki/Fixed-income_attribution see the section modelling the yield curve At the moment I have an optimisation process that finds the optimal set of parameters for the factors and the beta (weighting) coefficients. As can be seen if the factors tau1 and tau2 are the same the model matrix will be singular or if they are close it will be ill conditioned. The optimisation process works well for test data sets where the parameters used to produce the yeld curves are chosen to ensure there is no ill- conditioning. However when I apply it to real data best fit parameters seem to jump around alot from day to day even when the curve is hardly moving. I think this may well be because the optimal fits to the real data are based on parameter sets that produce ill-conditioned matrices meaning the mapping of the Beta weights is incorrect. I could just incorporate a penalty into the optimisation process that prevents badly conditioned matrices from being considered. However, there is nothing in principal (that I can think of) that means that the optimal mapping of these 50 observations onto the factors will result in a matrix that is not have at lest some ill conditioning. So I really need a way to handle these ill conditioned cases and to be able to compare the results from matrices where the not ill conditioned. I have reserched the topic and the main methods to deal with such situations is to use Ridge Regression, Principal Components Regressiona and Partial Least Squares. With Ridge Regression generally being seen as the best. (Do you agree). Even then there seems to be different kinds of Ridge Regression, unbiased ridge regression, mixture methods. eg http://www.utgjiu.ro/math/sma/v04/p08.pdf I don't want to remove any of the 4 factors ao anything like that indeed its imperative that all 4 factors remain. I just want to get the best fit to those factors. I just wondered what the best method would be for my problem? Ideally I also need the standard errors for the cofficients. Can anyone receommend code ideally in MATLAB or R that can implement the best method? Is cross validation the best way to compare the competing models? It is OK To compare ill conditioned and non ill conditioned models in this way?


title: Conditional expectation conditional on exponential random variable

text: Given $X_1$, $X_2$ are independent `EXP(1)` random variables, how do I compute the following expectations:   * $E[X_1 | X_1 + X_2]$   * $P(X_1 > 3 | X_1 + X_2)$   * $E[X_1 | min(X_1, t)]$


title: Penalized spline confidence intervals based on cluster-sandwich VCV

text: This is my first post here, but I've benefited a lot from this forum's results popping up in google search results. I've been teaching myself semi-parametric regression using penalized splines. Basically, $$ \hat{y} = X(X'X + \lambda D)^{-1}X'y $$ where X is full of knots in addition to my original variables, $\lambda$ is a penalty term applied to the knots. (See e.g. Wood 2006 or Ruppert, Wand & Carrol 2003 if you want to learn this stuff). My data has cluster-wise autocorrelation, and probably heteroskedasticity as well. My training is in econometrics, so naturally I want to represent the variance-covariance matrix as $$ X(X'X + \lambda D)^{-1}\left(\displaystyle\sum_j X'_j \epsilon_j \epsilon_j'X\right)(X'X + \lambda D)^{-1}X' $$ The above assumes independence between clusters but not within them, and doesn't assume homoskedasticity. This should give me correct se's. My problem is that I don't just care about se's. With semi-parametrics, I need confidence intervals on plots. **EDIT** some time later:   I've developed a quasi-brute-force way of doing this by hacking `mgcv`. The general approach is to calculate my own parameter covariance matrix and stuff it into `mgcv`'s. It is then used by other functions to calculate CI's and give CI's on plots This is very much a work in progress and I'm not confident in it yet, but I'm putting it up in hopes it might be a first step towards something that works, invite feedback, and perhaps eventually be useful to others. Here it is as it stands: **Another EDIT** : The function now tries to get the full penalty matrix by inverting the parameter covariance matrix and the subtracting off $X'X$. **Weirdly, the parameter covariance matrix supplied by`mgcv` is often not invertible. I have no idea why this would be the case,** after removing auto- dropped redundant dummies. **Anybody have any ideas here?** I'm using generalized inverses as a hack for the moment, and I'm not sure exactly what sort of distortion that will add to my results.               clusterGAM = function(m,clustvar){         require(MASS)         mm = predict(m,type='lpmatrix') #Get model X matrix         mm = mm[,colnames(mm) %in% names(coefficients(m))[rowSums(m$Vp)!=0]] #Cut out coefficients from model matrix that get auto-dropped due to collinearity          meat = matrix(rep(0,ncol(mm)*ncol(mm)),ncol(mm)) #Pre-fill a meat matrix          for (i in 1:length(unique(clustvar))){#This sums over the clusters           X = mm[clustvar == unique(clustvar)[i],]           e = residuals(m)[clustvar == unique(clustvar)[i]]           meat = meat+t(X) %*% e %*% t(e) %*% X           }          M <- length(unique(clustvar))          N <- length(clustvar)          K <- sum(m$edf) #use EDF instead of K         dfc <- (M/(M-1))*((N-1)/(N-K))#Degrees of freedom correction         Vp = vcov(m,dispersion=1)         Vp = Vp[rowSums(m$Vp)!=0,rowSums(m$Vp)!=0]              S = tryCatch(solve(Vp) - t(mm)%*%mm, error=function(e) e, finally=NA)         if(inherits(S, "error")) {             S = ginv(Vp) - t(mm)%*%mm             print("Warning:  Generalized inverse used to extract the penalty matrix")             }            vcv = tryCatch(dfc*solve(t(mm)%*%mm+S)%*%meat %*% solve(t(mm)%*%mm+S), error=function(e) e, finally=NA)         if(inherits(vcv, "error")) {             vcv = dfc*ginv(t(mm)%*%mm+S)%*%meat %*% ginv(t(mm)%*%mm+S)             print("Warning:  Generalized inverse used to calculate the Vp matrix")             } #Get the vcv         if (length(which(rowSums(m$Vp)==0)!=0)){ #If any zero rows got knocked out, put them back in           zeros = which(rowSums(m$Vp)==0)           for (i in 1:length(zeros)){            vcv = rbind(vcv[0:(zeros[i]-1),],0,vcv[(zeros[i]):nrow(vcv),])            vcv = cbind(vcv[,0:(zeros[i]-1)],0,vcv[,(zeros[i]):ncol(vcv)])            }           }          m$Vp = vcv #stuff the new vcv into the place where the plots and the summary gets made from         m         }      *edited to be able to handle tensors Some data to run it with               load(url('http://andrewcd.berkeley.edu/sdat'))     head(sdat)     N = nrow(sdat)     #set.seed(417)     sdat$x = with(sdat,log(y+1+runif(N)*T)+rexp(N)) #Making a fake smooth covariate             plot(sdat$x,sdat$y)             m = gam(y~ID+G:t+T+s(x)-1,data=sdat)#LSDV panel model with time x group FEs and one smooth function             par(mfrow=c(1,2))             plot(m,ylim=c(-.75,2))             summary(m)             m = clusterGAM(m,sdat$G)     plot(m,ylim=c(-.75,2)) #the confidence intervals shrink a bit.  This seems counterintuitive.       summary(m) #SE on T gets a little bit bigger, as it probably should if you account for repeated observations          #Now try it with a tensor     sdat$z = with(sdat,y+x+rnorm(N))         m = gam(y~ID+G:t+T+te(x,z)-1,data=sdat)         summary(m)         vis.gam(m,view=c('x','z'),se=2,theta=150,phi=20)         m = clusterGAM(m,sdat$G)     vis.gam(m,view=c('x','z'),se=2,theta=150,phi=20)      It could be broken either by my imperfect understanding of the math involved, or the internal workings of `mgcv`. As it stands, the standard errors on parametric dummies are inflated in ways that seem reasonable, but confidence intervals on smooth terms seem like they get shrunk, which seems wrong. Any comments towards improving this would be greatly appreciated. As a final note: I'm not a stats PhD. I'm a grad student doing applied work, and probably spending far more time than I should letting myself get bothered by functional form assumptions.


title: Possible subgroup differences in regression

text: I am running an OLS regression and it turned out that the effect of my main predictor X on the outcome Y depends on one of my control variables (age). The central focus of the paper is describing the multiple mediators that link X to Y, and splitting up the sample reduces my power drastically. What are my options? I have probed the interaction and calculated the regions of significance, splitting my sample based on these regions of significance eliminates most of my significant effects. Sigh.. Is it OK to simply include the interaction in all analyses: i.e, Does including the interaction adequately control for the subgroup differences?


title: Inverse prediction from binomial glm - how to obtain confidence intervals?

text: I am trying to use a glm with a binomial distribution (logit link) to analyse data from a dose response curve for the lethal effects of different bacterial strains. I now want to obtain estimates of the lethal dose 50 (LD50). My model is simply of the form:               survival ~ a + b*dose      Thus since I have a logit link, obtaining the estimates for the LD50 (i.e. p=0.5) simply means dividing -a/b. However, I am also wanting to know the confidence intervals for these estimates and that's where I am unsure about what to do. If anyone would have a suggestion for this, that would be great.


title: Convergence of K-means

text: I have a clustering algorithm which works iteratively like K-means, but there are some constraints on cluster sizes with lower and upper thresholds. Do you know any convergence proofs of K-means in general? Would they be applicable in the case of this alternative algorithm?


title: extrasolar planet puzzler

text: Based on the Kepler Space Observatory's findings, astronomer Seth Shostak estimated in 2011 that "within a thousand light-years of Earth," there are "at least 30,000" habitable planets (from this Wikipedia article). Using Dr. Shostak's estimate, within which radius R from Earth (in light years) can we be 95% certain there is at least one habitable planet?


title: Problem on Normal Distribution

text: My homework problem is as follows: I start with a string of 0's of length $n$. Then, with probability $p$, I change each bit to 1 identically and independently. After that, I count the number of 1's in the string (which, on average, will be $np$). Let the number of 1's be $P_1$. Continue this process 100 times, get integers $P_1, P_2,\ldots, P_{100}$. I want to estimate maximum in the set $\\{P_1, P_2, .., P_{100}\\}$.


title: Entering panel-data (cross-sectional time-series data) into SPSS

text: I don't know how to enter pooled data (I have data about 42 countries over 7 years per variable) into SPSS and then run a regression. So far, SPSS does not assign the observations to a certain country and year and rather treats each variable as independent. **How do I tell SPSS that I have 42 observations (one for each country) for 7 years per variable and not 294 (7*42) observations?** I've been looking at interaction variables and think this might be the key. Do I need an interaction term (Country*Year*Variable)? Once I manage to enter this correctly, 4 independent variables make up one latent variable which is supposed to have an effect on one dependent variable. I will use SPSS Amos for this purpose.


title: Density functions with ellipsoidal symmetry

text: I'm currently working on Hidden Markov models. I know that we must first train the observation sequence to obtain the model parameters such as the initial and transition probability. And this is done using the Baum-Welch (L.E Baum et al , 1970). The hidden states are continuous multivariate probability density functions which are assumed to have ellipsoid symmetry. But this algorithm encounters difficulties when the densities of the states are not log- concave.(Liporace, 1982). Please help me in understanding the relationship between the density which is non-log concave and the density which has an ellipsoid symmetry. Thanks


title: Extended Cox model with continuous time dependent covariate - how to structure data?

text: I need to run an extended Cox model with a time-varying covariate in R: let’s call it the "number of doses" (`X`). I am interested in the hazard ratio associated with **each level** of `X`, ie. how an additional dose affects the likelihood of recovery. `X` is a discrete, positive integer. It can either remain constant or increase by 1 during each observation time. `status=1` represents recovery and `d1`,`d2`,...`d6` are the number of doses received (created for modeling purposes) I am not sure how to set up the data. I’ve come up with two alternatives based on my reading and am attaching an excerpt of the data set-up and the code for each method, however I’m not sure which (if either) is correct for answering this question? They give fairly different results. The data has been artificially parsed at each observation time (`stop`):                   patient     start      stop doses d1 d2 d3 d4 d5 d6 status               1  0.000000  0.034000     1  1  0  0  0  0  0      0               1  0.034000 12.949399     2  1  1  0  0  0  0      0               1 12.949399 13.813242     3  1  1  1  0  0  0      0               1 13.813242 30.835070     4  1  1  1  1  0  0      0               2  0.000000  0.240000     1  1  0  0  0  0  0      0               2  0.240000  2.984179     2  1  1  0  0  0  0      0               2  2.984179  4.014723     3  1  1  1  0  0  0      0               2  4.014723  5.186506     4  1  1  1  1  0  0      0               2 20.869955 29.832999     4  1  1  1  1  0  0      0               2 29.832999 32.063887     5  1  1  1  1  1  0      0               2 32.063887 37.924743     6  1  1  1  1  1  1      1      * * * METHOD 1: **treat the number of doses as a factor**                   dt<-read.csv('data.csv',header=T)         surv<-Surv(start,stop,status,data=dt)         m<-coxph(surv ~ factor(doses),data=dt)      METHOD 2: **treat each dose as a binary variable**                   dt<-read.csv('data.csv',header=T)         surv<-Surv(start,stop,status,data=dt)         m<-coxph(surv ~ d1+d2+d3+d4+d5+d6, data=dt)      Does either method take into account that a patient had (`n-1`) doses in the previous period?


title: How to compute $f(R_1,R_2)$ from $f(R_1|X_1)$ and $f(R_2|X_2)$ given $X_1, X_2$ are dependent RVs?

text: Given $R_1=g(X_1,Y_1,Z_1)$ and $R_2=g(X_2,Y_2,Z_2)$ $Y_1,Y_2$ are iid jointly Normal $Z_1,Z_2$ are iid jointly Normal I have computed the conditional pdfs $f(R_1|X_1)$ and $f(R_2|X_2)$ and know that $X_1,X_2$ are dependent jointly Normal with correlation $\rho$ How can I compute $f(R_1,R_2)$ from $f(R_1|X_1)$ and $f(R_2|X_2)$ given $X_1, X_2$ are dependent RVs? **Edit:** To clarify based on the comments, I'm specifically looking at the case where $X_i$'s, $Y_i$'s and $Z_i$'s are all independent of each other with the only dependence being between $X_1$ and $X_2$


title: Means of groups A and B differ significantly. Can I build a Classifier that will classify values into A or B?

text: With regard to this question about classifications of groups which differ significantly, I have two somewhat newbie / general questions:   1. Lets say I'm doing a hypothesis test (the difference between two means) and I reject the null with 95% confidence (alpha = 0.05). Does this imply that I can predict/classify new values into A or B with a success rate of 95% ? If so, how this can be done ?   2. Not sure in what realm my question would fall, but does this somehow relate to the explanatory vs predictive modeling discussion ?


title: Persistent cluster IDs over similar inputs with k-means

text: I have multiple `kmeans` plots that I have generated in R. Specifically I have $5$ weeks and I generate $1$ `kmeans` plot per week. I am clustering on vectors. Most vectors in the $5$ `kmeans` plots will occur in each plot. What I am interested in determining is which vectors have changed cluster membership. To make this clear suppose I have a vector identified by the word "soccer" then I would like to see which cluster it belongs to in week1, week2, and so forth. Testing for membership change should be a simple task, for each `kmeans` clustering, each point is given an ID as to which cluster it belongs. I have 4 clusters so each week our example vector with name "soccer" could be tagged $1$, $2$, $3$, or $4$. The naive solution would be to just check the tag for a particular vector each week. However, this is not the case because R randomly selects the tags for each cluster. I know this is the case because each cluster represents some class of curves. You can visibly see that the kmeans algorithm has partitioned the vectors into 4 classes of curves. Are there ways to make the tag IDs for each cluster stay constant? That is if the cluster tagged in week1 with ID $2$ is the linear curve, then the clusters tagged $2$ in all remaining weeks will always be the linear cluster. Are there any initial conditions I can pass to `kmeans` to make this happen? I posted this question here instead of stackoverflow because I believe this questions requires more understanding of the `kmeans` algorithm.


title: Most relevant algorithms for Collaborative Filtering to test against

text: I am working on algorithms for collaborative filtering (CF). As part of this work, I want to compare a new algorithm to previous approaches to the problem. I am also surveying the most important methods and publications in publications. So far I have skimmed through around 50 articles, but I have difficulty in pin-pointing the most relevant/interesting. I want to test my CF-algorithm against a set of other CF-algorithms, both some simple baseline methods and some of the newer methods. Which methods should I compare my algorithm to? I am going to implement (or find implementations) as many of these as possible during the next couple of months, so it should be possible to get around a fair bit of different methods.


title: Find the limiting distribution of $W_n$

text: Having a little trouble with this one: Suppose $X_1, X_2, \ldots $ are iid standard normal random variables. Let $W_n = \sqrt{n} \frac{X_1 + \cdots + X_n}{X_1^2 + \cdots + X_n^2}$. Find the limiting distribution of $W_n$ as $n \to \infty$. Too bad convergence in distribution isn't closed under division. Can't get slutsky's to apply.


title: Pointers for understanding the derivation of inference in linear dynamic systems

text: I am trying to learn about the inference and maximization basically EM of the linear dynamic systems(Kalman filters for example) from Bishop's book of Pattern Recognition and Machine Learning. However, I am not being able to follow the derivations given in there. I have got the basic idea of what Kalman filters are and what they are used for. However, I am a bit confused with the learning steps(basically the derivation of the equations and all). They seem a bit complicated. I have spent lot of time trying to figure them out. But I still have some issues. Can anyone suggest me where I can get the idea. Because in the book, they haven't given the details of how it is derived(the equations. My question is how this is derived. It is a piece from the 2nd last image. I might be asking too much but I am really finding it difficult to get how this is derived. I would really appreciate if someone could give me some pointers ![enter image description here](http://i.stack.imgur.com/5mbZj.png) ![enter image description here](http://i.stack.imgur.com/yJL5g.png) ![enter image description here](http://i.stack.imgur.com/JuGNw.png) ![enter image description here](http://i.stack.imgur.com/XNsQy.png) ![enter image description here](http://i.stack.imgur.com/rIPGS.png) ![enter image description here](http://i.stack.imgur.com/FaQSN.png) ![enter image description here](http://i.stack.imgur.com/5RMci.png) ![enter image description here](http://i.stack.imgur.com/tnIaC.png) ![enter image description here](http://i.stack.imgur.com/GtpOL.png)


title: Obtaining resampling based estimates of prediction error in boosted regression tree model

text: I try to use the `gbm.fit()` function for a boosted regression tree model implemented in the R package gbm. To investigate e.g., the bootstrapped prediction error and all other functionalities I want to use the `errorest()` from the ipred package. I think `errorest()` does not accept the `gbm` output. Is there a workaround? Sorry, for the missing example. Please, see below               library(ipred)     library(gbm)     data(BostonHousing)     test <- gbm(medv ~ ., distribution = "gaussian",  data=BostonHousing)      I am not sure how to use the result in `errorest()`. Can someone give me a helping hand? Thanks!


title: Choosing optimal constraints for optimisation when the potential constraints are correlated

text: I am characterising some detectors. Each detector consists of a series of elements, thinking of a CCD in a camera is a reasonable approxiamtion. For each detector I have a distribution representing the response of each pixel. In the CCD analogy it's broadly equivalent to plotting a CDF showing the intensity required for a given pixel to fire (I don't know if CCDs work that way). While each detector is unique the plots are all relatively similar (and nothing like a parametric distribution). For I have a model that produces this curve as it's output. I wish to run a little optimisation task by setting constraints on the shape of the curve. For example I may want 20% of the voxels to achieve a minimum value and 80% to be less than some maximum. I am trying to think of an optimum way to decide which value to use as constraints. For example the value at which 20% of the pixels fire is closely correlated to many other values, whilst the value at which 40% of the pixles fire much less so. Is there a method by which I can use my correlation matrix to decide which would be the three best values to place a constraint on.


title: Determining sample size for uniform distribution

text: May you help me to decide what is the minimal sample size for a uniform distributed sample. Assume that I've find the sample average, standard deviation and the $\alpha$.


title: R: rugarch forecast function with skew-student distribution (sstd) returns "illegal" shape estimates

text: I sometimes get shape values below 4 when running rugarch::ugarchforecast(), but rugarch:::.sstdexkurt() doesn't accept shape below 4. rugarch:::.sstdexkurt() is supposed to return excess kurtosis of the skew- student distribution. Does this mean that rugarch::ugarchforecast() returns incorrect shape estimates? Can I use shape estimates below 4 to estimate the distribution?


title: Hypergeometric test

text: I have compared 1370 genes which are obtained from ChIP-sequencing analysis and 652 genes which are differentially regulated genes obtained by analyzing affymetrix 430.2.0 mouse array. When I intersect both these lists of genes, I get 37 genes in common. I would like to calculate the significance of this overlap. I was thinking to use `phyper` in R, but this requires the total number of genes. I am confused about which number to give. Should I give the total number of probes from affymetrix ChIP, or should I give the whole mouse genome number from ChIP- sequencing data? Can anybody suggest some ways to perform this significance test in R? Thanks.


title: Why is Lasso regression for high dimensional data better than Stepwise AIC?

text: I know Lasso eventually set some parameters to zero, acting like variable selection. I also read from paper talking about automated variable selection method like Stepwise AIC can be troublesome. So what are the advantages of using Lasso for variable selection over using the automated procedures like Stepwise AIC?


title: Generalized linear model with lasso regularization for continuous non-negative response

text: I have a big data problem with a large number of predictors and a non-negative response (time until inspection). For a full model I would use a glm with Gamma distributed response (link="log"). However I would like to find a small model. The "best subset glm" approach does not work for me as I run out of memory - it seems that it is not efficient enough for my setting (big data, weak computer). So I switched to the LASSO approach (using R packages `lars` or `glmnet`). `glmnet` even offers some distribution families besides the Gaussian but not the Gamma family. How can I do a lasso regularization for a glm with Gamma distributed response in R? Could it be a Cox-model (Cox net) for modelling some kind of waiting time? EDIT: As my data consists of all data points with the information about the time since the last inspection it really seems appropriate to apply a COX model. Putting data in the right format (as `Surv` does) and calling `glmnet` with `family="cox"` could do the job in my case of "waiting times" or survival analysis. In my data all data points "died" and the Cox model allows to analyse which ones "died" sooner. It seems as if in this case `family="gamma"` is not needed. Comments are very welcome.


title: Subtracting mahalanobis distance in one direction from total

text: Sincere appologies for what may be a stupid question, I'm new to pattern recognition and have exhausted all the textbooks which are aproachable to a beginner. I have substantially reworded this post to hopefully help explain myself better. Problem: Lets say I want to build a classification program which is designed to detect whether or not a patient is likely to have a genetic heart disorder based on their heart rate and blood pressure. While age does not influence the risk of the disease it is a huge determinant of HR and BP. Partial solution: I know that I can build a data set [AGE HR BP] of many healthy and diseased individuals and use the means and covariance matrices of each data set to perform a mahalanobis based classification. In this case I could then use the mahalanobis distance of a patients data to the class means of the healthy and unhealthy groups. However as I understand it the patients age would be partilly a determinant of whether they were healthy or not which is what we want to avoid. Question: Is there a way in which the age of a patient could be used to account for variances in HR and BP by including it in the covariance matrices, without it forming part of the decision making process. to avoid miss- classifications such as - [patient is over 50, majority of patients in unhealthy group are over 50 (by chance) therefore patient is unhealthy].


title: Avoiding numerical error propagation

text: I simulated a MA(1) process $y_t = \mu + \epsilon_t + \theta \epsilon_{t-1}$ in R and tried to recover the error terms $\epsilon_t$ just by computing them successively (supposing I know the parameters $\mu$ and $\theta$).               # Variables     n = 30     mu = 10     theta = 5     sigma = 3          # Simulation of the MA(1) process     y = rep(1, n)     epsilon1 = rep(1, n)          prev.epsilon = 0     for (t in 1:n)     {         epsilon1[t] = rnorm(1, mean = 0, sd = sigma)         y[t] = mu + epsilon1[t] + theta * prev.epsilon         prev.epsilon = epsilon1[t]     }          # Recovery of the error terms     epsilon2 = rep(1, n)          epsilon2[1] = y[1] - mu     for (t in 2:n)         epsilon2[t] = y[t] - mu - theta * epsilon2[t - 1]      The problem is that that recovery process involves numerical errors that propagate until the recovered values of $\epsilon_t$ radically diverge from the true ones. Yet, I fixed the sample size to 30, which is very small. The values start diverging from the 18th term and here are the last values:               > epsilon1[n]     [1] -0.6801702     > epsilon2[n]     [1] -595.0506      How can I deal with that issue? Note that the actual sample size is much larger than 30.


title: Gaussian Mixture Model - Density estimation

text: I am working on the estimation of density function and I am using the scikit- learn library. As example of code I am looking at the code at the following link. My doubt is at this following line of code:               Z = np.log(-clf.score_samples(XX)[0])      The function clf.score_samples returns the Log probabilities of each data point in X. Why is it applied the log for the second time? What the Log(Log P(X) ) stands for?


title: NP hard implementation optimisation using Monte Carlo method

text: I need to implement an algorithm ( or find an implementation) and optimise it using Monte Carlo method. This must be an NP hard such as the Travelling Salesman problem or the Knapsack problem. How can a problem be optimised exactly? Also I need to have the optimum value and compare with that, is there any site that gives the optimum for such problems? (since I can't find anything) Thanks


title: How do I do temporal correlations of matrices?

text: I have a some data like the following explaining the presence of a relation between various entities `(A,B,C,...)` in my system at time `t`:               At T=1             A B C D . . .     A  0 2 1 0 . . .     B  2 0 1 3 . . .     C  2 1 0 2 . . .     D  0 3 2 0 . . .     .  . . . . . . .           At T=2             A B C D     A  0 3 1 0 . . .     B  3 0 1 3 . . .     C  2 1 0 2 . . .      D  0 3 2 0 . . .     .  . . . . . . .           At T=3     ...      My goal is to understand how the matrix is changing temporally and characterize it in a meaningful way to be able to extract temporal points where the matrix underwent a significant change (again some meaningful metric) and also figure out the region inside the matrix where this change happened. I am just wondering if there is a standard way of doing this. Any suggestions? I am using Python and NumPy.


title: Assessing if set of sequences increase or decrease

text: I have a set of sequences of numbers, each sequence is independent from each other. I'd like to know if, "in general", these sequences increase, decrease or remain the same. What I have done so far is a fitted a linear model to each sequence, so that I can use the gradient to determine if the sequence is increasing or not. I am then using a Wilcoxson U test to test if to compare if the positive gradients are as large as the negative ones (in absolute value, of course). Is this a good solution to my problems? What are the threats of this solution? What would be a better one?


title: Looking for pattern of events in a time series

text: I have a time course experiment that follows 8 treatment groups of 12 fish for 24 hours with observations made at 5 second intervals. Among the measurements made is how far each fish travels (in mm) between observations. The 24 hours are divided into 1 dark period and 1 light period. Here is a plot of the movements of the 12 individual fish in treatment group H for the first hour of the dark period: ![control group during 1st hour of dark](http://i.stack.imgur.com/naIsV.png) You can see that some fish have long periods of inactivity, some short periods, and some have none _during this particular window_. I need to combine the data from all 12 fish in the treatment group in such a way as to identify the length and frequency of the rest periods during the entire dark period and the entire light period. I need to do this for each treatment group. Then I need to compare the differences between their rest period lengths and frequencies. I'm not a stats gal, and I'm completely at sea. The problem resembles sequence alignment to me (my bioinfomatics background), so I'm thinking Hidden Markov models, but this may be way off base. Could anyone suggest a good approach to this problem and perhaps a small example in R? Thanks!


title: Singular gradient error in nls with correct starting values

text: I'm trying to fit a line+exponential curve to some data. As a start, I tried to do this on some artificial data. The function is: $$y=a+b\cdot r^{(x-m)}+c\cdot x$$ It is effectively an exponential curve with a linear section, as well as an additional horizontal shift parameter ( _m_ ). However, when I use R's `nls()` function I get the dreaded " _singular gradient matrix at initial parameter estimates_ " error, even if I use the same parameters that I used to generate the data in the first place.   I've tried the different algorithms, different starting values and tried to use `optim` to minimise the residual sum of squares, all to no avail. I've read that a possible reason for this could be an over-parametrisation of the formula, but I don't think it is (is it?)   Does anyone have a suggestion for this problem? Or is this just an awkward model? A short example:               #parameters used to generate the data     reala=-3     realb=5     realc=0.5     realr=0.7     realm=1     x=1:11 #x values - I have 11 timepoint data     #linear+exponential function     y=reala + realb*realr^(x-realm) + realc*x     #add a bit of noise to avoid zero-residual data     jitter_y = jitter(y,amount=0.2)     testdat=data.frame(x,jitter_y)          #try the regression with similar starting values to the the real parameters     linexp=nls(jitter_y~a+b*r^(x-m)+c*x, data=testdat, start=list(a=-3, b=5, c=0.5, r=0.7, m=1), trace=T)      Thanks!


title: Problem on parametric learning when datasets are small

text: I'm currently writing a program to learn a TAN (Tree-Augmented Bayesian network) classifier from data, and I have almost finished it. I use the algorithm described in Friedman's paper 'Bayesian Network Classifiers', and for continuous variables, I assume a Gaussian distribution. I have tested my program on some dataset from UCI repo and basically it works fine. However my original task requires incremental training (and the data are all continuous). The problem is that initially the size of dataset is small, and there may be a variable whose values in the whole dataset are the same, which results in a variance of 0. The algorithm I use requires calculating the conditional mutual information of every to attribute variables give the class variable C, below is the equation I use for continuous variables, $$I(X_i; X_j|C) = -\frac{1}{2} \sum_{c=1}^{r} P(c) \log(1-\rho_c^2(X_i,X_j)),$$ $$\rho_c(X_i,X_j)=\frac{\sigma_{ij|c}}{\sqrt{\sigma_{i|c}^2\sigma_{j|c}^2}}$$ As you can see, if one variable's variance is 0, I will run into a zero division error. In order not to crash my program, in such cases I set the variable's variance to 1 hundredth of its expectation or 0.01 if the expectation is 0 (which means all values are 0). However, when there are two such variables and their values are all the same (in the case I came across, all 0), the program still crashes because then it will be calculating log of 0. How can I deal with such situations? As in an all discrete-variable case I use a Dirichlet prior to avoid dividing by zero, maybe I could similarly use some kind of prior for continuous variables? What are the common methods to deal with small datasets? I'm poor in statistics, hoping my question isn't too annoying.


title: Building a linear model for a ratio vs. percentage?

text: Suppose I want to build a model to predict some kind of ratio or percentage. For example, let's say I want to predict the number of boys vs. girls who will attend a party, and features of the party I can use in the model are things like amount of advertising for the party, size of the venue, whether there will be any alcohol at the party, etc. (This is just a made-up example; the features aren't really important.) My question is: what's the difference between predicting a ratio vs. a percentage, and how does my model change depending on which I choose? Is one better than the other? Is some other function better than either one? (I don't really care about the specific numbers of ratio vs. percentage; I just want to be able to identify which parties are more likely to be "boy parties" vs. "girl parties".) For example, I'm thinking:   * If I want to predict a percentage (say, `# boys / (# boys + # girls)`, then since my dependent feature is bounded between 0 and 1, I should probably use something like a logistic regression instead of a linear regression.   * If I want to predict a ratio (say, `# boys / # girls`, or `# boys / (1 + # girls)` to avoid dividing-by-zero errors), then my dependent feature is positive, so should I maybe apply some kind of (log?) transformation before using a linear regression? (Or some other model? What kind of regression models are used for positive, non-count data?)   * Is it better generally to predict (say) the percentage instead of the ratio, and if so, why?


title: Determining an optimal discretization of data from a continuous distribution

text: Suppose you have a data set $Y_{1}, ..., Y_{n}$ from a continuous distribution with density $p(y)$ supported on $[0,1]$ that is not known, but $n$ is pretty large so a kernel density (for example) estimate, $\hat{p}(y)$, is pretty accurate. For a particular application I need to transform the observed data to a finite number of categories to yield a new data set $Z_{1}, ..., Z_{n}$ with an implied mass function $g(z)$. A simple example would be $Z_{i} = 0$ when $Y_{i} \leq 1/2$ and $Z_{i} = 1$ when $Y_{i} > 1/2$. In this case the induced mass function would be $$ \hat{g}(0) = \int_{0}^{1/2} \hat{p}(y) dy, \ \ \ \hat{g}(1) = \int_{1/2}^{1} \hat{p}(y)dy$$ The two "tuning parameters" here are the number of groups, $m$, and the $(m-1)$ length vector of thresholds $\lambda$. Denote the induced mass function by $\hat{g}_{m,\lambda}(y)$. I'd like a procedure that answers, for example, "What is the best choice of $m, \lambda$ so that increasing the number of groups to $m+1$ (and choosing the optimal $\lambda$ there) would yield a negligible improvement?". I feel like perhaps a test statistic can be created (maybe with the difference in KL divergence or something similar) whose distribution can be derived. Any ideas or relevant literature? **Edit:** I have evenly spaced temporal measurements of a continous variable and am using an inhomogenous Markov chain to model the temporal dependence. Frankly, discrete state markov chains are much easier to handle and that is my motivation. The observed data are percentages. I'm currently using an ad hoc discretization that looks very good to me but I think this is an interesting problem where a formal (and general) solution is possible. **Edit 2:** Actually minimizing the KL divergence would be equivalent to not discretizing the data at all, so that idea is totally out. I've edited the body accordingly.


title: How to create a simulation for revenue distribution in some industry?

text: I'm doing some research on some industry, and I'm trying to predict some revenue stats. I know that the industry worth $220M and there are 360 businesses. I assume that the business income distribution would be similar to the income distribution in America: http://www.theglitteringeye.com/images/us-income- distribution.gif Can you please help me to come up with a function that will help me to simulate such distribution based on the amount of businesses and the industry worth? Thanks. Btw, I don't have any education in Stats, so sorry if the question is dumb.


title: How do we use logistic regression (scikit-learn) to predict values

text: Logistic regression can help to predict a value whether it would happen or no. I'd like to know how can I do that using sklearn. I'd like to know the probability if this event would happen or no. I have a huge dataset (20K lines and 20 columns). My data has 19 columns as predictors and last column as target (values between 0-10). To simplify work, I am using random data to understand how can I interpret data. A,B,C : Predictors target: as a target               from sklearn import linear_model     import pandas as pd          dataset = pd.DataFrame({'A':np.random.rand(100)*1000, 'B':np.random.rand(100)*100,  'C':np.random.rand(100)*10, 'target':np.random.rand(100)})          predictors= dataset.ix[:,['A','B','C']].values     target = dataset.ix[:,['target']].values          lr = linear_model.LogisticRegression()          lr.fit(predictors, target)          linear_model.LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,               intercept_scaling=1, penalty='l2', random_state=None, tol=0.0001)      Now, should I plot (lr.predict_proba) to get probability of every element ? what should I do in order to have probability of every line.


title: Trouble fitting a simple linear regression

text: I have been trying to do a simple linear regression of x3=Weeks Claimed against x4=Weeks compensated. I am having problems because my residual standard error is very high and also my residuals are not homoscedastic. If anyone could provide some guidance on how to fix these problems I would really appreciate it. Here is the output from R. ![IMG](http://i.imgur.com/q18e0zh.jpg) ![IMG2](http://i.imgur.com/DJsKegq.jpg) ![IMG](http://i.imgur.com/hYBiBwm.jpg) Thanks!!


title: Confidence interval calculated from values with different numbers of sample

text: I'm trying to calculate the estimated win probability between two players whose ratings (assume it's something like Elo rating in chess) are known. E.g. given player A's rating of 2000, player B's rating of 1800, A will have a higher win probability than B that can be determined by a function f, e.g. A has a 60% chance of beating B. The rating value itself is based on result of past games. Players with few games may have an unstable rating (i.e. not enough samples to correctly identify the right rating). Hence I'm treating the number of games equivalent to number of samples, and the more games the higher confidence I have for the player's ratings. My question is, given a different confidence level for player A and player B rating above (say player A has played 100 games while player B only played 1000 games), is there a way for me to derive the confidence level of the win probability?


title: Understanding stationarity with Inflation

text: I am looking at the link between inflation and insolvencies for an econometrics project. I have the raw quarterly insolvency data and raw quarterly CPI data for the UK (roughly 100 samples) from 1988-2013. Both of these are non-stationary at the 5%, and both are I(1) on checking the raw data. I am using 5% levels for all tests. However, my hypothesis is that **differences from the mean rate of inflation affect insolvencies** (so high or low values of inflation both cause higher insolvencies).For this I think I should use difference of logs of CPI, giving me the inflation rate (and then do the subtraction of the mean). Calculation in stata : abs(d.logCPI - mean(d.logCPI)). And, I am using the log of Insolvencies as I would expect there to be some exponential growth over time as the number of firms grows. I have a few questions on this:   1. Should I be testing for stationarity on the raw data or the variations on which my hypothesis is based (d.log CPI and logInsolvencies)   2. Difference of the log of CPI is stationary when I run dfuller, is this therefore I(0) or is it I(1) because the raw data is? My hypothesis is based on the difference of the log so is the percentage change CPI (inflation) technically my raw data?   3. The test on the Insolvencies raw data showed I(1) but the log tests as stationary without differencing so is that therefore I(0)? Using difference of logs of CPI and the log of Insolvencies will show me the impact (in number of insolvencies) inflation has on it. Basically, I’m trying to understand whether: 1\. I should check for stationarity on the raw data, the changed data, or both. 2\. How this is interpreted with order of integration. 3\. If the log of the data appears to be stationary, does that imply it is integrated to I(1)?


title: Exponential Distribution PDF: Why multiply by lambda?

text: Hoping you can help me understand the probability density function for the exponential distribution. Given that the distribution's PDF is described as follows when x > 0: $$\lambda e^{-\lambda x}$$ Why do we multiply by λ? Sure I'm missing something basic here, but seems that if we're modeling time to event, using e and the exponents would be enough. Is it because we might expect the event to happen multiple times in the time interval? Thank you in advance.


title: Covariance of logit transformed variables

text: I am having trouble evaluating the covariance of a non-linear transformation for two vectors $s_0^*$ and $s_1^*$ $$\left| \begin{array}{c} s_0^*\\\s_1^* \end{array} \right| \sim N\left( \left| \begin{array}{c}\mu_0^*\\\ \mu_1^* \end{array} \right|,\left| \begin{array}{c}\Sigma_s & \Gamma_{01} \\\ \Gamma_{10} & \Sigma_s \end{array} \right| \right)$$ I want to find the variance expression of $$ s_1 -s_0=\frac{\exp(s_1^*)}{1+\exp(s_1^*)}-\frac{\exp(s_0^*)}{1+\exp(s_0^*)}=\text{logit}(s_1^*)-\text{logit}(s_0^*)\\\$$ $\text{Var} \left(\frac{\exp(s_1^*)}{1+\exp(s_1^*)}-\frac{\exp(s_0^*)}{1+\exp(s_0^*)}\right)= \text{Var}( s_1)+\text{Var}(s_0)-2\text{cov}(s_1,s_0)$, I know how to find the variance of the two variance expressions from by the delta method, but I'm and having trouble finding an expression for the covariance between the two. I'm guessing that I should use Taylor expansions somehow, but it stops for me pretty quickly. Any tips on how I should go about finding this covariance? I am familiar with the relation $\text{cov}(s_0,s_1)=\text{E}(s_0s_1)+\text{E}(s_0)\text{E}(s_1)$ so is it the case that all I have to do is find $\text{E}(s_0 s_1)$? I'm pretty stuck here so any help is greatly appreciated.


title: Can covariates be categorial variables?

text: I have a 2x2 experimental design. In the experiment, I also collected the participants' professional qualifications (categorical variable- yes /no). I would like to test the effect of controlling this variable on the MSE. I understand that in ANCOVA, the covariate is either an interval or ordinal variable. So should I refer to the test as an ANCOVA if I include professional qualification as an additional variable? Or should it be a blocking variable since "prof. qualification" is categorical. Thanks in advance for your replies!


title: There is a load pdf and a capacity pdf, what is pdf of the surplus (the difference)?

text: I have a probability distribution function for expected available energy $\text{P}(t)$ and another pdf for expected load $\text{E}(t)$. How do I find a pdf for the set of random variables of $\text{P}(t)-\text{E}(t)$, i.e. the surplus energy of the system? Is it a stress-strength model?


title: Is there a limit to the SPSS random number generator?

text: Will the random number generator in SPSS be serviceable if I need 250,000 random numbers, or will the randomness start to degenerate? Asked another way, what practical limits are there to using the random number generator in SPSS to generate large numbers of random numbers?


title: Code to do Monte Carlo using SPSS?

text: Is SPSS command code available for doing Monte Carlo research? Thanks.


title: MA terms in arima

text: With the arima function I found some nice results, however now I have trouble interpreting them for use outside R. I am currently struggling with the MA terms, here is a short example:               ser=c(1, 14, 3, 9)        #Example series     mod=arima(ser,c(0,0,1))   #From {stats} library     mod          #Series: ser     #ARIMA(0,0,1) with non-zero mean     #     #Coefficients:     #          ma1  intercept     #      -0.9999     7.1000     #s.e.   0.5982     0.8762     #     #sigma^2 estimated as 7.676:  log likelihood = -10.56     #AIC = 27.11   AICc = Inf   BIC = 25.27          mod$resid          #Time Series:     #Start = 1     #End = 4     #Frequency = 1     #[1] -4.3136670  3.1436951 -1.3280435  0.6708065          predict(mod,n.ahead=5)          #$pred     #Time Series:     #Start = 5     #End = 9     #Frequency = 1     #[1] 6.500081 7.100027 7.100027 7.100027 7.100027     #     #$se     #Time Series:     #Start = 5     #End = 9     #Frequency = 1     #[1] 3.034798 3.917908 3.917908 3.917908 3.917908     ?arima      When looking at the specification this formula is presented: `X[t] = a[1]X[t-1] + … + a[p]X[t-p] + e[t] + b[1]e[t-1] + … + b[q]e[t-q]` Given my choice of AR and MA terms, and considering that I have included a constant this should reduce to: `X[t] = e[t] + b[1]e[t-1] + constant` However this does not hold up when i compare the results from R with manual calculations: `6.500081 != 6.429261 == -0.9999 * 0.6708065 + 7.1000` Furthermore I can also not succeed in reproducing the insample errors, assuming i know the first one this should be possible: `-4.3136670 * -0.9999 +7.1000 != 14 - 3.1436951` `3.1436951 * -0.9999 +7.1000 != 3 + 1.3280435` `-1.3280435 * -0.9999 +7.1000 != 9 - 0.6708065` I hope someone can shed some light on this matter so I will actually be able to use the nice results that I have obtained.


title: Is there a better alternative to a logit / probit regression when all dependent variables are dichotomous?

text: I'm working on a clinical trial dataset with binary response. All independent variables are also binary. My first impulse was to simply run a standard logit / probit regression and be done with it. But after thinking about it for a little longer, I kept thinking: is there something else I could do? My concerns arose from the fact that all variables on the RHS (right hand side, i.e., independent variables) are also binary (which is a first for me). Any suggestions?


title: Normalizing when doing count statistics of patient numbers on different arms of a clinical trial

text: I'm looking at clinical trial data where there are various numbers of patients on different arms (A, B, C, D). All arms receive a drug that is known to cause a toxicity, but some arms receive an additional drug which might contribute to the toxicity. I have been using 3 methods to identify the number of patients with the toxicity, and I was curious if the ratios of patients on the different arms differed between the methods and whether the ratios might have been higher on some arms (suggesting the additional toxicity of the added drug) than what is seen with the entire population. For example, number of patients identified with the toxicity with the different methods:                            A     B    C    D   Total     method 1:    8    18   26   16      68     method 2:   45   111   53  113     322       method 3:   68   158   69  153     448       Overall:   695  1331  441  925    3392      To compare, I calculated the percentage of patients on each arm for each of the methods (Method 1 Arm A: 8/68) and compared it to the overall population (695/3392). However, after doing this, I was told I need to normalize by the patients on each arm. I am not sure what is meant by this. I can divide the numbers of patients identified for each arm by the total number of patients on the arm, but then how does this help me do the comparison I am interested in? I don't have a stats background, so I am not sure how to tackle this question in general.


title: probablistic output for binary SVM classification

text: I'm using libSVM for a binary classification problem. After a test instance is assigned a label (1 or -1), I also want to know how likely it is assigned such a label. I'm thinking about calculating the distance from the instance to the hyperplane in the feature space. The larger the distance is, the more likely the label assignment is correct. Is my idea correct? Or is there already such an option in libSVM for my purpose?


title: Contradictory Results in autocorrelation tests

text: I have a time-series model, with stationary variables. Testing for outocorrelation hasn't been easy: \- a can't calculate DW stat because of my small number of observations; \- the GB test indicates that there is autocorrelation \- the Correlogram Q statistics tells me differently, indicating that there is no autocorrelation... what can cause this contradictory results? how do i manage this?


title: Legend and TraMiner seqplot.tentrop function

text: I have a problem with the legend of a TraMineR graph. I have computed the 4 different entropies for the 4 cohorts that I'm analyzing and I would like to plot them in the same graph using the `seqplot.tentrop` function. However, the default legend is too big and I would plot it separately reducing the dimensions, but I can't reproduce the colors of the legend in the new one I've created with the last command below.                   plot_colors <- c("blue", "red", "forestgreen", "purple")         seqplot.tentrop(myseq, group = job$cohort, ylim = c(0, 1),            col=plot_colors, withlegend=FALSE)         legend("bottomright", names(myseq),            c("1974-1978","1982-1986","1990-1994", "1998-2001"),           cex=0.8, col=plot_colors)      Is it possible to change colors/dotted lines?


title: Intrarater reliability

text: I would like to estimate the intrarater reliability in a MRI study. Different regions of interest were manually traced with a mouse-driven cursor in order to estimate their volumes. All data were traced by the same rater and a subset of the data was traced twice. How can I estimate the intra-raterreliability of the brain volumes? Specifically, which type of intraclass correlation coefficient would you suggest?


title: What is the appropriate machine learning algorithm for this problem?

text: I am working on a problem which looks like this: Input Variables   * Categorical     * a     * b     * c     * d   * Continuous     * e Output Variables   * Discrete(Integers)     * v     * x     * y   * Continuous     * z The major issue that I am facing is that Output Variables are not totally independent of each other and there is no relation that can be established between them. That is, there is a dependence but not due to the causality (one value being high doesn't imply that the other will be high too but the chances of other being higher will improve) An Example would be: v - Number of Ad Impressions x - Number of Ad Clicks y - Number of Conversions z - Revenue Now, for an Ad to be clicked, it has to first appear on a search, so Click is somewhat dependent on Impression. Again, for an Ad to be Converted, it has to be first clicked, so again Conversion is somewhat dependent on Click. So running 4 instances of the problem predicting each of the output variables doesn't make sense to me. Infact there should be some way to predict all 4 together taking care of their implicit dependencies. But as you can see, there won't be a direct relation, infact there would be a probability that is involved but which can't be worked out manually. Plus the output variables are not Categorical but are in fact Discrete and Continuous. Any inputs on how to go about solving this problem. Also guide me to existing implementations for the same and which toolkit to use to quickly implement the solution. Just a random guess - I think this problem can be targeted by Bayesian Networks. What do you think ?


title: Interval estimation of $\sigma^2$ with the reliability of $95\%$

text: **I got these values from the measurement by a telescope: 20.1, 20.2, 19.9, 20, 20.5, 20.5, 20, 19.8, 19.9, 20. I know that the actual distance is 20 km and the error of the measurement is not affected by the systematic error. What is the accuracy of the telescope? [Answer is: ${σ}^2$=0.061]** Now I am supposed to calculate ${σ}^2$ with a reliability of 95% using the interval estimation. So, ${σ}^2$ should be:   $ \frac{n* s^{2} }{ \chi^{2}_{0.975} \big(n-1\big) } ; \frac{n* s^{2} }{ \chi^{2}_{0.025} \big(n-1\big) } $ After using $n = 10$, $k = n-1 = 9$ and $s^{2} = 0.061$ I get this: $ \frac{10*0.061}{ \chi^{2}_{0.975} \big(9\big) }; \frac{10*0.061}{ \chi^{2}_{0.025} \big(9\big) } $ Is it right? In the numerator I should use n = 10 or n-1 = 9? In the denumerator I should use n-1 = 9 or n = 10 ? When I'm looking for the values in the table of $\chi^2$ quantiles, should I look for $0.95$ and $0.05$ or $0.975$ and $0.025$? When they say they need the reliability of $95\%$, I assume I will use $1- \frac{ \alpha }{2}$ and $\frac{ \alpha }{2}$. Since my $\alpha$ is $0.05$, I think $0.975$ and $0.025$ are the right values. [ Answer is: $σ$ = <0.173 , 0.433> ]


title: The values of data points near the start and the end of time series, when smoothed with simple moving average

text: I am making a task for a class, where I need to code simple moving average algorithm. In my copybook, I have this written: $$\bar{y_{i}}=\frac{y_{i-1}+y_{i}+y_{i+1}}{3}$$ for $n=3$ MA case, where near the start and end the values of smoothed time series are: $$\bar{y_{1}}=\frac{5y_{1}+2y_{2}-y_{3}}{6}$$ $$\bar{y_{n}}=\frac{5y_{n}+2y_{n-1}-y_{n-2}}{6}$$ For $n=5$ MA case it is: $$\bar{y_{i}}=\frac{y_{i-2}+y_{i-1}+y_{i}+y_{i+1}+y_{i+2}}{5}$$ with near-borders values of: $$\bar{y_{1}}=\frac{3y_{1}+2y_{2}+y_{3}-y_{4}}{5}$$ $$\bar{y_{2}}=\frac{4y_{1}+3y_{2}+2y_{3}+y_{4}}{10}$$ $$\bar{y_{n}}=\frac{3y_{n}+2y_{n-1}+y_{n-2}-y_{n-3}}{5}$$ $$\bar{y_{n-1}}=\frac{4y_{n}+3y_{n-1}+2y_{n-2}+y_{n-3}}{10}$$ But what would be the formulas to recalculate the near-border values in the general case? Thank you.


title: Mining association rules on relational data

text: I have recently started working on my master's thesis, which is a collaboration between my university and an IT company. The problem from the company's point of view is to identify correlations in their web marketing campaigns for various properties, with the purpose of applying this information to future ads in order to increase the relevant metrics (e.g. clicks or conversion rate). The area of the thesis is somewhat new to me, even though I have the basics in for example statistics, algorithms and AI. I am currently pursuing the idea of applying data mining techniques to the data set, and more specifically association rule mining as described by Agrawal et al. (1993). Association rule mining work on transactional databases rather than relational, though the latter can be converted into the former with the disadvantage that the new table will become incredibly large, especially considering the most natural transaction type in this case is an impression (i.e. one table row for every time an ad is shown) which easily leads to many, many millions of rows where only a fraction will actually lead to clicks. (A relevant question might be Find association rules for a given value, but it has not been answered yet.) My idea is thus to build a transaction table of impressions from the relational data, apply an association rule algorithm and finally filter out any rules that do not lead to a click, since such rules are at least for now of little interest. **The question I'd like to ask is if this is a reasonable approach?** I feel somewhat confident that this will work, however it doesn't feel like an obvious solution. This might not be a problem for the company who will use the result, however from an academic point of view I wish to cover my bases. **If this is not a good approach, which areas might be better suited for the purpose?** Don't hesitate to ask if there is anything you would like me to clarify!


title: Algorithms for 1D Gaussian mixture with equal variance and noise cluster

text: I would like to fit a Gaussian mixture model to some data. The data is 1D and I want to constrain all the Gaussians to have equal variance. I would also like to have a uniform background noise cluster to pick up points that are not in a real cluster. I don't know how many clusters there should be, and this can vary over a large range for different data sets (tens to thousands are possible), so ideally this should be selected automatically. It might also be useful to put some lower threshold on the mixing coefficient for each cluster to avoid fitting very tiny ones to the background noise. Depending on the data set there will be from hundreds to many thousands of data points. What would be a good way to pursue this problem? I am a bit overwhelmed by the possible options (EM, variational bayes, infinite GMM, dirichlet processes etc.), or if it is even possible.


title: some basic questions for process mining

text: Are there any good books and articles for introducing Process Mining? What are the state-of-the-art research results for process mining? What are the differences between process mining and time series prediction?


title: How can I test if given samples are taken from a Poisson distribution?

text: I know of normality tests, but how do I test for "Poisson-ness"? I have sample of ~1000 non-negative integers, which I suspect are taken from a Poisson distribution, and I would like to test that.


title: Alternative graphics to "handle bar" plots

text: In my area of research, a popular way of displaying data is to use a combination of a bar chart with "handle-bars". For example, ![enter image description here](http://i.stack.imgur.com/cZQ89.jpg) The "handle-bars" alternate between standard errors and standard deviations depending on the author. Typically, the sample sizes for each "bar" are fairly small - around six. These plots seem to be particularly popular in biological sciences - see the first few papers of BMC Biology, vol 3 for examples. So how would you present this data? **Why I dislike these plots** Personally I don't like these plots.   1. When the sample size is small, why not just display the individual data points.    2. Is it the sd or the se that is being displayed? No-one agrees which to use.    3. Why use bars at all. The data doesn't (usually) go from 0 but a first pass at the graph suggests it does.   4. The graphs don't give an idea about range or sample size of the data. **R script** This is the R code I used to generate the plot. That way you can (if you want) use the same data.                                                       #Generate the data     set.seed(1)     names = c("A1", "A2", "A3", "B1", "B2", "B3", "C1", "C2", "C3")     prevs = c(38, 37, 31, 31, 29, 26, 40, 32, 39)          n=6; se = numeric(length(prevs))     for(i in 1:length(prevs))       se[i] = sd(rnorm(n, prevs, 15))/n                                                  #Basic plot     par(fin=c(6,6), pin=c(6,6), mai=c(0.8,1.0,0.0,0.125), cex.axis=0.8)     barplot(prevs,space=c(0,0,0,3,0,0, 3,0,0), names.arg=NULL, horiz=FALSE,             axes=FALSE, ylab="Percent", col=c(2,3,4), width=5, ylim=range(0,50))                                                  #Add in the CIs     xx = c(2.5, 7.5, 12.5, 32.5, 37.5, 42.5,  62.5, 67.5, 72.5)     for (i in 1:length(prevs)) {       lines(rep(xx[i], 2), c(prevs[i], prevs[i]+se[i]))       lines(c(xx[i]+1/2, xx[i]-1/2), rep(prevs[i]+se[i], 2))     }                                                  #Add the axis     axis(2, tick=TRUE, xaxp=c(0, 50, 5))     axis(1, at=xx+0.1, labels=names, font=1,          tck=0, tcl=0, las=1, padj=0, col=0, cex=0.1)


title: Fast method for finding best metaparameters of SVM (that is faster than grid search)

text: I am using SVM models to do short term forecasting of air pollutants. To train a new model I need to find appropriate metaparameters for an SVM model (I mean C, gamma and so on). Libsvm documentation (and many other books I have read) suggests using grid search to find these parameters - so I basically train model for each combination of these parameters from a certain set and choose the best model. Is there any better way to find optimal (or near optimal) metaparameters? For me it is mainly a matter of computation time - one grid search of this problem takes about two hours (after I did some optimisations). Pros of grid search:   * It can be easily parallelized - if you have 20 CPUs it will run 20 times faster, parallelizing other methods is harder   * You check big parts of metaparameter space, so if there is a good solution you will find it.


title: Creating a maximum entropy Markov model from an existing multi-input maximum entropy classifier

text: I am intrigued by the concept of a Maximum Entropy Markov Model (MEMM), and I am thinking of using it for a Part of Speech (POS) tagger. At the moment, I am using a conventional Maximum Entropy (ME) classifier to tag each individual word. This uses a number of features, including the preceding two tags. MEMMs use the Viterbi algorithm to find the optimum path through the Markov Chain (ie. to find a complete optimum set of tags for the sentence rather than individual optimums for each word). Reading about it, this appears to be have a wonderful elegance and simplicity. However, each stage only relies on the "results" of the previous stage (i.e. as per a Markov Chain). However, my ME model uses the previous two stages (i.e. the tags for the preceding two words). It seems I have two possible approaches:   * As with a conventional Viterbi implementation, use a set of paths stored according to one (the previous) stage. My ME classifier would use this and a 'frozen' stage before this (frozen into the path under consideration) to produce the transfer function.   * Or I write the algorithm to keep track of two stages. This is more complicated and would no longer be a true Markov Model because each transfer function (i.e. from the ME Model) would depend on the preceding two stages and not one stage. It strikes me that the second will be more accurate, although it will be more complicated. I have yet to find any examples of this during my literature search. Has it been tried? Did the two stage approach give an improvement to the overall accuracy?


title: What is a good method for designing an HTM network in terms of depth and breadth?

text: Regarding Hawkins' HTM algorithm. Has anyone found a fruitful method to define the number of nodes in a layer and number of layers in a network?


title: Domain-agnostic feature engineering that retains semantic meaning?

text: Feature engineering is often an important component to machine learning (it was used heavily to win the KDD Cup in 2010). However, I find that most feature engineering techniques either   * destroy any intuitive meaning of the underlying features or    * are very specific to a particular domain or even particular types of features. A classic example of the former would be principal component analysis. It would seem to me that any knowledge that a subject-matter expert would have about the features would be destroyed by converting those features to principal components. Contrast that with a simple technique of converting a date into features for "day of month" and "day of week." The underlying meaning is still retained in the new features, but obviously this particular technique only applies for dates and not arbitrary features. Is there any standard body of feature engineering techniques that do not destroy the meaning of the underlying features while also being applicable to arbitrary domains (or at least a wide variety of domains)?


title: Bayesian and frequentist reasoning in plain English

text: How would you describe in plain English the characteristics that distinguish Bayesian from Frequentist reasoning?


title: How to calculate the boundary value for a random variable which is sum of variables with gamma and uniform distributions?

text: The variable is a sum of two random variable which obey gamma and uniform distributions, respectively. The parameters of the uniform distribution variable are determined, and the other's must be retrieved from real time calculation. I want to find the boundary of the summed variable, like using a three sigma rule would do for a normal distribution. Is there any good method to estimate that? Some inequality, like the Chebyshev inequality, might be too generic of a method so that it will not perform exactly as the three sigma rule would.


title: Fitting a particular Gaussian model

text: Using R or SAS, I want to fit the following Gaussian model: $$ \begin{pmatrix} y_{1j1} \\\ y_{1j2} \\\ y_{1j3} \\\ y_{2j1} \\\ y_{2j2} \\\ y_{2j3} \end{pmatrix} \sim_{\text{i.i.d.}} {\cal N} \left( \begin{pmatrix} \mu_1 \\\ \mu_1 \\\ \mu_1 \\\ \mu_2 \\\ \mu_2 \\\ \mu_2 \end{pmatrix} , \Sigma \right), j=1, \ldots n $$ with covariance matrix having the following structure: $$ \Sigma = \begin{pmatrix} \Sigma_0 & M \\\ M & \Sigma_0 \end{pmatrix} $$ with $\Sigma_0$ a "compound symmetry" (exchangeable) covariance matrix and $M=\delta \begin{pmatrix} 1 & 1 & 1 \\\ 1 & 1 & 1 \\\ 1 & 1 & 1 \end{pmatrix}$, $\delta \in \mathbb{R}$. Importantly, I need a general exchangeable matrix $\Sigma_0$, with possibly negative correlation. EDIT: In view of some comments (and even an answer) given below I should add a precision: I am not a beginner with PROC MIXED in SAS and nlme in R, and I know how to consult the documentations. But in spite of my experience I am not able to fit this model.


title: Logistic regression: Is it valid to code NA as 0?

text: Hi I have the following problem: I plan to do a logistic regression on web traffic data. One of my attributes is Number of days since last visit (type numeric in R). Now lots of visitors are firs-timers and I have no data for this attribute. Should I code them as NAs, or is is valid to code them as 0s? Thank you all


title: How to test for normality in a 2x2 ANOVA?

text: **Study Design:** I showed participants some information about sea-level rise, focusing the information in different ways, both in terms of the time-scale and the magnitude of potential rise. Thus I had a 2 (Time: 2050 or 2100) by 2 (Magnitude: Medium or High) design. There were also two control groups who received no information, only answering the questions for my DVs. **Questions:** I've always checked for normality within cells -- for the 2x2 portion of this design, it would mean looking for normality within 4 groups. However, reading some discussions here has made me second guess my methods. First, I've read that I should be looking at the normality of the residuals. **How can I check for normality of residuals (in SPSS or elsewhere)?** Do I have to do this for each of the 4 groups (6 including the controls)? I also read that **normality within groups implies normality of the residuals. Is this true?** (Literature references?) Again, does this mean looking at each of the 4 cells separately? In short, **what steps would _you_ take to determine whether your (2x2) data are not violating assumptions of normality?** References are always appreciated, even if just to point me in the right direction.


title: Calculate the Kullback-Leibler Divergence in practice?

text: I am using KL Divergence as a measure of dissimilarity between 2 $p.m.f.$ $P$ and $Q$. $$D_{KL}(P||Q) = \sum_{i=1}^N \ln \left( \frac{P_i}{Q_i} \right) P_i$$ $$=-\sum P(X_i)ln\left(Q(X_i)\right) + \sum P(X_i)ln\left(P(X_i)\right)$$ If $$P(X_i)=0$$ then we can easily calculate that $$P(X_i)ln\left(Q(X_i)\right)=0$$ $$P(X_i)ln\left(P(X_i)\right)=0$$ But if $$P(X_i)\ne0$$ and $$Q(X_i)=0$$ how to calculate the $$P(X_i)ln\left(Q(X_i)\right)$$


title: How many PCs should varimax rotation be applied to?

text: I have a list of 25 air pollutants many of which are strongly correlated. I was hoping to reduce down to a short list of eigenvectors which would each be composed of a small number of the pollutants. I mostly followed this tutorial video so far. When I tried varimax rotation on all of the PCs (as seems to be done at the end of this video) I get each of the 25 PCa having a score of 1 or -1 for only one variable, with every variable accounted for, i.e. all PCs and variables accounted for by unique pairings. Have I done something wrong or if not then why did the example in the video not also do the same? Could it be something to do with my data having many more dimensions? Depending on the selection criteria (Kaiser criterion; levelling-off of scree plot) the top 7 or top 13 are the important ones. Putting just these 7 or 13 components in gives more the sort of output I had expected. Does anyone have any suggestions which of these options (7 or 13 components) I should choose or is it a fairly arbitrary decision based on how much I want to compress the data?


title: What coefficient could I use to calculate the relative difficulty of a test in relation to others using only mean and population standard deviation?

text: I have a series of tests, all of them of different difficulty, and from each of them I get an average score and a population standard deviation; e.g:               86.94%  16.27%       84.17%   6.52%       99.55%   1.19%     81.70%   9.61%       99.82%   0.34%      Now, clearly the last test, that has an average of 99.82% and a standard deviation of 0.34% is a lot easier than the one with a mean of 81.7% and a standard deviation of 9.61% (supposing of course the variations of results are effect of the difficulty and not of the population that took the test). What I'm trying to figure out is a coefficient that uses only that data to get a number that properly reflects how much more difficult was one test than the others. I hope you can help. btw: sorry for the wrong tag, I'm not sure which would have been a good fit.


title: Estimating "Probability" in normal probability plot

text: I plotted normal probability plot in R using `qqnorm` and `qqline`. I want some help on:   1. How to estimate "probability" that a data has normal distribution? (I read in a paper that a probability of 0.10 is required to assume that a data is normally distributed).    2. Also, how to calculate correlation coefficient for a normal probability plot in R?   3. Is the normality test valid for nonlinear regression too? (This might be a silly question, excuse me for that!) Thanks!


title: Correlation coefficient in longitudinal data

text: We are performing a longitudinal analysis over different characteristics (variables) of a single subject of analysis. Our goal is not to make any predictions, but rather analyse correlations among these different variables. No data is ever missing, and all variables are guaranteed to have a value at a given point in time. My questions is: shall I stick to measuring correlation ignoring time (as no prediction is required for the study), modelling two variables' values at a given point in time as a mere observation (dot) in the scatter-plot, or should I use cross-correlation using time-series? Regards.


title: Out-of-sample vs. in-sample interpretation

text: I am running predictive regressions on stock returns and as expected less relations hold out-of-sample than in-sample, however in some cases I find a significant relation out-of-sample but not in-sample at the usual significance levels (90%, 95% or 99%) for a potential predictor series. Is this possible? How can this be interpreted? How could an iterative regression over the same dataset yield better results than using all the information at once... To evaluate the significance in-sample I simply check the t-statistic of the slope coefficient and out-of-sample I use a mean squared prediction error based method which compares the errors of my model containing the potential predictor to a model containing only the historical mean return.


title: Comparison of point-biserial and linear correlation coefficients

text: I have a continuous variable X that is associated with a continuous outcome Y and a dichotomous outcome Z. If I calculate Pearson's $r$ between X and Y, and the point-biserial coefficient $r_{pb}$ between X and Z, and find that the coefficient between X and Y is much higher than the coefficient between X and Z, can I say "The association between X and Y is much stronger than the association between X and Z," or are those coefficients non-comparable? If not, is there a reasonable way to compare the associations for the two pairs (X and Y, and X and Z)?


title: When do Post-hoc tests require an adjustment for multiple comparisons

text: What happens to the Type 1 error rate, the Type 2 error rate, and Power if you do not adjust for multiple comparisons when you should? What happens to these as a result of using something like Bonferroni’s adjustment?


title: How do I find the expected values for this type of Chi squared question?

text: > ![enter image description here](http://i.stack.imgur.com/tOPfV.jpg) I have been doing exam papers and preparing for my exam, but I always can't do these type of questions in statistics, after I get to know how to do one, I'll be able to do these, which are common on my papers. We have to use chi squared distribution, but I am unable to find the expected values; I don't know because they do not tell the expectation.


title: Does Normal distribution theory originate in the psychometric literature, especially reliability theory?

text: The basic statistical literature does not talk about the exact background of the normal distribution. Is the basis of this assumption in psychometry or it has an origin in pure statistics i.e. sampling theory.


title: How to quickly collect data (survey) from students in class?

text: What would be an easy and simple way to get online responses from students in class for use in class demonstrations. E.g. Height and Weight. A quick survey method in a computer class could make lab work more interesting. Any ideas? Thank you for reading.


title: Identifying structural breaks in regression with Chow test

text: I have some problems in using (and finding) the Chow test for structural breaks in a regression analysis using R. I want to find out if there are some structural changes including another variable (represents 3 spatial subregions). Namely, is the regression with the subregions better than the overall model. Therefore I need some statistical validation. I hope my problem is clear, isn't it? Kind regards   marco Toy example in R:               library(mlbench) # dataset     data("BostonHousing")          # data preparation     BostonHousing$region <- ifelse(BostonHousing$medv <=                                     quantile(BostonHousing$medv)[2], 1,                              ifelse(BostonHousing$medv <=                                     quantile(BostonHousing$medv)[3], 2,                             ifelse(BostonHousing$medv >                                     quantile(BostonHousing$medv)[4], 3, 1)))          BostonHousing$region <- as.factor(BostonHousing$region)          # regression without any subregion      reg1<- lm(medv ~ crim + indus + rm, data=BostonHousing)          summary(reg1)          # are there structural breaks using the factor "region" which     # indicates 3 spatial subregions     reg2<- lm(medv ~ crim + indus + rm + region, data=BostonHousing)      \------- subsequent entry I struggled with your suggested package "strucchange", not knowing how to use the "from" and "to" arguments correctly with my factor "region". Nevertheless, I found one hint to calculate it by hand (https://stat.ethz.ch/pipermail/r-help/2007-June/133540.html). This results in the following output, but now I am not sure if my interpetation is valid. The results from the example above below. Does this mean that region 3 is significant different from region 1? Contrary, region 2 is not? Further, each parameter (eg region1:crim) represents the beta for each regime and the model for this region respectively? Finally, the ANOVA states that there is a signif. difference between these models and that the consideration of regimes leads to a better model? Thank you for your advices! Best Marco               fm0 <- lm(medv ~ crim + indus + rm, data=BostonHousing)     summary(fm0)     fm1 <- lm(medv  ~ region / (crim + indus + rm), data=BostonHousing)     summary(fm1)     anova(fm0, fm1)      Results:               Call:     lm(formula = medv ~ region/(crim + indus + rm), data = BostonHousing)          Residuals:            Min         1Q     Median         3Q        Max      -21.079383  -1.899551   0.005642   1.745593  23.588334           Coefficients:                    Estimate Std. Error t value Pr(>|t|)         (Intercept)    12.40774    3.07656   4.033 6.38e-05 ***     region2         6.01111    7.25917   0.828 0.408030         region3       -34.65903    4.95836  -6.990 8.95e-12 ***     region1:crim   -0.19758    0.02415  -8.182 2.39e-15 ***     region2:crim   -0.03883    0.11787  -0.329 0.741954         region3:crim    0.78882    0.22454   3.513 0.000484 ***     region1:indus  -0.34420    0.04314  -7.978 1.04e-14 ***     region2:indus  -0.02127    0.06172  -0.345 0.730550         region3:indus   0.33876    0.09244   3.665 0.000275 ***     region1:rm      1.85877    0.47409   3.921 0.000101 ***     region2:rm      0.20768    1.10873   0.187 0.851491         region3:rm      7.78018    0.53402  14.569  < 2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 4.008 on 494 degrees of freedom     Multiple R-squared: 0.8142,     Adjusted R-squared: 0.8101      F-statistic: 196.8 on 11 and 494 DF,  p-value: < 2.2e-16          > anova(fm0, fm1)     Analysis of Variance Table          Model 1: medv ~ crim + indus + rm     Model 2: medv ~ region/(crim + indus + rm)       Res.Df     RSS Df Sum of Sq     F    Pr(>F)         1    502 18559.4                                      2    494  7936.6  8     10623 82.65 < 2.2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


title: Meta-analysis with zero values for mean and sd

text: I want to do a meta-analysis with count data for treatement/control cases (abundance of sea polichaetos). My problem is that I need to use zero values (an informative zero value) for the mean and standard deviation for one of the treatement, but R has a problem: "Studies with zero values for sd.e or sd.c get no weight in meta-analysis". I can agroup the case by Family (`byvar=Family`). There is an option in `metabin`, where you could add either a numerical value to each cell frequency for studies with a zero cell count or the character string "TA" which stands for treatment arm continuity correction, but I want to use `metacont` which don´t have this option. Can you help me? Thanks!               > datos<-read.table(file="datos_totales.csv",head=TRUE,sep=",",dec=",")      > datos[10:20,]     Study ControlN ControlMean ControlSD TreatedN TreatedMean TreatedSD       Family     10 Ali.Gmb.S1.Inv       3    893.3333 180.36999        3  1213.33333 410.52812 Capitellidae     11 Ali.Gmb.S2.Inv       3    640.0000 348.71192        3   666.66667 220.30282 Capitellidae     12 Ali.Gmb.S3.Inv       3    426.6667 184.75209        3   920.00000 628.64935 Capitellidae     13 Ali.Cul.S1.Ver       3    213.3333 115.47005        3     0.00000   0.00000 Cirratulidae     14 Ali.Cul.S2.Ver       3    160.0000 144.22205        3    40.00000  40.00000 Cirratulidae     15 Ali.Cul.S3.Ver       3    293.3333 234.37861        3    26.66667  46.18802 Cirratulidae     16 Ali.Cul.S1.Inv       3    653.3333 532.66625        3    13.33333  23.09401 Cirratulidae     17 Ali.Cul.S2.Inv       3    706.6667 335.45988        3    26.66667  23.09401 Cirratulidae     18 Ali.Cul.S3.In        3    666.6667 674.48746        3    80.00000  80.00000 Cirratulidae     19 Ali.Gmb.S1.Vr        3    280.0000  69.28203        3     0.00000   0.00000 Cirratulidae     20 Ali.Gmb.S2.Ver       3    226.6667 220.30282        3     0.00000   0.00000 Cirratulidae          >  library(meta)          > metaanalisis2<- metacont(TreatedN, TreatedMean, TreatedSD, ControlN, ControlMean,        ControlSD, data=datos, byvar=Family, label.c="Control", print.byvar=FALSE,        label.e="Treatment", studlab=Study)     **Mensajes de aviso perdidos     In metacont(TreatedN, TreatedMean, TreatedSD, ControlN, ControlMean,  :       Studies with zero values for sd.e or sd.c get no weight in meta-analysis**


title: Resources about Gibbs sampling in hybrid Bayesian networks

text: I've been trying to get my hands on a substantial resource for using Gibbs sampling in hybrid Bayesian networks, that is, networks with both continuous and discrete variables. So far I can't say I have succeeded. I'm interested in hybrid networks where there are no constraints regarding discrete children having continuous parents. Gibbs sampling is a very widely employed method in approximate inference in Bayesian methods, and yet, I can't seem to find detailed resources that focus on hybrid networks


title: Concentration bounds on a sequence of (0,1)-experiments

text: I have a $(0,1)$-experiment that returns $1$ with probability $p$ and $0$ with probability $1-p$. Let $X_i$ be the random variable describing the outcome of iteration $i$ of the experiment. I'm interested in the value $Y_t=\sum_{i=0}^t X_i$, which counts the number of successes since the beginning of the sequence. The sequence terminates as soon as $Y_t \geq (p+\delta) \cdot t$ (with $0 < \delta < 1-p$) is fulfilled, i.e., when the ratio between $0$ and $1$ is over the threshold $p+\delta$. Quite intuitively the sum will converge to $Y_t = p\cdot t$ (expected value) but the probability of reaching $Y_t = (p+\delta)\cdot t$ is non-zero, and will be reached in an infinite run. The question is when will it be reached? For this there are two interesting values I'd like to calculate: the probability of each step $t \rightarrow t+1$ being the last, i.e., $Pr[Y_{t+1} \geq (p+\delta) \cdot (t+1)|Y_t < (p+\delta) \cdot t]$ and the expected time until the sequence is terminated, $\mathrm{E}[\text{steps}]$.


title: Confusion related to HMM

text: I was reading this book related to Pattern Recognition and Machine learning by Bishop. However, I have a confusion related to a derivation ![enter image description here](http://i.stack.imgur.com/vJCY5.png) In the book when they calculate the marginal posterior of the hidden variables, I am ok up to. But in the page attached, I am confused what they mean by conditional probability of $z_{nk} = 1$ and they have $\gamma(z_{nk}) = \sum_{z}\gamma(z) z_{nk}$ what does z signify here? Lets say I have an example with z=[z1...zN] does that mean expectation over all possible combination of z lets say I have 3 possible values for each latent variable. So it means N^3 is it? I am confused any clarifications? Also I am confused about the maximization part During maximization they calculate the $\pi_{k}$ as $\pi_{k} = \frac{\gamma(z_{1k})}{\sum_{j=1}^{K} \gamma{z_{(1j)}}}$. Isn't the sum of the lower part equal to 1? So it doesn't change anything. Suggestions? ![enter image description here](http://i.stack.imgur.com/EQ7lm.png)


title: What is an ad hoc test?

text: What is an ad hoc test? What is opposite to an ad hoc test? How is it related to post hoc test? I found it in Goldfeld–Quandt test > The parametric Goldfeld–Quandt test offers a simple and intuitive diagnostic > for heteroskedastic errors in a univariate or multivariate regression model. > However some disadvantages arise under certain specifications or in > comparison to other diagnostics, namely the Breusch–Pagan test, as the > Goldfeld–Quandt test is somewhat of **an ad hoc test**. Primarily, the > Goldfeld–Quandt test requires that data be ordered along a known explanatory > variable. The parametric test orders along this explanatory variable from > lowest to highest. If the error structure depends on an unknown variable or > an unobserved variable the Goldfeld–Quandt test provides little guidance. > Also, error variance must be a monotonic function of the specified > explanatory variable. For example, when faced with a quadratic function > mapping the explanatory variable to error variance the Goldfeld–Quandt test > may improperly accept the null hypothesis of homoskedastic errors.


title: What is the formula of correlation coefficient against 3 or more variables?

text: This is the formula for calculating the coefficient of correlation between a single subject value ($z$) and two reference variables ($x$ and $y$). ![](http://i.imgur.com/Zv0ntSS.png) The question is how to expand it to cases when there are 3 or more variables in the reference set. Update: Adding some explanation as the notation was said to be unclear... $x$, $y$, $z$ etc are variables (lines). $R$ is the correlation coefficient, $R$zx is the coefficient of correlation between $z$ and $x$ (defined in the definition 2 here). $R$z,xy is the coefficient of correlation of a dependent variable $z$ to two independent variables: $x$ and $y$ (defined in the definition 1 here where it is also said that it "may also be expanded to more than two independent variables" but is not explained how).


title: What would you use as your default model for comparison?

text: I have a question. Say you are given with like 1500 instances as your dataset, in which instance will be categorized into 1 of the 3 classes. You're supposed to come up with the best model from the best learning algorithm. Before that, you're supposed to have a default model for comparison in order to come up with some analysis, i.e. how is this best model better than the default model. What would be the default model for comparison in this case? Thank you. Your help is very much appreciated :)


title: Normalization on frequently updated dataset

text: There are some normalization types like rescaling, standart score or modified standart score. I can apply these algorithms to large dataset. If the dataset that i am working on getting frequently updated like insert, update or delete then for every update i have to process whole data set again. It is costly to process whole data for only 1 update. Lets to be more clear If i normalize age of user so any user age update, user delete or user insert can affect normalization parameters like max age or min age etc. So in real life how we should handle this more efficiently from performance viewpoint ?


title: Is support vectors algorithms dependent?

text: In SVC problem, given all the coefficients fixed (C, gamma, etc), is it possible to get different decision functions and support vectors with different optimization strategies?


title: Dynamic conditional correlation reported in the table in Stata

text: The image below is an example result table of a dynamic conditional correlation (`dcc`) Multivariate GARCH model. I cannot figure out how the correlation numbers are computed. From what I understand, the `dcc garch` model gives you a set of correlations of variances over time. So, my questions are:   1. What do the correlations presented in the table stand for? Is it the average of the dynamic correlations?   2. I have generated my own result of two equations' dynamic conditional correlations series. How do you report them in the table like it shows in a table in Stata for a dynamic conditional correlation Garch model? Am I right that we just take a t-test between two series of correlations? ![enter image description here](http://i.stack.imgur.com/Pbcub.png)