title: Select the most probable learning algorithm from the testing set and output of a classifier?

text: Let's say I have a "mystery classifier". It is a "black box", I don't know exactly what it is doing, I know neither the training set nor the learning algorithm that was used. What I do have is a testing set, and the responses of the "mystery classifier" on this testing set. From this information, what can I know about the type of learning algorithm that was used? More specifically: (1) Can I know something about the type of supervision that was used ? For example, if i have several possible sets of labels for my test data, can I know which set of label was most likely used?) \--> For example, if I assume that the "mystery classifier" is a logistic regression, can I use the same "cost function" as the one from logistic regression to quantify the "distance" between the outputs of one "candidate classifier" and the "mystery classifier"? (2) Can I know something about the type of dimension reduction (supervised/unsupervised) that was used before classification, if any? \--> Is it true that a model which uses supervised dimension reduction + logistic regression (with the same label) would perform better than an unsupervised dimension reduction + logistic regression, because if dimension reduction is supervised then all the relevant information for the task will be there? Then, can I infer that my "mystery classifier" included some supervised dimensionality reduction related to the task if it performs way better than a simple PCA + logistic regression, or that it included some supervised dimension reduction unrelated to the task if it performs way worse? I assume there is the problem that algorithm can be better or worse for a lot of other reasons... can I know?


title: How to compute the total entropy of a system of binary outcomes

text: This is nowhere near my field of expertise, so if my question is poorly formatted for the context please feel free to edit. My question is at the end of the text below. Let's assume I flip a coin 5 times. It has an entropy of 5 bits, irrespective of if the coin is weighted. Equally if I roll a die 5 times the entropy would be log2(6^5) = 12.92 bits (would like clarification on this if possible - this is mostly my conjecture). Now lets assume the weighting of the die is based on the outcome of the coin flips. So would the following tell me the total entropy based on 5 coin flips and 3 die rolls: log2(6^(3xlog2(2^5))) = 248.156 bits So, more concisely how does one calculate the entropy (in bits) of a system where a second random variable is based on the results of a first?


title: Controlling for Type 1 Errors: Post-hoc testing on more than 1 ANOVA

text: I have conducted several repeated-measures ANOVA's with post-hoc testing via Tukey's HSD whenever the omnibus F-test was significant. I'm having some confusion about controlling for type 1 errors when presenting these post hoc tests. Q1: It's my understanding that Tukey's HSD will control the Type-1 error rate only for the set of post-hoc comparisons coming from one ANOVA. Is this correct? If so, do I need to take additional precautions to limit Type-1 errors with post-hoc comparisons from more than one ANOVA? Q2: If I do need to take further action, what should I do? Suppose I have 2 ANOVA's with 6 post-hoc comparisons each. If I did a Bonferonni correction, would I use the number of post-hoc tests (so, 12) or the number of ANOVA's (2), and would I still use Tukey's HSD at all? Is there a completely different approach/control I should consider?


title: Multivariate proportional data

text: I am looking for literature on what I call multivariate proportional data where a single observation is a vector of proportions that sum to 1. For example, each person weights their preferences for each of a fixed number of categories such that the total weight is fixed. More generally this could be multivariate data with a linear constraint on the Ys. Is this the proper term for such data and can anyone recommend references on the topic? Search engines (naturally) key in on the case where the response is a proportion, which is not quite what I am looking for.


title: Hypothesis testing to determine cluster outliers

text: I have a cluster of $p$-dimensional data from $n$ samples which is assumed to be normally distributed as a multivariate Gaussian with sample mean ${\bar{\mu}}$ and sample covariance matrix ${\bar{S}}$. Given a new sample $x_i$ I want to test whether or not the sample belongs to this cluster. I am doing the following: My null hypothesis, $H_o$, is: $x_i$ **does** belong to the cluster defined by the multivariate Gaussian with sample mean ${\bar\mu}$ and sample covariance $\bar{S}$ My alternate hypothesis, $H_a$, is: $x_i$ **does not** belong to the cluster My test statistic, $T(\bar{\mu},\bar{S})$, is defined as the Mahalanobis distance which, according to this post and the paper Gaussian mixture modeling by exploiting the Mahalanobis distance (page 13, Appendix), is distributed as a scaled ${\beta}$ distribution: $$ \frac{n}{(n-1)^2}T(\bar{\mu},\bar{S}) \sim \beta\left(\frac{p}{2}, \frac{(n-p-1)}{2}\right) $$ If I want to use a significance level ${\alpha} = 0.05$ then I believe I need to calculate my $p-value$ which, in this example, shows the probability of getting a sample as extreme (far from this distribution) or worse (further) from the underlying multivariate Gaussian as $x_i$. With this $p-value$ I will then check: **if** $p-value < {\alpha}$ : **then** reject $H_o$ i.e. this sample does not belong to this cluster **else** accept $H_o$ i.e. this sample does belong to the cluster Is this the correct way to perform this test? And if so, how do I calculate the $p$-value? **EDIT:(Adding more information and idea for solution)** It is my understanding that the $p-value$ represents the probability of obtaining a Mahalanobis distance ($d$) as extreme or worse than $d(x_i)$. This is calculated as $P(X{\geq}x_i)$, which is $1-P(X{\leq}x_i)$ which is $1 - $ the $cdf$ of the Beta function $\beta\left(\frac{p}{2}, \frac{(n-p-1)}{2}\right)$ So the $p-value = 1 - I_{x_i}\left(\frac{p}{2}, \frac{(n-p-1)}{2}\right)$ where $I_{x_i}\left(\frac{p}{2}, \frac{(n-p-1)}{2}\right)$ regularized incomplete beta function. Given that this is a two-tailed hypothesis test, with confidence level of $95%$ My test should be: $p-value = 2 {\times} \left(1 - I_{x_i}\left(\frac{p}{2}, \frac{(n-p-1)}{2}\right)\right)$ **if** $p-value {\leq} 0.05$ **then** : reject the null hypothesis (i.e. declare the sample as an outlier) **else** : accept the null hypothesis (i.e. sample is not an outlier) I think this is correct, but could someone more knowledgeable please confirm?


title: Calculating odds for a cricket match

text: I'm writing a program to help me gamble on cricket, but have run into a wee problem. I've got loads of data from previous games already and have compiled various stats from it. I have calculated that after 1 over of a 20 over match, the average score is 4.5 runs, and the average final score will be 145. Assume standard deviation is 2 and 20 respectively. So from there I can get a fairly decent idea of the odds for number of runs after 1 over, and at the end of the innings. (Assuming normal distributions N(4.5, 2) and N(145, 20)) My problem is, suppose after the first ball, they score 1 runs. How do I calculate the new odds? At first I thought I'd take 4.5/6 (number of balls in an over) = 0.75 as an "expected per ball" value and then say, 1 run scored is 1.33 times more than expected, therefore scale expected final score to 145*1.33 ~= 200 , and over score to 4.5*1.33 = 6. However this is going to lead to all sorts of problems: Firstly if they score 0 off the first ball, they would be expected to score 0 overall, which can't be right. Secondly, scaling 145 by 1.33 makes almost 200 runs, this jump is in no way justifiable after seeing just one ball. So I'm a bit stuck, not sure where to go from here. Sorry if this is a bit complicated, this one may be for cricket enthusiasts only!


title: Measurement error in dependent variable?

text: Suppose the cross-sectional model that I wish to estimate via OLS is $Y_{i} = a + b X_{i} + e_{i}$. Say that $Y_i$, $X_i$ and $e_i$ are normally distributed with mean 0 and constant variance. This specification adheres to all Gauss- Markov assumptions. But we don't observe $Y_i$. We instead observe $\hat{Y}_i = Y_i - \nu_i$. So I estimate $\hat{Y}_{i} = a + b X_{i} + e_{i}$, or equivalently $Y_i =a+bX_i+(e_i+\nu_i)$, or equivalently $Y_i=a+bX_i+\epsilon_i$, $\epsilon_i = e_i + \nu_i$. What happens in the following two cases: _**1)_** $\nu_i \sim N(q,\sigma_\nu)$, where $q$ is some positive constant. _**2)_** $\nu_i$ has mean 0 and variance $\sigma_\nu \forall i$ but is skewed (say, positively). _**3)_** $\nu_i \sim N(q_i,\sigma_\nu)$, where $q$ is a random variable whose value is conditional upon observation $i$. My inclination for 1) is that the $\hat{a} = a - q$. Not sure of 2) or 3) (preliminary ideas only). When I ask "what happens", I'm asking about (i) bias/consistency/efficiency of coefficient estimates (ii) accuracy of coefficient standard errors (iii) anything else. Undeveloped comments/suggestions are welcome as well!


title: What is the relationship between R-squared and p-value in a regression?

text: tl;dr - for OLS regression, does a higher R-squared also imply a higher P-value? Specifically for a single explanatory variable (Y = a + bX + e) but would also be interested to know for n multiple explanatory variables (Y = a + b1X + ... bnX + e). Context - I'm performing OLS regression on a range of variables and am trying to develop the best explanatory functional form by producing a table containing the R-squared values between the linear, logarithmic, etc., transformations of each explanatory (independent) variable and the response (dependent) variable. This looks a bit like: Variable name --linear form-- --ln(variable) --exp(variable)-- ...etc Variable 1 ------- R-squared ----R-squared ----R-squared --   ...etc... I'm wondering if R-squared is appropriate or if P-values would be better. Presumably there is some relationship, as a more significant relationship would imply higher explanatory power, but not sure if that is true in a rigorous way.


title: Cumulative counts or counts for Poisson regression

text: I have a set of data with measurements X1 and X2 across multiple time points, T1, T2 and T3. I would like to conduct a Poisson regression using X1 and X2 on the counts of a phenomenon. An example of how my dataset would like is as follows (pardon my inelegant code):               library(lme4)     set.seed(888)     df <- data.frame(SubID = sort(rep(1:10, 3)), Time = rep(1:3, 10), Count = sample(0:9, 30, replace = TRUE))     df <- data.frame(df, Cumulative.count = unlist(lapply(1:10, function(n) cumsum(df[cbind(c((1:10*3)-2), c(1:10*3))[n,1]:cbind(c((1:10*3)-2), c(1:10*3))[n,2],"Count"]))), x1 = rnorm(30), x2 = rnorm(30), x3 = rnorm(30))      As can be seen, I have 2 "counts" variable. One is cumulative, one is not. My question is should I be doing:               model.count <- glmer(Count ~ x1 + x2 + x3 + (Time|SubID), family = poisson(link = "log"), data = df)      OR:               model.cumulative.count <- glmer(Cumulative.count ~ x1 + x2 + x3 + (Time|SubID), family = poisson(link = "log"), data = df)


title: Cental limit theorem for average of iid variables

text: The Wikipedia entry on the CLT states at one point: "For fixed large $n$ one can also say that the distribution of $S_n$ is close to the normal distribution with mean $\mu$ and variance $\frac1n\sigma^2$." $S_n = (\sum_{i=1}^nX_n) / n$ and $X_i$ are iids with mean $\mu$ and variance $\sigma^2$. I don't quite see how this follows from the other, more formal definitions. Is this statement true and what is the source/proof?


title: How to implement credible 95% interval for median odds ratio using JAGS?

text: As described in Merlo et al (J Epidem Comm Health 2006), the 95% credible interval for MOR is calculated using MCMC. MOR is defined as $\exp(\sqrt{2\sigma^2}\times 0.675)$, where $\sigma$ is the level-2 variance of the random intercept $u$ from a null model of a hierarchical logistic regression. Does anyone have an idea of how to write a program for an Markov chain Monte Carlo to calculate the standard error of the median odds ratio (MOR) using rjags?   My dependent variable is outcome(alive/dead) and the clustering (level2)variable is Hospital. There are 140 hospitals and would like to see variations in outcome between hospitals. Other risk factors will be included later as independent level1 variables.


title: Simulating time-series given power and cross spectral densities

text: I am having trouble generating a set of stationary colored time-series, given the covariance matrix (their PSDs and CSDs). I know that, given two time-series $y_{I}(t)$ and $y_{J}(t)$, I can estimate their power spectral densities (PSDs) and cross spectral densities (CSDs) using many widely available routines, such as `psd()` and `csd()` functions in Matlab, etc. The PSDs and CSDs make up the covariance matrix: $$ \mathbf{C}(f) = \left( \begin{array}{cc} P_{II}(f) & P_{IJ}(f)\\\ P_{JJ}(f) & P_{JJ}(f) \end{array} \right)\;, $$ which is in general a function of frequency $f$. What happens if I want to do the reverse? **Given the covariance matrix, how do I generate a realisation of $y_{I}(t)$ and $y_{J}(t)$?** Please include any background theory, or point out any existing tools that do this (anything in Python would be great). ## My attempt Below is a description of what I have tried, and the problems I have noticed. It is a bit of a long read, and sorry if it contains terms that have been misused. If what is erroneous can be pointed out, that would be very helpful. But my question is the one in bold above.   1. The PSDs and CSDs can be written as the expectation value (or ensemble average) of the products of the Fourier transforms of the time-series. So, the covariance matrix can be written as: $$ \mathbf{C}(f) = \frac{2}{\tau} \langle \mathbf{Y}^{\dagger}(f) \mathbf{Y}(f) \rangle \;, $$ where $$ \mathbf{Y}(f) = \left( \begin{array}{cc} \tilde{y}_{I}(f) & \tilde{y}_{J}(f) \end{array} \right) \;. $$    2. A covariance matrix is a Hermitian matrix, having real eigenvalues that are either zero or positive. So, it can be decomposed into $$ \mathbf{C}(f) = \mathbf{X}(f) \boldsymbol\lambda^{\frac{1}{2}}(f) \: \mathbf{I} \: \boldsymbol\lambda^{\frac{1}{2}}(f) \mathbf{X}^{\dagger}(f) \;, $$ where $\lambda^{\frac{1}{2}}(f)$ is a diagonal matrix whose non-zero elements are the square-roots of $\mathbf{C}(f)$'s eigenvalues; $\mathbf{X}(f)$ is the matrix whose columns are the orthonormal eigenvectors of $\mathbf{C}(f)$; $\mathbf{I}$ is the identity matrix.    3. The identity matrix is written as $$ \mathbf{I} = \langle \mathbf{z}^{\dagger}(f) \mathbf{z}(f) \rangle \;, $$ where $$ \mathbf{z}(f) = \left( \begin{array}{cc} z_{I}(f) & z_{J}(f) \end{array} \right) \;, $$ and $\\{z_{i}(f)\\}_{i=I,J}$ are uncorrelated and complex frequency-series with zero mean and unit variance.    4. By using 3. in 2., and then compare with 1. The Fourier transforms of the time-series are: $$ \mathbf{Y}(f) = \sqrt{ \frac{\tau}{2} } \mathbf{z}(f) \: \boldsymbol\lambda^{\frac{1}{2}}(f) \: \mathbf{X}^{\dagger}(f) \;. $$   5. The time-series can then be obtained by using routines like the inverse fast Fourier transform. I have written a routine in Python for doing this:               def get_noise_freq_domain_CovarMatrix( comatrix , df , inittime , parityN , seed='none' , N_previous_draws=0 ) :         """                                                                                returns the noise time-series given their covariance matrix                                                           INPUT:                                                                                comatrix --- covariance matrix, Nts x Nts x Nf numpy array                                                              ( Nts = number of time-series. Nf number of positive and non-Nyquist frequencies )                              df --- frequency resolution         inittime --- initial time of the noise time-series                                                                    parityN --- is the length of the time-series 'Odd' or 'Even'                                                          seed --- seed for the random number generator                                                                         N_previous_draws --- number of random number draws to discard first                                                   OUPUT:                                                                                t --- time [s]                                                                                n --- noise time-series, Nts x N numpy array                                                                          """         if len( comatrix.shape ) != 3 :            raise InputError , 'Input Covariance matrices must be a 3-D numpy array!'         if comatrix.shape[0]  != comatrix.shape[1] :             raise InputError , 'Covariance matrix must be square at each frequency!'              Nts , Nf = comatrix.shape[0] , comatrix.shape[2]              if parityN == 'Odd' :             N = 2 * Nf + 1         elif parityN == 'Even' :             N = 2 * ( Nf + 1 )         else :             raise InputError , "parityN must be either 'Odd' or 'Even'!"         stime = 1 / ( N*df )         t = inittime + stime * np.arange( N )              if seed == 'none' :             print 'Not setting the seed for np.random.standard_normal()'             pass         elif seed == 'random' :             np.random.seed( None )         else :             np.random.seed( int( seed ) )         print N_previous_draws         np.random.standard_normal( N_previous_draws ) ;              zs = np.array( [ ( np.random.standard_normal((Nf,)) + 1j * np.random.standard_normal((Nf,)) ) / np.sqrt(2)                      for i in range( Nts ) ] )              ntilde_p = np.zeros( ( Nts , Nf ) , dtype=complex )         for k in range( Nf ) :             C = comatrix[ :,:,k ]             if not np.allclose( C , np.conj( np.transpose( C ) ) ) :                 print "Covariance matrix NOT Hermitian! Unphysical."             w , V = sp_linalg.eigh( C )             for m in range( w.shape[0] ) :                 w[m] = np.real( w[m] )                 if np.abs(w[m]) / np.max(w) < 1e-10 :                     w[m] = 0                 if w[m] < 0 :                     print 'Negative eigenvalue! Simulating unpysical signal...'                  ntilde_p[ :,k ] =  np.conj( np.sqrt( N / (2*stime) ) * np.dot( V , np.dot( np.sqrt( np.diag( w ) ) , zs[ :,k ] ) ) )              zerofill = np.zeros( ( Nts , 1 ) )         if N % 2 == 0 :             ntilde = np.concatenate( ( zerofill , ntilde_p , zerofill , np.conj(np.fliplr(ntilde_p)) ) , axis = 1 )         else :             ntilde = np.concatenate( ( zerofill , ntilde_p , np.conj(np.fliplr(ntilde_p)) ) , axis = 1 )         n = np.real( sp.ifft( ntilde , axis = 1 ) )         return t , n      I have applied this routine to PSDs and CSDs, the analytical expressions of which have been obtained from the modeling of some detector I'm working with. The important thing is that at all frequencies, they make up a covariance matrix (well at least they pass all those `if` statements in the routine). The covariance matrix is 3x3. The 3 time-series have been generated about 9000 times, and the estimated PSDs and CSDs, averaged over all these realisations are plotted below with the analytical ones. While the overall shapes agree, there are noticeable _noisy features_ at certain frequencies in the CSDs (Fig.2). After a close-up around the peaks in the PSDs (Fig.3), I noticed that the PSDs are actually _underestimated_ , and that the noisy features in the CSDs occur at just about the same frequencies as the peaks in the PSDs. I do not think that this is a coincidence, and that somehow power is leaking from the PSDs into the CSDs. I would have expected the curves to lie on top of one another, with this many realisations of the data. ![Figure 1: P11](http://i.stack.imgur.com/nPisE.png)   ![Figure 2: P12](http://i.stack.imgur.com/axOf2.png) ![Figure 2: P11 \(close- up\)](http://i.stack.imgur.com/v48Zi.png)


title: T-test / ANOVA on Box-Cox transformed non-normal data

text: Suppose I apply a Box-Cox transformation to my data and now it looks rather like a normal distribution. I then add another dataset, transform it by Box- Cox with the same lambda and run a t-test to compare the means. Would this approach make sense if my data is non-normal by its nature? In other words, is the fact that a Box-Cox transform produces a Gaussian-like distribution sufficient to then use standard methods for normally distributed data such as t-test and ANOVA? **Update** \- to formulate this question a bit more specifically: I want to test whether there are significant differences between the means of two samples. I can see that the distributions in each sample are very much non- normal. My question is: if I force them to look normal by using a transformation, will this be enough to essentially forget about their "original" non-normal nature for testing this hypothesis? **Update 2** \- I suppose my question is similar in spirit to this one, which asked the same thing about log-transformation.


title: Log, geometric and elasticity

text: Suppose I want to calculate the average annual growth in GDP per capita from year 1980 to 1988. I suppose this is done by ln(gdp per capita 1998) - ln(gdp per capita 1990) ? So for example: ln(900)-ln(100) / 8 = 0.2747 , so 27.47%. But if I take the start value times the growth increase i get: 100 x 1.2747^8 = 697 But if I use the geometric growth average (compound annual growth) to calculate growth I get: 8^√(900/100) = 1.316 , so 31.6%. and taking the start value times growth increase i get: 100 x 1.316^8 = 899.6 , so almost exact 900. Does the log difference method underestimate the growth increase? Why? Which is best to use? And if I have calculated the growth rate with the geometric approach and regress it against another percentage change in say increase of cars (growth of cars bought), would such a relationship be an elasticity between average growth in gdp per capita and average growth of cars?


title: Is delta (Δ) R² the same as adjusted R²?

text: Is delta (Δ) R² the same as adjusted R² ? I have delta R² in my notes somewhere and it didn't come up on a Google Search.


title: Minimum no. of observations required for statistical distribution fitting!

text: I have one conceptual (but sort of vague) question regarding distribution fitting. How many observations one would need to best fit any statistical distribution to given data. Like I am dealing with loss data and in one situation, I had only three observations still I could fit Normal distribution, which in reality I think is absurd. One more data and the scenario may change. Hence, is there any thumb rule for minimum number of observations in order to have reasonably stable distribution fitting? Thanking in advance Akshata


title: Interaction terms in poisson regression

text: I'm confused about the interpretation of interaction terms in poisson regressions. Here's a hypothetical dataset               Id  Group  Condition   Error rate     1   A      X           1     1   A      Y           2     2   B      X           3     2   B      Y           6      Imagine that individuals 1 and 2 are duplicated a number times (with random noise) to make two groups. If we analyse this data using a poisson regression (as it's count data) there will be no Group by Condition interaction (in both groups, error rates increase by 100% in condition Y, i.e. an Incidence Rate Ratio of 2). If we analyse this data using an ANOVA (and don't transform data) it is likely to yield a significant Group by Condition interaction (I've simulated this). While the poisson regression looks at proportional differences across conditions, the ANOVA is more sensitive to absolute differences. It seems to me that because the poisson regression uses a multiplicative scale (via log transformation), an additive effect under the poisson regression will become multiplicative when the same data is analysed using OLS (or ANOVA). So which interaction term is more "truthful"? Is it just a matter of interpretation? Or are there any a priori reasons for choosing one over the other. For example, could you argue that because count data exhibits ratio properties we should be analysing ratios and not intervals (according to Stevens)? Thanks and hope this makes sense.


title: Reference on factor analysis with categorical variables

text: I have been studying so many things on factor analysis with categorical variables. I am frustrated studying these many pdfs. I have 40 variables obtained from 40 questions. All of them are categorical. I can take them to be ordinal. The questions are different, so the answers are different too. It is not something like Likert scale where, for example, 1 means good for all questions, 2 means moderate for all questions, 3 means bad for all questions etc. It is like for one question:               "How often do senior management visit the wards to talk to staff?"          rarely or never ..................... 1         around once a year................... 2     around once a month.................. 3     around once a week................... 4      For another question:               "What is the average amount of training (per person) received by a management staff?"          Less than a day ..................... 1     Less than a week .................... 2     One to two weeks .................... 3      Etc. I want to factor analyze these variables: as can be seen here, 1, 2, 3 etc. differ in their meaning, and also the number of categories for each question differ. Another problem is that I have non-response (missing values) in the data. Here are my questions:   1. What is the actual and best method of factor analysis with this kind of data. Besides could you please give me one good reference with your suggested methodology? I will be grateful for your help.    2. If possible, please also give me indication on what to do with the missing values.   3. I need to calculate factor scores from the analysis. I have tried polychoric correlation, but can't get factor scores by this. Factor scores are very important for my analysis. I can't make further analysis without them.  * * * @this.is.not.a.nick: Thank you so much for your kind advice. I was suggested to use CATPCA too. But if using polychoric correlation can solve the problem of calculating factor scores then it's really great. But Andrea, chl and ttnphns, could you please confirm me whether `fa.poly()` or `fa.parallel.poly()` use principal component solution of the loadings and the specific factor or they just use maximum likelihood method for estimation of these parameters like `factanal()` command does? As I can not assume the data to be normally distributed so I guess using principal component method for estimating parameters will be good in this case. If these functions don't use principal component solution then I think I can do one thing here-   1. Calculate the polychoric correlation matrix r using package psych.   2. Compute loadings              f <- principal(r,nfactors=3,rotate="varimax",scores=T,residuals=T) #say, 3 factors taken     l <- print(f$loadings[c(1:ncol(data)),],cutoff=.0001) #data means original data        3. Compute scores              h <- t(l)%*%l #communality     s <- h%*%t(l) #as fhat_i=(L'L)^(-1)*L'*Z_i     data1 <- t(data[1,])     f1 <- s%*%data1      f1      And thus I can get $f_i$, $i=1,\dots,n$, manually from Bartlett's formula for scores for PC solution. (Reference: Applied Multivariate Statistical Analysis by Johnson and Wichern) Now, for this procedure to be implemented there should not be any missing value (NA). So, if I take the missing values to be 0 by which for each question means "no comment" is there any problem? I think it will then act as just another category of my categorical variables. As the variables are categorical I think I should not take mean or median. Am I right in my thoughts? Please suggest me. If possible please send me the paper regarding Polychoric versus Pearson correlations. I really had a lot of help from you and looking forward to contribute when I grow up.


title: Cluster on high dimensional categorical data (Images with keywords)

text: We're looking for clues to perform a Cluster Analysis in a DB with **+400K images** which have keywords associated to them.   Each image could have from 1 to 30 keywords.   Total **keywords count is +35K**.   Keywords are not assigned with uniform criteria, e.g. one dog image can have `[animal, canine, dog]` and other similar dog picture can have `[dog, friend]`. We will assume that keywords are strong related with images.   We are not expecting to classify all images, but we'll be happy to get about 50% of data in coherent groups.   We have used factor analysis with ordinal data to reduce variables and then perform cluster, but we are a little bit lost in this case.     1. How could we reduce the number of variables (keywords)?   2. How to filter cases (images) that are less likely to be grouped?   3. Is this problem already (publicly) solved? Thanks for your hints.


title: Mixfdr penalization parameter

text: I'm trying to use mixfdr in a simulation of p-values but sometimes the search for the penalization parameter spits out an error and shows some matrix entries to be infinite:               >m <- mixFdr(stdjk,plots=F,J=3,theonull=F,calibrate=T)     Fitting preliminary model      [1] 0.750 0.125 0.125     [1] -7412.522 -7412.522 -7412.522     [1] 3706.261 3706.261 3706.261     [,1]      [,2]     [,3]     [1,] 0.750 -7412.522 3706.261     [2,] 0.125 -7412.522 3706.261     [3,] 0.125 -7412.522 3706.261     [,1] [,2]     [1,] -Inf  Inf     [2,] -Inf  Inf     [3,] -Inf  Inf      [,1] [,2]     [1,]  0.1  Inf     [2,]  0.1  Inf     [3,]  0.1  Inf     Error in if (distance < tol) converged = TRUE :      missing value where TRUE/FALSE needed      Is this normal or am i missing some parameter?


title: LOESS and MA normalization in R?

text: Attempting to do loess on two variables `x` and `y` in R using MA normalization(see MA-plot; see also Bland-Altman or Tukey mean-difference plot) like this:               > x = rnorm(100) + 5     > y = x + 0.6 + rnorm(100)*0.8     > m = log(x/y)     > a = 0.5*log(x*y)      I want to normalize x and y in such a way that the average `m` is 0, as in standard MA normalization, and then back-calculate the correct x and y values. First running loess on MA:               > l = loess(m ~ a)      What is the way to get corrected `m` values then? Is this correct?               > mc <- predict(l, a)     # original MA plot     > plot(a,m)     # corrected MA plot     > plot(a,m-mc)      not clear to me what `predict` actually does in the case of `loess` objects and how it's different from using `l$residuals` in the object `l` returned by `loess` \- can someone explain? finally, how can I back calculate new `x` and `y` values based on this correction? Attempt to back-calculate `x` and `y`: the corrected `m` is,               new_m = m - mc      so `x,y` can be derived from the definition of `m`:               m = log(x/y) = log(x) - log(y)      therefore,               x = exp(new_m + log(y))     y = exp(-1*(new_m - log(x)))      but this is wrong; there's a missing scaling factor of a half and I'd like to see a derivation of where that comes from. Probably from the definition of `A` but I don't see why I can't get the same be rewriting `x,y` in terms of `A` **edit** if someone can explain where the half correction factor comes in when back calculating `x` \-- i.e. why the formula `x = exp(new_m + log(y))` is wrong -- i'd appreciate it.


title: Evaluating probability mass function for a discrete random variable

text: Consider the following pmf for a discrete random variable $X$: $f(x) = x/15$ when $x = 1, 2, 3, 4, 5$, and $f(x)= 0$ for all other values of $x$. What is $f(3)$? (I.e., the PMF evaluated at 3)


title: Interaction model significant in multiple linear regression

text: This question is not about interaction effect, but about interaction model in Hierarchical Linear Regression. I have 1 DV and 5 IVs. I want to see which of the IVs is significant predictor of DV. Although I am not controlling for any moderator variable here, but to do multiple regression I had to use hierarchical style where I put IVs in Block 1, and interaction terms (10 in number) in Block 2 of SPSS. I used non-centered DV but had centered all the IVs before forming their interaction terms. Now lets say I have IV1, IV2, IV3, IV4, and IV5 as the independent variables. In Model 1, IV1, IV2, IV3 are coming significant and Model 1 is significant too. In Model 2, IV3 and IV4 are coming significant while interaction terms are not significant, and Model itself is again significant. My question is that what do I have to do now? Which IVs should be considered significant predictors (I think Model 2 predictors)? Why did interaction terms altered the whole picture? By the way, this is first part of the whole statistical procedure. In second part, I have a moderator variable, M1. And I put it in the Block 1, IVs in Block 2, and interaction terms (this time 15 in number, IVs with M1 too) in Block 3. Here everything is fine. In Model 1, M1 is significant. In Model 2 with IVs and control predictor, only IV2 is significant, and in the Model 3, no interaction term was significant but IV4 was significant, IV2 became insignificant and Model itself was not significant. Therefore, since model is insignificant, so significance of IV4 does not matter. Similarly, does it mean that in first step of regression, without M1, second model should be accepted as it is significant? Is it ok if model with interaction terms comes significant but interaction terms are insignificant? From what I have been able to understand so far, interaction terms need to be insignificant but the model in which they are entered can be significant. Am I correct here?


title: How to integrate principal components with GLM?

text: How would I integrate the output of a principal components analysis with a GLM (assuming the PCA is used for variable selection for the GLM)?


title: Mixed linear models in R, help with nested terms and procedure

text: This is my first post, so sorry if it not optimally written. I have a paired samples at two time points in two groups, undergoing the same intervention. I want to test the effect of my intervention on weight. I'm using R. Here's some data to illustrate: my data frame called `df`:                     patients   timepoint        group          Weight             102            1            1            107.30             104            1            1             94.10             117            1            1            110.80             121            1            1            108.90             153            1            1             95.40             155            1            1            105.10             161            1            1             97.70             162            1            1             83.60             167            1            1             82.40             173            1            1             86.40             176            1            1             81.90             177            1            1             90.90             179            1            1             95.30             101            1            2             81.30             108            1            2             72.30             113            1            2             68.50             170            1            2             89.20             171            1            2             77.50             172            1            2             94.50             175            1            2             78.30             181            1            2             71.40             182            1            2             72.80             183            1            2             73.50             186            1            2             87.90             187            1            2             83.50             188            1            2             70.10             102            2            1            102.70             104            2            1             90.40             117            2            1            107.30             121            2            1            107.50             153            2            1             95.00             155            2            1            102.80             161            2            1             95.40             162            2            1             78.30             167            2            1             81.90             173            2            1             85.30             176            2            1             83.10             177            2            1             90.50             179            2            1             97.50             101            2            2             78.40             108            2            2             72.00             113            2            2             66.80             170            2            2             90.20             171            2            2             77.60             172            2            2             93.40             175            2            2             80.30             181            2            2             72.60             182            2            2             71.40             183            2            2             74.20             186            2            2             88.70             187            2            2             80.50             188            2            2             71.20      Since this is paired data (between time points) and unpaired (between groups), I guess I must use a mixed linear model. Im going for the lme4 package in R. "timepoints" and "group" will be my fixed effects (I exhaust both). "patients" will be my random effect, which also picks up that I have multiple responses per patient. I will use a random intercept model, but I actually expect that my patients differ with how they react to my experimental manipulation, so a random slopes model would be nice also. However, it seams I over-parametrize the model if doing so. I will use the likelighood ratio using anova() for a full model and a reduced model.               full = lmer(Weight ~ timepoint + group + (1|patients), data=df,             REML=FALSE)          reduced = lmer(Weight ~ group + (1|patients), data=df,                        REML=FALSE)      "timepoint" is the main factor in question. I now test using anova():               > anova(reduced,full)     Data: df     Models:     reduced: Weight ~ group + (1 | patients)     full: Weight ~ timepoint + group + (1 | patients)             Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)       reduced  4 309.72 317.52 -150.86   301.72                               full     5 306.02 315.78 -148.01   296.02 5.695      1    0.01701 *      Question is, Im I doing it correctly? And how do I interpret the result? BUT, what I really want is to test if the effect of time is different on the two groups. How do I do this? Thank you.


title: How to organize a series of data analysis results?

text: For several days, I collected a series of data analysis results for given data with different statistical methods, different sample size, different other analysis parameters, and so on. So right now, I have 100 or so analysis results sorted by test data error rate. Is there any way to group top of the results? As a data mining guy, it's temptating to run a clustering method to find a top cloud, or a regression tree to find important parameters to get top results. I would like to find a statistically valid method on this. **UPDATE:** Responding to Michelle's request on more details, my data is images and my analysis is to recognize certain object in the image that has something to do with a disease. But I doubt this information would help the answers to this question. Please correct me if I am wrong.


title: Does the link function in a GLM influence the plot of residuals?

text: Could an inappropriate link function contribute to non-uniform residuals in a GLM, or could that only be caused by an inappropriate error structure?


title: Basic Bayesian reasoning question -- what is wrong with this "equational" argument?

text: Let's say that we have a model: $y + x = \epsilon$, $\epsilon \sim N(0, 1)$. After observing a value for $y$, we can write $x$ as: $x = \epsilon - y$ Since $y$ is just a constant, and $\epsilon$ is standard normally distributed, it seems that we can claim: $\epsilon - y \sim N(-y, 1)$ and thus: $p(x \mid y) = N(x; -y, 1)$ where $N$ is the normal pdf. But this makes no sense as there is no mention of a prior for $x$ anywhere. So what went wrong with this "equational" reasoning argument?


title: Stratified sampling with multiple variables?

text: I don't know much about stats so I'm looking for a starting point here. Any resources or insights would be helpful. I'm conducting an e-learning experiment, in which students watch videos and then complete a survey which measures cognitive load and user satisfaction, and afterwards they complete a short assessment to test what they have learned. Because of the issues with pre-testing in learning, I want to stratify the experiment groups by three vairables:   1. Results of a previous examination   2. Attitude towards their course (measured by survey)   3. Attitude towards e-learning (measured by survey) By doing this I can minimise variance between each group. If I was stratifying by variable 1, I could ensure that each group had an equal amount of 'A' students, 'B' students etc. However, I'm not sure how to fairly divide students into groups given that I have three variables to control for. I could just play around with the groups until I get them roughly even, but I was wondering if there are statistical methods for stratifying with multiple variables? Thanks.


title: Which distribution to parametrize price histogram with exponential tail?

text: I have plotted price histograms which seem to follow and exponential tail very well. ![](http://i60.tinypic.com/24dr3vr.png) However, I would also like to parametrize the initial peak in the low price range. Can suggest what distribution would be appropriate for this? How can I research this in general?


title: Best way to visualize attrition using R?

text: Thru this site I've recently discovered Sankey Diagrams, a great way to visualize what is happening in a _traditional_ flow chart. Here is a good example of a Sankey Diagram by George M. Whitesides and George W. Crabtree, ![Don't Forget Long-Term Fundamental Research in Energy](http://i.stack.imgur.com/fKl53.png) Source; _Don't Forget Long-Term Fundamental Research in Energy_ , Science 9 February 2007:Vol. 315. no. 5813, pp. 796 - 798. After I realized that there was no Sankey R-package I found an R script online, unfortunately this script is quite raw and somewhat limited. With high hopes I asked for a Sankey R-package or a more mature function at stackoverflow, but to my surprise it seems as we do not have a mature function for building Sankey Diagrams in R. After I posted a bounty Geek On Acid was kind enough to suggest a small hack on the existing script which made it work more or less for my specific purpose. The improved R-script produced this diagram, ![Geek On Acid's R-Sankey Diagram](http://i.stack.imgur.com/lFahc.png) Source; stackoverflow.com. But, does the lack of a R package indicate that Sankey Diagrams isn't such an amazing way to visualize attrition using R in a data flow à la the one presented in the diagram above (see initial stackoverflow question for data and R code. Maybe there's a better way to visualize attrition. What do you think is the best way to visualize attrition in a data flow using R?


title: Chi-square with unbalanced design

text: I am a Ph.D. candidate of Applied Linguistics. Currently I am involved in a research study about English teachers' provision of corrective feedback on English L2 learners' errors. I have divided the teacher participants into two groups. I want to juxtapose the groups in terms of the number (frequency) of corrections they supply in response to the learners' errors. I suppose I should apply chi-square test. However, as the group sizes are different (10 teachers in group A and 15 teachers in group B), I am sort of uncertain about the applicability of chi-square test. I want to know if there is any statistical test that is capable of dealing with frequency data while not impacted by the lack of equality in group size.


title: Are degrees of freedom $n-1$ for both the sample standard deviation of the individual observations and for the standard error of the sample mean?

text: I collect $n$ ($<20$) i.i.d. observations from any distribution. In order to compute the sample variance, I take $$s^2=\sum_i \frac{(\bar{X}-X_i)^2}{n-1}$$ If I want to build a confidence interval for the mean, I need the standard error of the sample mean. Do I divide by $n$, or by $n-1$ again? $$S.E. \stackrel{?}{=}\sqrt{\frac{s^2}{n-1}}$$


title: Interpret Pearsons correlation output from R

text: I have conducted Pearson's correlation in R and need help interpreting the results.               Pearson's product-moment correlation          data:  A.C$Average.tortuosity and A.C$Area     t = -0.6168, df = 14, p-value = 0.5473     alternative hypothesis: true correlation is not equal to 0     95 percent confidence interval:      -0.6092369  0.3622606     sample estimates:            cor      -0.1626531      I am interested in reporting the _r_ and _p_ values. Would it be fair to say based on the the above output that: > There was a non-significant correlation between the two variables, _r_ (14) > = -.16, _p_ > .05.


title: How to test if all the levels have the same success/fail probability?

text: I have a csv data set like this.               Levels, Trial, Fail     1,100,1     2,200,3     3,300,4      How do I test if the probability of fail or success is the same across all the levels?


title: How to calculate joint probabilities from conditional probabilities in a Bayesian Network?

text: This question is about Bayesian Networks. I want to calculate the probability for certain events to be in a certain state knowing all conditional probabilities. Consider that I am totally new to Bayesian Networks. # The problem Consider a Bayesian Network. For all variables $A_i,B_i \dots Z_i$ in the network I know the potentials in the form $P\\{ A_i | par(A_i) \\}=P\\{A_i | B_1 \dots B_n\\}$ for all possible states $a_i,b_1 \dots b_n$. These potentials are calculated using subjective probabilities. # Calculating single probabilities What I know are just potentials. For example, for variable $A_i$ I do not know $P\\{A_i = a_i\\}$ but I know $P\\{A_i = a_i | B_1=b_1 \dots B_n=b_n\\}$. Where $B_1, B_2 \dots B_n$ are parents of $A_i$. I want to calculate $P\\{A_i = a_i\\}$. Consider that No restrictions are considered so this variable $A_i$ can have parents and children. As well as $B_1 \dots B_n$. It does not seem to complicated. Having potentials I can use the _total probability law_ and: $$P\\{ A_i = a_i \\} = \sum_{\forall b_1 \dots b_n} P\\{ A_i = a_i | B_1=b_1 \dots B_n=b_n \\} \cdot P\\{B_1=b_1 \dots B_n=b_n \\}$$ In the sum, the first part: $P\\{ A_i = a_i | B_1=b_1 \dots B_n=b_n \\}$ can be calculated using potentials. What about $P\\{B_1=b_1 \dots B_n=b_n \\}$. That is, How can I calculate the probability that some variables have some state contemporarly? Please note that $B_1 \dots B_n$ are not one parent of the other, they are on the same level of the network, they are peers, meaning they have no parent relationship. They are parents of $A_i$. # The question Generally speaking, how to calculate the probability that some variables are in some state contemporarly? When they have no parent relationships? Preconditions are: I know the conditional probabilites (potentials).


title: "Incorrect" p-value during probit regression when probability of a response is zero

text: I am running a probit regression with glm. I got a non-significant value (p~0.98) for an interaction that other statistical methods suggested should be significant. Looking into it, I think it's because 1 of the cells had zero positive ("YES") outcomes. When I artificially change a data point so that there's one positive outcome, I get a significant result of p < 0.02, even though this actually decreases the difference in the effect of the cue ("Cue") between the two conditions ("Condition"). The code in R is `glm(Response~Cue*Condition, family=binomial(link=probit), data=data)` The data look like this:               Condition = 0             Response:    0   1         Cue: 0            49   3              1            25  25            Condition = 1               Response:    0   1         Cue: 0            48   3                   1            47   0       The counts aren't quite exact but it shows the pattern, which is that the subject responds "appropriately" 50% of the time to Cue 1 in Condition 0 but does not do so in Condition 1 (and in fact, responds slightly more often in this condition to the incorrect cue, although this isn't significant). Changing the last row of counts to 46 1 gives the significant interaction I was expecting. My intuitions are hazy, but I imagine this has something to do with the probit function being undefined for probabilities of 0 or 1. Does this seem correct? If so, why does it return a non-significant value rather than a warning? More to the point, is there a way to deal with this? I can imagine building in a correction that gives a non-zero estimate for the probability based on the number of trials, but is there a standard, implemented solution?


title: Testing for a significant difference between ML estimates: Likelihood ratio or Wald test?

text: I am trying to test whether or not there is a significant difference between maximum likelihood estimates of two genetic parameters (selection and dominance) across two environments with genotype data from a cross. There are three genotype classes, C11, C12, and C22, and I am looking to maximize the likelihood below:               llhood2 <- function(theta, c, c11,c12,c22){          s<-theta[1]     d<-theta[2]          P11 = (1 - ((s * (1 - c) * (1 - c + (2 * d * c)))) / (4 - (s * (1 + (2 * d))))     P22 = (1 - (s * c * ((2 * d * (1 - c)) + c))) / (4 - (s * (1 + (2 * d))))     P12 = (1 - P11 - P22)     logL = c11*log(P11) + c12*log(P12) + c22*log(P22)     -logL        }      C is fixed depending on which genetic marker you are using, and c11, c12, and c22 are the data, representing genotypes AA, AB, and BB respectively. So I am jointly estimating s and d. Using the optim function in R and plugging in genotype data from each environment, I am able to get a point estimate for s and d in one environment (s1 and d1) and the variance of these estimates from the hessian, then I get parameter estimates using genotype data from the 2nd environment,(s2 and d2, and thier associated variances). So, two separate fittings of the model give me s1 var_s_1 and s2 var_s_2, etc. The Code for optim is below:               fit1 <- optim(inits, llhood, method="L-BFGS-B", hessian=T, lower=lower, upper=upper,     c=c,c11=c11,c12=c12,c22=c22, control=list(trace=1))          fit2 <- optim(inits, llhood, method="L-BFGS-B", hessian=T, lower=lower, upper=upper,     c=c,c11=n11,c12=n12,c22=n22, control=list(trace=1)      I then want to know if s1 is significantly different from s2 (likewise for d1 and d2). I thought this would be a straight forward test, but I have had a hard time finding documentation on how to do this properly, or most elegantly. As a crude approximation, I used a simple two-sample T test, but Im pretty sure that it is not quite right... Ive also considered a likelihood ratio test where I test whether the difference between the two estimates of s1 and s2, is not different from zero (reduced model) vs different from zero. But Im not sure how to set this up given my likelihood equations. Any specific help would be much appreciated. UPDATE: Here is some attempted code for the restricted and unrestricted likelihoods which will be used to perform a likelihood ratio test. First the restricted (null) where s and ss are equal- (c11-c22 are the data for env.1 , n11-n12 for environment 2)               llhoodnull <- function(theta,c,c11,c12,c22,n11,n12,n22){          s<-theta[1]      d<-theta[2]     ss=s     dd=d               P11 = (1 - (s * (1 - c) * (1 - c + (2 * d * c)))) / (4 - (s * (1 + (2 * d))))     P22 = (1 - (s * c * ((2 * d * (1 - c)) + c))) / (4 - (s * (1 + (2 * d))))     P12 = (1 - P11 - P22)     b11 = (1 - (ss * (1 - c) * (1 - c + (2 * dd * c)))) / (4 - (ss * (1 + (2 * dd))))     b22 = (1 - (ss * c * ((2 * dd * (1 - c)) + c))) / (4 - (ss * (1 + (2 * dd))))     b12 = (1 - b11 - b22)          logL = c11*log(P11) + c12*log(P12) + c22*log(P22)+n11*log(b11) + n12*log(b12) +          n22*log(b22)     -logL          }      I think that I may also need to fix the d so that I am only testing for the likelihood of s. And the unrestricted:               llhoodb <- function(theta,c,c11,c12,c22,n11,n12,n22){          s<-theta[1]      d<-theta[2]     ss<-theta[3]     dd<-theta[4]               P11 = (1 - (s * (1 - c) * (1 - c + (2 * d * c)))) / (4 - (s * (1 + (2 * d))))     P22 = (1 - (s * c * ((2 * d * (1 - c)) + c))) / (4 - (s * (1 + (2 * d))))     P12 = (1 - P11 - P22)     b11 = (1 - (ss * (1 - c) * (1 - c + (2 * dd * c)))) / (4 - (ss * (1 + (2 * dd))))     b22 = (1 - (ss * c * ((2 * dd * (1 - c)) + c))) / (4 - (ss * (1 + (2 * dd))))     b12 = (1 - b11 - b22)          logL = c11*log(P11) + c12*log(P12) + c22*log(P22)+n11*log(b11) + n12*log(b12) +          n22*log(b22)     -logL          }      This one gives me errors when I try to maximize it with optim()... Any thoughts on these?


title: Web Videos that can provide me with a better understanding of statistical concepts

text: I have an exam coming up and I would like to know are there youtube videos/playlist or any videos online that can explain the following concepts below in a simple manner. The concepts are: (Descriptive Statistics) Scales of measurement Histograms and other plots Empirical rule Summary statistics Percentiles (Probability) Expected values: Binomial distribution Poisson distribution Normal distribution Normal approximations to distributions Sampling Distributions Experimental and survey design (Inference) Mean from a single population Proportion from a single population Means from two independent samples Means from paired samples Proportions from two independent samples (Hypothesis Testing) Means from a single population Proportions from a single population Multinomial proportions from a single sample Two Means using independent samples Two Means using paired samples Two Proportions using independent samples Independence or homogeneity of proportions (same as chi-square tests) (Analysis of Variance:) CRD RCB (Trends, Regression, and Correlation) Trends Least squares Correlation Chi-square tests and contingency tables (Tests for independence) "I especially need help on Chi-square tests, analysis of variance, and Inference." Thanks Everyone for all your help!!


title: Marginal probability density of x (obs) obtained by integrating x(missing)

text: I am reading a text book on missing data, and a sentence below is slightly challenging for me to understand. > The marginal probability density of $ \left ( x_{obs}\ \right)$ is obtained > by integrating the missing data $ \left ( x_{miss}\ \right)$ $$ f\left ( > x_{obs}|\theta \right )=\int_{}^{}f\left ( x_{obs},x_{miss} \right > ){}dx_{miss}\\\\\\\\{}$$ I need help understanding this, with an example preferably.


title: Extraction of a decision boundary (LDA) after a systematic querying of the feature space and convolution with Sobel filter (examples in numpy)

text: I am doing some experiments with LDA ( _Linear Discriminant Analysis_ ), in python. Now I am at the point in which I would like to display the separation planes in the 3-dimensional feature space. I tough that a solution could be, initially, querying the feature space so that for each coordinate ("voxel"), the inferred category can be attached. This category is the _argmax_ of the posterior probabilities inferred for each category. Till here everything went fine, and the operation, although computationally expensive, was successful. At this point I would like a reliable method to extract a decision boundary from the queried feature space. In other words, I would like to get a 3D matrix filled with zeroes except for the occurrences of the separation plane, which should be ones. In order to use existing methods and not to re-invent the wheel, I decided to use the well-known Sobel filter in order to give some values at the points that belong to the decision boundary and zeroes elsewhere. I will illustrate what I am doing with a simple example, using _numpy_. Let's assume that the following forged matrix represents a 2D feature space that has been queried with a classifier:               >>> m = np.fliplr(np.triu(np.ones([10,10]))) + np.triu(np.ones([10,10]))     >>> m     array([[ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.],            [ 1.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  1.],            [ 1.,  1.,  2.,  2.,  2.,  2.,  2.,  2.,  1.,  1.],            [ 1.,  1.,  1.,  2.,  2.,  2.,  2.,  1.,  1.,  1.],            [ 1.,  1.,  1.,  1.,  2.,  2.,  1.,  1.,  1.,  1.],            [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.],            [ 1.,  1.,  1.,  1.,  0.,  0.,  1.,  1.,  1.,  1.],            [ 1.,  1.,  1.,  0.,  0.,  0.,  0.,  1.,  1.,  1.],            [ 1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.],            [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.]])      As you can see, the matrix `m` could well represent a bi-dimensional space partitioned in several categories by a classifier. Now, if I convolve that matrix with the Sobel filter, I obtain something that is promising:               >>> from scipy.ndimage.filters import sobel     >>> sobel(m)     array([[ 1.,  1.,  0.,  0.,  0.,  0.,  0.,  0., -1., -1.],            [ 2.,  3.,  1.,  0.,  0.,  0.,  0., -1., -3., -2.],            [ 1.,  3.,  3.,  1.,  0.,  0., -1., -3., -3., -1.],            [ 0.,  1.,  3.,  3.,  1., -1., -3., -3., -1.,  0.],            [ 0.,  0.,  1.,  3.,  2., -2., -3., -1.,  0.,  0.],            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],            [ 0.,  0., -1., -3., -2.,  2.,  3.,  1.,  0.,  0.],            [ 0., -1., -3., -3., -1.,  1.,  3.,  3.,  1.,  0.],            [-1., -3., -3., -1.,  0.,  0.,  1.,  3.,  3.,  1.],            [-3., -4., -1.,  0.,  0.,  0.,  0.,  1.,  4.,  3.]])      Since I want to use only zeroes and ones in the representation of the decision boundary, I can do the following:               >>> (np.abs(sobel(m)) > 0).astype(int)     array([[1, 1, 0, 0, 0, 0, 0, 0, 1, 1],            [1, 1, 1, 0, 0, 0, 0, 1, 1, 1],            [1, 1, 1, 1, 0, 0, 1, 1, 1, 1],            [0, 1, 1, 1, 1, 1, 1, 1, 1, 0],            [0, 0, 1, 1, 1, 1, 1, 1, 0, 0],            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],            [0, 0, 1, 1, 1, 1, 1, 1, 0, 0],            [0, 1, 1, 1, 1, 1, 1, 1, 1, 0],            [1, 1, 1, 1, 0, 0, 1, 1, 1, 1],            [1, 1, 1, 0, 0, 0, 0, 1, 1, 1]])      As you can see, now I have something that resembles what I would like to obtain, but still not: the decision boundary is too thick, almost three times that what is necessary. I could attempt to solve the problem by setting a different threshold value:               >>> (np.abs(sobel(m)) > 2).astype(int)     array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],            [0, 1, 0, 0, 0, 0, 0, 0, 1, 0],            [0, 1, 1, 0, 0, 0, 0, 1, 1, 0],            [0, 0, 1, 1, 0, 0, 1, 1, 0, 0],            [0, 0, 0, 1, 0, 0, 1, 0, 0, 0],            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],            [0, 0, 0, 1, 0, 0, 1, 0, 0, 0],            [0, 0, 1, 1, 0, 0, 1, 1, 0, 0],            [0, 1, 1, 0, 0, 0, 0, 1, 1, 0],            [1, 1, 0, 0, 0, 0, 0, 0, 1, 1]])      The problem is that this threshold value ($2$) is dependent on the context. What if the categories are in larger number, or differently numbered, or if I have an arbitrary number of dimensions? This is clearly not working. Do you have any ideas on how to solve this problem? Maybe another filter for the convolution or a different algorithm? thank you! :)


title: Two way scatter plot with time dimension

text: I am plotting GDP against exports for 134 countries. I have data for three years. For each year I can easily plot these scatter dots and see the relative country positions but I am interested in the change in positions in the graph for each country. In other words, I'd ideally like to plot GDP against exports for all countries for three years in one graph with smooth lines connecting three years for each country showing how their relative positions changed overtime. below is an example of such graph in 2D. Assume exports in on X axis and GDP is on Y axis then arrows will suggest how position of each country changed over time in the graph. countries on top right corner will show high level of exports and high level of GDP and vice versa. I am interested in figuring out how countries change their long term growth strategies over time w.r.t exports share in GDP. Any idea how to do this? ![enter image description here](http://i.stack.imgur.com/2HJUI.jpg)


title: If not minimizing SSEs, what am I doing?

text: I've got a problem that can be summarized as a linear regression. Thus it takes the following form: $$ Y=X \beta +\epsilon $$ where $Y$ and $\epsilon$ are vectors of size $N\times1$, $X$ is a matrix of size $N\times3$, and $\beta$ is a vector of size $3\times1$ (i.e., I am solving for three parameters). The proper way to solve this under an ordinary least squares framework is via the normal equations—that is, taking the partial derivative of the sum of squared errors with respect to each parameter, and solving for where that equals zero in order to minimize the sum of squared errors. However, this works rather poorly for me. On a whim, I tried just doing the following to solve for the parameters instead: $$ \beta = X^{-1}Y $$ This actually works out really, really nicely. The problem is, I can't justify it, and I don't know what it's doing! Note that it's basically ignoring the idea that there are any errors at all. Any ideas as to what's going on here? Thanks in advance.


title: Detecting subset of high sensitivity corresponding to genotype

text: ## Background There is a data set from a clinical trial, involving about 1000 patients received an behavior intervention. The intervention is supposed to affect a variable x1, and the change in x1 results the change in the response variable y. The data also contains the information of several(>10) snips for each patient. ## Analysis planned Given a subset of patients, fit y=beta1*x1+beta2*other variable+… I would like to find certain combinations of genotype, such that the patients with these genotypes have significantly high beta1. For each subset (i.e. one geno combination), the p-value of beta1 is calculated and I plan to perform some adjustment to them since all p's are correlated. What I did below is intended to see if the linear model, with selected variables, is appropriate. * * * ## Original question (before comments and edits) I performed thousands of tests and obtained p-values. Then I did two-sided uniform KS test of the p-values, the result claimed it's uniform. So does it mean that my model are wrong? Because I expect more small p-values near 0. This is a preliminary step before correcting the multiplicity. ![http://www.flickr.com/photos/96546690@N02/9352516904/](http://i.stack.imgur.com/lJtiO.png) The ks test:               One-sample Kolmogorov-Smirnov test          D = 0.0493, p-value = 1.388e-06          alternative hypothesis: two-sided      * * *


title: Under which conditions do gradient boosting machines outperform random forests?

text: Can Friedman's gradient boosting machine achieve better performance than random forests? If so, in which conditions or what kind of data set can make gbm better?


title: Negative binomial regression question - is it a poor model?

text: I am reading a very interesting article by Sellers and Shmueli on regression models for count data. Near the beginning (p. 944) they cite McCullaugh and Nelder (1989) saying that negative binomial regression is unpopular and has a problematic canonical link. I found the referred passage and it says (p. 374 of M and N) "Little use seems to have been made of the negative binomial distribution in applications; in particular, the use of the canonical link is problematical because it makes the linear predictor a function of a parameter of the variance function". On the previous page they give that link function as $$\eta = \log\left(\frac{\alpha}{1 + \alpha} \right) = \log\left( \frac{\mu}{\mu + k}\right)$$ and variance function $$V = \mu + \frac{\mu^2}{k}.$$ The distribution is given as $$Pr(Y = y; \alpha,k) = \frac{(y+k-1)!}{y!(k-1)!}\frac{\alpha^y}{(1+\alpha)^{y=k}}$$ I have found NB regression to be quite widely used (and recommended in several books). Are all these uses and recommendations in error? What are the consequences of this problematic link?


title: Weighting data sources in a bayesian model (BUGS)

text: I use a state space model to fit observations to a population dynamic model (using the BUGS language). In the "state" part, the dynamic model create a new "state" of the population (i.e. size and structure) for each year modeled. In the "observation" part, several outputs of the "state" part are fitted to observations given by two different sources. The size of the population, each year there was observation, is fitted using a lognormal.               for (i in 1:N) {        obsPop[i] ~ dlnorm(log(predictPop[yearObsPop[i]]),prec.obsPop[i])     }      The number of individual observed dead, each year there was observation, is fitted using a binomial               for (i in 1:N2){        obsDead[i] ~ dbin(pObsDead, deadPredict[yearObsDead[i]])     }      However, N (i.e. the number of year of observation of the population size) is smaller than N2 (i.e. the number of year of observation for mortalities) and the second source of data seems to have an excessive weight in the result.   * How can I do to weight differently each source of data in the global fitting of the model ? Thanks by advance !


title: SEM/CFA and data transformation

text: In order to linearize my model I use $\ln(x+1)$, so I come up with the following equation: $$ \ln(y+1) = b_0 + b_1 \ln(x_1 + 1) + b_2\ln(x_2+1) + \cdots + b_n\ln(x_n+1). $$ I hypothesize based on latent (unobserved variables), and as measures I use observed variables (or proxies) - the ones in the equation. My plan is to use SEM (Structural equation modeling) or/and CFA (confirmatory factor analysis) to validate my model. In doing my SEM/CFA should I use the transformed variables $\ln(x+1)$ or as they are observed (just $x$)?


title: Understanding Odds Ratios in Logistic Regression

text: I'm trying to understand how to interpret log odds ratios in logistic regression. Let's say I have the following output:               > mod1 = glm(factor(won) ~ bid, data=mydat, family=binomial(link="logit"))     > summary(mod1)          Call:     glm(formula = factor(won) ~ bid, family = binomial(link = "logit"),          data = mydat)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -1.5464  -0.6990  -0.6392  -0.5321   2.0124            Coefficients:                   Estimate Std. Error z value Pr(>|z|)         (Intercept) -2.133e+00  1.947e-02 -109.53   <2e-16 ***     bid          2.494e-03  5.058e-05   49.32   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for binomial family taken to be 1)              Null deviance: 83081  on 80337  degrees of freedom     Residual deviance: 80645  on 80336  degrees of freedom     AIC: 80649          Number of Fisher Scoring iterations: 4      So my equation would look like: $$\Pr(Y=1) = \frac{1}{1 + \exp\left(-[-2.13 + 0.002\times(\text{bid})]\right)}$$ From here I calculated probabilities from all bid levels. ![enter image description here](http://i.stack.imgur.com/5mLa9.png) I have been using this graph to say that at a 1000 bid, the probability of winning is x. At any given bid level, the probability of winning is x. I have a feeling that my interpretation is wrong because I'm not considering that these are log-odds. How should I really be interpreting this plot/these results?


title: how to compute minimum required vc dimension for a classifer to classify a specific data

text: Suppose we're given an N dimensional data to classify. To cope with this task we may choose a classifier that suits our desires more. However obviously not every classifier is capable of classifying every data for many reasons including its VCD. Now I want to know is there a way to compute the minimum required VC-dimension for a classifier(Hypothesis) to classify a specific data? As an example to make my question more clear, a classifier having VC-dimension of 2 is not capable of shattering this data set:               +++    ---    +++


title: Compare distributions with uneven and small sample sizes

text: _**I will be performing my analysis in R but any general mathematical-based answers will be appreciated._** I realize similar questions have been asked before but I am new to statistics so please excuse my ignorance. I do not know enough to be sure that approaches used in other threads I viewed are valid for my case so I would like to be sure. Here is the context, simplified for those unfamiliar with this kind of data: I am looking at gene expression data. Essentially, I have a number of categories that each have a set of genes in them. Each gene has an associated expression value (amount present in my sample, **log2 transformed** ). I also have three different test conditions (these conditions are completely independent of one another!). What this amounts to is a distribution. As an example: _Biological Process A:_   * Contains X **genes** and their **expression levels** for **Condition 1**   * Contains Y **genes** and their **expression levels** for **Condition 2**   * Contains Z **genes** and their **expression levels** for **Condition 3** ![Example showing my data set up](http://i.stack.imgur.com/T4xX8.png) I would like to **_compare_** these three distributions. By compare I mean that I want to see if there are any large shifts in either the shape or spread of the distribution. **EDIT: as Peter has pointed out, what I am really interested in is 1. Location, 2) Spread and 3) Shape with regards to the three distributions under comparison**. For example, if Condition 1 moved such that it was almost entirely below the x-axis and Condition 2 was entirely above it, that would be a very interesting shift. Additionally, if two Conditions went above the axis and one went down, this would also be an interesting shift. Essentially, I need a way to compare the distributions for the biggest changes in mean density above or below the x-axis. I have about 400+ "Biological Processes". The problem is that for some of them, the sample sizes may vary considerably (genes (n) = 3, 25, 10 for the three conditions respectively would be an extreme example in my case) so I do not think simply comparing the means is a fair approach. Visual comparison will take a very long time with 400 beanplots to look at and I would also like some statistics to back up the conclusions I may be making. I want to essentially say something like "Treating tissue with Condition 1 had the effect of increasing expression for most genes known to be involved in Biological Process A". Additionally, if this approach is invalid or there is a better way to look at this data I would be open to suggestions. Thank you.


title: Meta-analysis of prevalence at the country level

text: I'm working on a meta-analysis of prevalence data. The aim is to get estimates of prevalence at the country level. The main issue is that the disease is highly correlated with age, and the sample ages of included studies are highly heterogeneous. Only median age is available for most studies, so I can't use SMR-like tricks. I figured I could use meta-regression to solve this, including age as a fixed-effect and introducing study-level and country-level random-effects. The idea (that I took from Fowkes et al) was to use this model to make country-specific predictions of prevalence for each 5-year age group from 15 to 60 (using the median age of the group), and to apply these predictions to the actual population size of each of those groups in the selected country, in order to obtain total infected population and to calculate age-adjusted prevalence in the 15-60 population from that. I tried several ways to do this using R with packages `meta` and `mgcv`. I got some satisfying results, but I'm not that confident with my results and would appreciate some feedback. First is some simulated data, then the description of my different approaches:               data<-data.frame(id_study=c("UK1","UK2","UK3","FRA1","FRA2","BEL1","GER1","GER2","GER3"),                      country=c("UK","UK","UK","FRANCE","FRANCE","BELGIUM","GERMANY","GERMANY","GERMANY"),                      n_events=c(91,49,18,10,50,6,9,10,22),                      n_total=c(3041,580,252,480,887,256,400,206,300),                      study_median_age=c(25,50,58,30,42,26,27,28,36))      **Standard random-effect meta-analysis** with package `meta`. I used `metaprop()` to get a first estimate of the prevalence in each country without taking age into account, and to obtain weights. As expected, heterogeneity was very high, so I used weights from the random-effects model.                meta <- metaprop(event=n_events,n=n_total,byvar=country,sm="PLOGIT",method.tau="REML",data=data)      summary(meta)      data$weight<-meta$w.random      I used meta to get a first estimate of the prevalence without taking age into account, and to obtain weights. As expected, heterogeneity was very high, so I used weights from the random-effects model. **Generalized additive model** to include age with package `mgcv`. The `gam()` model parameters (k and sp) were chosen using BIC and GCV number (not shown here).                model <- gam( cbind(n_events,n_total-n_events) ~ s(study_median_age,bs="cr",k=4,sp=2) + s(country,bs="re"), weights=weight, data=data, family="binomial"(link=logit), method="REML")      plot(model,pages=1,residuals=T, all.terms=T, shade=T)      Predictions for each age group were obtained from this model as explained earlier. CI were obtained directly using `predict.gam()`, that uses the Bayesian posterior covariance matrix of the parameters. For exemple considering UK:                newdat<-data.frame(country="UK",study_median_age=seq(17,57,5))      link<-predict(model,newdat,type="link",se.fit=T)$fit      linkse<-predict(model,newdat,type="link",se.fit=T)$se      newdat$prev<-model$family$linkinv(link)      newdat$CIinf<-model$family$linkinv(link-1.96*linkse)      newdat$CIsup<-model$family$linkinv(link+1.96*linkse)      plot(newdat$prev~newdat$study_median_age, type="l",ylim=c(0,.12))      lines(newdat$CIinf~newdat$study_median_age, lty=2)      lines(newdat$CIsup~newdat$study_median_age, lty=2)      The results were satisfying, representing the augmentation of the prevalence with advanced age, with coherent confidence intervals. I obtained a total prevalence for the country using the country population structure (not shown, I hope it is clear enough). However, I figured I needed to include study-level random-effects since there was a high heterogeneity (even though I did not calculate heterogeneity after the meta-regression). **Introducing study-level random-effect** with package `gamm4`. Since `mgcv` models can't handle that much random-effect parameters, I had to switch to `gamm4`.                model2 <- gamm4(cbind(n_events,n_total-n_events) ~ s(study_median_age,bs="cr",k=4) + s(country,bs="re"), random=~(1|id_study), data=data, weights=weight, family="binomial"(link=logit))      plot(model2$gam,pages=1,residuals=T, all.terms=T, shade=T)           link<-predict(model2$gam,newdat,type="link",se.fit=T)$fit      linkse<-predict(model2$gam,newdat,type="link",se.fit=T)$se      newdat$prev2<-model$family$linkinv(link)      newdat$CIinf2<-model$family$linkinv(link-1.96*linkse)      newdat$CIsup2<-model$family$linkinv(link+1.96*linkse)      plot(newdat$prev2~newdat$study_median_age, type="l",col="red",ylim=c(0,0.11))      lines(newdat$CIinf2~newdat$study_median_age, lty=2,col="red")      lines(newdat$CIsup2~newdat$study_median_age, lty=2,col="red")      lines(newdat$prev~newdat$study_median_age, type="l",ylim=c(0,.12))      lines(newdat$CIinf~newdat$study_median_age, lty=2)      lines(newdat$CIsup~newdat$study_median_age, lty=2)      Since the study-level random effect was in the mer part of the fit, I didn't have to handle it. As you can see, I obtain rather different results, with a much smoother relation between age and prevalence, and quite different confidence intervals. It is even more different in the full-data analysis, where the CI are much wider in the model including study-level RE, to the point it is sometimes almost uninformative (prevalence between 0 and 15%, but if it is the way it is...). Moreover, the study-level RE model seems to be more stable when outliers are excluded. **So, my questions are:**   * Did I properly extract the weights from the metaprop() function and used them further?   * Did I properly built my `gam()` and `gamm4()` models? I read a lot about this, but I'm not used to this king of models.   * Which of these models should I use? I would really appreciate some help, since neither my teachers nor my colleagues could. It was a really harsh to conduct the systematic review, and very frustrating to struggle with the analysis... Thank you in advance!


title: Weak convergence implies uniform convergence in distribution?

text: Consider an empirical processes $\\{v_T(\theta), T\geq 1 \\}$ and assume it is weakly convergent to the stochastic process $\\{ v(\theta); \theta \in \Theta\\}$., i.e. $$ E^*(f(v_T(\theta)))\rightarrow E(f(v(\theta))) $$ as $T\rightarrow \infty$ $\forall f$ bounded uniformly continuous function in $\Theta$ and $\forall \theta \in \Theta$ where $E^*$ denote the outer expectation. I know that this implies convergence in distribution $\forall \theta \in \Theta$ as $T\rightarrow \infty$, i.e. $v_T(\theta)=O_p(1)$ $\forall \theta \in \Theta$, i.e. $$ \lim_{T \rightarrow \infty} F_T(a; \theta)=F(a;\theta) $$ $\forall a$ at which $F$ is continuous, where $F_T(\cdot;\theta)$ is the cdf of $v_T(\theta)$ and $F(\cdot;\theta)$ is the cdf of $v(\theta)$. Question: Does this imply also $$ \sup_{\theta \in \Theta}v_T(\theta)=O_p(1) $$ ?


title: Estimation of parameter when having two sources of errors

text: First of all sorry for the bad title. I picture this is a standard problem, i just don't know how to deal with it I have the following model \begin{equation} Y_{ij} = \mu + E_i + \varepsilon_{ij} \end{equation} Where $1 \leq i \leq N$ and $1 \leq j \leq M$. And $E_i \sim \mathcal{N}(0,\sigma^2_E$) and $\varepsilon_{ij}\sim \mathcal{N}(0,\sigma^2_\varepsilon)$ are error terms all independent among each other. Note that $E_i$ appear in many equations and $\varepsilon_{ij}$ in only one. I would like to know how is the standard way to estimate $\mu$ and the variance of the estimation, after observing all $Y_{ij}$. My first approach is to estimate by the mean \begin{equation} \hat{\mu} =\frac{1}{MN} \sum_{i=1}^{N}\sum_{j=1}^M Y_{ij} \end{equation} I am able to calculate the variance of $\hat{\mu}$, that yields \begin{equation} \text{Var}(\hat{\mu}) = \frac{\sigma^2_{\varepsilon}}{MN} + \frac{\sigma^2_E}{N} \end{equation} I can manage to estimate $\sigma^2_\varepsilon$ with \begin{equation} \hat{\sigma^2_\varepsilon} = \frac{1}{N(M-1)}\sum_{i=1}^M\sum_{j=1}^N(Y_{ij} - \bar{Y_i})^2 \end{equation} But I am not sure how to estimate $\sigma^2_E$.


title: Prediction error in least squares with a linear model

text: In the classical linear model with $$Y=X\beta +\epsilon,$$ where $Y \in \mathbb{R}^n$ is the observation, $X\in \mathbb{R}^{n\times p}$ is the known covariates, $\beta \in \mathbb{R}^p$ is the unknown parameter with, $p < n$, and $\epsilon \in \mathbb{R}^{n}$, $\epsilon \sim \mathcal{N}(0,\sigma^{2}I)$. The classical least squares estimator here would be $$\hat{\beta}= (X^TX)^{-1}X^TY.$$ If $\beta^{0}$ is the true parameter, then we have that the prediction error is given by $E=||X(\hat{\beta}- \beta^{0})||_2^2$. We have that, $$E=||X(\hat{\beta}- \beta^{0})||_2^2/\sigma^2=||X(X^TX)^{-1}X^T \epsilon||_2^2/\sigma^2=||\gamma||_2^2/\sigma^2.$$ It is easy to see that $\gamma \sim \mathcal{N}(0,\sigma^2 X(X^TX)^{-1}X^T)$. However it is claimed in `High Dimensional Statistics` by `Bühlmann and van de Geer` (on page 101) that $E$ is distributed according to a chi-square distribution with $p$ degrees of freedom. I can not see how this is true (It would be true if $\gamma \sim \mathcal{N}(0,D)$ for some diagonal matrix $D$ with non-negative diagonal entries.) Am I missing something here?


title: Regression model result interpretation

text: I am trying to learn few things about generalized least square regression modeling.   Here are the models that I am using., X1 ~ Y1 + Y2, X2 ~ Y1 + Y4, All predictor variables, Y1, Y2 and Y4 are significantly related with the response variables, but when I took the difference of two response variables and the difference of two predictor variables, something like this: X3 = X1 - X2, Y5 = Y2 - Y4, and use this model X3 ~ Y1 + Y5 Now, Y5 is not significant: I am confuse about it. As, Y5 is just a difference of Y2 and Y3 and X3 is the difference of X1 and X2. Can someone guide me?


title: Community stability measure in network analysis using iGraph in R

text: Does anyone know of a "community stability" measure that is already implemented in iGraph/R?? My approach is to permute the adjacency matrix and compute the association in the cross-frequency table between the original community membership and the membership in the new community structure from the perturbed network (where perturbations are single add/removal of edges). The problem is that I use Fisher-exact test with simulated p-value and all the results are the same.               GetStability <- function(ingraph,pv_parm=0.05,permutations= 1000,MC_Replicates=10000){            tempcomm <- fastgreedy.community(ingraph)         V(ingraph)$commMbr <-  tempcomm$membership            # get communities of initial graph       orig.communities.df <- data.frame(node=V(ingraph)$label, origComm = V(ingraph)$commMbr)            # get adjacency of graph       A <- get.adjacency(ingraph, type=c("upper"),edges=TRUE)       dimA <- dim(A)[1]       for(i in 1:(dimA-1)){          for(j in (i+1):dimA){             ij <- A[i,j]              if(ij > 0){                A[i,j] <- as.numeric(E(ingraph)$weight[ A[i,j] ])              }          }      }            numUpperDiagonals <- dimA*(dimA-1)/2            if(permutations > numUpperDiagonals){         permutations <- numUpperDiagonals       }            sampleIndices <- sample(which(upper.tri(A, diag = TRUE) == TRUE ), permutations)       all.p.values <- c()       for(sampleIndex in sampleIndices){               # sampleIndex <-  sampleIndices[1]          tempA <- A               tempA[sampleIndex] <- ifelse(tempA[sampleIndex] > 0, 0,     ifelse(tempA[sampleIndex]==0,1, 0))               rownames(tempA) <- V(ingraph)$label          #now set temp subgraph with new matrix          tempgraph <- graph.adjacency(tempA, mode = "upper", weighted=NULL, diag=FALSE)               V(tempgraph)$label <- rownames(tempA)              tempcomm <- fastgreedy.community(tempgraph)              V(tempgraph)$commMbr <-  tempcomm$membership               pert.communities.df <- data.frame(node=V(tempgraph)$label, pertComm = V(tempgraph)$commMbr)               merged.df <- merge(orig.communities.df, pert.communities.df, by = "node")               freqTable <- table(merged.df$pertComm,merged.df$origComm)               p.value <- NULL          try(p.value <- fisher.test(freqTable, simulate.p.value=TRUE, B=MC_Replicates)$p.value,silent = TRUE)               if(! is.null(p.value) ){             all.p.values <- c(all.p.values, p.value)           }             }            returnThisList <- list()       percentiles <- quantile(all.p.values,probs=c(0.025,0.5, 0.975))            returnThisList$all.p.values <- all.p.values           returnThisList$LB <- percentiles[1]       returnThisList$estimate <- percentiles[2]           returnThisList$UB <- percentiles[3]            returnThisList$stability <- sum(all.p.values < pv_parm) / length(all.p.values)       return(returnThisList)     }


title: Can a magician bias the result of a fair die?

text: Some time ago, Persi Diaconis was discussing probability and made the point that the level of information a person possessed about an event played a role in his accessing probabilities. To illustrate this point, he did a very nice demonstration showing how a magician can toss a fair coin in such a way that he will always get heads. **(See 2:10 in a video of Diaconishttp://www.youtube.com/watch?v=ZdA7TeoZD_g)** I am curious if this is possible using a fair die? TIA, Matt


title: How to dynamically calculate interval designations of data?

text: I am looking for an equation to calculate appropriate intervals to designate values as very low, low, medium, high, very high for dynamic data. for example if my dataset ranges from 0 to 100 very low= <=20, low=21-40, medium=41-60, high=61-80, very high= >81 However, if the range is small there should be fewer intervals. For example, if the range is 20 then there should be a low and high interval designation only. the intervals need to be logical intervals so they can be used for a graphical display. logical intervals would include intervals of 5, 10, 25, 50, 100 etc depending on the dataset. Is there an equation that can determine appropriate intervals for each unique dataset?


title: Outputting multiple maximum length sequences

text: So, I am mining patterns in a sequence like, ("a0cc0vaaaabaaadbaaabbaa00bvw"). Now, I came up with this. So, I am mining patterns in a sequence like, `("a0cc0vaaaabaaadbaaabbaa00bvw")`. Now, I came up with this.               reps <- function(s, n) paste(rep(s, n), collapse = "") # repeat s n times          find.string <- function(string, th = 3, len = floor(nchar(string)/th)) {         found <- FALSE         for(sublen in len:1) {             for(inlen in 0:sublen) {                 pat <- paste0("((.{", sublen-inlen, "})(.)(.{", inlen, "}))", reps("(\\2.\\4)", th-1))                 r <- regexpr(pat, string, perl = TRUE)                 if (attr(r, "capture.length")[1] > 0){                     found = TRUE                     break;                 }             }             if(found) break         }         if (r > 0) substring(string, r, r + attr(r, "capture.length")[1] - 1) else ""     }          find.string("a0cc0vaaaabaaadbaaabbaa00bvw"); # returns "aaaab"      All of this is put inside a function and when I run this on the input string, I get the output as `aaaab`. All is good uptill now. However, my function fails to work on input like this. `110110111abcabcabd`. Here, it outputs only one value, i.e. `110`. Now, I want all the maximum length patterns in the string. How can I modify this? A pattern is the maximum substring which appears consecutively for atleast three times. The substring can have noise of maximum one character. For eg: If the sequence is, 110110111, here the output will be 110


title: G/G/1 Simulation help

text: i am a complete nooby in statistics, therefore the question i am about to ask may sound very unprofessional **warning** :) so, we have a task to simulate a G/G/1 queueing system for particular parameters (Erlang distribution is assumed). I have been able (hopefully) to find out what mean and number of stages should be. We are using the QtsPlus for the simulation and supposed to do it for 10000 clients. QtsPlus, however, has following parameter fields:   * Number of Transactions to Simulate   * Checkpoint Message Frequency   * Maximum Probability Size and i am very confused about it. Which field corresponds to the number of clients and what values should i enter in the other fields? thanks a lot


title: Best way to characterize slowdown factor

text: I am currently measuring the execution time of a program I have written, but I am unsure how to represent the results. The program has two modes: normal mode and analysis mode. In analysis mode, it always runs slower than in normal mode. I am interested in the slowdown factor incurred by the analysis mode. I have measured the execution time for different program inputs, and noted the execution time for normal mode ($b_i$) and analysis mode ($a_i$). Basically, I can think of two ways to characterize the slowdown: Average slowdown = $AS$ = $\frac{1}{n}(\frac{a_1}{b_1} + ... + \frac{a_n}{b_n})$ Compound slowdown = $CS$ = $\frac{a_1 + ... + a_n}{b_1 + ... + b_n}$ Which method is preferrable and why? I have seen that either of $AS$ and $CS$ can be bigger, depending on the concrete values of the $a_i$ and $b_i$. Are there intuitive and precise conditions on $a_i, b_i$ for $AS > CS$ and vice- versa?


title: Differential Entropy of Gaussian Process

text: I have $N$ datapoints that have $d$ features in a GP and their covariance matrix $K$ and I want to calculate the differential entropy of that GP. Is this formula right? $E(I)= \frac{1}{2} \log((2πe)^d \det(K))$ Moreover, given $L$ the cholesky decomposition of $K$ is there a faster way? Thank you


title: Cross validating quantile regression

text: I applied quantile regression on some data and did it for tau = 0.25, 0.5, 0.75. After i got the estimates of each model, i did some cross validation on my hold out data. When i used the estimates for my hold out data i calculated in what quantile each of my independent variables lie, 25%, 50% or 75%. For instance if my x-values (independent variable) looks as follows, `x -> 1 2 3 4 5 6 7 8 9 11 10 12 13 14 15` I calculated my quantiles using the R function `cut` as `cut(x,4)` and got the following quantile ranges (0.986; 4.49], (4.49; 8], (8; 11.15], (11.15; 15]. I then divide the x vector into the correct quantile range and depending on what range it falls determine what quantile regression estimates i use, either 0.25, 0.5 or 0.75. So for the x-values 1 to 4 i will use the estimates of the 0.25, 5 to 8 i use the estimates of 0.5 and the rest i used 0.75. I do get very accurate results, i.t.o. forecast accuracy, but can someone please tell if this is the correct way to do my quantile cross validation.


title: Support Vectors Not Falling on Margin Lines for e1071 and kernlab packages in R?

text: In a previous thread, Computing the Decision Boundary of a linear SVM model(Computing the decision boundary of a linear SVM model), the following R code was given as a way to compute the formula of the hyperplane and it's margins, given an input set of hyperparameters:               library(kernlab)         set.seed(101)     x <- rbind(matrix(rnorm(120),,2),matrix(rnorm(120,mean=3),,2))     y <- matrix(c(rep(1,60),rep(-1,60)))     svp <- ksvm(x,y,type="C-svc")     plot(svp,data=x)     alpha(svp)  # support vectors whose indices may be found with alphaindex(svp)     b(svp)      # (negative) intercept      plot(scale(x), col=y+2, pch=y+2, xlab="", ylab="")     w <- colSums(coef(svp)[[1]] * x[unlist(alphaindex(svp)),])     b <- b(svp)     abline(b/w[1],-w[2]/w[1])     abline((b+1)/w[1],-w[2]/w[1],lty=2)     abline((b-1)/w[1],-w[2]/w[1],lty=2)      However, in the resulting plot, none of the support vectors actually fall on either margin line. If               abline(b/w[1],-w[2]/w[1]) #is supposed to be the maximum-marginal hyperplane and     abline((b+1)/w[1],-w[2]/w[1],lty=2) #is supposed to be the yi = +1 margin and     abline((b-1)/w[1],-w[2]/w[1],lty=2) #is supposed to be the yi = -1 margin      then why don't the support vectors, `alpha(svp)`, fall on either the +1 or the -1 margin? Now before you say it's because `kernelf = "rbfkernel"` by default and our hyperplane and margins have not been transformed according to the `rbfkernel`, I can explicitly call the linear kernel with `kernelf = "vanilladot"` (i.e `svp <- ksvm(x,y,type="C-svc",kernelf="vanilladot")`) and the hyperplane and margins are even farther from the support vectors defined by `alpha(svp)` then when `kernelf="rbfkernel"`. Note that I find similar results when I use the `e1071` package (which is perhaps not surprising as they both based on the `LIBSVM` library). If I am doing something wrong, or missing a crucial point (like perhaps `ksvm()` uses a soft margin by default and is hiding the slack variables from me), please let me know! Here is my code for the linear kernel with the `kernlab` package:               library(kernlab)     set.seed(101)     x <- rbind(matrix(rnorm(120),,2),matrix(rnorm(120,mean=3),,2))     y <- matrix(c(rep(1,60),rep(-1,60)))     svp <- ksvm(x,y,type="C-svc",kernelf="vanilladot")     plot(svp,data=x)     plot(x, col=y+2, pch=y+2, xlab="", ylab="")     w <- colSums(coef(svp)[[1]] * x[unlist(alphaindex(svp)),])     b <- b(svp)     plot(svp,data=x)     abline(b/w[1],-w[2]/w[1])     abline((b+1)/w[1],-w[2]/w[1],lty=2)     abline((b-1)/w[1],-w[2]/w[1],lty=2)


title: Can I use permutation tests for to avoid the multiple comparison problem in the context of proportions?

text: I am evaluating the effectiveness of 5 different methods to predict a particular binary outcome (call them 'Success' and 'Failure'). The data look like so:               Method    Sample_Size    Success    Percent_Success     1         28             4          0.14       2         19             4          0.21       3         24             7          0.29       4         21             13         0.61       5         22             9          0.40       I would like to perform a test among these 5 methods to assess the relative superiority of the methods. In other words, I want to order the methods in order of performance as method 1 > method 2 > ... method 5. To avoid the issue of multiple comparisons, I plan to do a permutation test along the following lines: Step 1: Pool all the data so that the overall sample size is 114 with overall 37 successes. Step 2: Randomly split the data into 5 groups with the corresponding sample sizes of 28, 19, 24, 21 and 22. Step 3: Increment a counter if the observed order of Percent_Success from step 2 is consistent with the ordering of my data. Step 4: Repeat steps 2 and 3 many times (say 10000). Desired p-value = Final Counter Value / 10000. Questions:   1. Is the above procedure ok?   2. Is there anything in R that would enable me to perform the above test?   3. Any suggestions for improvement or alternative methods would be helpful.


title: Inter-rater reliability: Nominal data, but *multiple ratings per unit/subject*, different raters

text: I have a data set, with each variable taking multiple values on a nominal scale. Separate raters could rate a given unit using more than one value per variable. That is, there are multiple ratings per unit/rater. How to proceed from here? I am aware that there are extensions of kappas (Fleiss) to support multiple ratings per unit/subject. Basically, the extension suggests computing separate Kappas per scale value, and then forming an aggregate. Are there alternative strategies? How is the combined kappa exactly computed, and the auxiliaries in this case: standard error, ...? (Unfortunately, right now, I do not have access to a copy of Fleiss' "Statistical methods for rates and proportions" because our library moves to a new location :/ Ch. 18 contains the necessary details) Besides: How would this apply to Krippendorff's alpha (which I am also considering)?


title: Probability of not drawing a word from a bag of letters in Scrabble

text: Suppose you had a bag with $n$ tiles, each with a letter on it. There are $n_A$ tiles with letter 'A', $n_B$ with 'B', and so on, and $n_*$ 'wildcard' tiles (we have $n = n_A + n_B + \ldots + n_Z + n_*$). Suppose you had a dictionary with a finite number of words. You pick $k$ tiles from the bag without replacement. How would you compute (or estimate) the probability that you can form zero words from the dictionary given the $k$ tiles selected? For those not familiar with Scrabble (TM), the wildcard character can be used to match any letter. Thus the word [ **BOOT** ] could be 'spelled' with the tiles 'B', '*', 'O', 'T'. To give some idea of the scale of the problem, $k$ is smallish, like 7, $n$ is around 100, and the dictionary contains about 100,000 words of size $k$ or smaller. **edit:** By 'form a word', I mean a word of length no greater than $k$. Thus, if the word [ **A** ] is in the dictionary, then by drawing even a single 'A' from the bag, one has 'formed a word'. The problem of wildcards is radically simplified if one can assume there are words of length 1 in the dictionary. For if there are, any draw of a wildcard automatically can match a length 1 word, and thus one can concentrate on the case where there are no wildcards. Thus the more slippery form of the problem has no 1-letter words in the dictionary. Also, I should explicitly state that the order in which the letters are drawn from the bag is immaterial. One does not have to draw the letters in the 'correct' order of the word.


title: Feature construction in R

text: I am wondering if there are any algorithms (perhaps genetic algorithms) in R for feature construction (deriving candidate predictors from existing predictors)? I am thinking of a routine to test higher order powers, interactions, ratios, and linear combinations and nonlinear functions of existing variables (sin, cos, atan etc). This could be a filter or wrapper routine (i.e. not using a learning algorithm or using one to define fitness of the feature). My goal is to "discover" potentially meaningful ratios and the like of existing predictors. Thanks!